Extracting isopycnals from multi-dimensional arrays

Technical
16

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.