Open tomvothecoder opened 1 year ago
tomvothecoder
commented
1 year ago 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
commented
1 year ago @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
commented
1 year ago Thanks @jasonb5!
dcherian
commented
1 year ago 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
commented
1 year ago 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
commented
1 year ago 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
commented
1 year ago 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?
Also MetPy has the log_interpolate_1d function:
pochedls
commented
11 months ago 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@pochedls @tomvothecoder Created PR in the cf_xarray repo https://github.com/xarray-contrib/cf-xarray/pull/528.
Is your feature request related to a problem?
Related to #45 and #388.
Discussed in https://github.com/xCDAT/xcdat/discussions/440
Describe the solution you'd like
First, check of any other xarray-based package includes this feature.
cdutil.vertical.reconstructPressureFromHybrid)Describe alternatives you've considered
No response
Additional context
No response