Find mean value in an area using emsarray

area_average.py

import emsarray
import pandas
import shapely

# Open the dataset
gbr4 = emsarray.open_dataset('https://dapds00.nci.org.au/thredds/dodsC/fx3/model_data/gbr4_bgc_GBR4_H2p0_B2p0_Chyd_Dcrt.ncml')

# Define the area of interest
point = shapely.Point(152.384000, -22.672241)
circle = shapely.buffer(point, 0.1)


# Coming soon: I'll add this as a function to dataset.ems
def selector_for_indices(dataset, indices):
    convention = dataset.ems
    selectors = map(convention.selector_for_index, indices)
    return pandas.DataFrame(selectors).to_xarray().drop_vars('index')


# Find all cells that intersect
# Coming soon: easier spatial queries
#   intersecting_cells = [
#       gbr4.ems.ravel_index(index)
#       for index in gbr4.ems.strtree.query(circle, predicate='intersects')
#   ]
intersecting_cells = [
    hit.index
    for polygon, hit in gbr4.ems.spatial_index.query(circle)
    if polygon.intersects(circle)
]


# Pick a variable. Trying to do this to the whole dataset will run out of memory.
# Coming soon: Convention.select_variables
#   ds = ds.ems.select_varibles(['eta', 'salt', 'temp'])
# https://emsarray.readthedocs.io/en/latest/api/conventions/interface/#emsarray.conventions.Convention.select_variables
temp = gbr4['temp']

# Select just the cells we want
temp = temp.isel(selector_for_indices(gbr4, intersecting_cells))

# Select the surface layer and the last 30 time steps
temp = temp.isel(k=-1, time=slice(-30, None))
print("Temperature in the area:", temp)
print("=" * 30)

# Compute the mean surface temperature in the area for each time step
mean_temp = temp.mean(dim='index')

print("Mean surface temperature:", mean_temp)

@mx-moth @frizwi I thought I would share my final code for working this problem out. This loop samples the area under a moving triangle then calculates the average of the sampled parameter, before looping through time.

There are a couple of bugs like why I have to do vor = ds but this works well for my purpose.

ds = ems.open_dataset('ereefdata_may_23_vor.nc')
vor = ds

t = 43
moving_tri_average = []

for i in range(22,288,4):
    
    triangle = shapely.Polygon([(Alon[i], Alat[i]),(Glon[i], Glat[i]),(Jlon[i], Jlat[i])]) #L
    
    def selector_for_indices(dataset, indices):
        convention = dataset.ems
        selectors = map(convention.selector_for_index, indices)
        return pandas.DataFrame(selectors).to_xarray().drop_vars('index')

    intersecting_cells = [
        hit.index
        for polygon, hit in ds.ems.spatial_index.query(triangle)
        if polygon.intersects(triangle)]
    
    tri_vor = vor.isel(time=t).isel(selector_for_indices(ds, intersecting_cells))
    tri_vor1 = tri_vor.to_array().to_numpy()
    tri_vor1 = tri_vor1.mean()
    moving_tri_average.append(tri_vor1)
    t=t+1