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:


Was that not fast enough for your use?

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


@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 :slight_smile:

@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 … :wink:

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})
# 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!


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

    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
    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 ()>
    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.

Screen Shot 2022-11-25 at 2.08.19 pm

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.