I would like to extract the location of isopycnals from a 3D xarray DataArray. The 3 dimensions are depth, latitude, and time and the data variable is density. I would like to know the coordinates (depth and latitude) of each isopycnal at each time step.

I can get something working by extracting the contours using plt.contour at each time step. However, this method seems inefficient.

Does anyone have a better solution?

Thank you, Kathy

aidanheerdegen
(Aidan Heerdegen, ACCESS-NRI Release Team Lead)
2

Hi Kathy,

I used the same method in the past to calculate isotherms:

aidanheerdegen
(Aidan Heerdegen, ACCESS-NRI Release Team Lead)
6

I sort of hoped that was what Iâ€™d provided in the initial post â€¦

The hackathon question was more along the lines of investigating and deciding on a â€śbest-practiceâ€ť example.

1 Like

aidanheerdegen
(Aidan Heerdegen, ACCESS-NRI Release Team Lead)
7

As @kathy.gunn only needs the value of the depth for a single density value Iâ€™d guess itâ€™s likely xhistogram or similar binning would do much more work than was required for this problem.

Usually binning is used to transform a diagnostic field from depth to density coordinates, so I believe wouldnâ€™t give the depth of a given density class?

At the moment, I have a code working that using for loops and plt.contour. This method is useful when I want to extract a specific isopycnal. It is slow if I want to extract every isopycnal, for example.

Iâ€™ve made a list of all the various ways I have found to do the above (or something similar), but I cannot post it as a new user. If anyone is interested, Iâ€™m happy to email it over. Iâ€™ve also pasted my code below, which I think is quite self-explanatory.

### Calculate the location of each isopycnal at every time step: N.B. the pressure and density bins must be the same size for this to work
isopycnal_layer_depths = xr.DataArray(np.zeros( (SA_section.shape) ) , dims=['isopycnal_surface','latitude','time'],coords={'isopycnal_surface': density_bins,'latitude': SA_section.latitude, 'time': SA_section.time})
isopycnal_layer_depths[:]=np.nan
# select a depth level to start at, otherwise this calculation takes ages
p_idx_to_start = (np.abs(presure_levels - 3000)).argmin()
print("This calculation starts at:", density_bins[p_idx_to_start], "kg/m^3")
for t in range(max_t_samples): # iterate over the time steps
print("progress:", np.round((t/max_t_samples)*100), "%")
for z in range(p_idx_to_start,max_z_samples): # iterate over the depth/density coordinates
g = density_bins[z] # density at this coordinate
# print("QC: time step =", t, "; density bin = ",g) # t = time step; g = density bin
cont = sigma3_section[:,:,t].plot.contour(levels=[g],yincrease=False) # this line will plot and extract contour info
for path in cont.collections[0].get_paths():
for (x,y),s in path.iter_segments(): # x = latitude and y = depth
isopycnal_layer_depths[z,:,t] = y;

Iâ€™m a little late to the party, but I think something like the following should return the depth of the deepest cell immediately above a crossing of the target value (warning untested code):

dzt is the time varying thickness of the model z levels. Because itâ€™s a z* vertical coordinate (as opposed to z), the thickness of the layers varies in time in line with the variations in sea surface height.

Hey everyone, I think this is a tricky problem to do efficiently and accurately, and probably good to add to the hackathon!

Iâ€™m also a little late to the party but have some thoughts:

Using xr.where() is fast, parallelised, and avoids the use of loops, e.g. @dougiesquire 's example or this notebook, using dzt (time-varying t-cell thickness). However, this means that the isopycnal depths will follow the grid depths (or slightly time-varying grid depths if you use dzt), when really we can probably do better. E.g. if you know the cells immediately lighter and denser than the target density, we can determine which is closer to the target density and add/remove some depth appropriately, assuming linearity. I think @Scott does this in the example that @aidanheerdegen posted, but I have a feeling the correction will add a lot of time to the computation.

An alternative is a binning strategy like @adele157 suggested. Something like this, which takes some isopycnal levels and bins the depth coordinate output dzt into isopycnal bins. You could then cumulatively sum the thicknesses of each isopycnal to get isopycnal depths. While this strategy does a similar fix of the discrete output by having a fraction where part of the cellâ€™s depth goes toward one isopycnal and part of it to the other isopycnal it is closest to, the results will be dependent on what densities you chose, so not accurate (unless you are binning everything else you are analysing into density space and treating it all as isopycnals).

Iâ€™d probably vote for the xr.where() option with a correction, but not sure how to do that efficiently!

2 Likes

aidanheerdegen
(Aidan Heerdegen, ACCESS-NRI Release Team Lead)
13

Scottâ€™s approach just does linear interpolation right? As I said, I think the approach I suggested could easily be extended to do the same. But Iâ€™m not sure whether this would be more performant than Scottâ€™s for loop. Iâ€™ll try to write a working example later today

aidanheerdegen
(Aidan Heerdegen, ACCESS-NRI Release Team Lead)
15

Yes youâ€™re correct, it does. I found something else he did that used a fit (maybe a spline), so was confused. It could be argued a more complex fit than linear interpolation is probably â€śover-fittingâ€ť.

Sorry I didnâ€™t get to this yesterday. This is the sort of thing I had in mind (again, untested code):

def isosurface(field, target, dim):
"""
Linearly interpolate a coordinate isosurface where a field
equals a target
Parameters
----------
field : xarray DataArray
The field in which to interpolate the target isosurface
target : float
The target isosurface value
dim : str
The field dimension to interpolate
Examples
--------
Calculate the depth of an isotherm with a value of 5.5:
>>> temp = xr.DataArray(
... range(10,0,-1),
... coords={"depth": range(10)}
... )
>>> isosurface(temp, 5.5, dim="depth")
<xarray.DataArray ()>
array(4.5)
"""
slice0 = {dim: slice(None, -1)}
slice1 = {dim: slice(1, None)}
field0 = field.isel(slice0).drop(dim)
field1 = field.isel(slice1).drop(dim)
crossing_mask_decr = (field0 > target) & (field1 <= target)
crossing_mask_incr = (field0 < target) & (field1 >= target)
crossing_mask = xr.where(
crossing_mask_decr | crossing_mask_incr, 1, np.nan
)
coords0 = crossing_mask * field[dim].isel(slice0).drop(dim)
coords1 = crossing_mask * field[dim].isel(slice1).drop(dim)
field0 = crossing_mask * field0
field1 = crossing_mask * field1
iso = (
coords0 + (target - field0) *
(coords1 - coords0) / (field1 - field0)
)
return iso.max(dim, skipna=True)

I think this does what Scottâ€™s function is trying to do â€“ linear interpolation of isosurface levels â€“ but itâ€™s about 10 times faster (at least for the data I was playing with) and allows for isosurface crossing from above and below. Incidentally, I think there might be an error in LINE 31 of Scottâ€™s function: df0 = target_field - field[i_lev1]
should be df0 = target_field - field[i_lev0]

More possible easy extensions to this function:

dealing with multiple isosurface crossing (i.e. donâ€™t just return the max)

returning one variable on an isosurface of a different variable (e.g. salinity on an isotherm) would be trivial

EDIT: added a docstring to the function with example usage, as per private request.

Note, if I make this change to Scottâ€™s function, then my and Scottâ€™s functions return the same result.

aidanheerdegen
(Aidan Heerdegen, ACCESS-NRI Release Team Lead)
18

Cool! FYI for all of us learning how discourse works, there is a copy symbol in the top right corner of the code snippet to allow for easy copying, which is a very nice feature.

@kathy.gunn if you think Dougieâ€™s answer solves your problem you can mark it as the solution (click on the three little dots to expand the menu under his post and click on the solution check box). That will link directly to the post in your original post to save people the hassle of searching through.