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-HEJSxI/python-xarray-2022.06.0/xarray/tutorial.py in open_dataset(name, cache, cache_dir, engine, **kws)
    123     try:
--> 124         import pooch
    125     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_3821363/3699872736.py in <module>
----> 1 ds = xr.tutorial.load_dataset("era5-2mt-2019-03-uk.grib", engine="cfgrib")

/build/python-xarray-HEJSxI/python-xarray-2022.06.0/xarray/tutorial.py in load_dataset(*args, **kwargs)
    261     load_dataset
    262     """
--> 263     with open_dataset(*args, **kwargs) as ds:
    264         return ds.load()
    265

/build/python-xarray-HEJSxI/python-xarray-2022.06.0/xarray/tutorial.py in open_dataset(name, cache, cache_dir, engine, **kws)
    124         import pooch
    125     except ImportError as e:
--> 126         raise ImportError(
    127             "tutorial.open_dataset depends on pooch to download and manage datasets."
    128             " 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_3821363/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_3821363/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)
   2835
   2836             self.patch.draw(renderer)
-> 2837             mimage._draw_list_compositing_images(
   2838                 renderer, self, artists, self.suppressComposite)
   2839

/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)
    536         self._done_img_factory = True
    537
--> 538         return super().draw(renderer=renderer, **kwargs)
    539
    540     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)
   3089             renderer.stop_rasterizing()
   3090
-> 3091         mimage._draw_list_compositing_images(
   3092             renderer, self, artists, self.figure.suppressComposite)
   3093

/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)
    149         except ValueError:
    150             warnings.warn('Unable to determine extent. Defaulting to global.')
--> 151         geoms = self._feature.intersecting_geometries(extent)
    152
    153         # Combine all the keyword args in priority order.

/usr/lib/python3/dist-packages/cartopy/feature/__init__.py in intersecting_geometries(self, extent)
    301         """
    302         self.scaler.scale_from_extent(extent)
--> 303         return super().intersecting_geometries(extent)
    304
    305     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)
    283         key = (self.name, self.category, self.scale)
    284         if key not in _NATURAL_EARTH_GEOM_CACHE:
--> 285             path = shapereader.natural_earth(resolution=self.scale,
    286                                              category=self.category,
    287                                              name=self.name)

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

/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)
    331         target_dir = os.path.dirname(target_path)
    332         if not os.path.isdir(target_dir):
--> 333             os.makedirs(target_dir)
    334
    335         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 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)
/tmp/ipykernel_3821363/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