Extract percentiles for each pixel in a dayofyear range and export it to xarray

[iuryt]

I am trying to export the percentile fields for a specific time range of all years from MERRA2 dataset.
The code I tried to run is here below, but this part of concatenation and saving just take so long that I didn't wait for it.

ds = xr.concat(ds, "percentile")
ds.to_netcdf("output.nc")

The code

# Load MERRA/slv/2 dataset
dataset = ee.ImageCollection('NASA/GSFC/MERRA/slv/2')

# Define the desired date range (Jan 1st to Jan 15th)
start_day = 1  
end_day = 15

# Filter by day-of-year range (using server-side filtering)
filtered_dataset = dataset.filter(ee.Filter.dayOfYear(start_day, end_day))


lonlim = [-(88+39/60), -(30+40/60)]
latlim = [-(41+23/60), 12+58/60]
bbox = ee.Geometry.BBox(lonlim[0], latlim[0], lonlim[1], latlim[1])
variables = ["PS", "T10M", "TQV", "U10M", "V10M"]

ds = []
for percentile in tqdm(np.arange(10, 100, 10)):
    percentile_image = filtered_dataset.select(variables).reduce(ee.Reducer.percentile([int(percentile)]))
    
    dsi = geemap.ee_to_xarray(
        percentile_image.clip(bbox)
    ).squeeze().transpose("lat","lon").sel(lon=slice(*lonlim), lat=slice(*latlim))
    
    dsi = dsi.rename({var: var.split("_")[0] for var in list(dsi)})
    dsi = dsi.assign_coords(percentile=percentile)
    ds.append(dsi)
ds = xr.concat(ds, "percentile")
ds.to_netcdf("output.nc")

The idea is to have a final data that each variable varies in (lat, lon, percentile) and the percentiles are computed for each pixel from all values in the specified time range of the years.

Anyone has any idea on how to perform this efficiently?

[giswqs]

Use the map function on the image collection. Never use a for loop on an image collection. It is computationally expensive.
https://developers.google.com/earth-engine/guides/ic_mapping

[iuryt]

Thanks for the quick response!

The problem is that even for a single image it takes too long. I waited for half-hour and just stopped the process.
Is this expected to be like this or am I doing something wrong?

# Load MERRA/slv/2 dataset
dataset = ee.ImageCollection('NASA/GSFC/MERRA/slv/2')

# Define the desired date range (Jan 1st to Jan 15th)
start_day = 1  
end_day = 15

# Filter by day-of-year range (using server-side filtering)
filtered_dataset = dataset.filter(ee.Filter.dayOfYear(start_day, end_day))


lonlim = [-(88+39/60), -(30+40/60)]
latlim = [-(41+23/60), 12+58/60]
bbox = ee.Geometry.BBox(lonlim[0], latlim[0], lonlim[1], latlim[1])
variables = ["PS", "T10M", "TQV", "U10M", "V10M"]

percentile = 90
percentile_image = filtered_dataset.select(variables).reduce(ee.Reducer.percentile([int(percentile)]))
    
dsi = geemap.ee_to_xarray(
        percentile_image.clip(bbox)
    ).squeeze().transpose("lat","lon").sel(lon=slice(*lonlim), lat=slice(*latlim))
    
dsi = dsi.rename({var: var.split("_")[0] for var in list(dsi)})
dsi = dsi.assign_coords(percentile=percentile)
dsi.to_netcdf("output.nc")

[giswqs]

ee_to_xarray is best suited for original EE dataset, not computed object. For long running process, you should export data to Google Drive

[iuryt]

Ok. So I export to Google Drive and download from there if I want to have it locally? I will test and post it here if I find a workable example.

[iuryt]

I noticed that it adds _p90 to the name of the variables instead of creating a new dimension. Is there a way to do all percentiles at once by creating a new dimension instead?

[giswqs]

Export the images to Google Drive, then use Xarray to deal with dimensions.

filtered_dataset = dataset.filter(ee.Filter.dayOfYear(start_day, end_day)).select(variables)
precentiles = ee.List.sequence(10, 100, 10)
collection = precentiles.map(lambda p: filtered_dataset.reduce(ee.Reducer.percentile([p])).clip(bbox))
image = ee.ImageCollection(collection).toBands()
image