Zarr multi-scales with pyramid_resample

This notebook demonstrates how to use pyramid_resample to generate pyramids using pyresample’s block based interpolators.

import dask  # noqa
import datatree
import xarray as xr

from ndpyramid import pyramid_resample

Setup output paths

You can set the output path to object storage or a memory store, after changing the S3 URI to a location that you have write access to.

S3 = False
if S3:
    output = "s3://carbonplan-scratch/pyramid_comparison/resampled.zarr"
else:
    import zarr

    output = zarr.storage.MemoryStore()

Open input dataset

In this example, we’ll use some of CarbonPlan’s downscaled CMIP6 data for the input dataset. It’s helpful to use chunked Zarr data as the input because we won’t need to manually rechunk the data, which can lead to a complicated Dask task graph if the rechunked data aren’t cached before pyramid generation.

ds = xr.open_zarr(
    'https://rice1.osn.mghpcc.org/carbonplan/cp-cmip/version1/rechunked_data/GARD-SV/ScenarioMIP.CCCma.CanESM5.ssp245.r1i1p1f1.month.GARD-SV.tasmax.zarr'
)
ds = ds.rio.write_crs("EPSG:4326")
ds

<xarray.Dataset> Size: 4GB
Dimensions:      (lat: 721, lon: 1440, time: 1020)
Coordinates:
  * lat          (lat) float32 3kB -90.0 -89.75 -89.5 -89.25 ... 89.5 89.75 90.0
  * lon          (lon) float32 6kB -180.0 -179.8 -179.5 ... 179.2 179.5 179.8
  * time         (time) datetime64[ns] 8kB 2015-01-01 2015-02-01 ... 2099-12-01
    spatial_ref  int64 8B 0
Data variables:
    tasmax       (time, lat, lon) float32 4GB dask.array<chunksize=(780, 72, 144), meta=np.ndarray>
Attributes: (12/17)
    Conventions:                     CF-1.8
    activity_id:                     ScenarioMIP
    cmip6_downscaling_contact:       hello@carbonplan.org
    cmip6_downscaling_explainer:     https://carbonplan.org/research/cmip6-do...
    cmip6_downscaling_institution:   CarbonPlan
    cmip6_downscaling_license:       CC-BY-4.0
    ...                              ...
    institution_id:                  CCCma
    member_id:                       r1i1p1f1
    references:                      Eyring, V., Bony, S., Meehl, G. A., Seni...
    source_id:                       CanESM5
    timescale:                       day
    variable_id:                     tasmax

xarray.Dataset

• Dimensions:
  • lat: 721
  • lon: 1440
  • time: 1020
• Coordinates: (4)
  • lat
    (lat)
    float32
    -90.0 -89.75 -89.5 ... 89.75 90.0

    long_name :

    latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([-90.  , -89.75, -89.5 , ...,  89.5 ,  89.75,  90.  ], dtype=float32)

  • lon
    (lon)
    float32
    -180.0 -179.8 ... 179.5 179.8

    long_name :

    longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([-180.  , -179.75, -179.5 , ...,  179.25,  179.5 ,  179.75],
          dtype=float32)

  • time
    (time)
    datetime64[ns]
    2015-01-01 ... 2099-12-01

    array(['2015-01-01T00:00:00.000000000', '2015-02-01T00:00:00.000000000',
           '2015-03-01T00:00:00.000000000', ..., '2099-10-01T00:00:00.000000000',
           '2099-11-01T00:00:00.000000000', '2099-12-01T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    array(0)

• Data variables: (1)
  • tasmax
    (time, lat, lon)
    float32
    dask.array<chunksize=(780, 72, 144), meta=np.ndarray>

    cell_measures :

    area: areacella

    cell_methods :

    area: mean time: maximum

    comment :

    maximum near-surface (usually, 2 meter) air temperature (add cell_method attribute 'time: max')

    history :

    2019-05-02T08:38:27Z altered by CMOR: Treated scalar dimension: 'height'. 2019-05-02T08:38:27Z altered by CMOR: Reordered dimensions, original order: lat lon time. 2019-05-02T08:38:27Z altered by CMOR: replaced missing value flag (1e+38) with standard missing value (1e+20).

    long_name :

    Daily Maximum Near-Surface Air Temperature

    original_name :

    STMX

    standard_name :

    air_temperature

    units :

    K

    Array Chunk Bytes 3.95 GiB 30.85 MiB Shape (1020, 721, 1440) (780, 72, 144) Dask graph 220 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([ -90.0, -89.75,  -89.5, -89.25,  -89.0, -88.75,  -88.5, -88.25,  -88.0,
           -87.75,
           ...
            87.75,   88.0,  88.25,   88.5,  88.75,   89.0,  89.25,   89.5,  89.75,
             90.0],
          dtype='float32', name='lat', length=721))

  • lon
    PandasIndex

    PandasIndex(Index([ -180.0, -179.75,  -179.5, -179.25,  -179.0, -178.75,  -178.5, -178.25,
            -178.0, -177.75,
           ...
             177.5,  177.75,   178.0,  178.25,   178.5,  178.75,   179.0,  179.25,
             179.5,  179.75],
          dtype='float32', name='lon', length=1440))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2015-01-01', '2015-02-01', '2015-03-01', '2015-04-01',
                   '2015-05-01', '2015-06-01', '2015-07-01', '2015-08-01',
                   '2015-09-01', '2015-10-01',
                   ...
                   '2099-03-01', '2099-04-01', '2099-05-01', '2099-06-01',
                   '2099-07-01', '2099-08-01', '2099-09-01', '2099-10-01',
                   '2099-11-01', '2099-12-01'],
                  dtype='datetime64[ns]', name='time', length=1020, freq=None))

• Attributes: (17)

  Conventions :

  CF-1.8

  activity_id :

  ScenarioMIP

  cmip6_downscaling_contact :

  hello@carbonplan.org

  cmip6_downscaling_explainer :

  https://carbonplan.org/research/cmip6-downscaling-explainer

  cmip6_downscaling_institution :

  CarbonPlan

  cmip6_downscaling_license :

  CC-BY-4.0

  cmip6_downscaling_license_url :

  https://creativecommons.org/licenses/by/4.0/

  cmip6_downscaling_method :

  GARD-SV

  cmip6_downscaling_repository :

  https://github.com/carbonplan/cmip6-downscaling

  cmip6_downscaling_web_map_tool :

  https://carbonplan.org/research/cmip6-downscaling

  experiment_id :

  ssp245

  institution_id :

  CCCma

  member_id :

  r1i1p1f1

  references :

  Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937-1958, doi:10.5194/gmd-9-1937-2016, 2016., Hersbach, H, Bell, B, Berrisford, P, et al. The ERA5 global reanalysis. Q J R Meteorol Soc. 2020; 146: 1999– 2049. https://doi.org/10.1002/qj.3803, Gutmann, E.D., J. Hamman, M.P. Clark, T. Eidhammer, A.W. Wood, and J.R. Arnold, 2021: A statistical downscaling framework for the production and testing of large ensembles of climate projections. In revision.

  source_id :

  CanESM5

  timescale :

  day

  variable_id :

  tasmax

Generate pyramids

We’ll use the pyramid_resample function to generate a pyramid with two levels. This function currently requires specifying which dimensions correspond to x and y. The data also must have longitude values between -180 and 180 and a dimension order of (time, y, x) for 3D datasets or (y, x) for 2D datasets. We’ll use the nearest and bilinear resampling algorithms for comparison with the other pyramid generation methods.

First, use the nearest neighbor method.

%%time
levels = 2
resampled_pyramid = pyramid_resample(ds, x="lon", y="lat", levels=levels, resampling="nearest")
resampled_pyramid.to_zarr(output, consolidated=True, mode='w')

CPU times: user 21.5 s, sys: 6.88 s, total: 28.4 s
Wall time: 24.2 s

Next, use the bilinear method.

%%time
levels = 2
resampled_pyramid = pyramid_resample(ds, x="lon", y="lat", levels=levels, resampling="bilinear")
resampled_pyramid.to_zarr(output, consolidated=True, mode='w')

CPU times: user 23.7 s, sys: 6.41 s, total: 30.1 s
Wall time: 21.5 s

Open and plot the result

dt = datatree.open_datatree(output, engine="zarr")
dt["1"].ds.isel(time=0).tasmax.plot()

<matplotlib.collections.QuadMesh at 0x7fe28fce15d0>