smooth_contour makes FilledContourPlot disappear

[sgdecker]

What went wrong?

As written, I get no contours in my frontogenesis plot, but if I comment out f.smooth_contour = 2 I do get the expected plot.

Oddly, in the second plot, I see the contours whether I have o.smooth_contour = 2 or not.

This is with matplotlib 3.8.2

Operating System

Linux

Version

Current main branch

Python Version

3.12.0

Code to Reproduce

import datetime
import numpy as np
import metpy.calc as mpcalc
from metpy.io import GempakGrid
from metpy.plots import FilledContourPlot, MapPanel, PanelContainer
from metpy.units import units

nam_file = '/nfs/wxdata1/ldmdata/gempak/model/nam/24103112_nam211.gem'
level = 850*units.hPa
plot_time = datetime.datetime(2024, 10, 31, 12)

p = level.m
gem_data = GempakGrid(nam_file)
u = gem_data.gdxarray(parameter='UREL', date_time=plot_time, level=p)[0] * units('m/s')
v = gem_data.gdxarray(parameter='VREL', date_time=plot_time, level=p)[0] * units('m/s')
T = gem_data.gdxarray(parameter='TMPK', date_time=plot_time, level=p)[0] * units('K')
th = mpcalc.potential_temperature(level, T)
frnt = mpcalc.frontogenesis(th, u, v)

f = FilledContourPlot()
f.data = frnt
f.scale = 1.08e9
f.contours = list(np.arange(-2.5, 2.6, .5))
f.colormap = 'PuOr'
f.colorbar = 'horizontal'
f.smooth_contour = 2

mp = MapPanel()
mp.area = (-105, -70, 30, 50)
mp.layers = 'coastline', 'states', 'borders'
mp.layers_edgecolor = 'brown'
mp.plots = f,

pc = PanelContainer()
pc.size = (11, 8.5)
pc.panels = mp,
pc.show()

omeg = gem_data.gdxarray(parameter='OMEG', date_time=plot_time, level=p)[0] * units('hPa/s')

o = FilledContourPlot()
o.data = omeg
o.contours = range(-9, 10, 1)
o.colormap = 'BrBG_r'
o.colorbar = 'horizontal'
o.smooth_contour = 2
o.plot_units = 'microbar/s'

mp = MapPanel()
mp.area = (-105,-70,30,50)
mp.layers = 'coastline', 'states', 'borders'
mp.layers_edgecolor = 'brown'
mp.plots = o,

pc = PanelContainer()
pc.size = (11, 8.5)
pc.panels = mp,
pc.show()

Errors, Traceback, and Logs

No response

[sgdecker]

This variation doesn't depend on a local file:

import datetime
import numpy as np
import xarray as xr
from xarray.backends import NetCDF4DataStore
import metpy.calc as mpcalc
from metpy.plots import FilledContourPlot, MapPanel, PanelContainer
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)
v = data['v-component_of_wind_isobaric'].metpy.sel(vertical=level)
T = data['Temperature_isobaric'].metpy.sel(vertical=level)
th = mpcalc.potential_temperature(level, T)
frnt = mpcalc.frontogenesis(th, u, v)

f = FilledContourPlot()
f.data = frnt
f.scale = 1.08e9
f.contours = list(np.arange(-2.5, 2.6, .5))
f.colormap = 'PuOr'
f.colorbar = 'horizontal'
f.smooth_contour = 2

mp = MapPanel()
mp.area = (-105, -70, 30, 50)
mp.layers = 'coastline', 'states', 'borders'
mp.layers_edgecolor = 'brown'
mp.plots = f,

pc = PanelContainer()
pc.size = (11, 8.5)
pc.panels = mp,
pc.show()

omeg = data['Vertical_velocity_pressure_isobaric'].metpy.sel(vertical=level)

o = FilledContourPlot()
o.data = omeg
o.contours = range(-9, 10, 1)
o.colormap = 'BrBG_r'
o.colorbar = 'horizontal'
o.smooth_contour = 2
o.plot_units = 'microbar/s'

mp = MapPanel()
mp.area = (-105, -70, 30, 50)
mp.layers = 'coastline', 'states', 'borders'
mp.layers_edgecolor = 'brown'
mp.plots = o,

pc = PanelContainer()
pc.size = (11, 8.5)
pc.panels = mp,
pc.show()

[kgoebber]

So I was able to look into this a bit and it appears that it is an issue with the fact that there are NAN values that result from the frontogenesis calculation. Underlying the smooth_contour parameter is our version of zoom_xarray, which is based on the SciPy ndimage zoom function.

The following will use the MetPy version of the function that makes it easier to work with an xarray DataArray

mpcalc.zoom_xarray(frnt[0]*1.08e9, 2, order=3)

The output ends up being all nan.

So the question is do we want to handle nan differently and is there a way to do that with the SciPy function that would work.

or

Work on our documentation and examples to highlight when something like this might fail and point users to other methods that should work for them. For example, using smooth_field will work - though the nan areas will increase slightly based on how the many iterations are used.

[sgdecker]

Ahh, good find! I do see one gridpoint in the Bahamas where the frontogenesis is NaN; the total deformation is very small there, and I did get a warning: /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)

While I agree that making smooth_contour work with NaNs is still a bug, I now consider there to be an additional bug in the frontogenesis function. It should be more careful with numerics and never return NaN. In those cases, it is the computation of beta that is going awry, but beta should always be computable with enough care (except in two cases: 1) the total deformation is exactly zero, but in that case the contribution of that term to the frontogenesis is zero anyway; 2) the temperature gradient magnitude is exactly zero, but in that case the frontogenesis itself is zero anyway).

I'll submit a separate bug report for this.