Rechunk the kerchunked dataset¶

In [1]:

import fsspec
import fsspec.implementations.reference
import zarr
import xarray as xr
from pathlib import Path

from rechunker import rechunk

In [21]:

import rechunker
rechunker.__version__

Out[21]:

'0.5.1'

Start a Dask Gateway cluster¶

Use a custom helper function ebd.start_dask_cluster to set options on this cluster. We don't have to use this helper, it just cuts down on lines of code in notebooks.

In [ ]:

import sys
import os
sys.path.append('/shared/users/rsignell/lib')
import ebdpy as ebd
os.environ['AWS_PROFILE'] = 'esip-qhub'  # use env vars for AWS credentials to write

client, cluster, gateway = ebd.start_dask_cluster(
    profile=os.environ['AWS_PROFILE'],
    worker_max=30,
    region='us-west-2', 
    use_existing_cluster=True,
    adaptive_scaling=False, 
    wait_for_cluster=True, 
    propagate_env=True)

Open Kerchunked Dataset JSON for one year¶

Open Entire Kerchunked Dataset, lazy loaded from parquet¶

In [3]:

s3_lazy_refs = 's3://esip-qhub-public/nwm/LDAS-1k/lazyrefs'

In [4]:

fs = fsspec.implementations.reference.DFReferenceFileSystem(s3_lazy_refs, lazy=True, target_options={"anon": True},
                                                            remote_protocol="s3", remote_options={"anon": True})
m = fs.get_mapper("")
ds = xr.open_dataset(m, engine="zarr", chunks={'time':1, 'y':3840, 'x':4608}, 
                     backend_kwargs=dict(consolidated=False))

In [5]:

ds

Out[5]:

<xarray.Dataset>
Dimensions:   (time: 116631, y: 3840, x: 4608, vis_nir: 2, soil_layers_stag: 4)
Coordinates:
  * time      (time) datetime64[ns] 1979-02-01T03:00:00 ... 2020-12-31T21:00:00
  * x         (x) float64 -2.303e+06 -2.302e+06 ... 2.303e+06 2.304e+06
  * y         (y) float64 -1.92e+06 -1.919e+06 ... 1.918e+06 1.919e+06
Dimensions without coordinates: vis_nir, soil_layers_stag
Data variables: (12/21)
    ACCET     (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    ACSNOM    (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    ALBEDO    (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    ALBSND    (time, y, vis_nir, x) float64 dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>
    ALBSNI    (time, y, vis_nir, x) float64 dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>
    COSZ      (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    ...        ...
    SNOWH     (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    SOIL_M    (time, y, soil_layers_stag, x) float64 dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>
    SOIL_W    (time, y, soil_layers_stag, x) float64 dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>
    TRAD      (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    UGDRNOFF  (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    crs       object ...
Attributes:
    Conventions:                CF-1.6
    GDAL_DataType:              Generic
    TITLE:                      OUTPUT FROM WRF-Hydro v5.2.0-beta2
    code_version:               v5.2.0-beta2
    model_configuration:        retrospective
    model_initialization_time:  1979-02-01_00:00:00
    model_output_type:          land
    model_output_valid_time:    1979-02-01_03:00:00
    model_total_valid_times:    472
    proj4:                      +proj=lcc +units=m +a=6370000.0 +b=6370000.0 ...

xarray.Dataset

Dimensions:
- time: 116631
- y: 3840
- x: 4608
- vis_nir: 2
- soil_layers_stag: 4
Coordinates: (3)
- time
  (time)
  datetime64[ns]
  1979-02-01T03:00:00 ... 2020-12-...
  long_name :
  valid output time
  standard_name :
  time
  valid_max :
  4862880
  valid_min :
  4778100
```
array(['1979-02-01T03:00:00.000000000', '1979-02-01T06:00:00.000000000',
       '1979-02-01T09:00:00.000000000', ..., '2020-12-31T15:00:00.000000000',
       '2020-12-31T18:00:00.000000000', '2020-12-31T21:00:00.000000000'],
      dtype='datetime64[ns]')
```
- x
  (x)
  float64
  -2.303e+06 -2.302e+06 ... 2.304e+06
  _CoordinateAxisType :
  GeoX
  long_name :
  x coordinate of projection
  resolution :
  1000.0
  standard_name :
  projection_x_coordinate
  units :
  m
```
array([-2303499.25, -2302499.25, -2301499.25, ...,  2301500.75,  2302500.75,
        2303500.75])
```
- y
  (y)
  float64
  -1.92e+06 -1.919e+06 ... 1.919e+06
  _CoordinateAxisType :
  GeoY
  long_name :
  y coordinate of projection
  resolution :
  1000.0
  standard_name :
  projection_y_coordinate
  units :
  m
```
array([-1919500.375, -1918500.375, -1917500.375, ...,  1917499.625,
        1918499.625,  1919499.625])
```

Data variables: (21)

ACCET

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Accumulated total ET
units :: mm
valid_range :: [-100000, 100000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

ACSNOM

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: accumulated melting water out of snow bottom
units :: mm
valid_range :: [0, 10000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

ALBEDO

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Surface albedo
units :: -
valid_range :: [0, 100]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

ALBSND

(time, y, vis_nir, x)

float64

dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: snowpack albedo, direct
units :: -
valid_range :: [0, 100]

	Array	Chunk
Bytes	30.03 TiB	135.00 MiB
Shape	(116631, 3840, 2, 4608)	(1, 3840, 1, 4608)
Dask graph	233262 chunks in 2 graph layers
Data type	float64 numpy.ndarray

ALBSNI

(time, y, vis_nir, x)

float64

dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: snowpack albedo, diffuse
units :: -
valid_range :: [0, 100]

	Array	Chunk
Bytes	30.03 TiB	135.00 MiB
Shape	(116631, 3840, 2, 4608)	(1, 3840, 1, 4608)
Dask graph	233262 chunks in 2 graph layers
Data type	float64 numpy.ndarray

COSZ

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Cosine of zenith angle
units :: -
valid_range :: [-100, 100]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

EDIR

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Direct from soil evaporation rate
units :: kg m-2 s-1
valid_range :: [-10000000, 10000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

FIRA

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Total net LW radiation to atmosphere
units :: W m-2
valid_range :: [-15000, 15000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

FSA

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Total absorbed SW radiation
units :: W m-2
valid_range :: [-15000, 15000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

FSNO

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Snow-cover fraction on the ground
units :: 1
valid_range :: [0, 1000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

HFX

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Total sensible heat to the atmosphere
units :: W m-2
valid_range :: [-15000, 15000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

LH

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Total latent heat to the atmosphere
units :: W m-2
valid_range :: [-15000, 15000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

QRAIN

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Rainfall rate on the ground
units :: mm s-1
valid_range :: [0, 10000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

QSNOW

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Snowfall rate on the ground
units :: mm s-1
valid_range :: [0, 10000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

SNEQV

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Snow water equivalent
units :: kg m-2
valid_range :: [0, 1000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

SNOWH

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Snow depth
units :: m
valid_range :: [0, 1000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

SOIL_M

(time, y, soil_layers_stag, x)

float64

dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: volumetric soil moisture, the dimensionless ratio of water volume (m3) to soil volume (m3)
units :: m3 m-3
valid_range :: [0, 100]

	Array	Chunk
Bytes	60.06 TiB	135.00 MiB
Shape	(116631, 3840, 4, 4608)	(1, 3840, 1, 4608)
Dask graph	466524 chunks in 2 graph layers
Data type	float64 numpy.ndarray

SOIL_W

(time, y, soil_layers_stag, x)

float64

dask.array<chunksize=(1, 3840, 1, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: liquid volumetric soil moisture
units :: m3 m-3
valid_range :: [0, 100]

	Array	Chunk
Bytes	60.06 TiB	135.00 MiB
Shape	(116631, 3840, 4, 4608)	(1, 3840, 1, 4608)
Dask graph	466524 chunks in 2 graph layers
Data type	float64 numpy.ndarray

TRAD

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Surface radiative temperature
units :: K
valid_range :: [0, 4000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

UGDRNOFF

(time, y, x)

float64

dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>

esri_pe_string :: PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
grid_mapping :: crs
long_name :: Accumulated underground runoff
units :: mm
valid_range :: [-10000, 10000000]

	Array	Chunk
Bytes	15.02 TiB	135.00 MiB
Shape	(116631, 3840, 4608)	(1, 3840, 4608)
Dask graph	116631 chunks in 2 graph layers
Data type	float64 numpy.ndarray

crs
()
object
...
GeoTransform :
-2303999.17655 1000.0 0 1919999.66329 0 -1000.0
_CoordinateAxes :
y x
_CoordinateTransformType :
Projection
earth_radius :
6370000.0
esri_pe_string :
PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
false_easting :
0.0
false_northing :
0.0
grid_mapping_name :
lambert_conformal_conic
inverse_flattening :
0.0
latitude_of_projection_origin :
40.0
long_name :
CRS definition
longitude_of_central_meridian :
-97.0
longitude_of_prime_meridian :
0.0
semi_major_axis :
6370000.0
spatial_ref :
PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DATUM["D_Sphere",SPHEROID["Sphere",6370000.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["false_easting",0.0],PARAMETER["false_northing",0.0],PARAMETER["central_meridian",-97.0],PARAMETER["standard_parallel_1",30.0],PARAMETER["standard_parallel_2",60.0],PARAMETER["latitude_of_origin",40.0],UNIT["Meter",1.0]];-35691800 -29075200 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision
standard_parallel :
[30.0, 60.0]
transform_name :
lambert_conformal_conic
```
[1 values with dtype=object]
```

Indexes: (3)

time

PandasIndex

PandasIndex(DatetimeIndex(['1979-02-01 03:00:00', '1979-02-01 06:00:00',
               '1979-02-01 09:00:00', '1979-02-01 12:00:00',
               '1979-02-01 15:00:00', '1979-02-01 18:00:00',
               '1979-02-01 21:00:00', '1979-02-02 00:00:00',
               '1979-02-02 03:00:00', '1979-02-02 06:00:00',
               ...
               '2020-12-30 18:00:00', '2020-12-30 21:00:00',
               '2020-12-31 00:00:00', '2020-12-31 03:00:00',
               '2020-12-31 06:00:00', '2020-12-31 09:00:00',
               '2020-12-31 12:00:00', '2020-12-31 15:00:00',
               '2020-12-31 18:00:00', '2020-12-31 21:00:00'],
              dtype='datetime64[ns]', name='time', length=116631, freq=None))

x

PandasIndex

PandasIndex(Float64Index([-2303499.25, -2302499.25, -2301499.25, -2300499.25, -2299499.25,
              -2298499.25, -2297499.25, -2296499.25, -2295499.25, -2294499.25,
              ...
               2294500.75,  2295500.75,  2296500.75,  2297500.75,  2298500.75,
               2299500.75,  2300500.75,  2301500.75,  2302500.75,  2303500.75],
             dtype='float64', name='x', length=4608))

y

PandasIndex

PandasIndex(Float64Index([-1919500.375, -1918500.375, -1917500.375, -1916500.375,
              -1915500.375, -1914500.375, -1913500.375, -1912500.375,
              -1911500.375, -1910500.375,
              ...
               1910499.625,  1911499.625,  1912499.625,  1913499.625,
               1914499.625,  1915499.625,  1916499.625,  1917499.625,
               1918499.625,  1919499.625],
             dtype='float64', name='y', length=3840))

Attributes: (10)
Conventions :
CF-1.6
GDAL_DataType :
Generic
TITLE :
OUTPUT FROM WRF-Hydro v5.2.0-beta2
code_version :
v5.2.0-beta2
model_configuration :
retrospective
model_initialization_time :
1979-02-01_00:00:00
model_output_type :
land
model_output_valid_time :
1979-02-01_03:00:00
model_total_valid_times :
472
proj4 :
+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@null +wktext +no_defs

Select high-priority vars only¶

In [6]:

ds = ds[['ACCET', 'SNEQV', 'FSNO', 'crs']]

In [7]:

ds

Out[7]:

<xarray.Dataset>
Dimensions:  (time: 116631, y: 3840, x: 4608)
Coordinates:
  * time     (time) datetime64[ns] 1979-02-01T03:00:00 ... 2020-12-31T21:00:00
  * x        (x) float64 -2.303e+06 -2.302e+06 ... 2.303e+06 2.304e+06
  * y        (y) float64 -1.92e+06 -1.919e+06 -1.918e+06 ... 1.918e+06 1.919e+06
Data variables:
    ACCET    (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    SNEQV    (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    FSNO     (time, y, x) float64 dask.array<chunksize=(1, 3840, 4608), meta=np.ndarray>
    crs      object ...
Attributes:
    Conventions:                CF-1.6
    GDAL_DataType:              Generic
    TITLE:                      OUTPUT FROM WRF-Hydro v5.2.0-beta2
    code_version:               v5.2.0-beta2
    model_configuration:        retrospective
    model_initialization_time:  1979-02-01_00:00:00
    model_output_type:          land
    model_output_valid_time:    1979-02-01_03:00:00
    model_total_valid_times:    472
    proj4:                      +proj=lcc +units=m +a=6370000.0 +b=6370000.0 ...

In [8]:

ds['ACCET'].isel(time=slice(0,144))

Out[8]:

<xarray.DataArray 'ACCET' (time: 144, y: 3840, x: 4608)>
dask.array<getitem, shape=(144, 3840, 4608), dtype=float64, chunksize=(1, 3840, 4608), chunktype=numpy.ndarray>
Coordinates:
  * time     (time) datetime64[ns] 1979-02-01T03:00:00 ... 1979-02-19
  * x        (x) float64 -2.303e+06 -2.302e+06 ... 2.303e+06 2.304e+06
  * y        (y) float64 -1.92e+06 -1.919e+06 -1.918e+06 ... 1.918e+06 1.919e+06
Attributes:
    esri_pe_string:  PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_Sphere",DAT...
    grid_mapping:    crs
    long_name:       Accumulated total ET
    units:           mm
    valid_range:     [-100000, 100000000]

In [9]:

ds.attrs

Out[9]:

{'Conventions': 'CF-1.6',
 'GDAL_DataType': 'Generic',
 'TITLE': 'OUTPUT FROM WRF-Hydro v5.2.0-beta2',
 'code_version': 'v5.2.0-beta2',
 'model_configuration': 'retrospective',
 'model_initialization_time': '1979-02-01_00:00:00',
 'model_output_type': 'land',
 'model_output_valid_time': '1979-02-01_03:00:00',
 'model_total_valid_times': 472,
 'proj4': '+proj=lcc +units=m +a=6370000.0 +b=6370000.0 +lat_1=30.0 +lat_2=60.0 +lat_0=40.0 +lon_0=-97.0 +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@null +wktext  +no_defs'}

set up zarr stores for temporary and final zarr stores on S3¶

In [10]:

fs_write = fsspec.filesystem('s3', anon=False, skip_instance_cache=True)

In [11]:

temp_name = 'esip-qhub/testing/usgs/nwm1km.tmp'
target_name = 'esip-qhub/testing/usgs/nwm1km.zarr'

In [ ]:

fs_write.rm(temp_name, recursive=True)

In [ ]:

fs_write.rm(target_name, recursive=True)

In [15]:

temp_store = fs_write.get_mapper(temp_name)
target_store = fs_write.get_mapper(target_name)

In [16]:

ds = ds.drop('crs')

In [17]:

a = len(ds.time)/(144/2)

In [18]:

b = (len(ds.x) * len(ds.y))/((96*2)*(132*2))

In [19]:

a/b

Out[19]:

4.640266927083334

Rechunk!¶

In [20]:

rechunked = rechunk(ds.isel(time=slice(0,144)), target_chunks={'y':96*2, 'x':132*2, 'time':144/2},
                    target_store=target_store, temp_store=temp_store, max_mem='3.5GiB')

---------------------------------------------------------------------------
ReadOnlyError                             Traceback (most recent call last)
Cell In[20], line 1
----> 1 rechunked = rechunk(ds.isel(time=slice(0,144)), target_chunks={'y':96*2, 'x':132*2, 'time':144/2},
      2                     target_store=target_store, temp_store=temp_store, max_mem='3.5GiB')

File /home/conda/users/6da996e8bac98e9914efe8740ae916e58e2ced6a618659b47fccce8f5bc9cd6a-20230314-182641-028005-117-pangeo/lib/python3.10/site-packages/rechunker/api.py:297, in rechunk(source, target_chunks, max_mem, target_store, target_options, temp_store, temp_options, executor)
    294 if isinstance(executor, str):
    295     executor = _get_executor(executor)
--> 297 copy_spec, intermediate, target = _setup_rechunk(
    298     source=source,
    299     target_chunks=target_chunks,
    300     max_mem=max_mem,
    301     target_store=target_store,
    302     target_options=target_options,
    303     temp_store=temp_store,
    304     temp_options=temp_options,
    305 )
    306 plan = executor.prepare_plan(copy_spec)
    307 return Rechunked(executor, plan, source, intermediate, target)

File /home/conda/users/6da996e8bac98e9914efe8740ae916e58e2ced6a618659b47fccce8f5bc9cd6a-20230314-182641-028005-117-pangeo/lib/python3.10/site-packages/rechunker/api.py:395, in _setup_rechunk(source, target_chunks, max_mem, target_store, target_options, temp_store, temp_options)
    392 attrs = _encode_zarr_attributes(attrs)
    394 if temp_store is not None:
--> 395     temp_group = zarr.group(temp_store)
    396 else:
    397     temp_group = None

File /home/conda/users/6da996e8bac98e9914efe8740ae916e58e2ced6a618659b47fccce8f5bc9cd6a-20230314-182641-028005-117-pangeo/lib/python3.10/site-packages/zarr/hierarchy.py:1355, in group(store, overwrite, chunk_store, cache_attrs, synchronizer, path, zarr_version)
   1352     requires_init = overwrite or not contains_group(store, path)
   1354 if requires_init:
-> 1355     init_group(store, overwrite=overwrite, chunk_store=chunk_store,
   1356                path=path)
   1358 return Group(store, read_only=False, chunk_store=chunk_store,
   1359              cache_attrs=cache_attrs, synchronizer=synchronizer, path=path,
   1360              zarr_version=zarr_version)

File /home/conda/users/6da996e8bac98e9914efe8740ae916e58e2ced6a618659b47fccce8f5bc9cd6a-20230314-182641-028005-117-pangeo/lib/python3.10/site-packages/zarr/storage.py:643, in init_group(store, overwrite, path, chunk_store)
    640     store['zarr.json'] = store._metadata_class.encode_hierarchy_metadata(None)  # type: ignore
    642 # initialise metadata
--> 643 _init_group_metadata(store=store, overwrite=overwrite, path=path,
    644                      chunk_store=chunk_store)
    646 if store_version == 3:
    647     # TODO: Should initializing a v3 group also create a corresponding
    648     #       empty folder under data/root/? I think probably not until there
    649     #       is actual data written there.
    650     pass

File /home/conda/users/6da996e8bac98e9914efe8740ae916e58e2ced6a618659b47fccce8f5bc9cd6a-20230314-182641-028005-117-pangeo/lib/python3.10/site-packages/zarr/storage.py:706, in _init_group_metadata(store, overwrite, path, chunk_store)
    704 key = _prefix_to_group_key(store, _path_to_prefix(path))
    705 if hasattr(store, '_metadata_class'):
--> 706     store[key] = store._metadata_class.encode_group_metadata(meta)  # type: ignore
    707 else:
    708     store[key] = encode_group_metadata(meta)

File /home/conda/users/6da996e8bac98e9914efe8740ae916e58e2ced6a618659b47fccce8f5bc9cd6a-20230314-182641-028005-117-pangeo/lib/python3.10/site-packages/zarr/storage.py:1398, in FSStore.__setitem__(self, key, value)
   1396 def __setitem__(self, key, value):
   1397     if self.mode == 'r':
-> 1398         raise ReadOnlyError()
   1399     key = self._normalize_key(key)
   1400     path = self.dir_path(key)

ReadOnlyError: object is read-only

In [ ]:

%%time
rechunked.execute(retries=10)

In [ ]:

zarr.convenience.consolidate_metadata(target_store)

Explore the rechunked dataset¶

In [ ]:

ds2 = xr.open_dataset(target_store, engine='zarr', chunks={})

In [ ]:

ds2

In [ ]:

ds2.ACCET

In [ ]:

import hvplot.xarray

In [ ]:

ds2.ACCET[:,2000,2000].hvplot(x='time')

In [ ]:

import rechunker

In [ ]:

rechunker.__version__

In [ ]: