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?