Why does calculating the mean values for each polygon take considerable time ?

[ysubash]

Hello everyone,

I've been using xarray for some time now and recently started using stackstac. I'm currently working on a problem where I need to extract the mean values of an index (such as NDVI) for each polygon (about 300 in number) and store them as individual CSV files for each polygon independently which is taking more than 24 hours.

However, the process of storing these values as individual CSV files for each polygon is taking more than 24 hours. I'm seeking assistance in optimizing this process to reduce the time required. Any suggestions would be greatly appreciated. Below are steps:

STEP1 :

I am fetching the Sentinel-2 Level 2A data from the AWS SpatioTemporal Asset Catalog (STAC) API hosted at "https://earth-search.aws.element84.com/v1" for a period of 1 year and the number of items retrieved are 135 and creating a stack.

bounding_box = array([75.80921939, 29.58215563, 75.85617824, 29.62902847])
start_date = '2023-01-01'
end_date = '2023-12-31'

sentinel_stack = stackstac.stack(
    items=s2_stac_items,
    assets = ['red','nir','scl']
    epsg= 4326,
    gdal_env=stackstac.DEFAULT_GDAL_ENV.updated(
    {'GDAL_HTTP_MAX_RETRY': 3,
    'GDAL_HTTP_RETRY_DELAY': 5,
    }),
    resampling=Resampling.nearest,
).where(lambda x: x > 0, other=np.nan).rename({'x': 'lon', 'y': 'lat'}).to_dataset(dim='band')  # sentinel-2 uses 0 as nodata

The dimensions of the sentinel_stack are 135 timestamps, with size of (11034, 10980, 3)

xarray.Dataset

Dimensions:
time: 135 lat: 11034 lon: 10980 band: 3

STEP2 : Computing ndvi using the bands.

sentinel_stack['ndvi'] = (sentinel_stack['nir'] - sentinel_stack['red'])/\
                        (sentinel_stack['nir'] + sentinel_stack['red'])
                        

STEP3: Extracting the mean ndvi for each polygon by clipping the xarray to polygon.

# Iterate over each polygon in the GeoDataFrame
for idx, row in gdf.iterrows():
    # Create a mask for the current polygon
    region = regionmask.Regions([row.geometry])
    
    bounds = row.geometry.bounds
    
    ndvi_temp = sentinel_stack_new['ndvi'].sel(
    lon=slice(bounds[0], bounds[2]),
    lat=slice(bounds[3], bounds[1]))
    
    mask = region.mask(ndvi_temp)
    masked_ndvi = ndvi_temp * mask
    
    spatial_mean = masked_ndvi.mean(dim=('lat', 'lon'))

    spatial_mean.load()

    ndvi_val=np.squeeze(spatial_mean).values

    
    # # Create a DataFrame to store the mean NDVI values along with timestamps
    df = pd.DataFrame({'Timestamp': spatial_mean.time, 'Mean_NDVI': ndvi_val})
    df.to_csv(f'polygon_{idx}_ndvi.csv', index=False)

[mathause]

You may get better answers asking the stackstac developpers. I can give you some tips but it will depend on how your data and mask actually looks:

• its generally advisable to slice (subset the array) before doing any computations (if your regions don't span the whole array (as you do in the loop; see total_bounds in geopandas))
• Aren't the mask values 0? (so for ndvi_temp * mask NaNs and 0s?)
• Instead of (ndvi_temp * mask).mean(dim=('lat', 'lon')) do ndvi_temp.groupby(mask).mean(im=('lat', 'lon')) (see also https://regionmask.readthedocs.io/en/stable/notebooks/mask_2D.html)
• You should test flox to do the groupby operations (you can potentially do them all at once, instead of looping over the regions; https://xarray.dev/blog/flox & https://flox.readthedocs.io/en/latest/)
• if your RAM allows calling load before the loop may give you the most speed up
• if not: the where is potentially called for every loop on the whole array (I am not sure it is clever enough to persist this step. Despite being counterintuitive, it might make sense do this in the loop after the sel operations)
• You can potentially avoid the where operation (sum the grid points and divide by the number of valid grid points)