You can run this notebook in a live session Binder or view it on Github.

GRIB Data Example

GRIB format is commonly used to disseminate atmospheric model data. With xarray and the cfgrib engine, GRIB data can easily be analyzed and visualized.

[1]:
import xarray as xr
import matplotlib.pyplot as plt

To read GRIB data, you can use xarray.load_dataset. The only extra code you need is to specify the engine as cfgrib.

[2]:
ds = xr.tutorial.load_dataset("era5-2mt-2019-03-uk.grib", engine="cfgrib")
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
File /build/python-xarray-iFXG89/python-xarray-2022.10.0/xarray/tutorial.py:131, in open_dataset(name, cache, cache_dir, engine, **kws)
    130 try:
--> 131     import pooch
    132 except ImportError as e:

ModuleNotFoundError: No module named 'pooch'

The above exception was the direct cause of the following exception:

ImportError                               Traceback (most recent call last)
Cell In [2], line 1
----> 1 ds = xr.tutorial.load_dataset("era5-2mt-2019-03-uk.grib", engine="cfgrib")

File /build/python-xarray-iFXG89/python-xarray-2022.10.0/xarray/tutorial.py:270, in load_dataset(*args, **kwargs)
    233 def load_dataset(*args, **kwargs) -> Dataset:
    234     """
    235     Open, load into memory, and close a dataset from the online repository
    236     (requires internet).
   (...)
    268     load_dataset
    269     """
--> 270     with open_dataset(*args, **kwargs) as ds:
    271         return ds.load()

File /build/python-xarray-iFXG89/python-xarray-2022.10.0/xarray/tutorial.py:133, in open_dataset(name, cache, cache_dir, engine, **kws)
    131     import pooch
    132 except ImportError as e:
--> 133     raise ImportError(
    134         "tutorial.open_dataset depends on pooch to download and manage datasets."
    135         " To proceed please install pooch."
    136     ) from e
    138 logger = pooch.get_logger()
    139 logger.setLevel("WARNING")

ImportError: tutorial.open_dataset depends on pooch to download and manage datasets. To proceed please install pooch.

Let’s create a simple plot of 2-m air temperature in degrees Celsius:

[3]:
ds = ds - 273.15
ds.t2m[0].plot(cmap=plt.cm.coolwarm)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [3], line 1
----> 1 ds = ds - 273.15
      2 ds.t2m[0].plot(cmap=plt.cm.coolwarm)

NameError: name 'ds' is not defined

With CartoPy, we can create a more detailed plot, using built-in shapefiles to help provide geographic context:

[4]:
import cartopy.crs as ccrs
import cartopy

fig = plt.figure(figsize=(10, 10))
ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines(resolution="10m")
plot = ds.t2m[0].plot(
    cmap=plt.cm.coolwarm, transform=ccrs.PlateCarree(), cbar_kwargs={"shrink": 0.6}
)
plt.title("ERA5 - 2m temperature British Isles March 2019")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [4], line 7
      5 ax = plt.axes(projection=ccrs.Robinson())
      6 ax.coastlines(resolution="10m")
----> 7 plot = ds.t2m[0].plot(
      8     cmap=plt.cm.coolwarm, transform=ccrs.PlateCarree(), cbar_kwargs={"shrink": 0.6}
      9 )
     10 plt.title("ERA5 - 2m temperature British Isles March 2019")

NameError: name 'ds' is not defined
---------------------------------------------------------------------------
PermissionError                           Traceback (most recent call last)
File /usr/lib/python3/dist-packages/IPython/core/formatters.py:339, in BaseFormatter.__call__(self, obj)
    337     pass
    338 else:
--> 339     return printer(obj)
    340 # Finally look for special method names
    341 method = get_real_method(obj, self.print_method)

File /usr/lib/python3/dist-packages/IPython/core/pylabtools.py:151, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    148     from matplotlib.backend_bases import FigureCanvasBase
    149     FigureCanvasBase(fig)
--> 151 fig.canvas.print_figure(bytes_io, **kw)
    152 data = bytes_io.getvalue()
    153 if fmt == 'svg':

File /usr/lib/python3/dist-packages/matplotlib/backend_bases.py:2295, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2289     renderer = _get_renderer(
   2290         self.figure,
   2291         functools.partial(
   2292             print_method, orientation=orientation)
   2293     )
   2294     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2295         self.figure.draw(renderer)
   2297 if bbox_inches:
   2298     if bbox_inches == "tight":

File /usr/lib/python3/dist-packages/matplotlib/artist.py:73, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
     71 @wraps(draw)
     72 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 73     result = draw(artist, renderer, *args, **kwargs)
     74     if renderer._rasterizing:
     75         renderer.stop_rasterizing()

File /usr/lib/python3/dist-packages/matplotlib/artist.py:50, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     47     if artist.get_agg_filter() is not None:
     48         renderer.start_filter()
---> 50     return draw(artist, renderer)
     51 finally:
     52     if artist.get_agg_filter() is not None:

File /usr/lib/python3/dist-packages/matplotlib/figure.py:2837, in Figure.draw(self, renderer)
   2834         # ValueError can occur when resizing a window.
   2836 self.patch.draw(renderer)
-> 2837 mimage._draw_list_compositing_images(
   2838     renderer, self, artists, self.suppressComposite)
   2840 for sfig in self.subfigs:
   2841     sfig.draw(renderer)

File /usr/lib/python3/dist-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    130 if not_composite or not has_images:
    131     for a in artists:
--> 132         a.draw(renderer)
    133 else:
    134     # Composite any adjacent images together
    135     image_group = []

File /usr/lib/python3/dist-packages/matplotlib/artist.py:50, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     47     if artist.get_agg_filter() is not None:
     48         renderer.start_filter()
---> 50     return draw(artist, renderer)
     51 finally:
     52     if artist.get_agg_filter() is not None:

File /usr/lib/python3/dist-packages/cartopy/mpl/geoaxes.py:558, in GeoAxes.draw(self, renderer, **kwargs)
    553         self.imshow(img, extent=extent, origin=origin,
    554                     transform=factory.crs, *factory_args[1:],
    555                     **factory_kwargs)
    556 self._done_img_factory = True
--> 558 return matplotlib.axes.Axes.draw(self, renderer=renderer, **kwargs)

File /usr/lib/python3/dist-packages/matplotlib/artist.py:50, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     47     if artist.get_agg_filter() is not None:
     48         renderer.start_filter()
---> 50     return draw(artist, renderer)
     51 finally:
     52     if artist.get_agg_filter() is not None:

File /usr/lib/python3/dist-packages/matplotlib/axes/_base.py:3091, in _AxesBase.draw(self, renderer)
   3088         a.draw(renderer)
   3089     renderer.stop_rasterizing()
-> 3091 mimage._draw_list_compositing_images(
   3092     renderer, self, artists, self.figure.suppressComposite)
   3094 renderer.close_group('axes')
   3095 self.stale = False

File /usr/lib/python3/dist-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    130 if not_composite or not has_images:
    131     for a in artists:
--> 132         a.draw(renderer)
    133 else:
    134     # Composite any adjacent images together
    135     image_group = []

File /usr/lib/python3/dist-packages/matplotlib/artist.py:50, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
     47     if artist.get_agg_filter() is not None:
     48         renderer.start_filter()
---> 50     return draw(artist, renderer)
     51 finally:
     52     if artist.get_agg_filter() is not None:

File /usr/lib/python3/dist-packages/cartopy/mpl/feature_artist.py:150, in FeatureArtist.draw(self, renderer, *args, **kwargs)
    148 except ValueError:
    149     warnings.warn('Unable to determine extent. Defaulting to global.')
--> 150 geoms = self._feature.intersecting_geometries(extent)
    152 # Combine all the keyword args in priority order.
    153 prepared_kwargs = style_merge(self._feature.kwargs,
    154                               self._kwargs,
    155                               kwargs)

File /usr/lib/python3/dist-packages/cartopy/feature/__init__.py:297, in NaturalEarthFeature.intersecting_geometries(self, extent)
    290 """
    291 Returns an iterator of shapely geometries that intersect with
    292 the given extent.
    293 The extent is assumed to be in the CRS of the feature.
    294 If extent is None, the method returns all geometries for this dataset.
    295 """
    296 self.scaler.scale_from_extent(extent)
--> 297 return super().intersecting_geometries(extent)

File /usr/lib/python3/dist-packages/cartopy/feature/__init__.py:106, in Feature.intersecting_geometries(self, extent)
    103 if extent is not None:
    104     extent_geom = sgeom.box(extent[0], extent[2],
    105                             extent[1], extent[3])
--> 106     return (geom for geom in self.geometries() if
    107             geom is not None and extent_geom.intersects(geom))
    108 else:
    109     return self.geometries()

File /usr/lib/python3/dist-packages/cartopy/feature/__init__.py:279, in NaturalEarthFeature.geometries(self)
    277 key = (self.name, self.category, self.scale)
    278 if key not in _NATURAL_EARTH_GEOM_CACHE:
--> 279     path = shapereader.natural_earth(resolution=self.scale,
    280                                      category=self.category,
    281                                      name=self.name)
    282     geometries = tuple(shapereader.Reader(path).geometries())
    283     _NATURAL_EARTH_GEOM_CACHE[key] = geometries

File /usr/lib/python3/dist-packages/cartopy/io/shapereader.py:279, in natural_earth(resolution, category, name)
    275 ne_downloader = Downloader.from_config(('shapefiles', 'natural_earth',
    276                                         resolution, category, name))
    277 format_dict = {'config': config, 'category': category,
    278                'name': name, 'resolution': resolution}
--> 279 return ne_downloader.path(format_dict)

File /usr/lib/python3/dist-packages/cartopy/io/__init__.py:203, in Downloader.path(self, format_dict)
    200     result_path = target_path
    201 else:
    202     # we need to download the file
--> 203     result_path = self.acquire_resource(target_path, format_dict)
    205 return result_path

File /usr/lib/python3/dist-packages/cartopy/io/shapereader.py:329, in NEShpDownloader.acquire_resource(self, target_path, format_dict)
    327 target_dir = os.path.dirname(target_path)
    328 if not os.path.isdir(target_dir):
--> 329     os.makedirs(target_dir)
    331 url = self.url(format_dict)
    333 shapefile_online = self._urlopen(url)

File /usr/lib/python3.10/os.py:215, in makedirs(name, mode, exist_ok)
    213 if head and tail and not path.exists(head):
    214     try:
--> 215         makedirs(head, exist_ok=exist_ok)
    216     except FileExistsError:
    217         # Defeats race condition when another thread created the path
    218         pass

File /usr/lib/python3.10/os.py:215, in makedirs(name, mode, exist_ok)
    213 if head and tail and not path.exists(head):
    214     try:
--> 215         makedirs(head, exist_ok=exist_ok)
    216     except FileExistsError:
    217         # Defeats race condition when another thread created the path
    218         pass

    [... skipping similar frames: makedirs at line 215 (3 times)]

File /usr/lib/python3.10/os.py:215, in makedirs(name, mode, exist_ok)
    213 if head and tail and not path.exists(head):
    214     try:
--> 215         makedirs(head, exist_ok=exist_ok)
    216     except FileExistsError:
    217         # Defeats race condition when another thread created the path
    218         pass

File /usr/lib/python3.10/os.py:225, in makedirs(name, mode, exist_ok)
    223         return
    224 try:
--> 225     mkdir(name, mode)
    226 except OSError:
    227     # Cannot rely on checking for EEXIST, since the operating system
    228     # could give priority to other errors like EACCES or EROFS
    229     if not exist_ok or not path.isdir(name):

PermissionError: [Errno 13] Permission denied: '/sbuild-nonexistent'
<Figure size 1000x1000 with 1 Axes>

Finally, we can also pull out a time series for a given location easily:

[5]:
ds.t2m.sel(longitude=0, latitude=51.5).plot()
plt.title("ERA5 - London 2m temperature March 2019")
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In [5], line 1
----> 1 ds.t2m.sel(longitude=0, latitude=51.5).plot()
      2 plt.title("ERA5 - London 2m temperature March 2019")

NameError: name 'ds' is not defined