As an example, I made 3 field polygons (labeled A,B,C) for which I want (field-averaged) time series of spectral index statistics (e.g. NDVI mean and stdDev) based on a (GEE) Sentinel-2 image collection (6 images in this example). I want to save these time series in a single Pandas dataframe, with data from the first field polygon (labeled 'A' in the field-name column) at the top, followed by the second field polygon (labeled 'B' in the field-name column), and with the data from the last field polygon at the bottom end of the dataframe (labeled 'C' in the field-name column). The final Pandas dataframe should have 3*6=18 rows and the following columns:

  • 'Field_name': with 'A' in the first 6 rows, followed by 6 rows with 'B', and then 6 rows with 'C'.
  • 'Date': 3 times the series with 6 image dates below one another.
  • 'NDVI_mean': first the 6 field-averaged NDVI values of field A, then for field B, then for field C.
  • 'NDVI_stdDev': first the 6 NDVI standard deviations of field A, then of field B, and then of field C.

The for-loop construction below sort of gives the dataframe I want (except from the column headers that are repeated). However, I guess it must be possible to use a 'nested function construction' instead (so mapping over both the field polygons and images), preferable on the server-side of Google Earth Engine. Can someone help me out? Thanking in advance.

# Libraries

import numpy as np
import pandas as pd
import ee
import geemap
ee.Initialize()

# Satellite image collection

start_date = '2018-06-14'
end_date = '2018-06-30'   

aoi = ee.Geometry.Point(5.995, 50.761)

S2 = (ee.ImageCollection('COPERNICUS/S2_SR')
    .filterBounds(aoi)
    .filterDate(start_date, end_date))

# Add NDVI as extra band to the images

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

S2 = S2.map(addNDVI)

# Field polygons - feature collection and dataframe

A = ee.Feature(ee.Geometry.Rectangle([[5.993, 50.763],[5.997, 50.766]]))
B = ee.Feature(ee.Geometry.Rectangle([[5.990, 50.769],[5.994, 50.772]]))
C = ee.Feature(ee.Geometry.Rectangle([[5.998, 50.767],[6.000, 50.770]]))

fields_fc = ee.FeatureCollection([A,B,C])

df_fields = geemap.ee_to_geopandas(fields_fc)
df_fields["field_id"] = [0,1,2]
df_fields["field_name"] = ['A','B','C']

# Visualise situation

Map = geemap.Map()
Map.centerObject(aoi,13)
Map.addLayer(S2.first(), {'min': 0, 'max': 4000,'bands': ['B4', 'B3', 'B2']}, "example image")
Map.addLayer(A,{'color':'red','strokeWidth': 5},"field A")
Map.addLayer(B,{'color':'green','strokeWidth': 5},"field B")
Map.addLayer(C,{'color':'blue','strokeWidth': 5},"field C")
Map

# For each field polygon, I want to calculate the mean and stdDev of the NDVI index

def Index_field_statistics(img):
    ndvi_mean = img.reduceRegion(reducer=ee.Reducer.mean(), geometry=fieldgeomx, scale=10).get('NDVI')
    ndvi_stdDev = img.reduceRegion(reducer=ee.Reducer.stdDev(), geometry=fieldgeomx, scale=10).get('NDVI')
    return img.set('field_name',namex).set('date',img.date().format('YYYY-MM-dd')).set('ndvi_mean',ndvi_mean).set('ndvi_stdDev',ndvi_stdDev)

# So here I want to loop over the field polygons (field_id, field_name) and extract the spectral statistics for each image 

field_id = df_fields['field_id'].values

field_name = df_fields['field_name'].values

FINAL_DATAFRAME = []

for x in range(0,len(field_id),1):

    idx = field_id[x]
    namex = field_name[x]

    fieldx = fields_fc.filterMetadata("system:index","equals",str(idx)).first()
    fieldgeomx = fieldx.geometry()

    S2_stats = S2.map(Index_field_statistics)
    S2_stats_list = S2_stats.reduceColumns(ee.Reducer.toList(4),['field_name','date','ndvi_mean','ndvi_stdDev']).values().get(0)

    df_stats = pd.DataFrame(S2_stats_list.getInfo(), columns=['field_name','date','ndvi_mean','ndvi_stdDev'])

    FINAL_DATAFRAME.append(df_stats)

FINAL_DATAFRAME
0

There are 0 best solutions below