You can run this notebook in a live session 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)
    /build/python-xarray-bL1QJU/python-xarray-0.21.1/xarray/tutorial.py in open_dataset(name, cache, cache_dir, engine, **kws)
        113     try:
    --> 114         import pooch
        115     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)
    /tmp/ipykernel_15500/3699872736.py in <module>
    ----> 1 ds = xr.tutorial.load_dataset("era5-2mt-2019-03-uk.grib", engine="cfgrib")

    /build/python-xarray-bL1QJU/python-xarray-0.21.1/xarray/tutorial.py in load_dataset(*args, **kwargs)
        223     open_dataset
        224     """
    --> 225     with open_dataset(*args, **kwargs) as ds:
        226         return ds.load()
        227

    /build/python-xarray-bL1QJU/python-xarray-0.21.1/xarray/tutorial.py in open_dataset(name, cache, cache_dir, engine, **kws)
        114         import pooch
        115     except ImportError as e:
    --> 116         raise ImportError(
        117             "tutorial.open_dataset depends on pooch to download and manage datasets."
        118             " To proceed please install pooch."

    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)
    /tmp/ipykernel_15500/1072882670.py in <module>
    ----> 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)
    /tmp/ipykernel_15500/2626464068.py in <module>
          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 )

    NameError: name 'ds' is not defined

    ---------------------------------------------------------------------------
    PermissionError                           Traceback (most recent call last)
    /usr/lib/python3/dist-packages/IPython/core/formatters.py in __call__(self, obj)
        339                 pass
        340             else:
    --> 341                 return printer(obj)
        342             # Finally look for special method names
        343             method = get_real_method(obj, self.print_method)

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

    /usr/lib/python3/dist-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
       2293                 )
       2294                 with getattr(renderer, "_draw_disabled", nullcontext)():
    -> 2295                     self.figure.draw(renderer)
       2296
       2297             if bbox_inches:

    /usr/lib/python3/dist-packages/matplotlib/artist.py in 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()

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

    /usr/lib/python3/dist-packages/matplotlib/figure.py in draw(self, renderer)
       2808
       2809             self.patch.draw(renderer)
    -> 2810             mimage._draw_list_compositing_images(
       2811                 renderer, self, artists, self.suppressComposite)
       2812

    /usr/lib/python3/dist-packages/matplotlib/image.py 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

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

    /usr/lib/python3/dist-packages/cartopy/mpl/geoaxes.py in draw(self, renderer, **kwargs)
        556         self._done_img_factory = True
        557
    --> 558         return matplotlib.axes.Axes.draw(self, renderer=renderer, **kwargs)
        559
        560     def _update_title_position(self, renderer):

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

    /usr/lib/python3/dist-packages/matplotlib/axes/_base.py in draw(self, renderer)
       3080             renderer.stop_rasterizing()
       3081
    -> 3082         mimage._draw_list_compositing_images(
       3083             renderer, self, artists, self.figure.suppressComposite)
       3084

    /usr/lib/python3/dist-packages/matplotlib/image.py 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

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

    /usr/lib/python3/dist-packages/cartopy/mpl/feature_artist.py in 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)
        151
        152         # Combine all the keyword args in priority order.

    /usr/lib/python3/dist-packages/cartopy/feature/__init__.py in intersecting_geometries(self, extent)
        295         """
        296         self.scaler.scale_from_extent(extent)
    --> 297         return super().intersecting_geometries(extent)
        298
        299     def with_scale(self, new_scale):

    /usr/lib/python3/dist-packages/cartopy/feature/__init__.py in intersecting_geometries(self, extent)
        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:

    /usr/lib/python3/dist-packages/cartopy/feature/__init__.py in 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)

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

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

    /usr/lib/python3/dist-packages/cartopy/io/shapereader.py in 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)
        330
        331         url = self.url(format_dict)

    /usr/lib/python3.10/os.py 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

    /usr/lib/python3.10/os.py 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

    /usr/lib/python3.10/os.py 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

    /usr/lib/python3.10/os.py 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

    /usr/lib/python3.10/os.py 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

    /usr/lib/python3.10/os.py 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

    /usr/lib/python3.10/os.py 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

    PermissionError: [Errno 13] Permission denied: '/sbuild-nonexistent'

    <Figure size 720x720 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)
    /tmp/ipykernel_15500/2485972357.py in <module>
    ----> 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