Creating a STAC Collection for a Virtual Icechunk Store

There is a virtual icechunk store that is publicly available at: s3://nasa-waterinsight/virtual-zarr-store/NLDAS-3-icechunk/

This notebook goes through the current thinking for how you would set up a STAC collection that points to that virtual icechunk store and provides all the information a user needs to interact with the virtual zarr store programmatically or via a web UI.

import json
import datetime

import icechunk
import pystac
import xstac
import zarr

import xarray as xr

Zarr can emit a lot of warnings about Numcodecs not being including in the Zarr version 3 specification yet – let’s suppress those.

import warnings

warnings.filterwarnings(
    "ignore",
    message="Numcodecs codecs are not in the Zarr version 3 specification*",
    category=UserWarning,
)

These are the PRs that need to land before you can open the virtual icechunk store with zarr directly:

• https://github.com/zarr-developers/zarr-python/pull/3369
• https://github.com/earth-mover/icechunk/pull/1161

Until then:

storage = icechunk.s3_storage(
    bucket="nasa-waterinsight",
    prefix="virtual-zarr-store/NLDAS-3-icechunk/",
    region="us-west-2",
    anonymous=True,
)

The bucket and prefix are from the icechunk href. The anonymous=True needs to come from somewhere else.

config = icechunk.RepositoryConfig.default()
config.set_virtual_chunk_container(
    icechunk.VirtualChunkContainer(
        "s3://nasa-waterinsight/NLDAS3/forcing/daily/",
        icechunk.s3_store(region="us-west-2")
    )
)
virtual_credentials = icechunk.containers_credentials(
    {
        "s3://nasa-waterinsight/NLDAS3/forcing/daily/": icechunk.s3_anonymous_credentials()
    }
)

Here we need the href for the internal storage bucket(s) (composed of bucket and prefix) and we need the region of that bucket. Then we need some way of providing credentials.

repo = icechunk.Repository.open(
    storage=storage,
    config=config,
    authorize_virtual_chunk_access=virtual_credentials,
)

session = repo.readonly_session(snapshot_id='YTNGFY4WY9189GEH1FNG')

Since icechunk manages versions (like git) we need some way of knowing which branch, tag or snapshot_id (similar to commit) to use.

ds = xr.open_zarr(session.store, consolidated=False, zarr_format=3)

Last of all we need a way of specifying that we are looking at icechunk here as well as the standard fields: consolidated, zarr_format that are already included in the STAC Zarr extension.

Note that it is possible that these last two are not actually required. Xarray should know that for icechunk stores consolidated is always false and similarly xarray should be able to infer the zarr format from the store itself.

ds

<xarray.Dataset> Size: 51TB
Dimensions:   (time: 8399, lat: 6500, lon: 11700)
Coordinates:
  * time      (time) datetime64[ns] 67kB 2001-01-02 2001-01-03 ... 2024-01-01
  * lat       (lat) float64 52kB 7.005 7.015 7.025 7.035 ... 71.97 71.98 71.99
  * lon       (lon) float64 94kB -169.0 -169.0 -169.0 ... -52.03 -52.01 -52.0
Data variables:
    PSurf     (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Qair      (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    SWdown    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Rainf     (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    LWdown    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Tair      (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Tair_min  (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Wind_E    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Wind_N    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Tair_max  (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
Attributes: (12/17)
    missing_value:          -9999.0
    time_definition:        daily
    shortname:              NLDAS_FOR0010_D_3.0
    title:                  NLDAS Forcing Data L4 Daily 0.01 x 0.01 degree V3...
    version:                3.0 beta
    institution:            NASA GSFC
    ...                     ...
    websites:               https://ldas.gsfc.nasa.gov/nldas/v3/ ; https://li...
    MAP_PROJECTION:         EQUIDISTANT CYLINDRICAL
    SOUTH_WEST_CORNER_LAT:  7.005000114440918
    SOUTH_WEST_CORNER_LON:  -168.9949951171875
    DX:                     0.009999999776482582
    DY:                     0.009999999776482582

xarray.Dataset

• Dimensions:
  • time: 8399
  • lat: 6500
  • lon: 11700
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2001-01-02 ... 2024-01-01

    long_name :

    time

    time_increment :

    one day

    begin_date :

    20010101

    begin_time :

    000000

    end_date :

    20010101

    end_time :

    235959

    array(['2001-01-02T00:00:00.000000000', '2001-01-03T00:00:00.000000000',
           '2001-01-04T00:00:00.000000000', ..., '2023-12-30T00:00:00.000000000',
           '2023-12-31T00:00:00.000000000', '2024-01-01T00:00:00.000000000'],
          shape=(8399,), dtype='datetime64[ns]')

  • lat
    (lat)
    float64
    7.005 7.015 7.025 ... 71.98 71.99

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    vmin :

    7.005000114440918

    vmax :

    71.9949951171875

    array([ 7.005   ,  7.015   ,  7.025   , ..., 71.974998, 71.984993, 71.994995],
          shape=(6500,))

  • lon
    (lon)
    float64
    -169.0 -169.0 ... -52.01 -52.0

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    vmin :

    -168.9949951171875

    vmax :

    -52.00499725341797

    array([-168.994995, -168.985001, -168.974991, ...,  -52.025002,  -52.014999,
            -52.004997], shape=(11700,))

• Data variables: (10)
  • PSurf
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    Pa

    standard_name :

    Surface pressure

    long_name :

    Surface pressure

    cell_methods :

    time: mean

    vmin :

    46024.41796875

    vmax :

    103161.1796875

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Qair
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    kg kg-1

    standard_name :

    Near surface specific humidity

    long_name :

    2-meter above ground Specific humidity

    cell_methods :

    time: mean

    vmin :

    5.8038440329255536e-05

    vmax :

    0.021006619557738304

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SWdown
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    W m-2

    standard_name :

    Surface incident shortwave radiation

    long_name :

    Shortwave radiation flux downwards (surface)

    cell_methods :

    time: mean

    vmin :

    0.0

    vmax :

    288.2958068847656

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Rainf
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    kg m-2

    standard_name :

    Total precipitation

    long_name :

    Total precipitation

    cell_methods :

    time: sum

    vmin :

    0.0

    vmax :

    60.617774963378906

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • LWdown
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    W m-2

    standard_name :

    Surface incident longwave radiation

    long_name :

    Longwave radiation flux downwards (surface)

    cell_methods :

    time: mean

    vmin :

    126.3478012084961

    vmax :

    450.3310546875

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Tair
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    K

    standard_name :

    Near surface air temperature

    long_name :

    2-meter above ground Temperature

    cell_methods :

    time: mean

    vmin :

    228.31163024902344

    vmax :

    304.30316162109375

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Tair_min
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    K

    standard_name :

    Daily minimum near surface air temperature

    long_name :

    Daily minimum 2-meter above ground Temperature

    cell_methods :

    time: minimum

    vmin :

    227.1617889404297

    vmax :

    300.4527282714844

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Wind_E
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    m s-1

    standard_name :

    Near surface eastward wind component

    long_name :

    10-meter above ground Zonal wind speed

    cell_methods :

    time: mean

    vmin :

    -12.775843620300293

    vmax :

    11.65600872039795

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Wind_N
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    m s-1

    standard_name :

    Near surface northward wind component

    long_name :

    10-meter above ground Meridional wind speed

    cell_methods :

    time: mean

    vmin :

    -10.617785453796387

    vmax :

    11.41434097290039

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Tair_max
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    K

    standard_name :

    Daily maximum near surface air temperature

    long_name :

    Daily maximum 2-meter above ground Temperature

    cell_methods :

    time: maximum

    vmin :

    230.1659393310547

    vmax :

    312.7441711425781

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2001-01-02', '2001-01-03', '2001-01-04', '2001-01-05',
                   '2001-01-06', '2001-01-07', '2001-01-08', '2001-01-09',
                   '2001-01-10', '2001-01-11',
                   ...
                   '2023-12-23', '2023-12-24', '2023-12-25', '2023-12-26',
                   '2023-12-27', '2023-12-28', '2023-12-29', '2023-12-30',
                   '2023-12-31', '2024-01-01'],
                  dtype='datetime64[ns]', name='time', length=8399, freq=None))

  • lat
    PandasIndex

    PandasIndex(Index([7.005000114440918, 7.015000343322754, 7.025000095367432,
           7.035000324249268, 7.045000076293945, 7.055000305175781,
           7.065000057220459, 7.075000286102295, 7.085000038146973,
           7.095000267028809,
           ...
           71.90499877929688, 71.91499328613281, 71.92499542236328,
           71.93499755859375, 71.94499206542969, 71.95499420166016,
           71.96499633789062,  71.9749984741211, 71.98499298095703,
            71.9949951171875],
          dtype='float64', name='lat', length=6500))

  • lon
    PandasIndex

    PandasIndex(Index([ -168.9949951171875, -168.98500061035156, -168.97499084472656,
           -168.96499633789062,  -168.9550018310547,  -168.9449920654297,
           -168.93499755859375, -168.92498779296875,  -168.9149932861328,
           -168.90499877929688,
           ...
           -52.095001220703125, -52.084999084472656,  -52.07499694824219,
            -52.06499481201172,  -52.05500030517578,  -52.04499816894531,
           -52.034996032714844, -52.025001525878906,  -52.01499938964844,
            -52.00499725341797],
          dtype='float64', name='lon', length=11700))

• Attributes: (17)

  missing_value :

  -9999.0

  time_definition :

  daily

  shortname :

  NLDAS_FOR0010_D_3.0

  title :

  NLDAS Forcing Data L4 Daily 0.01 x 0.01 degree V3.0 - *BETA*

  version :

  3.0 beta

  institution :

  NASA GSFC

  doi :

  TBD

  history :

  created on date: 2025-06-07T15:04:40

  reference :

  Maina-et-al_to-be-submitted; doi:TBD

  conventions :

  CF-1.6

  standard :

  ALMA v3 ; https://www.lmd.jussieu.fr/~polcher/ALMA/

  websites :

  https://ldas.gsfc.nasa.gov/nldas/v3/ ; https://lis.gsfc.nasa.gov/

  MAP_PROJECTION :

  EQUIDISTANT CYLINDRICAL

  SOUTH_WEST_CORNER_LAT :

  7.005000114440918

  SOUTH_WEST_CORNER_LON :

  -168.9949951171875

  DX :

  0.009999999776482582

  DY :

  0.009999999776482582

Ok! We have the xarray dataset lazily opened by accessing only the virtual icechunk store! Now the goal is to create a STAC collection that describes that virtual icechunk store and contains all the information we need for accessing it (all the inputs to the functions above).

Extract metadata for STAC collection

Now let’s see how xstac can help us extract the variables and represent them using the Datacube STAC Extension.

This section takes inspiration from https://github.com/stac-utils/xstac/blob/main/examples/nasa-nex-gddp-cmip6/generate.py

We’ll start with some of the hard-coded values that will need to be provided for any given dataset. As much as possible these should be lifted directly from the xr.Dataset attrs.

collection_id = "nldas-3"
description = (
    "NLDAS-3 provides a fine-scale (1 km) meteorological forcing (precipitation) in "
    "both retrospective and near real-time over North and Central America, including "
    "Alaska, Hawaii, and Puerto Rico, by leveraging high-quality gauge, satellite, "
    "and model datasets through advanced data assimilation methods. Read more: "
    "https://ldas.gsfc.nasa.gov/nldas/v3"
)
providers = [
    pystac.Provider(
        name="NLDAS",
        roles=["producer", "processor", "licensor"],
        url="https://ldas.gsfc.nasa.gov/nldas"
    )
]

I want to draw special attention to how we can use the Storage STAC Extension to capture that this particular bucket can be accessed anonymously. We can also capture the region within this blob.

storage_schemes = {
    "aws-s3-nasa-waterinsight": {
        "type": "aws-s3",
        "platform": "https://{bucket}.s3.{region}.amazonaws.com",
        "bucket": "nasa-waterinsight",
        "region": "us-west-2",
        "anonymous": True,
    }
}

Now let’s configure some metadata that can be gotten from the xr.Dataset itself.

title = ds.attrs["title"]
extents = pystac.Extent(
    spatial=pystac.SpatialExtent(bboxes=[list(ds.rio.bounds())]),
    temporal=pystac.TemporalExtent(
        intervals=[
            datetime.datetime.fromisoformat(str(ds.time.min().values)),
            datetime.datetime.fromisoformat(str(ds.time.max().values))
        ]
    ),
)

Now that we have all those values set, create a pystac.Collection:

template = pystac.Collection(
    collection_id,
    description=description,
    extent=extents,
    extra_fields={"storage:schemes": storage_schemes, "item_assets": {}},
    providers=providers,
    title=title,
    stac_extensions=[
        "https://stac-extensions.github.io/storage/v2.0.0/schema.json",
    ]
)

Now that we have a preliminary version of the STAC Collection we can pass it off to xstac to pull out the variables and dims using the Datacube STAC Extension.

collection = xstac.xarray_to_stac(
    ds,
    template,
    temporal_dimension="time",
    x_dimension="lon",
    y_dimension="lat",
    reference_system=4326,
    validate=False,
)

With that collection in hand we can create an asset that points to the virtual icechunk store and add it as a collection-level asset.

Add collection-level assets

The main concern of the asset is how to access the data. So we need the asset to contain all the information we need to pass into icechunk functions when opening the virtual icechunk store.

snapshot_id = "YTNGFY4WY9189GEH1FNG"
virtual_href = "s3://nasa-waterinsight/virtual-zarr-store/NLDAS-3-icechunk/"

name = "nldas-3"
virtual_key = f"{name}@{snapshot_id}"
storage_ref = "aws-s3-nasa-waterinsight"

collection.add_asset(
    virtual_key,
    pystac.Asset(
        virtual_href,
        title="NLDAS-3 Virtual Zarr Store",
        media_type="application/vnd.zarr+icechunk",  # I made this up: discussion https://earthmover-community.slack.com/archives/C07NQCBSTB7/p1756918042834049
        roles=["data", "references", "virtual", "latest-version"],
        extra_fields={
            "zarr:consolidated": False,
            "zarr:zarr_format": 3,
            "icechunk:snapshot_id": snapshot_id,
            "storage:refs": [storage_ref],
        }
    )
)

We also need to specify how to access the legacy files which potentially sit in their own bucket. We can use the Virtual Assets STAC Extension to capture that.

legacy_href = "s3://nasa-waterinsight/NLDAS3/forcing/daily/"
legacy_key = "nldas-3-legacy-bucket"

collection.add_asset(
    legacy_key,
    pystac.Asset(
        legacy_href,
        title="NLDAS-3 Legacy Bucket",
        media_type="application/x-netcdf",
        roles=["data"],
        extra_fields={
            "storage:refs": [storage_ref],
        }
    )
)

collection.assets[virtual_key].extra_fields["vrt:hrefs"] = [
    {
        "key": legacy_key,
        "href": f"https://raw.githubusercontent.com/NASA-IMPACT/dse-virtual-zarr-workshop/refs/heads/main/docs/examples/collection.json#/assets/{legacy_key}"
    },
]

We can also add information about how to render each variable. This uses the Render STAC Extension and specifies how applications (for instance titiler) should represent each variable visually.

renders = {}
for k in ds:
    if k.startswith("LW"):
        colormap_name = "inferno"
    elif k.startswith("PS"):
        colormap_name = "viridis"
    elif k.startswith("Q"):
        colormap_name = "plasma"
    elif k.startswith("Rain"):
        colormap_name = "cfastie"
    elif k.startswith("SW"):
        colormap_name = "magma"
    elif k.startswith("T"):
        colormap_name = "RdYlBu_r"
    elif k.startswith("Wind"):
        colormap_name = "PuOr"
        
    renders[k] = {
        "title": ds[k].attrs["long_name"],
        "assets": ["nldas-3"],
        "resampling": "average",
        "colormap_name": colormap_name,
        "rescale": [[ds[k].attrs["vmin"], ds[k].attrs["vmax"]]],
        "backend": "xarray",
    }

collection.extra_fields["renders"] = renders

Last of all we will add the Zarr STAC Extension and the Render STAC Extension to the list of stac_extensions. We need to do this after adding the asset because they relies on there being assets in the collection.

collection.stac_extensions = list(set(
    collection.stac_extensions + [
        "https://stac-extensions.github.io/render/v2.0.0/schema.json",
        "https://stac-extensions.github.io/virtual-assets/v1.0.0/schema.json",
        "https://stac-extensions.github.io/zarr/v1.1.0/schema.json",
        "https://stac-extensions.github.io/version/v1.2.0/schema.json",
    ]
))

Validate the collection:

collection.validate()

['https://schemas.stacspec.org/v1.1.0/collection-spec/json-schema/collection.json',
 'https://stac-extensions.github.io/zarr/v1.1.0/schema.json',
 'https://stac-extensions.github.io/datacube/v2.2.0/schema.json',
 'https://stac-extensions.github.io/version/v1.2.0/schema.json',
 'https://stac-extensions.github.io/storage/v2.0.0/schema.json',
 'https://stac-extensions.github.io/render/v2.0.0/schema.json',
 'https://stac-extensions.github.io/virtual-assets/v1.0.0/schema.json']

Dump the collection to json

with open("collection.json", "w") as f:
    json.dump(collection.to_dict(), f, indent=2)

Read the virtual icechunk using collection-level asset

There is functionality in xpystac to open the virtual icechunk store referenced by the virtual asset so with properly configured STAC, you don’t need to interact with it at all.

Now let’s read the collection in from where we stored it in the json, find the latest version of the icechunk store, and lazily open it as an xarray.Dataset

collection = pystac.Collection.from_file("collection.json")

assets = collection.get_assets(media_type="application/vnd.zarr+icechunk", role="latest-version")
assets

{'nldas-3@YTNGFY4WY9189GEH1FNG': <Asset href=s3://nasa-waterinsight/virtual-zarr-store/NLDAS-3-icechunk/>}

ds = xr.open_dataset(list(assets.values())[0], chunks={}, engine="stac")
ds

<xarray.Dataset> Size: 51TB
Dimensions:   (time: 8399, lat: 6500, lon: 11700)
Coordinates:
  * lat       (lat) float64 52kB 7.005 7.015 7.025 7.035 ... 71.97 71.98 71.99
  * lon       (lon) float64 94kB -169.0 -169.0 -169.0 ... -52.03 -52.01 -52.0
  * time      (time) datetime64[ns] 67kB 2001-01-02 2001-01-03 ... 2024-01-01
Data variables:
    Rainf     (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    PSurf     (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    LWdown    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Qair      (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    SWdown    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Wind_E    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Tair_max  (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Wind_N    (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Tair      (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
    Tair_min  (time, lat, lon) float64 5TB dask.array<chunksize=(1, 500, 900), meta=np.ndarray>
Attributes: (12/17)
    missing_value:          -9999.0
    time_definition:        daily
    shortname:              NLDAS_FOR0010_D_3.0
    title:                  NLDAS Forcing Data L4 Daily 0.01 x 0.01 degree V3...
    version:                3.0 beta
    institution:            NASA GSFC
    ...                     ...
    websites:               https://ldas.gsfc.nasa.gov/nldas/v3/ ; https://li...
    MAP_PROJECTION:         EQUIDISTANT CYLINDRICAL
    SOUTH_WEST_CORNER_LAT:  7.005000114440918
    SOUTH_WEST_CORNER_LON:  -168.9949951171875
    DX:                     0.009999999776482582
    DY:                     0.009999999776482582

xarray.Dataset

• Dimensions:
  • time: 8399
  • lat: 6500
  • lon: 11700
• Coordinates: (3)
  • lat
    (lat)
    float64
    7.005 7.015 7.025 ... 71.98 71.99

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    vmin :

    7.005000114440918

    vmax :

    71.9949951171875

    array([ 7.005   ,  7.015   ,  7.025   , ..., 71.974998, 71.984993, 71.994995],
          shape=(6500,))

  • lon
    (lon)
    float64
    -169.0 -169.0 ... -52.01 -52.0

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    vmin :

    -168.9949951171875

    vmax :

    -52.00499725341797

    array([-168.994995, -168.985001, -168.974991, ...,  -52.025002,  -52.014999,
            -52.004997], shape=(11700,))

  • time
    (time)
    datetime64[ns]
    2001-01-02 ... 2024-01-01

    long_name :

    time

    time_increment :

    one day

    begin_date :

    20010101

    begin_time :

    000000

    end_date :

    20010101

    end_time :

    235959

    array(['2001-01-02T00:00:00.000000000', '2001-01-03T00:00:00.000000000',
           '2001-01-04T00:00:00.000000000', ..., '2023-12-30T00:00:00.000000000',
           '2023-12-31T00:00:00.000000000', '2024-01-01T00:00:00.000000000'],
          shape=(8399,), dtype='datetime64[ns]')

• Data variables: (10)
  • Rainf
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    kg m-2

    standard_name :

    Total precipitation

    long_name :

    Total precipitation

    cell_methods :

    time: sum

    vmin :

    0.0

    vmax :

    60.617774963378906

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • PSurf
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    Pa

    standard_name :

    Surface pressure

    long_name :

    Surface pressure

    cell_methods :

    time: mean

    vmin :

    46024.41796875

    vmax :

    103161.1796875

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • LWdown
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    W m-2

    standard_name :

    Surface incident longwave radiation

    long_name :

    Longwave radiation flux downwards (surface)

    cell_methods :

    time: mean

    vmin :

    126.3478012084961

    vmax :

    450.3310546875

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Qair
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    kg kg-1

    standard_name :

    Near surface specific humidity

    long_name :

    2-meter above ground Specific humidity

    cell_methods :

    time: mean

    vmin :

    5.8038440329255536e-05

    vmax :

    0.021006619557738304

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SWdown
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    W m-2

    standard_name :

    Surface incident shortwave radiation

    long_name :

    Shortwave radiation flux downwards (surface)

    cell_methods :

    time: mean

    vmin :

    0.0

    vmax :

    288.2958068847656

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Wind_E
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    m s-1

    standard_name :

    Near surface eastward wind component

    long_name :

    10-meter above ground Zonal wind speed

    cell_methods :

    time: mean

    vmin :

    -12.775843620300293

    vmax :

    11.65600872039795

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Tair_max
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    K

    standard_name :

    Daily maximum near surface air temperature

    long_name :

    Daily maximum 2-meter above ground Temperature

    cell_methods :

    time: maximum

    vmin :

    230.1659393310547

    vmax :

    312.7441711425781

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Wind_N
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    m s-1

    standard_name :

    Near surface northward wind component

    long_name :

    10-meter above ground Meridional wind speed

    cell_methods :

    time: mean

    vmin :

    -10.617785453796387

    vmax :

    11.41434097290039

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Tair
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    K

    standard_name :

    Near surface air temperature

    long_name :

    2-meter above ground Temperature

    cell_methods :

    time: mean

    vmin :

    228.31163024902344

    vmax :

    304.30316162109375

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

  • Tair_min
    (time, lat, lon)
    float64
    dask.array<chunksize=(1, 500, 900), meta=np.ndarray>

    units :

    K

    standard_name :

    Daily minimum near surface air temperature

    long_name :

    Daily minimum 2-meter above ground Temperature

    cell_methods :

    time: minimum

    vmin :

    227.1617889404297

    vmax :

    300.4527282714844

    Array Chunk Bytes 4.65 TiB 3.43 MiB Shape (8399, 6500, 11700) (1, 500, 900) Dask graph 1419431 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([7.005000114440918, 7.015000343322754, 7.025000095367432,
           7.035000324249268, 7.045000076293945, 7.055000305175781,
           7.065000057220459, 7.075000286102295, 7.085000038146973,
           7.095000267028809,
           ...
           71.90499877929688, 71.91499328613281, 71.92499542236328,
           71.93499755859375, 71.94499206542969, 71.95499420166016,
           71.96499633789062,  71.9749984741211, 71.98499298095703,
            71.9949951171875],
          dtype='float64', name='lat', length=6500))

  • lon
    PandasIndex

    PandasIndex(Index([ -168.9949951171875, -168.98500061035156, -168.97499084472656,
           -168.96499633789062,  -168.9550018310547,  -168.9449920654297,
           -168.93499755859375, -168.92498779296875,  -168.9149932861328,
           -168.90499877929688,
           ...
           -52.095001220703125, -52.084999084472656,  -52.07499694824219,
            -52.06499481201172,  -52.05500030517578,  -52.04499816894531,
           -52.034996032714844, -52.025001525878906,  -52.01499938964844,
            -52.00499725341797],
          dtype='float64', name='lon', length=11700))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2001-01-02', '2001-01-03', '2001-01-04', '2001-01-05',
                   '2001-01-06', '2001-01-07', '2001-01-08', '2001-01-09',
                   '2001-01-10', '2001-01-11',
                   ...
                   '2023-12-23', '2023-12-24', '2023-12-25', '2023-12-26',
                   '2023-12-27', '2023-12-28', '2023-12-29', '2023-12-30',
                   '2023-12-31', '2024-01-01'],
                  dtype='datetime64[ns]', name='time', length=8399, freq=None))

• Attributes: (17)

  missing_value :

  -9999.0

  time_definition :

  daily

  shortname :

  NLDAS_FOR0010_D_3.0

  title :

  NLDAS Forcing Data L4 Daily 0.01 x 0.01 degree V3.0 - *BETA*

  version :

  3.0 beta

  institution :

  NASA GSFC

  doi :

  TBD

  history :

  created on date: 2025-06-07T15:04:40

  reference :

  Maina-et-al_to-be-submitted; doi:TBD

  conventions :

  CF-1.6

  standard :

  ALMA v3 ; https://www.lmd.jussieu.fr/~polcher/ALMA/

  websites :

  https://ldas.gsfc.nasa.gov/nldas/v3/ ; https://lis.gsfc.nasa.gov/

  MAP_PROJECTION :

  EQUIDISTANT CYLINDRICAL

  SOUTH_WEST_CORNER_LAT :

  7.005000114440918

  SOUTH_WEST_CORNER_LON :

  -168.9949951171875

  DX :

  0.009999999776482582

  DY :

  0.009999999776482582

Make sure you can actually access underlying data:

%%time

cape_rain = ds.Rainf.sel(
    time="2023-07-16",
    lat=slice(41.48, 42.10),
    lon=slice(-70.84, -69.77)
).compute()

CPU times: user 3.05 s, sys: 596 ms, total: 3.65 s
Wall time: 5.89 s

cape_rain

<xarray.DataArray 'Rainf' (lat: 62, lon: 107)> Size: 53kB
array([[ 5.16139793,  4.73033953,  4.79954004, ...,  2.80910897,
         2.8472836 ,  2.64635181],
       [ 4.856112  ,  4.68789768,  4.71580458, ...,  2.80869484,
         2.87434697,  2.80177665],
       [ 4.72029495,  4.83154964,  4.46945238, ...,  2.90991831,
         3.02153993,  3.10612035],
       ...,
       [ 7.01298904,  7.19748259,  7.40456963, ..., 18.69691277,
        18.66467857, 18.25347519],
       [ 7.89184809,  8.40285873,  7.87267876, ..., 17.93091774,
        19.15982437, 18.79140854],
       [ 8.36418152,  8.91619682,  8.63145733, ..., 18.75998878,
        20.46075821, 20.56634521]], shape=(62, 107))
Coordinates:
  * lat      (lat) float64 496B 41.49 41.49 41.51 41.51 ... 42.08 42.08 42.1
  * lon      (lon) float64 856B -70.83 -70.82 -70.81 ... -69.79 -69.78 -69.77
    time     datetime64[ns] 8B 2023-07-16
Attributes:
    units:          kg m-2
    standard_name:  Total precipitation
    long_name:      Total precipitation
    cell_methods:   time: sum
    vmin:           0.0
    vmax:           60.617774963378906

xarray.DataArray

'Rainf'

• lat: 62
• lon: 107

• 5.161 4.73 4.8 4.876 4.935 5.637 ... 17.67 17.78 18.76 20.46 20.57

  array([[ 5.16139793,  4.73033953,  4.79954004, ...,  2.80910897,
           2.8472836 ,  2.64635181],
         [ 4.856112  ,  4.68789768,  4.71580458, ...,  2.80869484,
           2.87434697,  2.80177665],
         [ 4.72029495,  4.83154964,  4.46945238, ...,  2.90991831,
           3.02153993,  3.10612035],
         ...,
         [ 7.01298904,  7.19748259,  7.40456963, ..., 18.69691277,
          18.66467857, 18.25347519],
         [ 7.89184809,  8.40285873,  7.87267876, ..., 17.93091774,
          19.15982437, 18.79140854],
         [ 8.36418152,  8.91619682,  8.63145733, ..., 18.75998878,
          20.46075821, 20.56634521]], shape=(62, 107))

• Coordinates: (3)
  • lat
    (lat)
    float64
    41.49 41.49 41.51 ... 42.08 42.1

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    vmin :

    7.005000114440918

    vmax :

    71.9949951171875

    array([41.485001, 41.494999, 41.505001, 41.514999, 41.525002, 41.535   ,
           41.545002, 41.555   , 41.564999, 41.575001, 41.584999, 41.595001,
           41.605   , 41.615002, 41.625   , 41.635002, 41.645   , 41.654999,
           41.665001, 41.674999, 41.685001, 41.695   , 41.705002, 41.715   ,
           41.724998, 41.735001, 41.744999, 41.755001, 41.764999, 41.775002,
           41.785   , 41.795002, 41.805   , 41.814999, 41.825001, 41.834999,
           41.845001, 41.855   , 41.865002, 41.875   , 41.885002, 41.895   ,
           41.904999, 41.915001, 41.924999, 41.935001, 41.945   , 41.955002,
           41.965   , 41.974998, 41.985001, 41.994999, 42.005001, 42.014999,
           42.025002, 42.035   , 42.045002, 42.055   , 42.064999, 42.075001,
           42.084999, 42.095001])

  • lon
    (lon)
    float64
    -70.83 -70.82 ... -69.78 -69.77

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    vmin :

    -168.9949951171875

    vmax :

    -52.00499725341797

    array([-70.834999, -70.824997, -70.814995, -70.805   , -70.794998, -70.784996,
           -70.774994, -70.764999, -70.754997, -70.744995, -70.735001, -70.724998,
           -70.714996, -70.704994, -70.695   , -70.684998, -70.674995, -70.665001,
           -70.654999, -70.644997, -70.634995, -70.625   , -70.614998, -70.604996,
           -70.594994, -70.584999, -70.574997, -70.564995, -70.555   , -70.544998,
           -70.534996, -70.524994, -70.514999, -70.504997, -70.494995, -70.485001,
           -70.474998, -70.464996, -70.454994, -70.445   , -70.434998, -70.424995,
           -70.415001, -70.404999, -70.394997, -70.384995, -70.375   , -70.364998,
           -70.354996, -70.344994, -70.334999, -70.324997, -70.314995, -70.305   ,
           -70.294998, -70.284996, -70.274994, -70.264999, -70.254997, -70.244995,
           -70.235001, -70.224998, -70.214996, -70.204994, -70.195   , -70.184998,
           -70.174995, -70.165001, -70.154999, -70.144997, -70.134995, -70.125   ,
           -70.114998, -70.104996, -70.094994, -70.084999, -70.074997, -70.064995,
           -70.055   , -70.044998, -70.034996, -70.024994, -70.014999, -70.004997,
           -69.994995, -69.985001, -69.974998, -69.964996, -69.954994, -69.945   ,
           -69.934998, -69.924995, -69.915001, -69.904999, -69.894997, -69.884995,
           -69.875   , -69.864998, -69.854996, -69.844994, -69.834999, -69.824997,
           -69.814995, -69.805   , -69.794998, -69.784996, -69.774994])

  • time
    ()
    datetime64[ns]
    2023-07-16

    long_name :

    time

    time_increment :

    one day

    begin_date :

    20010101

    begin_time :

    000000

    end_date :

    20010101

    end_time :

    235959

    array('2023-07-16T00:00:00.000000000', dtype='datetime64[ns]')

• Indexes: (2)
  • lat
    PandasIndex

    PandasIndex(Index([ 41.48500061035156, 41.494998931884766, 41.505001068115234,
            41.51499938964844, 41.525001525878906,  41.53499984741211,
            41.54500198364258,  41.55500030517578, 41.564998626708984,
            41.57500076293945, 41.584999084472656, 41.595001220703125,
            41.60499954223633,   41.6150016784668,             41.625,
            41.63500213623047,  41.64500045776367, 41.654998779296875,
           41.665000915527344,  41.67499923706055, 41.685001373291016,
            41.69499969482422,  41.70500183105469,  41.71500015258789,
           41.724998474121094,  41.73500061035156, 41.744998931884766,
           41.755001068115234,  41.76499938964844, 41.775001525878906,
            41.78499984741211,  41.79500198364258,  41.80500030517578,
           41.814998626708984,  41.82500076293945, 41.834999084472656,
           41.845001220703125,  41.85499954223633,   41.8650016784668,
                       41.875,  41.88500213623047,  41.89500045776367,
           41.904998779296875, 41.915000915527344,  41.92499923706055,
           41.935001373291016,  41.94499969482422,  41.95500183105469,
            41.96500015258789, 41.974998474121094,  41.98500061035156,
           41.994998931884766, 42.005001068115234,  42.01499938964844,
           42.025001525878906,  42.03499984741211,  42.04500198364258,
            42.05500030517578, 42.064998626708984,  42.07500076293945,
           42.084999084472656, 42.095001220703125],
          dtype='float64', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Index([-70.83499908447266, -70.82499694824219, -70.81499481201172,
           -70.80500030517578, -70.79499816894531, -70.78499603271484,
           -70.77499389648438, -70.76499938964844, -70.75499725341797,
            -70.7449951171875,
           ...
           -69.86499786376953, -69.85499572753906,  -69.8449935913086,
           -69.83499908447266, -69.82499694824219, -69.81499481201172,
           -69.80500030517578, -69.79499816894531, -69.78499603271484,
           -69.77499389648438],
          dtype='float64', name='lon', length=107))

• Attributes: (6)

  units :

  kg m-2

  standard_name :

  Total precipitation

  long_name :

  Total precipitation

  cell_methods :

  time: sum

  vmin :

  0.0

  vmax :

  60.617774963378906