Image workflows

Image workflows#

Packages#

import ee
import geemap

try:
    ee.Initialize(project="ee-pflugmacher")
except Exception as e:
    ee.Authenticate()
    ee.Initialize(project="ee-pflugmacher")

Set up map window#

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

Image properties#

Images contain several types of metadata information such as band names, projection information, properties, and other metadata. To understand what methods are available for an ee.Image object, you may consult the API documentation (see Client libraries).

We can retrieve the band names of the ee.Image using the bandNames method.

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.

l8_img.bandNames()

List (19 elements)
- 0:SR_B1
- 1:SR_B2
- 2:SR_B3
- 3:SR_B4
- 4:SR_B5
- 5:SR_B6
- 6:SR_B7
- 7:SR_QA_AEROSOL
- 8:ST_B10
- 9:ST_ATRAN
- 10:ST_CDIST
- 11:ST_DRAD
- 12:ST_EMIS
- 13:ST_EMSD
- 14:ST_QA
- 15:ST_TRAD
- 16:ST_URAD
- 17:QA_PIXEL
- 18:QA_RADSAT

Get a list of all metadata properties.

properties = l8_img.propertyNames()
properties

List (98 elements)
- 0:DATA_SOURCE_ELEVATION
- 1:WRS_TYPE
- 2:system:id
- 3:REFLECTANCE_ADD_BAND_1
- 4:REFLECTANCE_ADD_BAND_2
- 5:DATUM
- 6:REFLECTANCE_ADD_BAND_3
- 7:REFLECTANCE_ADD_BAND_4
- 8:REFLECTANCE_ADD_BAND_5
- 9:REFLECTANCE_ADD_BAND_6
- 10:REFLECTANCE_ADD_BAND_7
- 11:system:footprint
- 12:REFLECTIVE_SAMPLES
- 13:system:version
- 14:GROUND_CONTROL_POINTS_VERSION
- 15:SUN_AZIMUTH
- 16:DATA_SOURCE_TIRS_STRAY_LIGHT_CORRECTION
- 17:UTM_ZONE
- 18:DATE_ACQUIRED
- 19:ELLIPSOID
- 20:system:time_end
- 21:DATA_SOURCE_PRESSURE
- 22:LANDSAT_PRODUCT_ID
- 23:STATION_ID
- 24:TEMPERATURE_ADD_BAND_ST_B10
- 25:DATA_SOURCE_REANALYSIS
- 26:REFLECTANCE_MULT_BAND_7
- 27:system:time_start
- 28:REFLECTANCE_MULT_BAND_6
- 29:L1_PROCESSING_LEVEL
- 30:PROCESSING_SOFTWARE_VERSION
- 31:L1_DATE_PRODUCT_GENERATED
- 32:ORIENTATION
- 33:REFLECTANCE_MULT_BAND_1
- 34:WRS_ROW
- 35:REFLECTANCE_MULT_BAND_3
- 36:REFLECTANCE_MULT_BAND_2
- 37:TARGET_WRS_ROW
- 38:REFLECTANCE_MULT_BAND_5
- 39:REFLECTANCE_MULT_BAND_4
- 40:THERMAL_LINES
- 41:TIRS_SSM_POSITION_STATUS
- 42:GRID_CELL_SIZE_THERMAL
- 43:IMAGE_QUALITY_TIRS
- 44:TRUNCATION_OLI
- 45:NADIR_OFFNADIR
- 46:CLOUD_COVER
- 47:REQUEST_ID
- 48:EARTH_SUN_DISTANCE
- 49:GEOMETRIC_RMSE_VERIFY
- 50:TIRS_SSM_MODEL
- 51:COLLECTION_CATEGORY
- 52:SCENE_CENTER_TIME
- 53:GRID_CELL_SIZE_REFLECTIVE
- 54:SUN_ELEVATION
- 55:ALGORITHM_SOURCE_SURFACE_TEMPERATURE
- 56:TEMPERATURE_MAXIMUM_BAND_ST_B10
- 57:CLOUD_COVER_LAND
- 58:GEOMETRIC_RMSE_MODEL
- 59:ROLL_ANGLE
- 60:COLLECTION_NUMBER
- 61:DATE_PRODUCT_GENERATED
- 62:L1_REQUEST_ID
- 63:DATA_SOURCE_OZONE
- 64:SATURATION_BAND_1
- 65:DATA_SOURCE_WATER_VAPOR
- 66:SATURATION_BAND_2
- 67:SATURATION_BAND_3
- 68:IMAGE_QUALITY_OLI
- 69:SATURATION_BAND_4
- 70:LANDSAT_SCENE_ID
- 71:SATURATION_BAND_5
- 72:MAP_PROJECTION
- 73:SATURATION_BAND_6
- 74:SENSOR_ID
- 75:SATURATION_BAND_7
- 76:SATURATION_BAND_8
- 77:WRS_PATH
- 78:SATURATION_BAND_9
- 79:TARGET_WRS_PATH
- 80:L1_PROCESSING_SOFTWARE_VERSION
- 81:TEMPERATURE_MULT_BAND_ST_B10
- 82:L1_LANDSAT_PRODUCT_ID
- 83:PROCESSING_LEVEL
- 84:ALGORITHM_SOURCE_SURFACE_REFLECTANCE
- 85:GROUND_CONTROL_POINTS_MODEL
- 86:SPACECRAFT_ID
- 87:TEMPERATURE_MINIMUM_BAND_ST_B10
- 88:GEOMETRIC_RMSE_MODEL_Y
- 89:REFLECTIVE_LINES
- 90:GEOMETRIC_RMSE_MODEL_X
- 91:THERMAL_SAMPLES
- 92:system:asset_size
- 93:DATA_SOURCE_AIR_TEMPERATURE
- 94:GROUND_CONTROL_POINTS_VERIFY
- 95:system:index
- 96:system:bands
- 97:system:band_names

Get a specific metadata property:

band_info = l8_img.get('CLOUD_COVER')
band_info

0.04

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 calculate the NDVI from the Landsat 8 image. But first, we need to rescale the digital numbers to reflectance values using the gains and offset described in the data catalog.

l8_sr = l8_img.multiply(0.0000275).add(-0.2)

We can select individual bands:

nir = l8_sr.select("SR_B5")
red = l8_sr.select("SR_B4")

… and calculate the NDVI:

ndvi = nir.subtract(red) \
              .divide(nir.add(red)) \
              .rename('NDVI')
              

ndviVis = {"min":-1, "max":1, "palette": ['red', 'white', 'green']}
Map.addLayer(ndvi, ndviVis, 'NDVI')
Map

Logical operations#

We can create a binary mask based on the NDVI. For example, we might want to separate areas with low NDVI (bare or sparse vegetation - but not exact) from vegetation areas.

veg_mask = ndvi.gt(0.4)
Map.addLayer(veg_mask, {}, 'High NDVI')

Masking#

We can apply the mask using ee.Image.mask() or ee.Image.updateMask(). For example, we can do the following to mask out bare areas:

veg_image = ndvi.updateMask(veg_mask)
Map.addLayer(veg_image, ndviVis, 'NDVI vegetation')

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

Batch environment#

Recall, Earth Engine data types are server-side objects that are created and processed entirely on cloud infrastructure, and that there are two different processing environments for you to request information: 1) an interactive environment optimized for speedy and small requests, and 2) a batch environment optimized for high-latency parallel processing of large amounts of data.

Batch processing is triggered when you use an import or export function (e.g., Export.* or ee.batch.). The batch environment let’s you export large datasets to your Assets Folder or Google Cloud Storage. You can manage (view and cancel) tasks using the task manager webpage, the task tab on the online code editor, and programmatically in Python. Look here for more information.

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.

bounds = ee.Geometry.Polygon(
        [[[12.232347553388268, 53.243183673553276],
          [12.232347553388268, 52.0636255084833],
          [14.737230365888268, 52.0636255084833],
          [14.737230365888268, 53.243183673553276]]])

task = ee.batch.Export.image.toAsset(image=maxNdvi,
                                     description='maxNdvi', assetId='projects/ee-pflugmacher/assets/maxNdvi',
                                     region=bounds, crs= 'EPSG:3035', 
                                     crsTransform= [30, 0, 100, 0, -30, 100]
)

# start task
task.start()

Export via geemap#

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

Export to Google Drive#

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

# task.start()

Task management in Python#

tasks = ee.batch.Task.list()
for task in tasks:
    print(f'Task ID: {task.id}, State: {task.state}, Description: {task.config.get("description", "No description")}')

tasks = ee.batch.Task.list()
for task in tasks:
    print(f'Task ID: {task.id}, State: {task.state}, Description: {task.config.get("description", "No description")}')
    
for task in tasks:
    print(f'Task ID: {task.id}, State: {task.state}, Description: {task.config.get("description", "No description")}')
    if task.state == 'FAILED':
        print(f'Error: {task.error_message}')