Using Xarray for Data read and selection

Use Xarray module to read in model data from nomads server.

This example uses the xarray module to access data from the nomads server for archive NAM analysis data via OPeNDAP. Xarray makes it easier to select times and levels, although you still have to know the coordinate variable name. A simple 500 hPa plot is created after selecting with xarray.

Import all of our needed modules

from datetime import datetime

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.units import units
import numpy as np
import xarray as xr

Accessing data using Xarray¶

# Specify our date/time of product desired
dt = datetime(2016, 4, 16, 18)

# Construct our OPeNDAP access URL
base_url = 'https://www.ncei.noaa.gov/thredds/dodsC/model-namanl-old/'
data = xr.open_dataset(f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/'
                       f'namanl_218_{dt:%Y%m%d}_{dt:%H}00_000.grb').metpy.parse_cf()

syntax error, unexpected WORD_WORD, expecting SCAN_ATTR or SCAN_DATASET or SCAN_ERROR
context: <!DOCTYPE^ HTML PUBLIC "-//IETF//DTD HTML 2.0//EN"><html><head><title>503 Service Unavailable</title></head><body><h1>Service Unavailable</h1><p>The server is temporarily unable to service yourrequest due to maintenance downtime or capacityproblems. Please try again later.</p><p>Additionally, a 404 Not Founderror was encountered while trying to use an ErrorDocument to handle the request.</p></body></html>

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/file_manager.py:219, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
    218 try:
--> 219     file = self._cache[self._key]
    220 except KeyError:

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/lru_cache.py:56, in LRUCache.__getitem__(self, key)
     55 with self._lock:
---> 56     value = self._cache[key]
     57     self._cache.move_to_end(key)

KeyError: [<class 'netCDF4._netCDF4.Dataset'>, ('https://www.ncei.noaa.gov/thredds/dodsC/model-namanl-old/201604/20160416/namanl_218_20160416_1800_000.grb',), 'r', (('clobber', True), ('diskless', False), ('format', 'NETCDF4'), ('persist', False)), 'ff51a862-0ff3-4da5-8b8c-ab008eef2b5f']

During handling of the above exception, another exception occurred:

OSError                                   Traceback (most recent call last)
Cell In[2], line 6
      2 dt = datetime(2016, 4, 16, 18)
      3 
      4 # Construct our OPeNDAP access URL
      5 base_url = 'https://www.ncei.noaa.gov/thredds/dodsC/model-namanl-old/'
----> 6 data = xr.open_dataset(f'{base_url}{dt:%Y%m}/{dt:%Y%m%d}/'
      7                        f'namanl_218_{dt:%Y%m%d}_{dt:%H}00_000.grb').metpy.parse_cf()

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/api.py:613, 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)
    601 decoders = _resolve_decoders_kwargs(
    602     decode_cf,
    603     open_backend_dataset_parameters=backend.open_dataset_parameters,
   (...)    609     decode_coords=decode_coords,
    610 )
    612 overwrite_encoded_chunks = kwargs.pop("overwrite_encoded_chunks", None)
--> 613 backend_ds = backend.open_dataset(
    614     filename_or_obj,
    615     drop_variables=drop_variables,
    616     **decoders,
    617     **kwargs,
    618 )
    619 ds = _dataset_from_backend_dataset(
    620     backend_ds,
    621     filename_or_obj,
   (...)    632     **kwargs,
    633 )
    634 return ds

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/netCDF4_.py:771, in NetCDF4BackendEntrypoint.open_dataset(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, group, mode, format, clobber, diskless, persist, auto_complex, lock, autoclose)
    749 def open_dataset(
    750     self,
    751     filename_or_obj: T_PathFileOrDataStore,
   (...)    768     autoclose=False,
    769 ) -> Dataset:
    770     filename_or_obj = _normalize_path(filename_or_obj)
--> 771     store = NetCDF4DataStore.open(
    772         filename_or_obj,
    773         mode=mode,
    774         format=format,
    775         group=group,
    776         clobber=clobber,
    777         diskless=diskless,
    778         persist=persist,
    779         auto_complex=auto_complex,
    780         lock=lock,
    781         autoclose=autoclose,
    782     )
    784     store_entrypoint = StoreBackendEntrypoint()
    785     with close_on_error(store):

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/netCDF4_.py:529, in NetCDF4DataStore.open(cls, filename, mode, format, group, clobber, diskless, persist, auto_complex, lock, lock_maker, autoclose)
    525 else:
    526     manager = CachingFileManager(
    527         netCDF4.Dataset, filename, mode=mode, kwargs=kwargs, lock=lock
    528     )
--> 529 return cls(manager, group=group, mode=mode, lock=lock, autoclose=autoclose)

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/netCDF4_.py:429, in NetCDF4DataStore.__init__(self, manager, group, mode, lock, autoclose)
    427 self._group = group
    428 self._mode = mode
--> 429 self.format = self.ds.data_model
    430 self._filename = self.ds.filepath()
    431 self.is_remote = is_remote_uri(self._filename)

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/netCDF4_.py:538, in NetCDF4DataStore.ds(self)
    536 @property
    537 def ds(self):
--> 538     return self._acquire()

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/netCDF4_.py:532, in NetCDF4DataStore._acquire(self, needs_lock)
    531 def _acquire(self, needs_lock=True):
--> 532     with self._manager.acquire_context(needs_lock) as root:
    533         ds = _nc4_require_group(root, self._group, self._mode)
    534     return ds

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/contextlib.py:141, in _GeneratorContextManager.__enter__(self)
    139 del self.args, self.kwds, self.func
    140 try:
--> 141     return next(self.gen)
    142 except StopIteration:
    143     raise RuntimeError("generator didn't yield") from None

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/file_manager.py:207, in CachingFileManager.acquire_context(self, needs_lock)
    204 @contextmanager
    205 def acquire_context(self, needs_lock: bool = True) -> Iterator[T_File]:
    206     """Context manager for acquiring a file."""
--> 207     file, cached = self._acquire_with_cache_info(needs_lock)
    208     try:
    209         yield file

File ~/micromamba/envs/metpy-cookbook/lib/python3.14/site-packages/xarray/backends/file_manager.py:225, in CachingFileManager._acquire_with_cache_info(self, needs_lock)
    223     kwargs = kwargs.copy()
    224     kwargs["mode"] = self._mode
--> 225 file = self._opener(*self._args, **kwargs)
    226 if self._mode == "w":
    227     # ensure file doesn't get overridden when opened again
    228     self._mode = "a"

File src/netCDF4/_netCDF4.pyx:2522
-> 2522 'Could not get source, probably due dynamically evaluated source code.'

File src/netCDF4/_netCDF4.pyx:2159
-> 2159 'Could not get source, probably due dynamically evaluated source code.'

OSError: [Errno -70] NetCDF: DAP server error: 'https://www.ncei.noaa.gov/thredds/dodsC/model-namanl-old/201604/20160416/namanl_218_20160416_1800_000.grb'

NAM data is in a projected coordinate and you get back the projection X and Y values in km

# Create a 2-d meshgrid of our x, y coordinates
# manually converted to meters (km * 1000)
#x, y = np.meshgrid(data['x'].values * 1000, data['y'].values * 1000)
x = data.Geopotential_height_isobaric.metpy.x.metpy.convert_units('meter').values
y = data.Geopotential_height_isobaric.metpy.y.metpy.convert_units('meter').values

Getting the valid times in a more useable format

# Get the valid times from the file
vtimes = data.Geopotential_height_isobaric.metpy.time.data.astype('datetime64[ms]').astype('O')
print(vtimes)

Xarray has some nice functionality to choose the time and level that you specifically want to use. In this example the time variable is ‘time’ and the level variable is ‘isobaric1’. Unfortunately, these can be different with each file you use, so you’ll always need to check what they are by listing the coordinate variable names

# print(data.Geopotential_height.coords)
hght_500 = data.Geopotential_height_isobaric.metpy.sel(time1=vtimes[0], vertical=500*units.hPa)
uwnd_500 = data['u-component_of_wind_isobaric'].metpy.sel(time1=vtimes[0], vertical=500*units.hPa)
vwnd_500 = data['v-component_of_wind_isobaric'].metpy.sel(time1=vtimes[0], vertical=500*units.hPa)

Now make the 500-hPa map¶

# Must set data projection, NAM is LCC projection
datacrs = data.Geopotential_height_isobaric.metpy.cartopy_crs

# A different LCC projection for the plot.
plotcrs = ccrs.LambertConformal(central_latitude=45., central_longitude=-100.,
                                standard_parallels=[30, 60])

fig = plt.figure(figsize=(17., 11.))
ax = plt.axes(projection=plotcrs)
ax.coastlines('50m', edgecolor='black')
ax.add_feature(cfeature.STATES, linewidth=0.5)
ax.set_extent([-130, -67, 20, 50], ccrs.PlateCarree())

clev500 = np.arange(5100, 6000, 60)
cs = ax.contour(x, y, mpcalc.smooth_n_point(hght_500, 9, 5), clev500,
                colors='k', linewidths=2.5, linestyles='solid', transform=datacrs)
ax.clabel(cs, fontsize=12, colors='k', inline=1, inline_spacing=8,
          fmt='%i', rightside_up=True, use_clabeltext=True)

# Here we put boxes around the clabels with a black boarder white facecolor
# `labelTexts` necessary as ~cartopy.mpl.contour.GeoContourSet.clabel
# does not return list of texts as of 0.18
for t in cs.labelTexts:
    t.set_bbox({'fc': 'w'})

# Transform Vectors before plotting, then plot wind barbs.
wind_slice = slice(None, None, 16)
ax.barbs(x[wind_slice], y[wind_slice],
         uwnd_500.data[wind_slice, wind_slice], vwnd_500.data[wind_slice, wind_slice],
         length=7, transform=datacrs)

# Add some titles to make the plot readable by someone else
plt.title('500-hPa Geopotential Heights (m)', loc='left')
plt.title(f'VALID: {vtimes[0]}', loc='right');