Frontogenesis should not return nan

[sgdecker]

What went wrong?

The frontogenesis function sometimes returns nans, but it should be more careful and never do so (assuming the input data doesn't contain nans). Given the warning, I assume overflow is occurring during the computation of the angle beta, but I haven't investigated further yet.

Operating System

Linux

Version

latest commit on main branch

Python Version

3.12.1

Code to Reproduce

import datetime
import xarray as xr
from xarray.backends import NetCDF4DataStore
import metpy.calc as mpcalc
from metpy.units import units
from siphon.catalog import TDSCatalog

level = 850*units.hPa
plot_time = datetime.datetime(2024, 10, 31, 12)

ds_name = 'NAM_CONUS_80km_20241031_1200.grib1'
cat_name = ('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/NAM/'
            'CONUS_80km/' + ds_name + '/catalog.xml')

cat = TDSCatalog(cat_name)
ds = cat.datasets[ds_name]
ncss = ds.subset()
query = ncss.query()
query.time(plot_time).variables('all')
nc = ncss.get_data(query)
data = xr.open_dataset(NetCDF4DataStore(nc)).metpy.parse_cf()

u = data['u-component_of_wind_isobaric'].metpy.sel(vertical=level).metpy.assign_latitude_longitude()
v = data['v-component_of_wind_isobaric'].metpy.sel(vertical=level).metpy.assign_latitude_longitude()
T = data['Temperature_isobaric'].metpy.sel(vertical=level).metpy.assign_latitude_longitude()
th = mpcalc.potential_temperature(level, T)
frnt = mpcalc.frontogenesis(th, u, v)

print(frnt.isel(x=3,y=1).values[0])

Errors, Traceback, and Logs

/home/decker/local/miniconda3/envs/met_course/lib/python3.12/site-packages/pint/facets/plain/quantity.py:998: RuntimeWarning: invalid value encountered in divide
magnitude = magnitude_op(new_self._magnitude, other._magnitude)

[sgdecker]

Note that in my original case (in which I used smooth_gaussian on the data before computing frontogenesis), the warning was slightly different, as was the location of the nan:

/usr/local/python/3.8/envs/met_course/lib/python3.11/site-packages/pint/facets/numpy/numpy_func.py:307: RuntimeWarning: invalid value encountered in arcsin
  result_magnitude = func(*stripped_args, **stripped_kwargs)

The code above is based on simplifying my original code to its essence, but there could be two ways the frontogenesis calculation is going bad here.

[DWesl]

Assuming the warnings are coming from frontogenesis rather than one of the functions it calls, the only division is here:

MetPy/src/metpy/calc/kinematics.py

Line 570 in b00b5d0

beta = np.arcsin((-ddx_theta * np.cos(psi) - ddy_theta * np.sin(psi)) / mag_theta)

where the denominator is the magnitude of the gradient of the potential temperature. Assuming you've checked the grid spacing to ensure that doesn't go to zero anywhere, that one would be due to some 3x3 box having a constant potential temperature. One way to avoid follow-on effects would be to set the angle to zero at locations where the magnitude of the gradient is zero, that is:

beta[mag_theta == 0] = 0

Gaussian smoothing would almost certainly make the potential temperature field non-constant, so the calculation can move to the next function call, the arcsine, where it complains of input values outside of $[-1.1]$. I would guess this to be roundoff error, so the problem is arcsin(1+1e-6) or arcsin(1+1e-15). One option would be np.clip to ensure the values are between -1 and 1, another would be to use trigonometric identities (cosine of difference, co-function, and cosine double-angle are the most visible) to avoid that problem.

dir_theta = np.arctan2(ddy_theta/mag_theta, ddx_theta/mag_theta)
beta = np.arcsin(-np.cos(dir_theta - psi))

or

beta = -(np.pi/2 - (dir_theta - psi))

with some logic to wrap the result into $[-\frac{\pi}{2}, \frac{\pi}{2}]$, though symmetry considerations could avoid the wrapping logic, given the next line:

MetPy/src/metpy/calc/kinematics.py

Line 572 in b00b5d0

return 0.5 * mag_theta * (tdef * np.cos(2 * beta) - div)

I suspect using the cosine double-angle identity below might be the best way to do things, especially if you also use sin_beta as a variable instead of beta, dropping the arcsine in the first line. You may want to add a line or two of comments to explain things if you go that route.

$$\cos(2\beta) = \cos^2(\beta) - \sin^2(\beta) = 1 - 2 \sin^2(\beta) = 1 - 2 \sin^2(\arcsin(-\cos(\alpha_{\nabla\theta} - \psi))) = 1 - 2 \cos^2(\alpha_{\nabla\theta} - \psi)$$

where cos(dir_theta - psi) = (ddx_theta * cos(psi) + ddy_theta * sin(psi)) / mag_theta, as line 570 currently has it.

Does one of these changes fix things for you?

[sgdecker]

Great sleuthing @DWesl ! There are two different ways I am getting a nan, so I am setting up two test cases to step through with pdb to see exactly what is going wrong.

For test case 1, it is indeed division by zero in the calculation of beta, because with the grid packing done in producing the dataset, you can have neighboring grid points where the temperature is exactly the same value, and that is the case here. Your suggested fix:

beta[mag_theta == 0] = 0

takes care of this one.

For test case 2, where I had done the gaussian smoothing (so the temperatures are no longer identical), something else is happening, but I am out of time at the moment to investigate. I'll return to this when I get the chance.

[sgdecker]

In my second test case, the problem is that the arcsin argument is ever so slightly greater than 1 due to round-off error. Using np.clip as suggested fixes that one!