Image workflows

Image workflows#

Packages#

import ee
import geemap

try:
    ee.Initialize()
except Exception as e:
    ee.Authenticate()
    ee.Initialize()

Set up map window#

ndviVis = {"min":0.2, "max":1, "palette": ['white', 'green']}

Map = geemap.Map(center=[53, 12], zoom=8)

Create an image#

Let’s create two image objects from the GEE data catalog, specifically two NDVI composites from August and April.

ndviAugust = ee.Image('LANDSAT/LC8_L1T_32DAY_NDVI/20150813')
ndviApril = ee.Image('LANDSAT/LC8_L1T_32DAY_NDVI/20150407')

Map.addLayer(ndviAugust , ndviVis, 'NDVI August')
Map.addLayer(ndviApril , ndviVis, 'NDVI April')

To understand what methods are available for an ee.Image object, you may consult the API documentation (see Client libraries). For example, the addBands() method adds a band to an image object. You must access the method from your image object not from ee.Image(), i.e., ee.Image.addBands() does not work.

ndviStack = ndviAugust.addBands(ndviApril)

Image math#

Image-based math operators are accessed via the created image object, e.g., image.add(), image.subtract, image.multiply(), image.divide(). Let’s subtract the April NDVI from the August NDVI to calculate the change in NDVI.

ndviChange = ndviAugust.subtract(ndviApril)
Map.addLayer(ndviChange, {}, 'delta NDVI')

Logical operations#

We can create a binary mask based on the NDVI change. For example, we might want to separate areas that had an increase in NDVI by 0.2 (or greater) from areas with stable or declining NDVI (NDVI < 0.2).

ndviIncreaseMask = ndviChange.gte(0.2)
Map.addLayer(ndviIncreaseMask, {}, 'Where NDVI Increased')

Masking#

We can apply the mask using ee.Image.mask() or ee.Image.updateMask(). For example, if we want to mask out stable NDVI areas from the NDVI change map, we can do the following:

ndviNoIncreaseMask = ndviChange.lte(0.2)
ndviChangeMasked = ndviChange.updateMask(ndviNoIncreaseMask)
Map.addLayer(ndviChangeMasked, {}, 'NDVI Increase')

Image properties#

Images conatain several types of metadata information such as band names, projection information, properties, and other metadata.

Get the band names as a list.

bandNames = ndviAugust.bandNames()
print('Band names: ', bandNames.getInfo())

Band names:  ['NDVI']

Get a list of all metadata properties.

properties = ndviAugust.propertyNames()
print('Metadata properties: ', properties.getInfo())

Metadata properties:  ['system:time_start', 'system:time_end', 'system:version', 'system:id', 'system:index', 'system:bands', 'system:band_names']

Get a specific metadata property.

band_info = ndviAugust.get('system:bands')
print(band_info.getInfo())

{'NDVI': {'data_type': {'type': 'PixelType', 'precision': 'float', 'min': -1, 'max': 1}, 'crs': 'EPSG:4326', 'crs_transform': [1, 0, 0, 0, 1, 0]}}

Filter collections#

Usually, we want to limit our operations to a subset of a collection. We can do that by filtering the collection by date, location, and metadata attributes.

l8sr = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
geometry = ee.Geometry.Point([13.38, 52.52])

l8c = l8sr.filterDate('2016-07-01', '2017-07-08') \
          .filter(ee.Filter.lt('CLOUD_COVER', 70)) \
          .filterBounds(geometry)

Map through a collection#

To apply a function to every Image in an ImageCollection use the imageCollection.map() method. The only argument to map() is a function which takes one parameter: an ee.Image. For example, the following code calculates the NDVI for each image. The resulting collection holds NDVI images corresponding to the Landsat image collection.

def calcNdvi(image):
  return image.normalizedDifference(["SR_B5", "SR_B4"]).rename("ndvi")

ndviCollection = l8c.map(calcNdvi)

# adds first image from collection to map canvas
Map.addLayer(ndviCollection.first(), ndviVis, "max NDVI")

Reduce a collection#

From the collection of NDVI images, we can now calculate aggregation statistics such as the mean, max, etc. Aggregating a collection to a single image is called reducing. Reducer’s (reducing functions) are available to you in two ways: 1) via the ee.Image object and 2) via the ee.Reducer object (help).

maxNdvi = ndviCollection.reduce(ee.Reducer.max())

medianNdvi = ndviCollection.median()

Map.addLayer(maxNdvi, {"min":0, "max":1, "palette": ['white', 'green']}, "max NDVI")

Map.addLayer(medianNdvi, {"min":0, "max":1, "palette": ['white', 'green']}, "median NDVI")
Map.centerObject(geometry, 10)

Map

Projections and scale#

The Earth Engine Data Catalog contains a vast amount of satellite data and geospatial data that all come in different projections and spatial resolutions. Fortunately, you do not have to worry about harmonizing between different coordinate reference systems and pixel sizes, the Earth Engine takes care of this in the background, but you need to understand some principles. The projection and scale (pixel size) is determined by the final output rather than the input, but the actual projection is performed right in the beginning of the analysis (see here). For example, if the output is the interactive Map display all input data are automatically projected to Web Mercator. The scale of the output in the Map window depends on the zoom factor. If you zoom out to see the entire globe, the scale is going to be very large, i.e., the computations are performed with pixels that might be hundreds of kilometers in size. This way you can quickly generate visual results though at the expense of spatial precision. Conversely, the scale will decrease the further you zoom in the Map display. In batch mode, the projection and scale is set in the export function using the arguments crs and scale or crs and crsTransform. If you do not specify crs and scale the projection will default to the geographic reference system WGS84 (EPSG=4326) and a potentially large scale. For many land cover and land use applications it is recommended to specify an equal area projection. An alternative to specifying scale is the crsTransform parameter which combines information on scale and the origin of the output pixel grid.

This is useful for aligning the pixel grids of multiple outputs. The crsTransform is a list of the following format [a, b, c, d, e, f], where

a is width of a pixel in map coordinates
b is row rotation (typically zero)
c is the x-coordinate of the upper-left corner of a pixel in the pixel grid
d is the column rotation (typically zero)
e is the height of a pixel (typically negative)
f is the y-coordinate of the upper-left corner of a pixel in the pixel grid

Export data#

The Earth Engine is very good in filtering and extracting information from huge datasets, e.g., creating maps and calculating some sort of statistic, but these outputs often require further statistical analysis or visualizations. The Earth Engine can do that to some degree, but sometimes it is more convenient to export the outputs and analyze them further on a local machine.

Export to asset#

The export functions allow you to compute results in batch mode as well as save to intermediate outputs to your Earth Engine account as an asset. The code below exports the image maxNdvi the asset. Note, it also specifies the projection and pixel size of 30 x 30 m (Projections and scale) using crs and crsTransform.

task = ee.batch.Export.image.toAsset(image=maxNdvi,
                                     description='maxNdvi', assetId='maxNdvi',
                                     region=Map.get_bounds(), crs= 'EPSG:3035', 
                                     crsTransform= [30, 0, 100, 0, -30, 100], 
)

# uncomment to export
# task.start()

Export via geemap#

geemap.ee_export_image(
    maxNdvi, filename='maxNDVI.tif', scale=90, 
    region=Map.get_bounds(), file_per_band=True
)

Export to Google Drive#

task = ee.batch.Export.image.toDrive(image=maxNdvi, folder='export', 
                                     description= 'landsat',
                                     region=Map.get_bounds(), scale=30)

# task.start()