-
-
Save mx-moth/4b0b2e2c99e8a35968446e1d6adf4ffb to your computer and use it in GitHub Desktop.
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) |
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
Bad formating, but that is my example