Search code examples
pythongoogle-earth-engine

Get results in an Earth Engine python script


I'm trying to get NDVI mean in every polygon in a feature collection with earth engine python API. I think that I succeeded getting the result (a feature collection in a feature collection), but then I don't know how to get data from it. The data I want is IDs from features and ndvi mean in each feature.

import datetime
import ee
ee.Initialize()

#Feature collection
fc = ee.FeatureCollection("ft:1s57dkY_Sg_E_COTe3sy1tIR_U-5Gw-BQNwHh4Xel");
fc_filtered = fc.filter(ee.Filter.equals('NUM_DECS', 1))
#Image collection
Sentinel_collection1 = (ee.ImageCollection('COPERNICUS/S2')).filterBounds(fc_filtered)
Sentinel_collection2 = Sentinel_collection1.filterDate(datetime.datetime(2017, 1, 1),datetime.datetime(2017, 8, 1))


# NDVI function to use with ee map
def NDVIcalc (image):
  red = image.select('B4')
  nir = image.select('B8')
  ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI')

  #NDVI mean calculation with reduceRegions
  MeansFeatures = ndvi.reduceRegions(reducer= ee.Reducer.mean(),collection= fc_filtered,scale= 10)

  return (MeansFeatures)

#Result that I don't know to get the information: Features ID and NDVI mean
result = Sentinel_collection2.map(NDVIcalc)

Solution

  • If the result is small, you pull them into python using result.getInfo(). That will give you a python dictionary containing a list of FeatureCollection (which are more dictionaries). However, if the results are large or the polygons cover large regions, you'll have to Export the collection instead.

    That said, there are probably some other things you'll want to do first:

    1) You might want to flatten() the collection, so it's not nested collections. It'll be easier to handle that way.

    2) You might want to add a date to each result so you know what time the result came from. You can do that with a map on the result, inside your NDVIcalc function

    return MeansFeatures.map(lambda f : f.set('date', image.date().format())
    

    3) If what you really want is a time-series of NDVI over time for each polygon (most common), then restructuring your code to map over polygons first will be easier:

    Sentinel_collection = (ee.ImageCollection('COPERNICUS/S2')
        .filterBounds(fc_filtered)
        .filterDate(ee.Date('2017-01-01'),ee.Date('2017-08-01')))
    
    def GetSeries(feature):
      def NDVIcalc(img):
        red = img.select('B4')
        nir = img.select('B8')
        ndvi = nir.subtract(red).divide(nir.add(red)).rename(['NDVI'])
        return (feature
                .set(ndvi.reduceRegion(ee.Reducer.mean(), feature.geometry(), 10))
                .set('date', img.date().format("YYYYMMdd")))
    
      series = Sentinel_collection.map(NDVIcalc)
      // Get the time-series of values as two lists.
      list = series.reduceColumns(ee.Reducer.toList(2), ['date', 'NDVI']).get('list')
      return feature.set(ee.Dictionary(ee.List(list).flatten()))
    
    result = fc_filtered.map(GetSeries)
    print(result.getInfo())
    

    4) And finally, if you're going to try to Export the result, you're likely to run into an issue where the columns of the exported table are selected from whatever columns the first feature has, so it's good to provide a "header" feature that has all columns (times), that you can merge() with the result as the first feature:

    # Get all possible dates.
    dates = ee.List(Sentinel_collection.map(function(img) {
          return ee.Feature(null, {'date': img.date().format("YYYYMMdd") })
    }).aggregate_array('date'))
    
    # Make a default value for every date.
    header = ee.Feature(null, ee.Dictionary(dates, ee.List.repeat(-1, dates.size())))
    output = header.merge(result)
    ee.batch.Export.table.toDrive(...)