ERA5_ARCO Pressure Level Data Exploration: Matplotlib¶

Overview¶

A team at Google Research & Cloud are making parts of the ECMWF Reanalysis version 5 (aka ERA-5) accessible in a Analysis Ready, Cloud Optimized (aka ARCO) format.

In this notebook, we will do the following:

Access the ERA-5 ARCO catalog
Select a particular dataset and variable from the catalog

Prerequisites¶

Concepts	Importance	Notes
Cartopy	Necessary
Xarray	Necessary

Time to learn: 30 minutes

Imports¶

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from metpy import calc as mpcalc
from metpy.units import units
import metpy
import numpy as np
from datetime import datetime, timedelta

Access the ARCO ERA-5 catalog on Google Cloud¶

Let’s open the 37-level isobaric surfaces reanalysis Zarr file

reanalysis = xr.open_zarr(
    'gs://gcp-public-data-arco-era5/ar/1959-2022-full_37-1h-0p25deg-chunk-1.zarr-v2', 
    consolidated=True
)

reanalysis

<xarray.Dataset>
Dimensions:                                           (time: 552264,
                                                       latitude: 721,
                                                       longitude: 1440,
                                                       level: 37)
Coordinates:
  * latitude                                          (latitude) float32 90.0...
  * level                                             (level) int64 1 2 ... 1000
  * longitude                                         (longitude) float32 0.0...
  * time                                              (time) datetime64[ns] 1...
Data variables: (12/31)
    10m_u_component_of_wind                           (time, latitude, longitude) float32 dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>
    10m_v_component_of_wind                           (time, latitude, longitude) float32 dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>
    2m_temperature                                    (time, latitude, longitude) float32 dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>
    angle_of_sub_gridscale_orography                  (latitude, longitude) float32 dask.array<chunksize=(721, 1440), meta=np.ndarray>
    anisotropy_of_sub_gridscale_orography             (latitude, longitude) float32 dask.array<chunksize=(721, 1440), meta=np.ndarray>
    geopotential                                      (time, level, latitude, longitude) float32 dask.array<chunksize=(1, 37, 721, 1440), meta=np.ndarray>
    ...                                                ...
    total_precipitation                               (time, latitude, longitude) float32 dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>
    type_of_high_vegetation                           (latitude, longitude) float32 dask.array<chunksize=(721, 1440), meta=np.ndarray>
    type_of_low_vegetation                            (latitude, longitude) float32 dask.array<chunksize=(721, 1440), meta=np.ndarray>
    u_component_of_wind                               (time, level, latitude, longitude) float32 dask.array<chunksize=(1, 37, 721, 1440), meta=np.ndarray>
    v_component_of_wind                               (time, level, latitude, longitude) float32 dask.array<chunksize=(1, 37, 721, 1440), meta=np.ndarray>
    vertical_velocity                                 (time, level, latitude, longitude) float32 dask.array<chunksize=(1, 37, 721, 1440), meta=np.ndarray>

geop = reanalysis.geopotential
temp = reanalysis.temperature

geop

<xarray.DataArray 'geopotential' (time: 552264, level: 37, latitude: 721,
                                  longitude: 1440)>
dask.array<open_dataset-geopotential, shape=(552264, 37, 721, 1440), dtype=float32, chunksize=(1, 37, 721, 1440), chunktype=numpy.ndarray>
Coordinates:
  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
  * level      (level) int64 1 2 3 5 7 10 20 30 ... 850 875 900 925 950 975 1000
  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.0 359.2 359.5 359.8
  * time       (time) datetime64[ns] 1959-01-01 ... 2021-12-31T23:00:00
Attributes:
    long_name:      Geopotential
    short_name:     z
    standard_name:  geopotential
    units:          m**2 s**-2

Select time and level ranges from the dataset.

nHours = 6 
startTime = datetime(1993,3,13,18)
endTime = startTime + timedelta(hours=nHours)

RESTRICT YOUR TIME/LEVEL RANGES!
The full dataset is over 70 TB ... if you try to load in too large a range (or forget to specify the start/stop points), you will likely run out of system memory!

geopSub1 = geop.sel(time=slice(startTime, endTime), level=500)
tempSub1 = temp.sel(time=slice(startTime, endTime), level=850)

Convert to dam and deg C.

Z = mpcalc.geopotential_to_height(geopSub1).metpy.convert_units('dam')

T = tempSub1.metpy.convert_units('degC')

<xarray.DataArray 'truediv-2f0b3ac192d450f946fff895443b747d' (time: 7,
                                                              latitude: 721,
                                                              longitude: 1440)>
<Quantity(dask.array<mul, shape=(7, 721, 1440), dtype=float32, chunksize=(1, 721, 1440), chunktype=numpy.ndarray>, 'decameter')>
Coordinates:
  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
    level      int64 500
  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.0 359.2 359.5 359.8
  * time       (time) datetime64[ns] 1993-03-13T18:00:00 ... 1993-03-14

<xarray.DataArray 'temperature' (time: 7, latitude: 721, longitude: 1440)>
<Quantity(dask.array<truediv, shape=(7, 721, 1440), dtype=float32, chunksize=(1, 721, 1440), chunktype=numpy.ndarray>, 'degree_Celsius')>
Coordinates:
  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
    level      int64 850
  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.0 359.2 359.5 359.8
  * time       (time) datetime64[ns] 1993-03-13T18:00:00 ... 1993-03-14
Attributes:
    long_name:      Temperature
    short_name:     t
    standard_name:  air_temperature