# Extracting isopycnals from multi-dimensional arrays

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

Hi Kathy,

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

https://gist.github.com/aidanheerdegen/11451b280d7e07a8f4c078820908ab4a

Was that not fast enough for your use?

@Scott used a different (arguably better) method here:

https://gist.github.com/ScottWales/8c2b37e4919a02d5930a89780cb9b714

@navidcy @rmholmes @adele157 Is this a potential hackathon candidate?

Could be! But I also wish that @kathy.gunn get help with this before the hackathon

@claireyungâ€™s Binning transformation example does this. You can bin `dzt` using these methods to get isopycnal depths.

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

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?

Adele, what is dzt? I assume it is a variable from the model, but I do not know what it means.

Hello everyone. Thank you for your suggestions.

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):

``````def isosurface(field, target, dim):
field1 = field.shift({dim: -1})

crossing_mask_decr = (field > target) & (field1 <= target)
crossing_mask_incr = (field < target) & (field1 >= target)
crossing_mask = xr.where(
crossing_mask_decr | crossing_mask_incr, 1, np.nan
)

coords_masked = crossing_mask * field[dim]

return coords_masked.max(dim, skipna=True)
``````

This approach could be pretty easily extended to interpolate the isosurface levels.

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

This is similar to the Cosima-recipe for pulling out the age at the bottom of the ocean, but isnâ€™t as precise as a fit like Scottâ€™s code linked above.

It is an embarrassingly parallel problem, so if it isnâ€™t fast enough the best approach would be to run it in parallel with dask.

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

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.

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.

Claire - thank you for the detailed reply. I am going to try each of the different approaches (if I can get them working).

Cheers, Kathy

Dougie - Thank you! My plan is to try this code this week, Iâ€™ll let you know how it goes.