[Feature]: Consider adding vertical regridding feature to reconstruct pressure from hybrid

[tomvothecoder]

Is your feature request related to a problem?

Related to #45 and #388.

Discussed in #440

^{Originally posted by tomvothecoder March 27, 2023}
Hi @xCDAT/core-developers or anybody in the xCDAT community. Have any of you come across a package that implements an xarray-based function similar to cdutil.vertical.reconstructPressureFromHybrid? I am working on refactoring e3sm_diags and e3sm_to_cmip to use xarray/xCDAT and there are references to this function that I need to replace. Thanks!

Docstring:

def reconstructPressureFromHybrid(ps, A, B, Po):
    """
    Reconstruct the Pressure field on sigma levels, from the surface pressure
    :param Ps: Surface pressure
    :param A: Hybrid Conversion Coefficient, such as: p=B.ps+A.Po.
    :param B: Hybrid Conversion Coefficient, such as: p=B.ps+A.Po.
    :param Po: Hybrid Conversion Coefficient, such as: p=B.ps+A.Po
    :param Ps: surface pressure
    .. note::
        A and B are 1d sigma levels.
        Po and Ps must have same units.
    :returns: Pressure field, such as P=B*Ps+A*Po.
    :Example:
        .. doctest:: vertical_reconstructPressureFromHybrid
            >>> P=reconstructPressureFromHybrid(ps,A,B,Po)
    """
	...

Examples Usages:

• https://github.com/E3SM-Project/e3sm_diags/blob/4a82b0ec3b2a81233b294fa505c6aa8228ecbd27/e3sm_diags/driver/utils/general.py#L132-L146
• https://github.com/E3SM-Project/e3sm_to_cmip/blob/70f7dad381ef304b022826b540df0e7b46ba9912/e3sm_to_cmip/cmor_handlers/vars/pfull.py#L38-L47
• https://github.com/E3SM-Project/e3sm_to_cmip/blob/70f7dad381ef304b022826b540df0e7b46ba9912/e3sm_to_cmip/cmor_handlers/vars/phalf.py#L39-L50

Describe the solution you'd like

First, check of any other xarray-based package includes this feature.

• If exist - check if the implementation meets our requirements (based on cdutil.vertical.reconstructPressureFromHybrid)
  • If meets our requirements -- just use this library API instead
  • If does not meet our requirements -- open a new PR after Add API extending xgcm vertical regridding #388 is merged
• If none exist -- open a new PR after Add API extending xgcm vertical regridding #388 is merged

Describe alternatives you've considered

No response

Additional context

No response

[tomvothecoder]

Hey @jasonb5, can you follow the "Describe the solution you'd like" and see what we should do?
If no API exists in another xarray-based package, I'm hoping to have this available by the next release around end of May.

e3sm_to_cmip and e3sm_diags are being refactored to use xarray/xcdat and both use cdutil.vertical.reconstructPressureFromHybrid() which blocks progress for these efforts.

[jasonb5]

@tomvothecoder I'll check this out and see if I can find anything existing otherwise I'll go off the original plan I had for this in #388.

[tomvothecoder]

Thanks @jasonb5!

[dcherian]

I would happily merge this in cf-xarray. We already support a few ocean parametric coordinates. The code would go here https://github.com/xarray-contrib/cf-xarray/blob/c36a3687b3820176c394580f86a67ca53b9cadbd/cf_xarray/accessor.py#L2507

[tomvothecoder]

Hey @dcherian, yes @jasonb5 came across that feature in cf-xarray! He mentioned it might be a good opportunity to make a PR to cf-xarray to support decoding other parametric vertical coordinates.

Ideally, we'd like to have this feature out relatively soon so we can start refactoring some other packages with it.
If cf-xarray has a quick turnaround time for releases, we can definitely contribute with a PR. Alternatively, we can implement this feature in xCDAT in the interim and then port it over to cf-xarray at a later date.

[dcherian]

If cf-xarray has a quick turnaround time for releases

Yes, at the moment I just release when there's some nice new feature.

[tomvothecoder]

Hey @jasonb5 and @chengzhuzhang, I found these features from geocat-comp. Can you take a look and see if it fits our requirements so we don't end up duplicating functionalities?

• https://geocat-comp.readthedocs.io/en/latest/user_api/generated/geocat.comp.interpolation.interp_hybrid_to_pressure.html
• https://geocat-comp.readthedocs.io/en/latest/user_api/generated/geocat.comp.interpolation.interp_sigma_to_hybrid.html

Also MetPy has the log_interpolate_1d function:

• https://unidata.github.io/MetPy/latest/examples/sigma_to_pressure_interpolation.html#sphx-glr-examples-sigma-to-pressure-interpolation-py
• Usuage request: Sigma to pressure conversion using ERA-I Unidata/MetPy#1037

[pochedls]

In testing xcdat v0.6.0, I wanted to be able to create the pressure coordinate from hybrid (generically across a bunch of models). I ended up creating the following function, which takes the embedded formula (e.g., ds.lev.formula) and converts it into a statement that can be executed by xarray (e.g., interprete_formula(ds.lev.formula) yields p = ds['a']*ds['p0']+ds['b']*ds['ps'] where ds.lev.formula = p = a*p0 + b*ps). So you can then do:

f = ds.plev.formula
fstring = interprete_formula(f)
exec(f)
print(p)

<xarray.DataArray (lev: 91, time: 1980, lat: 128, lon: 256)>
...
Coordinates:

• lev (lev) float64 0.9988 0.9959 0.992 ... 5.61e-05 2.951e-05 9.869e-06
• lat (lat) float64 -88.93 -87.54 -86.14 -84.74 ... 86.14 87.54 88.93
• lon (lon) float64 0.0 1.406 2.812 4.219 ... 354.4 355.8 357.2 358.6
• time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00

This is a clunky function (though it so far appears to be working on CMIP data...). I am absolutely sure you developers can come up with a better method for interpreting the hybrid formula on the fly, but I thought I'd put it here in case it is useful. I also recognize using exec/eval is not a best practice...maybe this can be avoided somehow (I have ideas if anyone ends up working on this).

def interprete_formula(f):
    operator_chars = ['*', '(', ')', '+', '-', '/']
    # first find/remove all parenthesis that are specifying axes
    # p(n,k,j,i) = ap(k) + b(k)*ps(n,j,i) --> p = ap + b*ps
    # get formula indices corresponding to a left hand parenthesis (
    pinds = [i for i in range(len(f)) if f[i] == '(']
    # initialize list to replace sections of formula (e.g., '(i,j,k)' --> '')
    replace = []
    # loop over left hand parenthesis instances
    for p in pinds:
        pref = ''  # initialize parenthesis string
        stopsearch = False  # set stop flag to false
        # loop from left hand parentheis to end of formula
        for i in range(p, len(f)):
            # if stopsearch, continue
            if stopsearch:
                continue
            # get indexed character
            c = f[i]
            # check for close parenthesis )
            # if closed, stopsearch
            if c == ')':
                pref += c # add character to parenthesis string
                stopsearch = True
            else:
                pref += c # add character to parenthesis string
        # check for arithmetic operators inside parenthesis string
        ocs = [c for c in pref[1:-1] if c in operator_chars]
        # if no arithmetic operators, queue parenthesis string for removal
        if len(ocs) == 0:
            replace.append(pref)
    # remove all parenthesis strings with no arithmetic operators
    for r in replace:
        f = f.replace(r, '')

    # further process formula
    fout = ''  # initialize blank formula
    vid = ''  # initialize blank variable
    lhs = f.split('=')[0] # get lhs
    f = f.split('=')[1] # get rhs
    # loop over each char
    for s in f:
        # ignore blank spaces
        if s == ' ':
            continue
        # if arithmetic operator, need to add to formula
        if s in operator_chars:
            # if a variable is not empty, write it out first
            if len(vid) > 0:
                fout = fout + 'ds["' + vid + '"]'
                vid = ''  # reset variable accumulator
            # write out arithmetic operator
            fout = fout + s
        else:
            # char must be part of a variable - add to variable accumulator
            vid += s
    # finished looping over all characters
    # write out variable to formula
    fout = fout + 'ds["' + vid + '"]'
    # add lhs
    fout = lhs + ' = ' + fout
    # return formula
    return fout