Water Mission Area STAC Catalog

Water Mission Area STAC Catalog#

The USGS Water Resources Mission Area (WMA) has created a STAC Catalog to help users find and access data related to our water modeling projects. STAC Catalogs provide a mechanism to expose spatiotemporal datasets using a machine-readable format, allowing users to discover and access datasets in a standardized way. The WMA Catalog hosts a variety of datasets related to hydro-terrestrial modeling. The scope of the WMA STAC Catalog datasets may include:

  • common hydrologic model inputs such as climate or other forcing datasets

  • hydrologic model outputs

  • observational datasets related to hydrology and/or water budgets

We have exposed our STAC Catalog through a pygeoapi endpoint that is compliant with the OGC API suite of standards. This endpoint allows both API access, as well as a user interface for browsing our data assets.

This STAC Catalog is part of a modernized replacement for the legacy Geo Data Portal.

Finding the data you need and reading it into your workflow#

You have two ways to discover the data assets in the WMA STAC Catalog:

  1. Using the pygeoapi user interface to browse our datasets through a user interface.

  2. Reading the pygeoapi API endpoint into a workflow using a library designed for reading STAC catalogs, such as PySTAC.

We will describe/demonstrate both methods below.

# first import the packages we will need
import pystac
import xarray as xr
import zarr
from packaging.version import Version

Option 1: Find and access a dataset using the PyGEOAPI or Radiant Earth user interface#

Step 1: Explore the catalog through either our pygeoapi user interface or the Radiant Earth STAC Browser and find the dataset you want to use.

  • This catalog contains both datasets and collections of datasets in the top level of the catalog. Therefore, some assets that you click may open up directly to a dataset explore page (with a map viewer and metadata about the dataset), while others will open up to another catalog listing. As you drill down through the collections, you will eventually land on a dataset explore page.

  • You can learn more about each dataset at its source publication point, which is linked in the cite-as property displayed on each dataset page (if available).

  • Click through the catalog until you find a dataset that you want to use. For this example, we will use the gridMET dataset.

Step 2: Identify which asset you want to use to access the dataset you’ve chosen.

  • Review the Title and Description of each asset listed in the Assets section. Each asset is a different copy of the dataset that may be stored in a different format or location.

  • Typically, we recommend using the asset with the title “Free access to zarr via S3 API”.

  • Copy the URL, Open Keywords, and Storage Options for the asset you have chosen, and store them as python variables, demonstrated below:

# URL
zarr_url = 's3://hytest/conus404/conus404_hourly.zarr/'

# Open Keywords
# note that you will need to capitalize the True or False in "consolidated" - these keys are stored as lowercase boolens in the STAC catalog json
# but should be capitalized for python
open_keywords = {
  "chunks": {},
  "consolidated": True,
  "engine": "zarr"
}

# Storage Options
# note that you will need to capitalize the True or False in "anon" - these keys are stored as lowercase boolens in the STAC catalog json
# but should be capitalized for python
storage_options = {
  "anon": True,
  "client_kwargs": {
    "endpoint_url": "https://usgs.osn.mghpcc.org/"
  }
}

Open the dataset using xarray#

Now we will open and view the dataset using xarray.

Please note that there are two major zarr format specifications: 2 and 3. If you are using the python package zarr>=3.0.0, you must specify the format of the zarr store you are trying to open in the xarray.open_dataset function. If you are using zarr<3.0.0, you do not need to specify the format, as it will default to version 2.

if Version(zarr.__version__) < Version("3.0.0"):
    ds = xr.open_dataset(
        zarr_url,
        storage_options=storage_options,
        **open_keywords
    )
else:
    ds = xr.open_dataset(
    zarr_url,
    storage_options=storage_options,
    **open_keywords,
    zarr_format=2
    )
ds

<xarray.Dataset> Size: 222TB
Dimensions:         (time: 376945, y: 1015, x: 1367, bottom_top_stag: 51,
                     bottom_top: 50, soil_layers_stag: 4, x_stag: 1368,
                     y_stag: 1016, snow_layers_stag: 3, snso_layers_stag: 7)
Coordinates:
    lat             (y, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lat_u           (y, x_stag) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lat_v           (y_stag, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lon             (y, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lon_u           (y, x_stag) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lon_v           (y_stag, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
  * time            (time) datetime64[ns] 3MB 1979-10-01 ... 2022-10-01
  * x               (x) float64 11kB -2.732e+06 -2.728e+06 ... 2.732e+06
  * y               (y) float64 8kB -2.028e+06 -2.024e+06 ... 2.028e+06
Dimensions without coordinates: bottom_top_stag, bottom_top, soil_layers_stag,
                                x_stag, y_stag, snow_layers_stag,
                                snso_layers_stag
Data variables: (12/153)
    ACDEWC          (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACDRIPR         (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACDRIPS         (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACECAN          (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACEDIR          (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACETLSM         (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ...              ...
    ZNU             (bottom_top) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    ZNW             (bottom_top_stag) float32 204B dask.array<chunksize=(51,), meta=np.ndarray>
    ZS              (soil_layers_stag) float32 16B dask.array<chunksize=(4,), meta=np.ndarray>
    ZSNSO           (time, snso_layers_stag, y, x) float32 15TB dask.array<chunksize=(144, 7, 175, 175), meta=np.ndarray>
    ZWT             (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    crs             int64 8B ...
Attributes: (12/148)
    AER_ANGEXP_OPT:                  1
    AER_ANGEXP_VAL:                  1.2999999523162842
    AER_AOD550_OPT:                  1
    AER_AOD550_VAL:                  0.11999999731779099
    AER_ASY_OPT:                     1
    AER_ASY_VAL:                     0.8999999761581421
    ...                              ...
    WEST-EAST_PATCH_START_STAG:      1
    WEST-EAST_PATCH_START_UNSTAG:    1
    W_DAMPING:                       1
    YSU_TOPDOWN_PBLMIX:              0
    history:                         Tue Mar 29 16:35:22 2022: ncrcat -A -vW ...
    history_of_appended_files:       Tue Mar 29 16:35:22 2022: Appended file ...

Option 2: Find and access a dataset using the PySTAC python library#

Before we begin, we will define a helper function that can be used to drill down through the STAC Catalog and extract key metadata.

def get_children(catalog, collection_id=None):
    """
    This function retrieves a specified collection from a STAC catalog/collection and prints key metadata 
    for exploring/accessing the datasets contained within it.
    If there is no collection ID provided, the collections in the top level of the catalog will be printed.
    If a collection ID is provided, it will retrieve the collection with that ID from the input catalog/collection.
    If the collection ID points to a dataset, it will print the assets available for the dataset.
    If the collection ID points to another collection, it will list the child collections in the IDed collection.

    Args:
        catalog (pystac.Catalog | pystac.Collection): The STAC catalog/collection object.
        collection_id (str): The ID of the collection or dataset to retrieve from catalog.
    
    Returns:
        collection (pystac.Catalog | pystac.Collection): The collection object corresponding to the provided ID
                                                         or the top-level catalog if no ID is provided.
    """
    dataset = False
    if collection_id:
        collection = catalog.get_child(collection_id)
        if collection.assets:
            dataset = True
            print(f"{collection_id} is a dataset. Please review the assets below and select one to open.")

        else:
            print(f"{collection_id} is a collection. Please review the child items and select one to open in the next cell.")
    else:
        collection = catalog
    if dataset==True:
        # List the assets
        for asset in collection.assets:
            print(f"Asset ID: {asset}")
            print(f"    Title: {collection.assets[asset].title}")
            print(f"    Description: {collection.assets[asset].description}")
    else:
        collections = list(collection.get_collections())
        print(f"Number of collections: {len(collections)}")
        print("Collections IDs:")
        for child_collection in collections:
            id = child_collection.id
            cite_as = "Not available"
            for link in child_collection.links:
                if link.rel == "cite-as":
                    cite_as = link.target
            print(f"- {id}, Source: {cite_as}")
    return collection

Step 1: Explore the catalog using the PySTAC library and find the dataset you want to use.

  • The following code will read the WMA STAC Catalog and print the items in the top level of the catalog.

  • This catalog contains both datasets and collections of datasets in the top level of the catalog. Therefore, some collections may point to a dataset endpoint, while others will point to another collection of several datasets. As you drill down through the collections, you will eventually land on a dataset.

  • You can learn more about each dataset at its source publication point, which the helper function defined above prints out (if available).

  • Identify a dataset that you want to use. For this example, we will use the gridMET dataset.

Note: We are working on enabling stac-search to our pygeoapi endpoint, which will allow you to search for datasets using keywords and other metadata which will help facilitate this exploration process.

# url for the WMA STAC Catalog
catalog_url = "https://api.water.usgs.gov/gdp/pygeoapi/stac/stac-collection/"

# use pystac to read the catalog
catalog = pystac.Catalog.from_file(catalog_url)

# list the collections in the catalog
catalog = get_children(catalog)

Number of collections: 43
Collections IDs:
- AIEM_permafrost, Source: Not available
- CA-BCM-2014, Source: https://doi.org/10.21429/dye8-h568
- FLET, Source: Not available
- GMO, Source: https://doi.org/10.21429/v7ys-6n72
- GMO_New, Source: https://doi.org/10.21429/c6s4-ve70
- LOCA2, Source: Not available
- PRISM_v2, Source: https://prism.oregonstate.edu/
- PuertoRico, Source: Not available
- RedRiver, Source: Not available
- SPEI, Source: Not available
- TTU_2019, Source: Not available
- TopoWx2017, Source: Not available
- WUS_HSP, Source: Not available
- alaska_et_2020, Source: Not available
- bcca, Source: Not available
- bcsd_mon_vic, Source: http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/
- bcsd_obs, Source: http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/
- cmip5_bcsd, Source: https://doi.org/10.21429/sbxv-1n90
- conus404-pgw, Source: Not available
- conus404, Source: Not available
- cooper, Source: Not available
- cprep, Source: Not available
- dcp_compressed, Source: https://doi.org/https%3A//doi.org/10.21429/j9f1-b218
- gridMET, Source: https://www.climatologylab.org/gridmet.html
- hawaii_2018, Source: Not available
- iclus, Source: Not available
- loca, Source: Not available
- maca-vic, Source: Not available
- macav2, Source: Not available
- maurer, Source: https://doi.org/10.21429/m7y0-xy02
- mows, Source: Not available
- nc-casc-snow, Source: Not available
- nlcd, Source: Not available
- notaro_2018, Source: Not available
- pacis, Source: Not available
- puerto_rico, Source: Not available
- red_river_2018, Source: https://doi.org/10.21429/em59-hn43
- serap, Source: Not available
- slr2d, Source: https://doi.org/10.21429/66gt-dm26
- snodas, Source: https://doi.org/10.7265/N5TB14TC
- ssebopeta, Source: Not available
- stageiv_combined, Source: https://pubs.usgs.gov/fs/2013/3035/
- wicci, Source: https://doi.org/10.21429/dtp5-z505

# select a collection from the catalog, replace the collection ID with the one you want to use:
collection = get_children(catalog, collection_id="conus404")

conus404 is a collection. Please review the child items and select one to open in the next cell.

Number of collections: 4
Collections IDs:
- conus404_daily, Source: https://doi.org/10.5066/P9PHPK4F
- conus404_hourly, Source: https://doi.org/10.5066/P9PHPK4F
- conus404_monthly, Source: https://doi.org/10.5066/P9PHPK4F
- conus404_xtrm_daily, Source: https://doi.org/10.5066/P9PHPK4F

# if the above was a collection, uncomment the line below and enter the collection ID 
# you want to use from the parent collection you selected:
collection = get_children(collection, collection_id="conus404_hourly")

conus404_hourly is a dataset. Please review the assets below and select one to open.
Asset ID: zarr-s3-osn
    Title: Free access to zarr via S3 API
    Description: Free, public access to zarr data store via the S3 API. This data is stored on an Open Storage Network Pod.
Asset ID: zarr-disk-hovenweep
    Title: USGS internal access to zarr via Hovenweep
    Description: USGS internal access to zarr data store via disk storage attached to the Hovenweep supercomputer.
Asset ID: default
    Title: None
    Description: None

Step 2: Identify which asset you want to use to access the dataset you’ve chosen.

  • Review the Title and Description of each asset printed above. Each asset is a different copy of the dataset that may be stored in a different format or location.

  • Typically, we recommend people use the zarr-s3-osn asset unless they have a reason to use a different asset.

  • Copy the Asset ID for the asset you have chosen, and paste in the code below to read the asset metadata.

# replace with the asset ID you want to use:
selected_asset_id = "zarr-s3-osn"

# read the asset metadata
asset = collection.assets[selected_asset_id]

Open the dataset using xarray#

Now we will open and view the dataset using xarray.

Please note that there are two major zarr format specifications: 2 and 3. If you are using the python package zarr>=3.0.0, you must specify the format of the zarr store you are trying to open in the xarray.open_dataset function. If you are using zarr<3.0.0, you do not need to specify the format, as it will default to version 2.

if Version(zarr.__version__) < Version("3.0.0"):
    ds = xr.open_dataset(
        asset.href,
        storage_options=asset.extra_fields['xarray:storage_options'],
        **asset.extra_fields['xarray:open_kwargs']
    )
else:
    ds = xr.open_dataset(
    asset.href,
    storage_options=asset.extra_fields['xarray:storage_options'],
    **asset.extra_fields['xarray:open_kwargs'],
    zarr_format=2
    )
ds

<xarray.Dataset> Size: 222TB
Dimensions:         (time: 376945, y: 1015, x: 1367, bottom_top_stag: 51,
                     bottom_top: 50, soil_layers_stag: 4, x_stag: 1368,
                     y_stag: 1016, snow_layers_stag: 3, snso_layers_stag: 7)
Coordinates:
    lat             (y, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lat_u           (y, x_stag) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lat_v           (y_stag, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lon             (y, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lon_u           (y, x_stag) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
    lon_v           (y_stag, x) float32 6MB dask.array<chunksize=(175, 175), meta=np.ndarray>
  * time            (time) datetime64[ns] 3MB 1979-10-01 ... 2022-10-01
  * x               (x) float64 11kB -2.732e+06 -2.728e+06 ... 2.732e+06
  * y               (y) float64 8kB -2.028e+06 -2.024e+06 ... 2.028e+06
Dimensions without coordinates: bottom_top_stag, bottom_top, soil_layers_stag,
                                x_stag, y_stag, snow_layers_stag,
                                snso_layers_stag
Data variables: (12/153)
    ACDEWC          (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACDRIPR         (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACDRIPS         (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACECAN          (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACEDIR          (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ACETLSM         (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    ...              ...
    ZNU             (bottom_top) float32 200B dask.array<chunksize=(50,), meta=np.ndarray>
    ZNW             (bottom_top_stag) float32 204B dask.array<chunksize=(51,), meta=np.ndarray>
    ZS              (soil_layers_stag) float32 16B dask.array<chunksize=(4,), meta=np.ndarray>
    ZSNSO           (time, snso_layers_stag, y, x) float32 15TB dask.array<chunksize=(144, 7, 175, 175), meta=np.ndarray>
    ZWT             (time, y, x) float32 2TB dask.array<chunksize=(144, 175, 175), meta=np.ndarray>
    crs             int64 8B ...
Attributes: (12/148)
    AER_ANGEXP_OPT:                  1
    AER_ANGEXP_VAL:                  1.2999999523162842
    AER_AOD550_OPT:                  1
    AER_AOD550_VAL:                  0.11999999731779099
    AER_ASY_OPT:                     1
    AER_ASY_VAL:                     0.8999999761581421
    ...                              ...
    WEST-EAST_PATCH_START_STAG:      1
    WEST-EAST_PATCH_START_UNSTAG:    1
    W_DAMPING:                       1
    YSU_TOPDOWN_PBLMIX:              0
    history:                         Tue Mar 29 16:35:22 2022: ncrcat -A -vW ...
    history_of_appended_files:       Tue Mar 29 16:35:22 2022: Appended file ...

You can use whichever of the two methods described above to find and access datasets from our STAC Catalog. Now that you have the dataset open in xarray, you can build on the workflow to visualize, analyze, or do any other data processing you might need to do. You don’t need to download your own copy of the dataset - you can perform your analysis directly on the dataset from the source we provide. If your workflow exceeds your computer’s memory, you can use the dask library to parallelize your analysis and take advantage of HPC or scalable cloud computing resources.

Reporting Issues with the Catalog#

If you find any issues with our STAC catalog or the datasets it contains, please reach out to us in one of the following ways:

  • If you have a code.usgs.gov account, you can open an issue

  • E-mail us at mdmf@usgs.gov

We will do our best to address any issues you find in a timely manner.

Contributing to the Catalog#

We are very interested in expanding the datasets available in our STAC Catalog to include more datasets that are relevant to water modeling. If you have a dataset that you think should be included, we would love to hear from you! We are particularly interested in datasets that are very large or difficult to access. Please open an issue if you think you have a dataset that should be included in the catalog and tell us a bit about the dataset, who is using it, and why you would like to add it to our STAC Catalog.

Contents