Smoothing contours

[sgdecker]

In GEMPAK, if I have jagged contours (which can happen on relatively coarse grids especially), I would use the CONTUR parameter to smooth them. In the Python world, scipy.ndimage.zoom seems to fill this role, but:

1. It returns a plain ndarray, losing all the metadata
2. The result has n times as many points in each direction, so the metadata that was lost wouldn't match the shape of the new data anyway.

Is there a way to overcome these issues and plot smooth contours?

[dopplershift]

This example is what @kgoebber put together for that. It might make sense to turn this into an example in the metpy docs...

[sgdecker]

Thanks for the feedback, but I'm going to be a party pooper in my response to Kevin (which is really a response to both of you).

[kgoebber]

I've got on my to do list to add a trait that can be set to do smoothing. I was planning on implementing the nine point smoother for the declarative syntax (and allow the number of passes to be input), but could be persuaded to use one of the other (e.g., gaussian).

I also have an example of smoothing a parameter, then adding it to an xarray dataset that will work as a substitute until we get that code in the declarative framework: https://github.com/kgoebber/valpo_courses/blob/master/wxtech/notebooks/metpy_calculations_wind_speed.ipynb

[sgdecker]

Thanks for the feedback. I should have been more clear, since "smoothing" is somewhat vague. I am looking to do what GEMPAK's CONTUR parameter does, which is something different than either a Gaussian filter or n-point smoothing. Here is a figure showing what I mean

It's the lower-right panel that I'm trying to plot. I can reproduce the other panels (aside from a projection difference that I didn't bother to match) with this code:

from datetime import datetime

import xarray as xr
from xarray.backends import NetCDF4DataStore
from siphon.catalog import TDSCatalog
from metpy.plots.declarative import ContourPlot, MapPanel, PanelContainer
import metpy.calc as mpcalc

def get_nam211(init_time, valid_time):
    ymd = init_time.strftime('%Y%m%d')
    hr = init_time.strftime('%H')
    filename = f'{ymd}_{hr}00.grib1'
    cat_name = ('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/NAM/'
                'CONUS_80km/NAM_CONUS_80km_' + filename +
                '/catalog.xml')
    ds_name = 'NAM_CONUS_80km_' + filename
    
    cat = TDSCatalog(cat_name)
    ds = cat.datasets[ds_name]
    ncss = ds.subset()
    query = ncss.query()
    query.time(valid_time).variables('all')
    nc = ncss.get_data(query)
    data = xr.open_dataset(NetCDF4DataStore(nc)).metpy.parse_cf()
    return data

init_time = datetime(2021, 9, 23, 12)
plot_time = datetime(2021, 9, 24, 0)

nam = get_nam211(init_time, plot_time)
tmpc = nam['Temperature_isobaric'][0,15,:,:] - 273.15
tmpc_gauss = mpcalc.smooth_gaussian(tmpc, 8)
tmpc_9pt = mpcalc.smooth_n_point(tmpc, 9)

raw = ContourPlot()
raw.data = tmpc
raw.contours = list(range(0,20))
raw.linecolor = 'green'
raw.clabels = True

gauss = ContourPlot()
gauss.data = tmpc_gauss
gauss.contours = list(range(0,20))
gauss.linecolor = 'green'
gauss.clabels = True

ninept = ContourPlot()
ninept.data = tmpc_9pt
ninept.contours = list(range(0,20))
ninept.linecolor = 'green'
ninept.clabels = True

mp1 = MapPanel()
mp1.layout = (2, 2, 1)
mp1.area = [-87, -78, 40.5, 45]
mp1.projection = 'lcc'
mp1.layers = ['states', 'coastline', 'borders']
mp1.title = 'Raw Data'
mp1.plots = [raw]

mp2 = MapPanel()
mp2.layout = (2, 2, 2)
mp2.area = [-87, -78, 40.5, 45]
mp2.projection = 'lcc'
mp2.layers = ['states', 'coastline', 'borders']
mp2.title = 'Gaussian Filter'
mp2.plots = [gauss]

mp3 = MapPanel()
mp3.layout = (2, 2, 3)
mp3.area = [-87, -78, 40.5, 45]
mp3.projection = 'lcc'
mp3.layers = ['states', 'coastline', 'borders']
mp3.title = '9-pt Smoother'
mp3.plots = [ninept]


pc = PanelContainer()
pc.size = (15,10)
pc.panels = [mp1, mp2, mp3]
pc.show()

But what code will give me the lower-right panel?

Here's what GEMPAK has to say about CONTUR:

      CONTUR sets attributes for the contour algorithms:
  
 		subbox factor / number of smoothing passes
  
      Each grid box will be subdivided into the number of subboxes
      specified by the subbox factor.  Increasing this value will
      produce smoother contours, but will increase contouring times
      and make metafiles larger.  The new default value in GEMPAK 5.2
      is 0.
  
      A simple three-point smoothing filter is used if the number
      of smoothing passes is greater than 0.  The default is 0.

The subbox factor is what I'm looking for. This Stack Overflow question pointed me to scipy.ndimage.zoom, and that indeed works... if I'm not trying to draw contours on a map. This is what triggered my question.

[jthielen]

I'm unaware of anything that already does this in the xarray realm, and due to the change in coordinates, you can't just use apply_ufunc or the like with scipy.ndimage.zoom as-is. So, I think this (with some more careful wording/documentation to clarify this "spline interpolation super-sampling" smoothing vs. the existing "noise filtering" smoothing) would be a worthwhile feature request for MetPy. Below is a very quick mockup of a scipy.ndimage.zoom wrapper:

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from scipy.ndimage import zoom as scipy_zoom

def xarray_zoom(
    input, zoom, output=None, order=3, mode='constant', cval=0.0, prefilter=True, *,
    grid_mode=False,
):
    # Zoom data
    zoomed_data = scipy_zoom(
        input.data, zoom, output=output, order=order, mode=mode, cval=cval,
        prefilter=prefilter, grid_mode=grid_mode
    )
    
    # Zoom dimension coordinates
    if not np.iterable(zoom):
        zoom = tuple(zoom for _ in input.dims)
    zoomed_dim_coords = {}
    for dim_name, dim_zoom in zip(input.dims, zoom):
        if dim_name in input.coords:
            zoomed_dim_coords[dim_name] = scipy_zoom(
                input[dim_name].data, dim_zoom, order=order, mode=mode, cval=cval,
                prefilter=prefilter, grid_mode=grid_mode
            )
    
    # Reconstruct (ignoring non-dimension coordinates)
    return xr.DataArray(
        zoomed_data, dims=input.dims, coords=zoomed_dim_coords, attrs=input.attrs
    )
    
# Test with quarter circle
size = 6
x, y = np.mgrid[:size, :size]
distance = np.sqrt((x) ** 2 + (y) ** 2)
da = xr.DataArray(
    distance / distance.max(),
    coords={'y': np.linspace(260, 270, size), 'x': np.linspace(30, 40, size)},
    dims=('y', 'x')
)
da_zoomed = xarray_zoom(da, 3)

# Plot
fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))
da.plot.contour(ax=ax0)
da_zoomed = da_zoomed.plot.contour(ax=ax1)
plt.show()

[sgdecker]

Thanks, and very impressive! My sad attempt to populate the fourth panel:

tmpc_zoom = zoom(tmpc, 3)
x_zoom = zoom(tmpc['x'], 3)
y_zoom = zoom(tmpc['y'], 3)
tmpc_contur = xr.DataArray(tmpc_zoom, dims=['y', 'x'],
                         coords={'x': x_zoom, 'y': y_zoom})
tmpc_contur = tmpc_contur.metpy.assign_crs(grid_mapping_name='lambert_conformal_conic',
                                           standard_parallel=[25, 25],
                                           longitude_of_central_meridian=-95)

gave me an incorrect plot:

I am getting "crs attribute is not available" when I try to plot the result of applying your function to tmpc, but it certainly seems like many steps in the right direction!

[kgoebber]

Okay, so here is what I've got...

A couple of things. There is an easier way to get the data you want without all of the excess code. The structure of the data is uniform and really only requires knowledge of the date to change which dataset you are accessing, then the reading can be done with the OPeNDaP link directly. I use your init_date variable and common formatting aspects for dates with Python f-strings.

Then I use the methods available on the dataset to pull the vertical level and desired plot time and converting to Celsius using MetPy's unit aware capabilities. Doing so keeps the attributes coming along with your variables (they are lost when doing simple calculations like subtracting 273.15), which then allows @jthielen function to work better.

Note: I did have to add one piece to the xarray_zoom function to bring along the metpy_crs variable, but just added that attribute after the work of figuring out the new x and y coords from the zoomed in aspect.

You'll also be able to see how I played a little with the figure size and then the arguments in the pc.save() call to get rid of excess white space.

from datetime import datetime

import metpy.calc as mpcalc
from metpy.plots.declarative import ContourPlot, MapPanel, PanelContainer
from metpy.units import units
from scipy.ndimage import zoom as scipy_zoom
import numpy as np
import xarray as xr

def xarray_zoom(
    input, zoom, output=None, order=3, mode='constant', cval=0.0, prefilter=True, *,
    grid_mode=False,
):
    # Zoom data
    zoomed_data = scipy_zoom(
        input.data, zoom, output=output, order=order, mode=mode, cval=cval,
        prefilter=prefilter, grid_mode=grid_mode
    )
    
    # Zoom dimension coordinates
    if not np.iterable(zoom):
        zoom = tuple(zoom for _ in input.dims)
    zoomed_dim_coords = {}
    for dim_name, dim_zoom in zip(input.dims, zoom):
        if dim_name in input.coords:
            zoomed_dim_coords[dim_name] = scipy_zoom(
                input[dim_name].data, dim_zoom, order=order, mode=mode, cval=cval,
                prefilter=prefilter, grid_mode=grid_mode
            )
    zoomed_dim_coords['metpy_crs'] = input.metpy_crs
    # Reconstruct (ignoring non-dimension coordinates)
    return xr.DataArray(
        zoomed_data, dims=input.dims, coords=zoomed_dim_coords, attrs=input.attrs
    )

init_time = datetime(2021, 9, 23, 12)

nam = xr.open_dataset('https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/NAM/'
                      f'CONUS_80km/NAM_CONUS_80km_{init_time:%Y%m%d}_{init_time:%H}00.grib1').metpy.parse_cf()


plot_time = datetime(2021, 9, 24, 0)

tmpc = nam['Temperature_isobaric'].metpy.sel(vertical=850*units.hPa, time=plot_time).metpy.convert_units('degC')
tmpc_gauss = mpcalc.smooth_gaussian(tmpc, 8)
tmpc_9pt = mpcalc.smooth_n_point(tmpc, 9)
tmpc_zoom = xarray_zoom(tmpc, 10)

raw = ContourPlot()
raw.data = tmpc
raw.contours = list(range(0,20))
raw.linecolor = 'green'
raw.clabels = True

gauss = ContourPlot()
gauss.data = tmpc_gauss
gauss.contours = list(range(0,20))
gauss.linecolor = 'green'
gauss.clabels = True

ninept = ContourPlot()
ninept.data = tmpc_9pt
ninept.contours = list(range(0,20))
ninept.linecolor = 'green'
ninept.clabels = True

zoom = ContourPlot()
zoom.data = tmpc_zoom
zoom.contours = list(range(0,20))
zoom.linecolor = 'green'
zoom.clabels = True

mp1 = MapPanel()
mp1.layout = (2, 2, 1)
mp1.area = [-87, -78, 40.5, 45]
mp1.projection = 'lcc'
mp1.layers = ['states', 'coastline', 'borders']
mp1.title = 'Raw Data'
mp1.plots = [raw]

mp2 = MapPanel()
mp2.layout = (2, 2, 2)
mp2.area = [-87, -78, 40.5, 45]
mp2.projection = 'lcc'
mp2.layers = ['states', 'coastline', 'borders']
mp2.title = 'Gaussian Filter'
mp2.plots = [gauss]

mp3 = MapPanel()
mp3.layout = (2, 2, 3)
mp3.area = [-87, -78, 40.5, 45]
mp3.projection = 'lcc'
mp3.layers = ['states', 'coastline', 'borders']
mp3.title = '9-pt Smoother'
mp3.plots = [ninept]

mp4 = MapPanel()
mp4.layout = (2, 2, 4)
mp4.area = [-87, -78, 40.5, 45]
mp4.projection = 'lcc'
mp4.layers = ['states', 'coastline', 'borders']
mp4.title = 'Scipy Zoom'
mp4.plots = [zoom]

pc = PanelContainer()
pc.size = (15,12)
pc.panels = [mp1, mp2, mp3, mp4]
pc.show()
pc.save('test_zoom.png', dpi=150, bbox_inches='tight')

[sgdecker]

Thanks for putting this together! Because of a bug in metpy.quantify() I just punted on working with the units in my example. I guess the interpolation algorithm isn't quite the same as what GEMPAK does (those 1 contours in Michigan are not in the GEMPAK plot), but it seems reasonably close!

[kgoebber]

If you add the key word argument order=1 to the xarray_zoom function call you get closer to what GEMPAK produced. Simpler interpolation than the default of order=3.