Zarr multi-scales with pyramid_reproject

This notebook demonstrates how to use pyramid_reproject to generate pyramids using rasterios’s reprojection functionality, which wraps GDAL.

import dask  # noqa
import datatree
import xarray as xr

from ndpyramid import pyramid_reproject

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.

# Open the dataset
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_reproject function to generate a pyramid with two levels. We’ll use the nearest and bilinear resampling algorithms for comparison with the other pyramid generation methods.

%%time
levels = 2
reprojected_pyramid = pyramid_reproject(ds, levels=levels, resampling="nearest")
reprojected_pyramid.to_zarr(output, consolidated=True, mode='w')

CPU times: user 40 s, sys: 19.2 s, total: 59.2 s
Wall time: 56.8 s

%%time
levels = 2
reprojected_pyramid = pyramid_reproject(ds, levels=levels, resampling="bilinear")
reprojected_pyramid.to_zarr(output, consolidated=True, mode='w')

CPU times: user 1min 28s, sys: 31.5 s, total: 2min
Wall time: 1min 57s

Open and plot the result

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

<matplotlib.collections.QuadMesh at 0x7fcfd0306500>