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)
/build/python-xarray-G8o5tX/python-xarray-2022.06.0~rc1/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_10394/3699872736.py in <module>
----> 1 ds = xr.tutorial.load_dataset("era5-2mt-2019-03-uk.grib", engine="cfgrib")

/build/python-xarray-G8o5tX/python-xarray-2022.06.0~rc1/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-G8o5tX/python-xarray-2022.06.0~rc1/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_10394/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_10394/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_10394/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