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-VFGSTu/python-xarray-2022.12.0/xarray/tutorial.py:132, in open_dataset(name, cache, cache_dir, engine, **kws)
    131 try:
--> 132     import pooch
    133 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-VFGSTu/python-xarray-2022.12.0/xarray/tutorial.py:284, in load_dataset(*args, **kwargs)
    247 def load_dataset(*args, **kwargs) -> Dataset:
    248     """
    249     Open, load into memory, and close a dataset from the online repository
    250     (requires internet).
   (...)
    282     load_dataset
    283     """
--> 284     with open_dataset(*args, **kwargs) as ds:
    285         return ds.load()

File /build/python-xarray-VFGSTu/python-xarray-2022.12.0/xarray/tutorial.py:134, in open_dataset(name, cache, cache_dir, engine, **kws)
    132     import pooch
    133 except ImportError as e:
--> 134     raise ImportError(
    135         "tutorial.open_dataset depends on pooch to download and manage datasets."
    136         " To proceed please install pooch."
    137     ) from e
    139 logger = pooch.get_logger()
    140 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
Error in callback <function _draw_all_if_interactive at 0x7f0d8df82200> (for post_execute):
---------------------------------------------------------------------------
PermissionError                           Traceback (most recent call last)
File /usr/lib/python3/dist-packages/matplotlib/pyplot.py:119, in _draw_all_if_interactive()
    117 def _draw_all_if_interactive():
    118     if matplotlib.is_interactive():
--> 119         draw_all()

File /usr/lib/python3/dist-packages/matplotlib/_pylab_helpers.py:132, in Gcf.draw_all(cls, force)
    130 for manager in cls.get_all_fig_managers():
    131     if force or manager.canvas.figure.stale:
--> 132         manager.canvas.draw_idle()

File /usr/lib/python3/dist-packages/matplotlib/backend_bases.py:2054, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
   2052 if not self._is_idle_drawing:
   2053     with self._idle_draw_cntx():
-> 2054         self.draw(*args, **kwargs)

File /usr/lib/python3/dist-packages/matplotlib/backends/backend_agg.py:405, in FigureCanvasAgg.draw(self)
    401 # Acquire a lock on the shared font cache.
    402 with RendererAgg.lock, \
    403      (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
    404       else nullcontext()):
--> 405     self.figure.draw(self.renderer)
    406     # A GUI class may be need to update a window using this draw, so
    407     # don't forget to call the superclass.
    408     super().draw()

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

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

File /usr/lib/python3/dist-packages/matplotlib/figure.py:3071, in Figure.draw(self, renderer)
   3068         # ValueError can occur when resizing a window.
   3070 self.patch.draw(renderer)
-> 3071 mimage._draw_list_compositing_images(
   3072     renderer, self, artists, self.suppressComposite)
   3074 for sfig in self.subfigs:
   3075     sfig.draw(renderer)

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

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

File /usr/lib/python3/dist-packages/cartopy/mpl/geoaxes.py:538, in GeoAxes.draw(self, renderer, **kwargs)
    533         self.imshow(img, extent=extent, origin=origin,
    534                     transform=factory.crs, *factory_args[1:],
    535                     **factory_kwargs)
    536 self._done_img_factory = True
--> 538 return super().draw(renderer=renderer, **kwargs)

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

File /usr/lib/python3/dist-packages/matplotlib/axes/_base.py:3107, in _AxesBase.draw(self, renderer)
   3104         a.draw(renderer)
   3105     renderer.stop_rasterizing()
-> 3107 mimage._draw_list_compositing_images(
   3108     renderer, self, artists, self.figure.suppressComposite)
   3110 renderer.close_group('axes')
   3111 self.stale = False

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

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

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

File /usr/lib/python3/dist-packages/cartopy/feature/__init__.py:305, in NaturalEarthFeature.intersecting_geometries(self, extent)
    298 """
    299 Returns an iterator of shapely geometries that intersect with
    300 the given extent.
    301 The extent is assumed to be in the CRS of the feature.
    302 If extent is None, the method returns all geometries for this dataset.
    303 """
    304 self.scaler.scale_from_extent(extent)
--> 305 return super().intersecting_geometries(extent)

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

File /usr/lib/python3/dist-packages/cartopy/feature/__init__.py:287, in NaturalEarthFeature.geometries(self)
    285 key = (self.name, self.category, self.scale)
    286 if key not in _NATURAL_EARTH_GEOM_CACHE:
--> 287     path = shapereader.natural_earth(resolution=self.scale,
    288                                      category=self.category,
    289                                      name=self.name)
    290     geometries = tuple(shapereader.Reader(path).geometries())
    291     _NATURAL_EARTH_GEOM_CACHE[key] = geometries

File /usr/lib/python3/dist-packages/cartopy/io/shapereader.py:283, in natural_earth(resolution, category, name)
    279 ne_downloader = Downloader.from_config(('shapefiles', 'natural_earth',
    280                                         resolution, category, name))
    281 format_dict = {'config': config, 'category': category,
    282                'name': name, 'resolution': resolution}
--> 283 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:333, in NEShpDownloader.acquire_resource(self, target_path, format_dict)
    331 target_dir = os.path.dirname(target_path)
    332 if not os.path.isdir(target_dir):
--> 333     os.makedirs(target_dir)
    335 url = self.url(format_dict)
    337 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'
---------------------------------------------------------------------------
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:2314, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2308     renderer = _get_renderer(
   2309         self.figure,
   2310         functools.partial(
   2311             print_method, orientation=orientation)
   2312     )
   2313     with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2314         self.figure.draw(renderer)
   2316 if bbox_inches:
   2317     if bbox_inches == "tight":

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

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

File /usr/lib/python3/dist-packages/matplotlib/figure.py:3071, in Figure.draw(self, renderer)
   3068         # ValueError can occur when resizing a window.
   3070 self.patch.draw(renderer)
-> 3071 mimage._draw_list_compositing_images(
   3072     renderer, self, artists, self.suppressComposite)
   3074 for sfig in self.subfigs:
   3075     sfig.draw(renderer)

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

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

File /usr/lib/python3/dist-packages/cartopy/mpl/geoaxes.py:538, in GeoAxes.draw(self, renderer, **kwargs)
    533         self.imshow(img, extent=extent, origin=origin,
    534                     transform=factory.crs, *factory_args[1:],
    535                     **factory_kwargs)
    536 self._done_img_factory = True
--> 538 return super().draw(renderer=renderer, **kwargs)

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

File /usr/lib/python3/dist-packages/matplotlib/axes/_base.py:3107, in _AxesBase.draw(self, renderer)
   3104         a.draw(renderer)
   3105     renderer.stop_rasterizing()
-> 3107 mimage._draw_list_compositing_images(
   3108     renderer, self, artists, self.figure.suppressComposite)
   3110 renderer.close_group('axes')
   3111 self.stale = False

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

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

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

File /usr/lib/python3/dist-packages/cartopy/feature/__init__.py:305, in NaturalEarthFeature.intersecting_geometries(self, extent)
    298 """
    299 Returns an iterator of shapely geometries that intersect with
    300 the given extent.
    301 The extent is assumed to be in the CRS of the feature.
    302 If extent is None, the method returns all geometries for this dataset.
    303 """
    304 self.scaler.scale_from_extent(extent)
--> 305 return super().intersecting_geometries(extent)

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

File /usr/lib/python3/dist-packages/cartopy/feature/__init__.py:287, in NaturalEarthFeature.geometries(self)
    285 key = (self.name, self.category, self.scale)
    286 if key not in _NATURAL_EARTH_GEOM_CACHE:
--> 287     path = shapereader.natural_earth(resolution=self.scale,
    288                                      category=self.category,
    289                                      name=self.name)
    290     geometries = tuple(shapereader.Reader(path).geometries())
    291     _NATURAL_EARTH_GEOM_CACHE[key] = geometries

File /usr/lib/python3/dist-packages/cartopy/io/shapereader.py:283, in natural_earth(resolution, category, name)
    279 ne_downloader = Downloader.from_config(('shapefiles', 'natural_earth',
    280                                         resolution, category, name))
    281 format_dict = {'config': config, 'category': category,
    282                'name': name, 'resolution': resolution}
--> 283 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:333, in NEShpDownloader.acquire_resource(self, target_path, format_dict)
    331 target_dir = os.path.dirname(target_path)
    332 if not os.path.isdir(target_dir):
--> 333     os.makedirs(target_dir)
    335 url = self.url(format_dict)
    337 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