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
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.
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 @Aidan 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!
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:
max)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.