Slicing/isosurfaces methods for 3D grids.

[lesommer]

Eventually, 3D grids should come with methods for projecting 3D fields onto 2D isosurfaces (ex : temperature on an isopycnal).

This issue gathers informations / ideas as to how to proceed.

[lesommer]

an example of code from @ocefpaf , including a discussion on how to accelerate the code with numba.

https://ocefpaf.github.io/python4oceanographers/blog/2015/10/05/isosurfaces/

[ocefpaf]

BTW that code is now a lightweight package to calculate isosurfaces. Sugestions and new use cases are welcomed! https://github.com/ioos/ciso

…

On Dec 2, 2016 12:19 PM, "Julien Le Sommer" ***@***.***> wrote: an example of code from @ocefpaf <https://github.com/ocefpaf> , including a discussion on how to accelerate the code with numba. https://ocefpaf.github.io/python4oceanographers/blog/ 2015/10/05/isosurfaces/ — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#42 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AA6BL6Te9inJWtgOIIVlXFHVEeVinK2nks5rECj9gaJpZM4LCmKv> .

[rabernat]

ciso looks very nice and simple. (would be even nicer if it could work directly with xarray data! ;)

I have implemented a somewhat different approach in xgcm:
https://github.com/xgcm/xgcm/blob/master/xgcm/regridding.py
(this code has not been updated in a while...)

The approach uses bins and is conservative--it transforms from height coordinates to isopycnal coordinates, preserving the total water column thickness (or other integral properties). To calculate an isopycnal projection of a tracer, you would do something like

thickness = regrid_vertical(dz, rho, rho_levs, dim='z')
tracer_thickness_weighted = regrid_vertical(dz * tracer, rho, rho_levs, dim='z')
tracer_rho_mean = tracer_thickness_weighted / thickness

# check that vertically integrated tracer is conserved
assert (tracer * dz).sum(dim='z') == (tracer_rho_mean * thickness).sum(dim='rho')

Because it uses binning, this approach also does not require the isosurface field to be monotonic in the vertical.

This can be useful in certain cases.

[hetland]

I just created an isosurface calculator in xarray using a similar algorithm to ciso. I stole the idea from my octant package. This method also requires monotonic data on which to define the isosurface. However, there are lots of use cases where this is true, like constant depth (not trivial for stretched coordinate ocean models) and since the algorithm is coordinate agnostic, you can get values along a constant latitude or longitude as well (again, not trivial for curvilinear grids).

def xisoslice(ds, coord, iso_array, iso_value, projected_array):
    Nm = len(ds[coord]) - 1

    lslice = {coord: slice(None, -1)}
    uslice = {coord: slice(1, None)}

    prop = ds[iso_array] - iso_value

    propl = prop.isel(**lslice)
    propl.coords[coord] = np.arange(Nm)
    propu = prop.isel(**uslice)
    propu.coords[coord] = np.arange(Nm)

    zc = xr.where((propu*propl)<0.0, 1.0, 0.0)

    varl = ds[projected_array].isel(**lslice)
    varl.coords[coord] = np.arange(Nm)
    varu = ds[projected_array].isel(**uslice)
    varu.coords[coord] = np.arange(Nm)

    propl = (propl*zc).sum(coord)
    propu = (propu*zc).sum(coord)
    varl = (varl*zc).sum(coord)
    varu = (varu*zc).sum(coord)

    return varu - propl*(varu-varl)/(propu-propl)

Use cases:

Getting salinity at a specified depth:

result = xisoslice(ds, 's_rho', 'z_rho', -20, 'salt')

Getting the depth of a salinity isosurface:

result = xisoslice(ds, 's_rho', 'salt', 28.0, 'z_rho')

Getting a cross-section along a constant longitude:

salt_slice = xisoslice(ds, 'xi_rho', 'lon_rho', -92.0, 'salt')
salt_slice.coords['z_rho'] = xisoslice(ds, 'xi_rho', 'lon_rho', -92.0, 'z_rho')
salt_slice.coords['lat_rho'] = xisoslice(ds, 'xi_rho', 'lon_rho', -92.0, 'lat_rho')

In this last case, the coordinates also need to be sliced in a similar way, hence the last two lines.

I think it would be nice if this code could live somewhere easy to find and use. I don't think xarray is the right place. Perhaps xgcm?

[ejhyer]

@hetland this was a great find. I would certainly support seeing it integrated into xarray. I would point out that it requires monotonic values in the iso_array, as you noted, and actually as written in your example requires monotonically ascending values.