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)

import emsarray
import pandas
import shapely
import numpy as np

ds = emsarray.open_dataset('https://dapds00.nci.org.au/thredds/dodsC/fx3/gbr1_2.0/gbr1_simple_2023-06-01.nc')

#point = shapely.Point(152.384000, -22.672241)
#circle = shapely.buffer(point, 0.1)

triangle = shapely.Polygon([(147.783257, -19.414204),
                    (147.836966, -19.344435),
                    (147.739891, -19.339585)])
# 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')

intersecting_cells = [
    hit.index
    for polygon, hit in ds.ems.spatial_index.query(triangle)
    if polygon.intersects(triangle)]

#u
u = ds['u']
# Select just the cells we want
u = u.isel(selector_for_indices(ds, intersecting_cells))
# Select the surface layer
u = u.isel(k=-1)
# Compute the mean u velocity in the area for each time step
du = u.diff('u')
du = du.rename('du')
mean_du = du.mean(dim='index')

#u
v = ds['v']
# Select just the cells we want
v = v.isel(selector_for_indices(ds, intersecting_cells))
# Select the surface layer
v = v.isel(k=-1)
# Compute the mean u velocity in the area for each time step
dv = v.diff('v')
dv = dv.rename('dv')
mean_dv = dv.mean(dim='index')

# Compute the mean dx in the area
x = (u.longitude).to_numpy()
dx = (np.diff(x))*111139 #converts degrees to meters
mean_dx = np.median(dx) #median chosen 

# Compute the mean dy in the area
y = (v.latitude).to_numpy()
dy = (np.diff(y))*111139 #converts degrees to meters
mean_dy = np.median(dy) #median chosen 

#Compute the du/dx
du_dx = mean_du/mean_dx
du_dx = du_dx.rename('du_dx')

#Compute the dv/dy
dv_dy = mean_dv/1500
dv_dy = dv_dy.rename('dv_dy')

#Issue I'm an orders of magnitude out as it is expected the dy is approx dx in magnitude. But maybe I'm right? 

#compute Divergence
div = du_dx + dv_dy

mean_dx
mean_dy
div.plot()

Bad formating, but that is my example

If you want to include large blocks of code in Markdown use triple backticks to surround the code. I've edited your comment to correct the formatting.

```python
triangle = shapely.Polygon([
    (147.783257, -19.414204),
    (147.836966, -19.344435),
    (147.739891, -19.339585),
])
```

Hey Matt,

I think the issue might be with the diff method on the u and v variables, I haven't run the code to verify but suspect that after the selection is done, it's a 1D array, and hence the differences are not necessarily of the adjacent cells that you are probably looking for. It's possible that they are in "row" major order but a) there's no guarantee that the i/j dimensions line up east/west and b) in any case, there will be wrap arounds

Tim suggests to do the differences before the selection is made as that will operate on the original 2D array, but you'll still need to be careful wrt to the i/j indices. First you need to establish which one goes east/west and n/s and then pass that index to the diff operator. After you've got the du and dv sorted out, you can then do the selection and go from there. Also note that the degree to metres factor for lat/lon is only the same at the equator i.e. for latitude it is longitude dependent as the world narrows at the poles.

Finally, it looks like you're trying to calculate the divergence so you might also want to think about the type of numerical difference (backward/forward euler etc) that is the most appropriate.

Hope this helps!

Cheers,
-Farhan

@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