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

GRIB Data Example

GRIB format is commonly used to disemminate 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')
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-2-783584127f97> in <module>
----> 1 ds = xr.tutorial.load_dataset('era5-2mt-2019-03-uk.grib', engine='cfgrib')

/build/python-xarray-6gtuKj/python-xarray-0.16.2/xarray/tutorial.py in load_dataset(*args, **kwargs)
    111     open_dataset
    112     """
--> 113     with open_dataset(*args, **kwargs) as ds:
    114         return ds.load()
    115

/build/python-xarray-6gtuKj/python-xarray-0.16.2/xarray/tutorial.py in open_dataset(name, cache, cache_dir, github_url, branch, **kws)
     76         # May want to add an option to remove it.
     77         if not _os.path.isdir(longdir):
---> 78             _os.mkdir(longdir)
     79
     80         url = "/".join((github_url, "raw", branch, fullname))

FileNotFoundError: [Errno 2] No such file or directory: '/sbuild-nonexistent/.xarray_tutorial_data'

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)
<ipython-input-3-1c96aded89da> 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)
<ipython-input-4-246f01288f90> in <module>
      4 ax = plt.axes(projection=ccrs.Robinson())
      5 ax.coastlines(resolution='10m')
----> 6 plot = ds.t2m[0].plot(cmap=plt.cm.coolwarm, transform=ccrs.PlateCarree(), cbar_kwargs={'shrink':0.6})
      7 plt.title('ERA5 - 2m temperature British Isles March 2019')

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 <lambda>(fig)
    246
    247     if 'png' in formats:
--> 248         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    249     if 'retina' in formats or 'png2x' in formats:
    250         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/usr/lib/python3/dist-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, **kwargs)
    130         FigureCanvasBase(fig)
    131
--> 132     fig.canvas.print_figure(bytes_io, **kw)
    133     data = bytes_io.getvalue()
    134     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)
   2191                            else suppress())
   2192                     with ctx:
-> 2193                         self.figure.draw(renderer)
   2194
   2195                     bbox_inches = self.figure.get_tightbbox(

/usr/lib/python3/dist-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     39                 renderer.start_filter()
     40
---> 41             return draw(artist, renderer, *args, **kwargs)
     42         finally:
     43             if artist.get_agg_filter() is not None:

/usr/lib/python3/dist-packages/matplotlib/figure.py in draw(self, renderer)
   1861
   1862             self.patch.draw(renderer)
-> 1863             mimage._draw_list_compositing_images(
   1864                 renderer, self, artists, self.suppressComposite)
   1865

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

/usr/lib/python3/dist-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     39                 renderer.start_filter()
     40
---> 41             return draw(artist, renderer, *args, **kwargs)
     42         finally:
     43             if artist.get_agg_filter() is not None:

/usr/lib/python3/dist-packages/cartopy/mpl/geoaxes.py in draw(self, renderer, **kwargs)
    477         self._done_img_factory = True
    478
--> 479         return matplotlib.axes.Axes.draw(self, renderer=renderer, **kwargs)
    480
    481     def _update_title_position(self, renderer):

/usr/lib/python3/dist-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     39                 renderer.start_filter()
     40
---> 41             return draw(artist, renderer, *args, **kwargs)
     42         finally:
     43             if artist.get_agg_filter() is not None:

/usr/lib/python3/dist-packages/matplotlib/cbook/deprecation.py in wrapper(*inner_args, **inner_kwargs)
    409                          else deprecation_addendum,
    410                 **kwargs)
--> 411         return func(*inner_args, **inner_kwargs)
    412
    413     return wrapper

/usr/lib/python3/dist-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
   2745             renderer.stop_rasterizing()
   2746
-> 2747         mimage._draw_list_compositing_images(renderer, self, artists)
   2748
   2749         renderer.close_group('axes')

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

/usr/lib/python3/dist-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     39                 renderer.start_filter()
     40
---> 41             return draw(artist, renderer, *args, **kwargs)
     42         finally:
     43             if artist.get_agg_filter() is not None:

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

/usr/lib/python3/dist-packages/cartopy/feature/__init__.py in intersecting_geometries(self, extent)
    300         """
    301         self.scaler.scale_from_extent(extent)
--> 302         return super(NaturalEarthFeature, self).intersecting_geometries(extent)
    303
    304     def with_scale(self, new_scale):

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

/usr/lib/python3/dist-packages/cartopy/feature/__init__.py in geometries(self)
    282         key = (self.name, self.category, self.scale)
    283         if key not in _NATURAL_EARTH_GEOM_CACHE:
--> 284             path = shapereader.natural_earth(resolution=self.scale,
    285                                              category=self.category,
    286                                              name=self.name)

/usr/lib/python3/dist-packages/cartopy/io/shapereader.py in natural_earth(resolution, category, name)
    293     format_dict = {'config': config, 'category': category,
    294                    'name': name, 'resolution': resolution}
--> 295     return ne_downloader.path(format_dict)
    296
    297

/usr/lib/python3/dist-packages/cartopy/io/__init__.py in path(self, format_dict)
    220         else:
    221             # we need to download the file
--> 222             result_path = self.acquire_resource(target_path, format_dict)
    223
    224         return result_path

/usr/lib/python3/dist-packages/cartopy/io/shapereader.py in acquire_resource(self, target_path, format_dict)
    344         target_dir = os.path.dirname(target_path)
    345         if not os.path.isdir(target_dir):
--> 346             os.makedirs(target_dir)
    347
    348         url = self.url(format_dict)

/usr/lib/python3.9/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.9/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.9/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.9/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.9/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.9/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.9/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)
<ipython-input-5-d6577a89d149> 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