gradient_richardson_number

[Marston]

I'm trying to generate the Richardson number, but found in lieu, this gradient function in metpy.calc.
gri = mc.gradient_richardson_number(gmh, theta, u, v, vertical_dim=1) but the results are filled with values ranging from inf to -7.0 x 1e32. There is a warning:

/home/jyn/.local/lib/python3.10/site-packages/pint/facets/plain/quantity.py:1119: 
RuntimeWarning: divide by zero encountered in true_divide
magnitude = magnitude_op(new_self._magnitude, other._magnitude)

I actuality, I just want the Ri number and not the gradient. There is a gradient function in Metpy, so I'm wondering why the Ri gradient given instead. Is there a way to retrieve the Ri number?

[dopplershift]

That function is not calculating the gradient of the Richardson number, it's calculating what's known as the "gradient Richardson number". See the docs to see the details of what the equation looks like. This is trying to do full derivatives in the vertical direction. What are the dimensions of your data?

Is it possible that "Bulk Richardson number" is what you're really looking for? This is similar, but explicitly looking across a similar layer. We don't have that yet (#628), though we'd be happy to see someone contribute an implementation.

[Marston]

I do not have my Holton handy (still buried in some box in the attic). Just went with the name.
I'm looking into how to calculate it but I'm not a great programmer. The difficulty is the finding a module for PDE. I looked into how Metpy calculated derivatives, but I didn't find one for a PDE. Is there such a function in Metpy?
I can try and calculate the BRN. I'll look into it next week. If I'm successful, I would gladly contribute to metpy.

[dopplershift]

In general, you're always approximating partial derivatives using finite differencing, which is what MetPy does. So it sounds like gradient_richardson_number might be what you're looking for--we just need to figure out what's causing the error. So we really need to know what your data look like? What are the dimension sizes? What are the corresponding coordinates? The metadata?

[Marston]

import pandas as pd
import numpy as np
import xarray as xr
from metpy.units import units
import metpy.calc as mc

def create_da(arr, time, lev, lat, lon, units, description):
    ts = xr.DataArray(arr,
        coords={
            "time": time,
            "lev": xr.DataArray(lev, dims=('lev',), attrs={'units': 'Pa', 'axis': 'Z'}),
            "lat": xr.DataArray(lat, dims=('lat',), attrs={'units': 'degrees', 'axis': 'Y'}),
            "lon": xr.DataArray(lon, dims=('lon',), attrs={'units': 'degrees', 'axis': 'X'}),
        },
    )
    ts.attrs['units'] = units
    ts.attrs['description'] = description
    return ts

tgmh = np.random.uniform(low=-615, high=31700, size=(1, 31, 721, 1440)) 
tmxr = np.random.uniform(low=0, high=0.02766, size=(1, 31, 721, 1440))  
ttheta = np.random.uniform(low=3.01000366, high=640.92000732, size=(1, 31, 721, 1440))  
vWind = np.random.uniform(low=-63.43000031, high=67.23999786, size=(1, 31, 721, 1440))
uWind = np.random.uniform(low=-48.02999878, high=115.94999695, size=(1, 31, 721, 1440))  
time = pd.date_range(start='1/1/2018', periods=1, freq='H')
lev = np.array([101300, 100000,  97500,  95000,  92500,  90000,  87500,  85000,
                82500,  80000,  77500,  75000,  72500,  70000,  65000,  60000,
                55000,  50000,  45000,  40000,  35000,  30000,  25000,  20000,
                15000,  10000,   7000,   5000,   3000,   2000,   1000])
res = 0.25
lon = np.arange(-180, 179.75+res, res)
lat = np.arange(90, -90-res, -res)
tgmh = create_da(tgmh, time, lev, lat, lon, 'm', "Geometric height")
tmxr = create_da(tmxr, time, lev, lat, lon, 'kg/kg', "Mixing Ratio")
ttheta = create_da(ttheta, time, lev, lat, lon, 'K', "gradient Ri")
tv = create_da(vWind, time, lev, lat, lon, 'm/s', 'V Wind')
tu = create_da(uWind, time, lev, lat, lon, 'm/s', 'U Wind')
# taken from NCL 
tthv = ttheta * ( 1 + 0.61 * tmxr )
tthv.attrs['units'] = 'K'

tgri = mc.gradient_richardson_number(tgmh, tthv, tu, tv, vertical_dim=1)
#-----------------------
print(tgri.max(), tgri.min())

<xarray.DataArray ()>
<Quantity(9.77938957e+10, 'dimensionless')> <xarray.DataArray ()>
<Quantity(-6.09264569e+08, 'dimensionless')>

This test using canned data from my model (but taking the maxima and minima from my data to generate some random data) gives reasonable results (see max and min at the end of the code block.
However, when I put in my model data instead, I get inf and nan. Not sure where this is coming from since I use the maxima and minima to generate the input data.

u max = 115.94999694824219 | min = -48.029998779296875
v max = 67.23999786376953 | min = -63.43000030517578
theta max = 640.9200073242188 | min = 3.0100036621093977
gmh max = 31700.0 | min = -615.0
mxr max = 0.027659999206662178 | min = 0.0
thv max = 640.9200073242188 | min = 3.010058745175023

# My output using raw model output data:
# gri max = inf | min = -6.934701267702697e+32  
# I noticed that if I drop all the indices where the results are not finite then the dimension shows a systematic pattern i.e., some levs and 
# latitudes fall away:
# After masking: dims = (1, 16, 695, 1440)
# Before masking: dims = (1, 31, 721, 1440)

mask = gri.where(~np.isfinite(gri), drop=True)
mask.coords

Coordinates:
  * time     (time) datetime64[ns] 2021-09-01
  * lev      (lev) int64 101300 100000 97500 95000 ... 72500 70000 65000 60000
  * lat      (lat) float64 90.0 83.5 83.25 83.0 ... -89.25 -89.5 -89.75 -90.0
  * lon      (lon) float64 -180.0 -179.8 -179.5 -179.2 ... 179.2 179.5 179.8

# Notice how all the levels above 600 hPa are not returned as finite values. 
# Also notice how the only these lats are missing:
# [89.75,  89.5,  89.25,  89.0,  88.75,  88.5,  88.25,  88.0,  87.75,  87.5,  87.25,  87.0,  86.75,  86.5,  86.25,  86.0,  85.75, 85.5,  85.25,  85.0,  
# 84.75, 84.5, 84.25, 84.0,  83.75,  -46.5]

Finally, I switched out the data points in the canned data and the results were stable and clean. Therefore I can only conclude that the error is coming from the dataarrays that are being passed to metpy. I can try and upload this data in a pickle file but it's about 1.5 GB. Is there a way to upload the input data of that size?

Late addition:

I copied the model data to the canned dataarrays (as in example) above to see if the inf values were coming from the data structures, but the results showed inf, which means it is the values in the data that is causing the problem, not something in the dataarrays.

[sgdecker]

Here's my guess. The formula for the gradient Richardson number has vertical shear (squared) in the denominator. Your dataset is on isobaric coordinates going as high as 1013 hPa. There are many areas of the world where that is underground. If the dataset is using a constant wind for points that are underground, the shear will be zero, and division by zero will happen, leading to inf or nan depending on if the shear comes out to be exactly zero or just something really close to zero (because of roundoff error). So I'd check to see if you have winds that are constant with height at underground gridpoints.

[Marston]

That makes a lot of sense. I didn't realize that - how blind of me. I'll have a deeper look at this.