Remapping from depth to density space and back again (xgcm, xarray)

[kathy.gunn]

Hello everyone. I need to re-map some data (temperature and salinity) from pressure-latitude space into density-latitude space. In density-latitude space, I do some calculations. Then, I need to re-map the results back into pressure-latitude space, so I can do some final calculations in pressure-latitude space.

I can do the first step, which is re-mapping temperature and salinity from pressure-latitude space into density-latitude space. Below is my working code

from xgcm import Grid
# Input Data must be 3D with depth, x-coord (lat), time.
def recast_property_from_depth_to_gamma_space(max_x_samples,max_z_samples,max_t_samples, density_bins, gridded_density, gridded_temperature, gridded_salinity):
    if gridded_density.shape[0] != max_z_samples:
        raise ValueError('Gridded fields must have max_z_samples rows')
    elif gridded_density.shape[1] != max_x_samples:
        raise ValueError('Gridded fields must have max_x_samples cols')
    elif gridded_density.shape[2] != max_t_samples:
        raise ValueError('Gridded fields must have max_t_samples in 3rd dim')
    elif gridded_density.ndim != 3:
        raise ValueError('Gridded fields must have 3 dimensions (pressure, latitude, time)')
      
    # xgcm has to work on an xarray.Dataset, so first step is to convert and combine the xarray dataArrays into data sets
    ds = gridded_salinity.to_dataset().rename({'salt':'SA'})
    # add tempearture to the data set as a variable. 
    ds = ds.assign({'CT': gridded_temperature})
    # xgcm also needs to know how density varies alongside temperature and salinity, so add that to the data set as a data variable. 
    ds = ds.assign({'density': gridded_density})
    # Now, do the conversion
    # First create an xgcm grid object
    grid = Grid(ds, coords={'Z': {'center':'pressure'}}, periodic=False)
    # and transform into density coordinates
    gridded_salinity_gamma_space = grid.transform(ds.SA, 'Z', density_bins, target_data=ds.density)
    gridded_temperature_gamma_space = grid.transform(ds.CT, 'Z', density_bins, target_data=ds.density)
    # transpose that the the order of density, latitude, time is re-instated
    gridded_salinity_gamma_space = gridded_salinity_gamma_space.transpose('density','latitude','time')
    gridded_temperature_gamma_space = gridded_temperature_gamma_space.transpose('density','latitude','time')    
    return gridded_temperature_gamma_space, gridded_salinity_gamma_space

# define the target values in density, linearly spaced
[gridded_temperature_gamma_space, gridded_salinity_gamma_space] = recast_property_from_depth_to_gamma_space(max_x_samples, max_z_samples, max_t_samples, density_bins, gridded_gamma, gridded_temperature, gridded_salinity);
gridded_salinity_gamma_space

But, I cannot figure out how to re-map that back into pressure-latitude space.

There may well be a more sensible way to do this calculation. Basically, I need to know how temperature and salinity map into the spaces of pressure-latitude, which is fixed over time, and density-latitude, which is not fixed over time.

Below is my failing code and a screenshot out of what the error is

gridded_field = dT_dgamma
gridded_field # this is in density space, I want to know where  dT_dgamma sits in depth/pressure space
ds = gridded_field.to_dataset() # make sure it is a dataSet, not a data Array
ds = ds.assign({'pressure': p_levels_to_interpolate_to}) # tell it what pressure is (pressure does not vary with time or space)
grid = Grid(ds, coords={'Z': {'center':'density'}}, periodic=False)
grid.transform(ds.CT, 'Z', p_levels_to_interpolate_to, target_data=ds.pressure)

[edoddridge (Ed Doddridge)]

Possibly you’re after a smarter solution, but to convert back into depth space, I have always just stacked up the thicknesses of the density layers.

E.g.

    Y_layered_plotting = -hbar.mean(dim=['longitude']).cumsum(dim='sigma').compute()

[Aidan (Aidan Heerdegen, ACCESS-NRI Release Team Lead)]

Did @edoddridge’s reply answer your question @kathy.gunn?

If not can I suggest you try uploading a simple example to gist.github.com and link to it here, as I think the source of your error is probably something fairly mundane to do with coordinate names or similar, but it is hard to know for sure without some code to play with.

[kathy.gunn]

@edoddridge

the

Hi Ed,

Thank you for the reply - I’ve just managed to come back to this after the Christmas break. Your solution seems elegant and simple. Just thinking about how I can incorporate this into my code…So hbar is the thickness of each layer? Does that have a standard name in the model output?

Best wishes, Kathy

[edoddridge (Ed Doddridge)]

Yep, hbar is the time-mean thickness of the layers.

It’s not part of the standard model output for any of the ACCESS-OM2 models, so you’ll need to compute it yourself. The easiest way is to regrid the vertical coordinate from depth to density - the calculation is analogous to what you’ve already done with salinity and temperature.

I just had a look and couldn’t see a cookbook example that does it. I can dig up my code if you need it.

[Aidan (Aidan Heerdegen, ACCESS-NRI Release Team Lead)]

@edoddridge

I just had a look and couldn’t see a cookbook example that does it. I can dig up my code if you need it.

Sounds like something that might be useful to add to the hackathon @edoddridge?

[dougiesquire (Dougie Squire, ACCESS-NRI Ocean Modelling Team)]

Hi @kathy.gunn. It’s hard to know exactly what’s wrong with your code without being able to run it, but the following might be helpful.

@kathy.gunn

gridded_field = dT_dgamma
gridded_field # this is in density space, I want to know where  dT_dgamma sits in depth/pressure space
ds = gridded_field.to_dataset() # make sure it is a dataSet, not a data Array
ds = ds.assign({'pressure': p_levels_to_interpolate_to}) # tell it what pressure is (pressure does not vary with time or space)
grid = Grid(ds, coords={'Z': {'center':'density'}}, periodic=False)
grid.transform(ds.CT, 'Z', p_levels_to_interpolate_to, target_data=ds.pressure)

To do this transformation, target_data should be the pressure levels corresponding to the density level dimension in ds. I think you are passing the pressure levels you want to interpolate to.

I don’t know whether it’s good practice, but one way to get the pressure levels in density space is to also interpolate them when you first transform from pressure to density. E.g., you could do a full round trip (pressure → density → pressure) as follows:

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from xgcm import Grid


# Some example test data
pressure = np.arange(2, 12)
density = xr.DataArray(np.log(pressure), coords={'pressure': pressure})
salinity = xr.DataArray(
    np.flip(np.log(pressure) * 0.5 + np.random.rand(len(pressure))), 
    coords={'pressure':pressure})

ds = xr.Dataset({'salinity': salinity, 'density': density})

# Set up xgcm Grid
grid = Grid(ds, coords={'Z': {'center':'pressure'}}, periodic=False)

# Interpolate salinity and pressure onto specified density levels
density_target = np.linspace(0, 3, 50)
salinity_density_space = grid.transform(
    ds.salinity, 'Z', density_target, target_data=ds.density
)
pressure_density_space = grid.transform(
    ds.pressure, 'Z', density_target, target_data=ds.density
)

ds_density_space = xr.Dataset(
    {"salinity": salinity_density_space, "pressure": pressure_density_space}
)

# Set up new xgcm Grid
grid_density_space = Grid(
    ds_density_space, coords={'Z': {'center':'density'}}, periodic=False
)

# Interpolate back to pressure levels
pressure_target = ds.pressure.values
salinity_roundtrip = grid_density_space.transform(
    ds_density_space.salinity, 
    'Z',
    pressure_target, 
    target_data=ds_density_space.pressure
)

Note, because of the interpolation, you don’t get back the answer you started with (though it gets closer the more levels you add to density_target):

salinity.plot(label="original salinity")
salinity_roundtrip.plot(label="roundtripped salinity")
plt.legend()

[kathy.gunn]

dougie - thank you, and once again you were right. I was passing the wrong thing to xgcm.

In my original question this line:
ds = ds.assign({'pressure': p_levels_to_interpolate_to})
should have been
ds = ds.assign({'pressure': pressure_gamma_space})

Where pressure_gamma_space is how pressure varies in density space, that is the pressure of discretized isopycnals.

I’ve pasted my updated and code below in case it is useful for anyone.

# isopycnal layer depths is really pressure in density space, it is simply discretized. Therefore, we can use the time-varying pressure field (in density space) to convert back to depth. 
high_res_p_levels_to_interpolate_to = np.arange(shallow_depth_limit,deep_depth_limit,10);
# xgcm has to work on an xarray.Dataset, so first step is to convert and combine the xarray dataArrays into data sets
ds = dH_dt.to_dataset(name='dH_dt')
# assign any other variables to be converted back into depth space
ds = ds.assign({'CT': gridded_temperature_gamma_space})
# xgcm needs to know how pressure (i.e. the depths of the isopycnals) varies in density space, so we add as a variable
ds = ds.assign({'pressure': pressure_gamma_space})
# Now, do the conversion
# First create an xgcm grid object
grid = Grid(ds, coords={'Z': {'center':'density'}}, periodic=False)
# and transform into density coordinates
gridded_temperature_depth_space = grid.transform(ds.CT, 'Z', high_res_p_levels_to_interpolate_to, target_data=ds.pressure)