Zarr in Cloud Object Storage

In this tutorial, we’ll cover the following:

• Finding a cloud hosted Zarr archive of CMIP6 dataset(s)

• Remote data access to a single CMIP6 dataset (sea surface height)

• Calculate future predicted sea level change in 2100 compared to 2015

import gcsfs
import pandas as pd
import xarray as xr

Finding cloud native data

Cloud-native data means data that is structured for efficient querying across the network. Typically, this means having metadata that describes the entire file in the header of the file, or having a a separate pointer file (so that there is no need to download everything first).

Quite commonly, you’ll see cloud-native datasets stored on these three object storage providers, though there are many other ones too.

• Amazon Simple Storage Service (S3)

• Azure Blob Storage

• Google Cloud Storage

Getting cloud hosted CMIP6 data

The Coupled Model Intercomparison Project Phase 6 (CMIP6) dataset is a rich archive of modelling experiments carried out to predict the climate change impacts. The datasets are stored using the Zarr format, and we’ll go over how to access it.

Sources:

• https://esgf-node.llnl.gov/search/cmip6/

• CMIP6 data hosted on Google Cloud - https://console.cloud.google.com/marketplace/details/noaa-public/cmip6

• Pangeo/ESGF Cloud Data Access tutorial - https://pangeo-data.github.io/pangeo-cmip6-cloud/accessing_data.html

First, let’s open a CSV containing the list of CMIP6 datasets available

df = pd.read_csv("https://cmip6.storage.googleapis.com/pangeo-cmip6.csv")
print(f"Number of rows: {len(df)}")
df.head()

Number of rows: 514818

	activity_id	institution_id	source_id	experiment_id	member_id	table_id	variable_id	grid_label	zstore	dcpp_init_year	version
0	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	ps	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
1	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	rsds	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
2	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	rlus	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
3	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	rlds	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706
4	HighResMIP	CMCC	CMCC-CM2-HR4	highresSST-present	r1i1p1f1	Amon	psl	gn	gs://cmip6/CMIP6/HighResMIP/CMCC/CMCC-CM2-HR4/...	NaN	20170706

Over 5 million rows! Let’s filter it down to the variable and experiment we’re interested in, e.g. sea surface height.

For the variable_id, you can look it up given some keyword at https://docs.google.com/spreadsheets/d/1UUtoz6Ofyjlpx5LdqhKcwHFz2SGoTQV2_yekHyMfL9Y

For the experiment_id, download the spreadsheet from ES-DOC/esdoc-docs, go to the ‘experiment’ tab, and find the one you’re interested in.

Another good place to find the right model runs is https://esgf-node.llnl.gov/search/cmip6 (once you get your head around the acronyms and short names).

Below, we’ll filter to CMIP6 experiments matching:

• Sea Surface Height Above Geoid [m] (variable_id: zos)

• Shared Socioeconomic Pathway 5 (experiment_id: ssp585)

df_zos = df.query("variable_id == 'zos' & experiment_id == 'ssp585'")
df_zos

	activity_id	institution_id	source_id	experiment_id	member_id	table_id	variable_id	grid_label	zstore	dcpp_init_year	version
12081	ScenarioMIP	NOAA-GFDL	GFDL-ESM4	ssp585	r1i1p1f1	Omon	zos	gn	gs://cmip6/CMIP6/ScenarioMIP/NOAA-GFDL/GFDL-ES...	NaN	20180701
12082	ScenarioMIP	NOAA-GFDL	GFDL-ESM4	ssp585	r1i1p1f1	Omon	zos	gr	gs://cmip6/CMIP6/ScenarioMIP/NOAA-GFDL/GFDL-ES...	NaN	20180701
45176	ScenarioMIP	AWI	AWI-CM-1-1-MR	ssp585	r1i1p1f1	Omon	zos	gn	gs://cmip6/CMIP6/ScenarioMIP/AWI/AWI-CM-1-1-MR...	NaN	20181218
54265	ScenarioMIP	CNRM-CERFACS	CNRM-CM6-1	ssp585	r1i1p1f2	Omon	zos	gn	gs://cmip6/CMIP6/ScenarioMIP/CNRM-CERFACS/CNRM...	NaN	20190219
69190	ScenarioMIP	CNRM-CERFACS	CNRM-CM6-1	ssp585	r4i1p1f2	Omon	zos	gn	gs://cmip6/CMIP6/ScenarioMIP/CNRM-CERFACS/CNRM...	NaN	20190410
...	...	...	...	...	...	...	...	...	...	...	...
497517	ScenarioMIP	MIROC	MIROC-ES2L	ssp585	r2i1p1f2	Omon	zos	gr1	gs://cmip6/CMIP6/ScenarioMIP/MIROC/MIROC-ES2L/...	NaN	20201222
497957	ScenarioMIP	MIROC	MIROC-ES2L	ssp585	r10i1p1f2	Omon	zos	gn	gs://cmip6/CMIP6/ScenarioMIP/MIROC/MIROC-ES2L/...	NaN	20201222
497958	ScenarioMIP	MIROC	MIROC-ES2L	ssp585	r10i1p1f2	Omon	zos	gr1	gs://cmip6/CMIP6/ScenarioMIP/MIROC/MIROC-ES2L/...	NaN	20201222
502853	ScenarioMIP	EC-Earth-Consortium	EC-Earth3-CC	ssp585	r1i1p1f1	Omon	zos	gn	gs://cmip6/CMIP6/ScenarioMIP/EC-Earth-Consorti...	NaN	20210113
505073	ScenarioMIP	CMCC	CMCC-ESM2	ssp585	r1i1p1f1	Omon	zos	gn	gs://cmip6/CMIP6/ScenarioMIP/CMCC/CMCC-ESM2/ss...	NaN	20210126

272 rows × 11 columns

There’s 272 modelled scenarios for SSP5. Let’s just get the URL to the first one in the list for now.

print(df_zos.zstore.iloc[0])

gs://cmip6/CMIP6/ScenarioMIP/NOAA-GFDL/GFDL-ESM4/ssp585/r1i1p1f1/Omon/zos/gn/v20180701/

Reading from the remote Zarr storage

In many cases, you’ll need to first connect to the cloud provider. The CMIP6 dataset allows anonymous access, but for some cases, you may need to authenticate.

fs = gcsfs.GCSFileSystem(token="anon")

Next, we’ll need a mapping to the Google Storage object. This can be done using fs.get_mapper.

A more generic way (for other cloud providers) is to use fsspec.get_mapper instead.

store = fs.get_mapper(
    "gs://cmip6/CMIP6/ScenarioMIP/NOAA-GFDL/GFDL-ESM4/ssp585/r1i1p1f1/Omon/zos/gn/v20180701/"
)

With that, we can open the Zarr store like so.

ds = xr.open_zarr(store=store, consolidated=True)
ds

<xarray.Dataset> Size: 2GB
Dimensions:    (bnds: 2, y: 576, x: 720, vertex: 4, time: 1032)
Coordinates:
  * bnds       (bnds) float64 16B 1.0 2.0
    lat        (y, x) float32 2MB dask.array<chunksize=(576, 720), meta=np.ndarray>
    lat_bnds   (y, x, vertex) float32 7MB dask.array<chunksize=(576, 720, 4), meta=np.ndarray>
    lon        (y, x) float32 2MB dask.array<chunksize=(576, 720), meta=np.ndarray>
    lon_bnds   (y, x, vertex) float32 7MB dask.array<chunksize=(576, 720, 4), meta=np.ndarray>
  * time       (time) object 8kB 2015-01-16 12:00:00 ... 2100-12-16 12:00:00
    time_bnds  (time, bnds) object 17kB dask.array<chunksize=(1032, 2), meta=np.ndarray>
  * x          (x) float64 6kB -299.8 -299.2 -298.8 -298.2 ... 58.75 59.25 59.75
  * y          (y) float64 5kB -77.91 -77.72 -77.54 -77.36 ... 89.47 89.68 89.89
Dimensions without coordinates: vertex
Data variables:
    zos        (time, y, x) float32 2GB dask.array<chunksize=(61, 576, 720), meta=np.ndarray>
Attributes: (12/49)
    Conventions:            CF-1.7 CMIP-6.0 UGRID-1.0
    activity_id:            ScenarioMIP
    branch_method:          standard
    branch_time_in_child:   60225.0
    branch_time_in_parent:  60225.0
    comment:                <null ref>
    ...                     ...
    tracking_id:            hdl:21.14100/328fd441-31bb-474f-8fd1-7e3420b72524...
    variable_id:            zos
    variant_info:           N/A
    variant_label:          r1i1p1f1
    netcdf_tracking_ids:    hdl:21.14100/328fd441-31bb-474f-8fd1-7e3420b72524...
    version_id:             v20180701

xarray.Dataset

• Dimensions:
  • bnds: 2
  • y: 576
  • x: 720
  • vertex: 4
  • time: 1032
• Coordinates: (9)
  • bnds
    (bnds)
    float64
    1.0 2.0

    long_name :

    vertex number

    array([1., 2.])

  • lat
    (y, x)
    float32
    dask.array<chunksize=(576, 720), meta=np.ndarray>

    bounds :

    lat_bnds

    cell_methods :

    y: point x: point

    long_name :

    latitude

    standard_name :

    latitude

    units :

    degrees_north

    Array Chunk Bytes 1.58 MiB 1.58 MiB Shape (576, 720) (576, 720) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

  • lat_bnds
    (y, x, vertex)
    float32
    dask.array<chunksize=(576, 720, 4), meta=np.ndarray>

    axis :

    Y

    long_name :

    latitude bounds

    units :

    degrees_north

    Array Chunk Bytes 6.33 MiB 6.33 MiB Shape (576, 720, 4) (576, 720, 4) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

  • lon
    (y, x)
    float32
    dask.array<chunksize=(576, 720), meta=np.ndarray>

    bounds :

    lon_bnds

    cell_methods :

    y: point x: point

    long_name :

    longitude

    standard_name :

    longitude

    units :

    degrees_east

    Array Chunk Bytes 1.58 MiB 1.58 MiB Shape (576, 720) (576, 720) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

  • lon_bnds
    (y, x, vertex)
    float32
    dask.array<chunksize=(576, 720, 4), meta=np.ndarray>

    axis :

    X

    long_name :

    longitude bounds

    units :

    degrees_east

    Array Chunk Bytes 6.33 MiB 6.33 MiB Shape (576, 720, 4) (576, 720, 4) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

  • time
    (time)
    object
    2015-01-16 12:00:00 ... 2100-12-...

    axis :

    T

    bounds :

    time_bnds

    calendar_type :

    noleap

    description :

    Temporal mean

    long_name :

    time

    standard_name :

    time

    array([cftime.DatetimeNoLeap(2015, 1, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2015, 2, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2015, 3, 16, 12, 0, 0, 0, has_year_zero=True),
           ...,
           cftime.DatetimeNoLeap(2100, 10, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2100, 11, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2100, 12, 16, 12, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • time_bnds
    (time, bnds)
    object
    dask.array<chunksize=(1032, 2), meta=np.ndarray>

    long_name :

    time axis boundaries

    Array Chunk Bytes 16.12 kiB 16.12 kiB Shape (1032, 2) (1032, 2) Dask graph 1 chunks in 2 graph layers Data type object numpy.ndarray

  • x
    (x)
    float64
    -299.8 -299.2 ... 59.25 59.75

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    degrees

    array([-299.75, -299.25, -298.75, ...,   58.75,   59.25,   59.75])

  • y
    (y)
    float64
    -77.91 -77.72 ... 89.68 89.89

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    degrees

    array([-77.907938, -77.723813, -77.539688, ...,  89.472   ,  89.6832  ,
            89.8944  ])

• Data variables: (1)
  • zos
    (time, y, x)
    float32
    dask.array<chunksize=(61, 576, 720), meta=np.ndarray>

    cell_measures :

    area: areacello

    cell_methods :

    area: mean where sea time: mean

    long_name :

    Sea Surface Height Above Geoid

    original_name :

    zos

    standard_name :

    sea_surface_height_above_geoid

    units :

    m

    Array Chunk Bytes 1.59 GiB 96.50 MiB Shape (1032, 576, 720) (61, 576, 720) Dask graph 17 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (4)
  • bnds
    PandasIndex

    PandasIndex(Index([1.0, 2.0], dtype='float64', name='bnds'))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([2015-01-16 12:00:00, 2015-02-15 00:00:00, 2015-03-16 12:00:00,
                 2015-04-16 00:00:00, 2015-05-16 12:00:00, 2015-06-16 00:00:00,
                 2015-07-16 12:00:00, 2015-08-16 12:00:00, 2015-09-16 00:00:00,
                 2015-10-16 12:00:00,
                 ...
                 2100-03-16 12:00:00, 2100-04-16 00:00:00, 2100-05-16 12:00:00,
                 2100-06-16 00:00:00, 2100-07-16 12:00:00, 2100-08-16 12:00:00,
                 2100-09-16 00:00:00, 2100-10-16 12:00:00, 2100-11-16 00:00:00,
                 2100-12-16 12:00:00],
                dtype='object', length=1032, calendar='noleap', freq=None))

  • x
    PandasIndex

    PandasIndex(Index([-299.75, -299.25, -298.75, -298.25, -297.75, -297.25, -296.75, -296.25,
           -295.75, -295.25,
           ...
             55.25,   55.75,   56.25,   56.75,   57.25,   57.75,   58.25,   58.75,
             59.25,   59.75],
          dtype='float64', name='x', length=720))

  • y
    PandasIndex

    PandasIndex(Index([-77.90793753487047,  -77.7238126046114, -77.53968767435234,
           -77.35556274409328,  -77.1714378138342, -76.98731288357514,
           -76.80318795331607, -76.61906302305701, -76.43493809279795,
           -76.25081316253888,
           ...
            87.99359944564549,   88.2047995039986,   88.4159995623517,
            88.62719962070481,  88.83839967905793,  89.04959973741103,
            89.26079979576413,  89.47199985411724,  89.68319991247033,
            89.89439997082347],
          dtype='float64', name='y', length=576))

• Attributes: (49)

  Conventions :

  CF-1.7 CMIP-6.0 UGRID-1.0

  activity_id :

  ScenarioMIP

  branch_method :

  standard

  branch_time_in_child :

  60225.0

  branch_time_in_parent :

  60225.0

  comment :

  <null ref>

  contact :

  gfdl.climate.model.info@noaa.gov

  creation_date :

  2019-06-19T11:11:02Z

  data_specs_version :

  01.00.27

  experiment :

  update of RCP8.5 based on SSP5

  experiment_id :

  ssp585

  external_variables :

  areacello

  forcing_index :

  1

  frequency :

  mon

  further_info_url :

  https://furtherinfo.es-doc.org/CMIP6.NOAA-GFDL.GFDL-ESM4.ssp585.none.r1i1p1f1

  grid :

  ocean data on native grid tripolar - nominal 0.5 deg latitude/longitude

  grid_label :

  gn

  history :

  File was processed by fremetar (GFDL analog of CMOR). TripleID: [exper_id_FlJGh4Wo6W,realiz_id_kt2pvOSbWt,run_id_1S546GMKbs]

  initialization_index :

  1

  institution :

  National Oceanic and Atmospheric Administration, Geophysical Fluid Dynamics Laboratory, Princeton, NJ 08540, USA

  institution_id :

  NOAA-GFDL

  license :

  CMIP6 model data produced by NOAA-GFDL is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License (https://creativecommons.org/licenses/). Consult https://pcmdi.llnl.gov/CMIP6/TermsOfUse for terms of use governing CMIP6 output, including citation requirements and proper acknowledgment. Further information about this data, including some limitations, can be found via the further_info_url (recorded as a global attribute in this file). The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  mip_era :

  CMIP6

  nominal_resolution :

  50 km

  parent_activity_id :

  CMIP

  parent_experiment_id :

  historical

  parent_mip_era :

  CMIP6

  parent_source_id :

  GFDL-ESM4

  parent_time_units :

  days since 1850-1-1

  parent_variant_label :

  r1i1p1f1

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  realm :

  ocean

  references :

  see further_info_url attribute

  source :

  GFDL-ESM4 (2018): atmos: GFDL-AM4.1 (Cubed-sphere (c96) - 1 degree nominal horizontal resolution; 360 x 180 longitude/latitude; 49 levels; top level 1 Pa) ocean: GFDL-OM4p5 (GFDL-MOM6, tripolar - nominal 0.5 deg; 720 x 576 longitude/latitude; 75 levels; top grid cell 0-2 m) seaIce: GFDL-SIM4p5 (GFDL-SIS2.0, tripolar - nominal 0.5 deg; 720 x 576 longitude/latitude; 5 layers; 5 thickness categories) land: GFDL-LM4.1 aerosol: interactive atmosChem: GFDL-ATMCHEM4.1 (full atmospheric chemistry) ocnBgchem: GFDL-COBALTv2 landIce: GFDL-LM4.1 (GFDL ID: 2019_0301)

  source_id :

  GFDL-ESM4

  source_type :

  AOGCM AER CHEM BGC

  status :

  2020-10-03;created; by gcs.cmip6.ldeo@gmail.com

  sub_experiment :

  none

  sub_experiment_id :

  none

  table_id :

  Omon

  title :

  NOAA GFDL GFDL-ESM4 model output prepared for CMIP6 update of RCP8.5 based on SSP5

  tracking_id :

  hdl:21.14100/328fd441-31bb-474f-8fd1-7e3420b72524 hdl:21.14100/5060ff3e-ac5b-41a7-b8be-6e00d5aca243 hdl:21.14100/37918073-5c8f-4fc2-8b5b-0edd752c0637 hdl:21.14100/f51f7134-2e58-4a9c-84e9-5a742fa818c6 hdl:21.14100/a39ee598-7962-4e3a-8943-ac5a8b091bcd

  variable_id :

  zos

  variant_info :

  N/A

  variant_label :

  r1i1p1f1

  netcdf_tracking_ids :

  hdl:21.14100/328fd441-31bb-474f-8fd1-7e3420b72524 hdl:21.14100/5060ff3e-ac5b-41a7-b8be-6e00d5aca243 hdl:21.14100/37918073-5c8f-4fc2-8b5b-0edd752c0637 hdl:21.14100/f51f7134-2e58-4a9c-84e9-5a742fa818c6 hdl:21.14100/a39ee598-7962-4e3a-8943-ac5a8b091bcd

  version_id :

  v20180701

Selecting time slices

Let’s say we want to calculate sea level change between 2015 and 2100. We can access just the specific time points needed using xr.Dataset.sel.

zos_2015jan = ds.zos.sel(time="2015-01-16").squeeze()
zos_2100dec = ds.zos.sel(time="2100-12-16").squeeze()

Sea level change would just be 2100 minus 2015.

sealevelchange = zos_2100dec - zos_2015jan

Note that up to this point, we have not actually downloaded any (big) data yet from the cloud. This is all working based on metadata only.

To bring the data from the cloud to your local computer, call .compute. This will take a while depending on your connection speed.

sealevelchange = sealevelchange.compute()

We can do a quick plot to show how Sea Level is predicted to change between 2015-2100 (from one modelled experiment).

sealevelchange.plot.imshow()

<matplotlib.image.AxesImage at 0x7f6781db2ff0>

Notice the blue parts between -40 and -60 South where sea level has dropped? That’s to do with the Antarctic ice sheet losing mass and resulting in a lower gravitational pull, resulting in a relative decrease in sea level. Over most of the Northern Hemisphere though, sea level rise has increased between 2015 and 2100.

That’s all! Hopefully this will get you started on accessing more cloud-native datasets!