By: Kevin Goebbert
This example uses the declarative syntax available through the MetPy package to allow a more convenient method for creating simple maps of atmospheric data. To plot aboslute vorticity, the data is scaled and reassigned to the xarray object for use in the declarative plotting interface.
from datetime import datetime
import xarray as xr
from metpy.plots import declarative
from metpy.units import units# Set date for desired dataset
dt = datetime(2012, 10, 31, 12)
# Open dataset from NCEI
ds = xr.open_dataset('https://www.ncei.noaa.gov/thredds/dodsC/'
f'model-gfs-g4-anl-files-old/{dt:%Y%m}/{dt:%Y%m%d}/'
f'gfsanl_4_{dt:%Y%m%d}_{dt:%H}00_000.grb2'
).metpy.parse_cf()
# Subset Data to be just over CONUS
ds_us = ds.sel(lon=slice(360-160, 360-40), lat=slice(65, 10))---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[2], line 5
2 dt = datetime(2012, 10, 31, 12)
4 # Open dataset from NCEI
----> 5 ds = xr.open_dataset('https://www.ncei.noaa.gov/thredds/dodsC/'
6 f'model-gfs-g4-anl-files-old/{dt:%Y%m}/{dt:%Y%m%d}/'
7 f'gfsanl_4_{dt:%Y%m%d}_{dt:%H}00_000.grb2'
8 ).metpy.parse_cf()
10 # Subset Data to be just over CONUS
11 ds_us = ds.sel(lon=slice(360-160, 360-40), lat=slice(65, 10))
File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/api.py:587, in open_dataset(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, create_default_indexes, inline_array, chunked_array_type, from_array_kwargs, backend_kwargs, **kwargs)
584 kwargs.update(backend_kwargs)
586 if engine is None:
--> 587 engine = plugins.guess_engine(filename_or_obj)
589 if from_array_kwargs is None:
590 from_array_kwargs = {}
File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/plugins.py:212, in guess_engine(store_spec, must_support_groups)
204 else:
205 error_msg = (
206 "found the following matches with the input file in xarray's IO "
207 f"backends: {compatible_engines}. But their dependencies may not be installed, see:\n"
208 "https://docs.xarray.dev/en/stable/user-guide/io.html \n"
209 "https://docs.xarray.dev/en/stable/getting-started-guide/installing.html"
210 )
--> 212 raise ValueError(error_msg)
ValueError: did not find a match in any of xarray's currently installed IO backends ['netcdf4', 'scipy', 'gini']. Consider explicitly selecting one of the installed engines via the ``engine`` parameter, or installing additional IO dependencies, see:
https://docs.xarray.dev/en/stable/getting-started-guide/installing.html
https://docs.xarray.dev/en/stable/user-guide/io.htmlContour Intervals¶
Since absolute vorticity rarely goes below zero in the Northern Hemisphere, we can set up a list of contour levels that doesn’t include values near but greater than zero. The following code yields a list containing: [-8, -7, -6, -5, -4, -3, -2, -1, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]
# Absolute Vorticity colors
# Use two different colormaps from matplotlib and combine into one color set
clevs_500_avor = list(range(-8, 1, 1))+list(range(8, 46, 1))The Plot¶
Using the declarative interface in MetPy to plot the 500-hPa Geopotential Heights and Absolute Vorticity.
# Set Contour Plot Parameters
contour = declarative.ContourPlot()
contour.data = ds_us
contour.time = dt
contour.field = 'Geopotential_height_isobaric'
contour.level = 500 * units.hPa
contour.linecolor = 'black'
contour.linestyle = '-'
contour.linewidth = 2
contour.clabels = True
contour.contours = list(range(0, 20000, 60))
# Set Color-filled Contour Parameters
cfill = declarative.FilledContourPlot()
cfill.data = ds_us
cfill.time = dt
cfill.field = 'Absolute_vorticity_isobaric'
cfill.level = 500 * units.hPa
cfill.contours = clevs_500_avor
cfill.colormap = 'PuOr_r'
cfill.image_range = (-45, 45)
cfill.colorbar = 'horizontal'
cfill.scale = 1e5
# Panel for plot with Map features
panel = declarative.MapPanel()
panel.layout = (1, 1, 1)
panel.area = 'uslcc'
panel.projection = 'area'
panel.layers = ['coastline', 'states', 'borders']
panel.layers_edgecolor = ['black', 'grey', 'black']
panel.layers_linewidth = [1.25, .75, 1]
panel.title = (f'{cfill.level} GFS Geopotential Heights'
f'and Absolute Vorticity at {dt}')
panel.plots = [cfill, contour]
# Bringing it all together
pc = declarative.PanelContainer()
pc.size = (15, 14)
pc.panels = [panel]
pc.show()