!pip install pyearthtools

Loading a Remote Zarr file - ERA5 ARCO example

Zarr is an increasingly popular format for storing gridded environmental science data, especially when working on cloud environments. The nature of a zarr, with consolidated metadat which enables interaction with a dataset as a unified whole, data divided into separate files called chunks to enable parallel read and write, which suits a massively parallel distributed file system like the object stores used by cloud providers, makes it efficient and intuive for many machine-learning platforms and workflows.

PyEarthTools provides a data accessor class for zarr files. As PyEarthTools uses xarray for data I/O, which through fsspec supports loading data directly from remote sources, it can load data from such sources.

Set up imports for demo

import xarray

import matplotlib.pyplot
import cartopy

import pyearthtools.data

import pyearthtools.pipeline

Set up the accessor

Fore this demo we will be using the copy of the ECMWF ERA5 Reanalsyis dataset hosted by Google Research on the Google Cloud Storage system.

remote_url = 'gs://gcp-public-data-arco-era5/ar/full_37-1h-0p25deg-chunk-1.zarr-v3'

zarr_opts = {'chunks': None, 'storage_options': {'token': 'anon'}}

Usually we are only interested in a subset of the data, so we can select a particular time, pressure levels and set of variables that we want to us

select_variables = ['temperature', 'specific_humidity']
select_pressure_levels = [200, 500, 700, 850, 1000]
select_time = ('2024-01-01','2025-01-01')
select_region = {'latitude': (-40,-5),
                 'longitude': (110,155)
                }

select_dict = {
    'level': select_pressure_levels,
    'time': slice(*select_time),
}

era5_accessor = pyearthtools.data.archive.ZarrTimeIndex(
       remote_url,
       variables=select_variables,
       open_kwargs=zarr_opts,
    remote=True,
    transforms=pyearthtools.data.transforms.TransformCollection([
        pyearthtools.data.transform.coordinates.Select(select_dict),
        pyearthtools.data.transform.region.Bounding(select_region['latitude'][0],select_region['latitude'][1], select_region['longitude'][0], select_region['longitude'][1])
    ]),
    )

Using the data of interest

We can now look at the xarray dataset through the time index function. If we do this, the xarray lazy loading paradigm means only the metadata will have been downloaded form the remote location, until we try to do some thing which requires the data. For example after exmaing the metadata, we will then plot some of the data which will take longer as the data will be transfered from the remote location. This will be quite efficient as only the chunks pertaining to the data selected (variables, time period, spatial region) will be transferred.

era5_accessor['2024-01-23 17:00']

<xarray.Dataset> Size: 1MB
Dimensions:            (time: 1, level: 5, latitude: 141, longitude: 181)
Coordinates:
  * time               (time) datetime64[ns] 8B 2024-01-23T17:00:00
  * level              (level) int64 40B 200 500 700 850 1000
  * latitude           (latitude) float32 564B -5.0 -5.25 -5.5 ... -39.75 -40.0
  * longitude          (longitude) float32 724B 110.0 110.2 ... 154.8 155.0
Data variables:
    specific_humidity  (time, level, latitude, longitude) float32 510kB ...
    temperature        (time, level, latitude, longitude) float32 510kB ...
Attributes:
    valid_time_start:       1940-01-01
    last_updated:           2026-06-01 03:27:53.232444+00:00
    valid_time_stop:        2025-12-31
    valid_time_stop_era5t:  2026-05-26

xarray.Dataset

• Dimensions:
  • time: 1
  • level: 5
  • latitude: 141
  • longitude: 181
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2024-01-23T17:00:00

    array(['2024-01-23T17:00:00.000000000'], dtype='datetime64[ns]')

  • level
    (level)
    int64
    200 500 700 850 1000

    long_name :

    level

    units :

    Hectopascal(hPa)

    array([ 200,  500,  700,  850, 1000])

  • latitude
    (latitude)
    float32
    -5.0 -5.25 -5.5 ... -39.75 -40.0

    long_name :

    latitude

    units :

    degrees_north

    array([ -5.  ,  -5.25,  -5.5 ,  -5.75,  -6.  ,  -6.25,  -6.5 ,  -6.75,  -7.  ,
            -7.25,  -7.5 ,  -7.75,  -8.  ,  -8.25,  -8.5 ,  -8.75,  -9.  ,  -9.25,
            -9.5 ,  -9.75, -10.  , -10.25, -10.5 , -10.75, -11.  , -11.25, -11.5 ,
           -11.75, -12.  , -12.25, -12.5 , -12.75, -13.  , -13.25, -13.5 , -13.75,
           -14.  , -14.25, -14.5 , -14.75, -15.  , -15.25, -15.5 , -15.75, -16.  ,
           -16.25, -16.5 , -16.75, -17.  , -17.25, -17.5 , -17.75, -18.  , -18.25,
           -18.5 , -18.75, -19.  , -19.25, -19.5 , -19.75, -20.  , -20.25, -20.5 ,
           -20.75, -21.  , -21.25, -21.5 , -21.75, -22.  , -22.25, -22.5 , -22.75,
           -23.  , -23.25, -23.5 , -23.75, -24.  , -24.25, -24.5 , -24.75, -25.  ,
           -25.25, -25.5 , -25.75, -26.  , -26.25, -26.5 , -26.75, -27.  , -27.25,
           -27.5 , -27.75, -28.  , -28.25, -28.5 , -28.75, -29.  , -29.25, -29.5 ,
           -29.75, -30.  , -30.25, -30.5 , -30.75, -31.  , -31.25, -31.5 , -31.75,
           -32.  , -32.25, -32.5 , -32.75, -33.  , -33.25, -33.5 , -33.75, -34.  ,
           -34.25, -34.5 , -34.75, -35.  , -35.25, -35.5 , -35.75, -36.  , -36.25,
           -36.5 , -36.75, -37.  , -37.25, -37.5 , -37.75, -38.  , -38.25, -38.5 ,
           -38.75, -39.  , -39.25, -39.5 , -39.75, -40.  ], dtype=float32)

  • longitude
    (longitude)
    float32
    110.0 110.2 110.5 ... 154.8 155.0

    long_name :

    longitude

    units :

    degrees_east

    array([110.  , 110.25, 110.5 , 110.75, 111.  , 111.25, 111.5 , 111.75, 112.  ,
           112.25, 112.5 , 112.75, 113.  , 113.25, 113.5 , 113.75, 114.  , 114.25,
           114.5 , 114.75, 115.  , 115.25, 115.5 , 115.75, 116.  , 116.25, 116.5 ,
           116.75, 117.  , 117.25, 117.5 , 117.75, 118.  , 118.25, 118.5 , 118.75,
           119.  , 119.25, 119.5 , 119.75, 120.  , 120.25, 120.5 , 120.75, 121.  ,
           121.25, 121.5 , 121.75, 122.  , 122.25, 122.5 , 122.75, 123.  , 123.25,
           123.5 , 123.75, 124.  , 124.25, 124.5 , 124.75, 125.  , 125.25, 125.5 ,
           125.75, 126.  , 126.25, 126.5 , 126.75, 127.  , 127.25, 127.5 , 127.75,
           128.  , 128.25, 128.5 , 128.75, 129.  , 129.25, 129.5 , 129.75, 130.  ,
           130.25, 130.5 , 130.75, 131.  , 131.25, 131.5 , 131.75, 132.  , 132.25,
           132.5 , 132.75, 133.  , 133.25, 133.5 , 133.75, 134.  , 134.25, 134.5 ,
           134.75, 135.  , 135.25, 135.5 , 135.75, 136.  , 136.25, 136.5 , 136.75,
           137.  , 137.25, 137.5 , 137.75, 138.  , 138.25, 138.5 , 138.75, 139.  ,
           139.25, 139.5 , 139.75, 140.  , 140.25, 140.5 , 140.75, 141.  , 141.25,
           141.5 , 141.75, 142.  , 142.25, 142.5 , 142.75, 143.  , 143.25, 143.5 ,
           143.75, 144.  , 144.25, 144.5 , 144.75, 145.  , 145.25, 145.5 , 145.75,
           146.  , 146.25, 146.5 , 146.75, 147.  , 147.25, 147.5 , 147.75, 148.  ,
           148.25, 148.5 , 148.75, 149.  , 149.25, 149.5 , 149.75, 150.  , 150.25,
           150.5 , 150.75, 151.  , 151.25, 151.5 , 151.75, 152.  , 152.25, 152.5 ,
           152.75, 153.  , 153.25, 153.5 , 153.75, 154.  , 154.25, 154.5 , 154.75,
           155.  ], dtype=float32)

• Data variables: (2)
  • specific_humidity
    (time, level, latitude, longitude)
    float32
    ...

    long_name :

    Specific humidity

    short_name :

    q

    standard_name :

    specific_humidity

    units :

    kg kg**-1

    [127605 values with dtype=float32]

  • temperature
    (time, level, latitude, longitude)
    float32
    ...

    long_name :

    Temperature

    short_name :

    t

    standard_name :

    air_temperature

    units :

    K

    [127605 values with dtype=float32]

• Attributes: (4)

  valid_time_start :

  1940-01-01

  last_updated :

  2026-06-01 03:27:53.232444+00:00

  valid_time_stop :

  2025-12-31

  valid_time_stop_era5t :

  2026-05-26

fig1 = matplotlib.pyplot.figure(figsize=(16,10))
ax1=fig1.add_subplot(1,1,1,projection=cartopy.crs.PlateCarree())
era5_accessor['2024-01-23 17:00']['temperature'].sel(level=850)[0].plot.contourf(ax=ax1)
ax1.coastlines(color='w')

<cartopy.mpl.feature_artist.FeatureArtist at 0x7f105bf97590>