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
Fontconfig error: No writable cache directories
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-VrETFG/python-xarray-2023.01.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-VrETFG/python-xarray-2023.01.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-VrETFG/python-xarray-2023.01.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 0x7f641def2480> (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:3082, in Figure.draw(self, renderer)
3079 # ValueError can occur when resizing a window.
3081 self.patch.draw(renderer)
-> 3082 mimage._draw_list_compositing_images(
3083 renderer, self, artists, self.suppressComposite)
3085 for sfig in self.subfigs:
3086 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:3100, in _AxesBase.draw(self, renderer)
3097 a.draw(renderer)
3098 renderer.stop_rasterizing()
-> 3100 mimage._draw_list_compositing_images(
3101 renderer, self, artists, self.figure.suppressComposite)
3103 renderer.close_group('axes')
3104 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 <frozen os>:215, in makedirs(name, mode, exist_ok)
File <frozen os>:215, in makedirs(name, mode, exist_ok)
[... skipping similar frames: makedirs at line 215 (3 times)]
File <frozen os>:215, in makedirs(name, mode, exist_ok)
File <frozen os>:225, in makedirs(name, mode, exist_ok)
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:3082, in Figure.draw(self, renderer)
3079 # ValueError can occur when resizing a window.
3081 self.patch.draw(renderer)
-> 3082 mimage._draw_list_compositing_images(
3083 renderer, self, artists, self.suppressComposite)
3085 for sfig in self.subfigs:
3086 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:3100, in _AxesBase.draw(self, renderer)
3097 a.draw(renderer)
3098 renderer.stop_rasterizing()
-> 3100 mimage._draw_list_compositing_images(
3101 renderer, self, artists, self.figure.suppressComposite)
3103 renderer.close_group('axes')
3104 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 <frozen os>:215, in makedirs(name, mode, exist_ok)
File <frozen os>:215, in makedirs(name, mode, exist_ok)
[... skipping similar frames: makedirs at line 215 (3 times)]
File <frozen os>:215, in makedirs(name, mode, exist_ok)
File <frozen os>:225, in makedirs(name, mode, exist_ok)
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