Creating Maps with Google Earth Engine

Google Earth Engine (GEE) is a powerful cloud-based system for analysing massive amounts of remote sensing data. One area where Google Earth Engine shines is the ability to calculate time series of values extracted from a deep stack of imagery. While GEE is great at crunching numbers, it has limited cartographic capabilities. That’s where QGIS comes in. Using the Google Earth Engine Plugin for QGIS and Python, you can combine the computing power of GEE with the cartographic capabilities of QGIS. In this post, I will show how to write PyQGIS code to programmatically fetch time-series data, and render a map template to create an animated maps like below.

In QGIS, go to Plugins → Manage and Install Plugins and install the ‘Google Earth Engine’ plugin. After installation, you will be prompted to authenticate using your Google Earth Engine account.

I prepared a map layout using the QGIS Print Layout and saved it as a template file. Download the template ndvi_map.qpt and save it to your Desktop. Note that for each frame of the animation, we will have to extract the image data and render the template with values for date_label, ndvi_label, and season_label.

We will use the Sentinel-2 Surface Reflectance (SR) data to compute the NDVI values for the year 2019 over a farm in northern India. The concept behind extracting the time series from an image collection is nicely demonstrated in this tutorial by Nicholas Clinton. This example adapts this code for Python and adds a few enhancements to pick images where 100% of farm area is cloud-free.

Earth Engine code integrates seamlessly with PyQGIS code. You can use the Earth Engine Python API just the way you would use it elsewhere in the Python Console. Open Plugins → Python Console. Click ‘Show Editor’ button to open the built-in editor. Copy/Paste the code below and click ‘Run’ to execute the code.

import ee
from ee_plugin import Map
import os
from datetime import datetime

# Script assumes you put the ndvi_map.qpt file on your Desktop
home_dir = os.path.join(os.path.expanduser('~'))
template = os.path.join(home_dir, 'Desktop', 'ndvi_map.qpt')

geometry = ee.Geometry.Polygon([[
    [79.38757620268325, 27.45829434648896],
    [79.38834214903852, 27.459092793050313],
    [79.38789690234205, 27.459397436737895],
    [79.38718343474409, 27.458592985177017]
    ]])
farm = ee.Feature(geometry, {'name': 'farm'})
fc = ee.FeatureCollection([farm])
Map.centerObject(fc)
empty = ee.Image().byte();
outline = empty.paint(**{
  'featureCollection': fc,
  'color': 1,
  'width': 1
});
viz_params = {'bands': ['B4', 'B3', 'B2'], 'min': 0, 'max': 2000}

def maskCloudAndShadows(image):
    cloudProb = image.select('MSK_CLDPRB')
    snowProb = image.select('MSK_SNWPRB')
    cloud = cloudProb.lt(5)
    snow = snowProb.lt(5)
    scl = image.select('SCL');
    shadow = scl.eq(3) #3 = cloud shadow
    cirrus = scl.eq(10) # 10 = cirrus
    # Cloud and Snow probability less than 5% or cloud shadow classification
    mask = (cloud.And(snow)).And(cirrus.neq(1)).And(shadow.neq(1))
    return image.updateMask(mask);

def addNDVI(image):
    ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi')
    return image.addBands([ndvi])

start_date = '2019-01-01'
end_date = '2019-12-31'

collection = ee.ImageCollection('COPERNICUS/S2_SR')\
    .filterDate(start_date, end_date)\
    .map(maskCloudAndShadows)\
    .map(addNDVI)\
    .filter(ee.Filter.intersects('.geo', farm.geometry()))

def get_ndvi(image):
    stats = image.select('ndvi').reduceRegion(**{
        'geometry': farm.geometry(), 
        'reducer': ee.Reducer.mean().combine(**{
            'reducer2': ee.Reducer.count(),
            'sharedInputs': True}
            ).setOutputs(['mean', 'pixelcount']), 
        'scale': 10
        })
    ndvi = stats.get('ndvi_mean')
    pixelcount = stats.get('ndvi_pixelcount')
    return ee.Feature(None, {
        'ndvi': ndvi,
        'validpixels': pixelcount,
        'id': image.id(),
        'date': ee.Date(image.get('system:time_start')).format('YYYY-MM-dd')
      })
      

with_ndvi = collection.map(get_ndvi)
# Find how many pixels in the farm extent
max_validpixels = with_ndvi.aggregate_max('validpixels').getInfo()

def select_color(ndvi):
    ndvi_colors = {
        0.3: QColor('#dfc27d'),
        0.5: QColor('#c2e699'),
        1: QColor('#31a354')}

    for max_value, color in ndvi_colors.items():
        if ndvi < max_value:
            return color

def select_season(date_str):
    seasons = {
        'Kharif': [7, 8, 9,10],
        'Rabi': [11, 12, 1, 2, 3],
        'Summer': [4, 5, 6]
    }
    date = datetime.strptime(date_str, '%Y-%m-%d')
    month = date.month
    for season, months in seasons.items():
        if month in months:
            return season
            
features = with_ndvi.getInfo()['features']
for feature in features:
    ndvi = feature['properties']['ndvi']
    validpixels = feature['properties']['validpixels']
    date = feature['properties']['date']
    id = feature['properties']['id']
    # The following condition ensures we pick images where 
    # all pixels in the farm geometry are unmasked
    if ndvi and (validpixels == max_validpixels):
        image = ee.Image(collection.filter(
            ee.Filter.eq('system:index', id)).first())
        Map.addLayer(image, viz_params, 'Image')
        Map.addLayer(outline, {'palette': '0000FF'}, 'farm')
        project = QgsProject.instance()
        layout = QgsLayout(project)
        layout.initializeDefaults()
        with open(template) as f:
            template_content = f.read()
        doc = QDomDocument()
        doc.setContent(template_content)
        # adding to existing items
        items, ok = layout.loadFromTemplate(doc, QgsReadWriteContext(), False)
        ndvi_label = layout.itemById('ndvi_label')
        ndvi_label.setText('{:.2f}'.format(ndvi))
        ndvi_label.setFontColor(select_color(ndvi))
        
        date_label = layout.itemById('date_label')
        date_label.setText(date)
        
        season = select_season(date)
        season_label = layout.itemById('season_label')
        season_label.setText(season)

        exporter = QgsLayoutExporter(layout)
        output_image = os.path.join(home_dir, 'Desktop', '{}.png'.format(date))
        exporter.exportToImage(output_image, QgsLayoutExporter.ImageExportSettings())

The script would iterate through each image and render the template with appropriate image and data. All of the processing is done in the cloud using Google Earth Engine and only the results are streamed to QGIS. If all went well, in a few minutes, the system would have crunched through Gigabytes of data and you will have a bunch of images on your desktop. You can animate them using a program like ezgif.com

7 Comments

Leave a Comment

  1. I try this code it’s work with this coordinates when I try it in my study area it run but there’s a problem with export images in desktop
    I sent mail and attached my code
    can you help me to solve the problem
    thanks in advance

  2. Thank you so much for sharing this work! It’s really amazing.

    I have a doubt about the code: in this example a temporary monitoring is done based on a single polygon defined in line 10.
    Is it possible to iterate on the different polygons within a .shp (being this uploaded as an asset to EE?

    I think the logic in pseudocode would be to iterate over the polygon or ROI, and then apply this code to get the images and load them into Qgis.

    Then, to export each image, you would have to center the zoom and extension on it in Qgis.

    Could you please help me with this? Thanks!
    Jorge

    • You don’t even need to upload the shapefile to Earth Engine. You can open the shapefile in QGIS, select the layer and run the following code to iterate over the features.

      View formatted code

      import json
      import ee
      from ee_plugin import Map

      layer = iface.activeLayer()
      for feature in layer.getFeatures():
      geom = feature.geometry()
      coordinates = json.loads(geom.asJson())[‘coordinates’][0]
      geometry = ee.Geometry.Polygon(coordinates)
      farm = ee.Feature(geometry, {‘name’: ‘farm’})
      fc = ee.FeatureCollection([farm])
      ……

  3. Hi Ujaval. I tried to run your code in qgis 3.10LTR. It worked but nvdi outputs (PGN files) were blank in my desktop, only your template but no NDVI image :(. There is no error in running the code, is it possible?

  4. Thanks Ujaval for sharing such application. As you described the Plugin would ask for authentication but the the version 0.0.2 not asking such authentication after installation rather showing a error message such as ‘Couldn’t load plugin ‘ee_plugin’ due to an error when calling its classFactory() method…’. Do you have any idea how to resolve it?

  5. Same problem as Nestor and faten above, the program works and NDVI images are properly mappes in canvas but the output is blank. Any indication to output the images is highly welcome! Thank you in advance!

Leave a Reply