Tip

For an interactive online version click here:

Example: Reading raster data

This example illustrates the how to read raster data using the HydroMT DataCatalog with the raster, netcdf and raster_tindex drivers.

# import hydromt and setup logging
import hydromt
from hydromt.log import setuplog

logger = setuplog("read raster data", log_level=10)

2024-02-26 15:12:09,517 - read raster data - log - INFO - HydroMT version: 0.9.4

# Download artifacts for the Piave basin to `~/.hydromt_data/`.
data_catalog = hydromt.DataCatalog(logger=logger, data_libs=["artifact_data"])

2024-02-26 15:12:09,543 - read raster data - data_catalog - INFO - Reading data catalog archive artifact_data v0.0.8

2024-02-26 15:12:09,544 - read raster data - data_catalog - INFO - Parsing data catalog from /home/runner/.hydromt_data/artifact_data/v0.0.8/data_catalog.yml

Raster driver

To read raster data and parse it into a xarray Dataset or DataArray we use the open_mfraster() method. All driver_kwargs in the data catalog yaml file will be passed to this method. The raster driver supports all GDAL data formats, including the often used GeoTiff of Cloud Optimized GeoTiff (COG) formats. Tiled datasets can also be passed as a virtual raster tileset (vrt) file.

As an example we will use the MERIT Hydro dataset which is a set of GeoTiff files with identical grids, one for each variable of the datasets.

# inspect data source entry in data catalog yaml file
data_catalog["merit_hydro"]

crs: 4326
data_type: RasterDataset
driver: raster
meta:
  category: topography
  paper_doi: 10.1029/2019WR024873
  paper_ref: Yamazaki et al. (2019)
  source_license: CC-BY-NC 4.0 or ODbL 1.0
  source_url: http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro
  source_version: 1.0
path: /home/runner/.hydromt_data/artifact_data/v0.0.8/merit_hydro/{variable}.tif

We can load any RasterDataset using DataCatalog.get_rasterdataset(). Note that if we don’t provide any arguments it returns the full dataset with nine data variables and for the full spatial domain. Only the data coordinates are actually read, the data variables are still lazy Dask arrays.

ds = data_catalog.get_rasterdataset("merit_hydro")
ds

2024-02-26 15:12:09,605 - read raster data - rasterdataset - INFO - Reading merit_hydro raster data from /home/runner/.hydromt_data/artifact_data/v0.0.8/merit_hydro/{variable}.tif

<xarray.Dataset> Size: 97MB
Dimensions:      (x: 1680, y: 1920)
Coordinates:
  * x            (x) float64 13kB 11.6 11.6 11.6 11.6 ... 13.0 13.0 13.0 13.0
  * y            (y) float64 15kB 46.8 46.8 46.8 46.8 ... 45.2 45.2 45.2 45.2
    spatial_ref  int64 8B 0
Data variables:
    basins       (y, x) uint32 13MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    elevtn       (y, x) float32 13MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    flwdir       (y, x) uint8 3MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    hnd          (y, x) float32 13MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    lndslp       (y, x) float32 13MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    rivwth       (y, x) float32 13MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    strord       (y, x) uint8 3MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    uparea       (y, x) float32 13MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
    upgrid       (y, x) int32 13MB dask.array<chunksize=(1920, 1680), meta=np.ndarray>
Attributes:
    category:        topography
    paper_doi:       10.1029/2019WR024873
    paper_ref:       Yamazaki et al. (2019)
    source_license:  CC-BY-NC 4.0 or ODbL 1.0
    source_url:      http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro
    source_version:  1.0

xarray.Dataset

• Dimensions:
  • x: 1680
  • y: 1920
• Coordinates: (3)
  • x
    (x)
    float64
    11.6 11.6 11.6 ... 13.0 13.0 13.0

    array([11.600417, 11.60125 , 11.602083, ..., 12.997917, 12.99875 , 12.999583])

  • y
    (y)
    float64
    46.8 46.8 46.8 ... 45.2 45.2 45.2

    array([46.799583, 46.79875 , 46.797917, ..., 45.202083, 45.20125 , 45.200417])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    GeoTransform :

    11.600000000000023 0.0008333333333333203 0.0 46.8 0.0 -0.0008333333333333345

    x_dim :

    x

    y_dim :

    y

    array(0)

• Data variables: (9)
  • basins
    (y, x)
    uint32
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    0

    source_file :

    basins.tif

    Array Chunk Bytes 12.30 MiB 12.30 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type uint32 numpy.ndarray

  • elevtn
    (y, x)
    float32
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    -9999.0

    source_file :

    elevtn.tif

    Array Chunk Bytes 12.30 MiB 12.30 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type float32 numpy.ndarray

  • flwdir
    (y, x)
    uint8
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    247

    source_file :

    flwdir.tif

    Array Chunk Bytes 3.08 MiB 3.08 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type uint8 numpy.ndarray

  • hnd
    (y, x)
    float32
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    -9999.0

    source_file :

    hnd.tif

    Array Chunk Bytes 12.30 MiB 12.30 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type float32 numpy.ndarray

  • lndslp
    (y, x)
    float32
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    -9999.0

    source_file :

    lndslp.tif

    Array Chunk Bytes 12.30 MiB 12.30 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type float32 numpy.ndarray

  • rivwth
    (y, x)
    float32
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    -9999.0

    source_file :

    rivwth.tif

    Array Chunk Bytes 12.30 MiB 12.30 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type float32 numpy.ndarray

  • strord
    (y, x)
    uint8
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    255

    source_file :

    strord.tif

    Array Chunk Bytes 3.08 MiB 3.08 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type uint8 numpy.ndarray

  • uparea
    (y, x)
    float32
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    -9999.0

    source_file :

    uparea.tif

    Array Chunk Bytes 12.30 MiB 12.30 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type float32 numpy.ndarray

  • upgrid
    (y, x)
    int32
    dask.array<chunksize=(1920, 1680), meta=np.ndarray>

    AREA_OR_POINT :

    Area

    _FillValue :

    -9999

    source_file :

    upgrid.tif

    Array Chunk Bytes 12.30 MiB 12.30 MiB Shape (1920, 1680) (1920, 1680) Dask graph 1 chunks in 3 graph layers Data type int32 numpy.ndarray

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([11.600416666666689, 11.601250000000022, 11.602083333333356,
           11.602916666666689, 11.603750000000023, 11.604583333333355,
            11.60541666666669, 11.606250000000022, 11.607083333333355,
            11.60791666666669,
           ...
           12.992083333333333, 12.992916666666666,           12.99375,
           12.994583333333333, 12.995416666666667,           12.99625,
           12.997083333333332, 12.997916666666667,           12.99875,
           12.999583333333334],
          dtype='float64', name='x', length=1680))

  • y
    PandasIndex

    PandasIndex(Index([ 46.79958333333333,           46.79875, 46.797916666666666,
            46.79708333333333,           46.79625,  46.79541666666666,
            46.79458333333333, 46.793749999999996,  46.79291666666666,
            46.79208333333333,
           ...
            45.20791666666666,  45.20708333333333,           45.20625,
           45.205416666666665,  45.20458333333333,           45.20375,
            45.20291666666666,  45.20208333333333, 45.201249999999995,
            45.20041666666666],
          dtype='float64', name='y', length=1920))

• Attributes: (6)

  category :

  topography

  paper_doi :

  10.1029/2019WR024873

  paper_ref :

  Yamazaki et al. (2019)

  source_license :

  CC-BY-NC 4.0 or ODbL 1.0

  source_url :

  http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro

  source_version :

  1.0

The data can be visualized with the .plot() xarray method. We replace all nodata values with NaNs with .raster.mask_nodata().

ds["elevtn"].raster.mask_nodata().plot(cmap="terrain")

<matplotlib.collections.QuadMesh at 0x7f1be96d5110>

We can request a (spatial) subset data by providing additional variables and bbox / geom arguments. Note that these return a smaller spatial extent and just two data variables. The variables argument is especially useful if each variable of the dataset is saved in a separate file and the {variable} key is used in the path argument of the data source (see above) to limit which files are actually read. If a single variable is requested a DataArray instead of a Dataset is returned unless the single_var_as_array argument is set to False (True by default).

bbox = [11.70, 45.35, 12.95, 46.70]
ds = data_catalog.get_rasterdataset(
    "merit_hydro", bbox=bbox, variables=["elevtn"], single_var_as_array=True
)
ds

2024-02-26 15:12:12,172 - read raster data - rasterdataset - INFO - Reading merit_hydro raster data from /home/runner/.hydromt_data/artifact_data/v0.0.8/merit_hydro/{variable}.tif

2024-02-26 15:12:12,188 - read raster data - rasterdataset - DEBUG - Clip to [11.700, 45.350, 12.950, 46.700] (epsg:4326))

<xarray.DataArray 'elevtn' (y: 1620, x: 1500)> Size: 10MB
dask.array<getitem, shape=(1620, 1500), dtype=float32, chunksize=(1620, 1500), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 12kB 11.7 11.7 11.7 11.7 ... 12.95 12.95 12.95
  * y            (y) float64 13kB 46.7 46.7 46.7 46.7 ... 45.35 45.35 45.35
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     -9999.0
    source_file:    elevtn.tif

xarray.DataArray

'elevtn'

• y: 1620
• x: 1500

• dask.array<chunksize=(1620, 1500), meta=np.ndarray>

  Array Chunk Bytes 9.27 MiB 9.27 MiB Shape (1620, 1500) (1620, 1500) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

• Coordinates: (3)
  • x
    (x)
    float64
    11.7 11.7 11.7 ... 12.95 12.95

    array([11.700417, 11.70125 , 11.702083, ..., 12.947917, 12.94875 , 12.949583])

  • y
    (y)
    float64
    46.7 46.7 46.7 ... 45.35 45.35

    array([46.699583, 46.69875 , 46.697917, ..., 45.352083, 45.35125 , 45.350417])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    GeoTransform :

    11.600000000000023 0.0008333333333333203 0.0 46.8 0.0 -0.0008333333333333345

    x_dim :

    x

    y_dim :

    y

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([11.700416666666687, 11.701250000000021, 11.702083333333354,
           11.702916666666688,  11.70375000000002, 11.704583333333353,
           11.705416666666688,  11.70625000000002, 11.707083333333355,
           11.707916666666687,
           ...
           12.942083333333334, 12.942916666666667, 12.943750000000001,
           12.944583333333334, 12.945416666666667, 12.946250000000001,
           12.947083333333333, 12.947916666666668,           12.94875,
           12.949583333333333],
          dtype='float64', name='x', length=1500))

  • y
    PandasIndex

    PandasIndex(Index([ 46.69958333333333,           46.69875, 46.697916666666664,
            46.69708333333333,           46.69625,  46.69541666666667,
           46.694583333333334, 46.693749999999994,  46.69291666666666,
            46.69208333333333,
           ...
            45.35791666666666,  45.35708333333333, 45.356249999999996,
            45.35541666666666,  45.35458333333333,           45.35375,
           45.352916666666665,  45.35208333333333,  45.35124999999999,
            45.35041666666666],
          dtype='float64', name='y', length=1620))

• Attributes: (3)

  AREA_OR_POINT :

  Area

  _FillValue :

  -9999.0

  source_file :

  elevtn.tif

As mentioned earlier, all driver_kwargs in the data entry of the catalog yaml file will be passed to the open_mfraster() method. Here we show how these arguments can be used to concatenate multiple raster files along a new dimension.

# update data source entry (this is normally done manually before initializing the data catalog!)
data_catalog["merit_hydro"].driver_kwargs.update(concat_dim="variable", concat=True)
data_catalog["merit_hydro"]

crs: 4326
data_type: RasterDataset
driver: raster
driver_kwargs:
  concat: true
  concat_dim: variable
meta:
  category: topography
  paper_doi: 10.1029/2019WR024873
  paper_ref: Yamazaki et al. (2019)
  source_license: CC-BY-NC 4.0 or ODbL 1.0
  source_url: http://hydro.iis.u-tokyo.ac.jp/~yamadai/MERIT_Hydro
  source_version: 1.0
path: /home/runner/.hydromt_data/artifact_data/v0.0.8/merit_hydro/{variable}.tif

# this returns a DataArray (single variable) with a new 'variable' dimension
da = data_catalog.get_rasterdataset("merit_hydro")
da

2024-02-26 15:12:12,207 - read raster data - rasterdataset - INFO - Reading merit_hydro raster data from /home/runner/.hydromt_data/artifact_data/v0.0.8/merit_hydro/{variable}.tif

<xarray.DataArray 'basins' (variable: 9, y: 1920, x: 1680)> Size: 232MB
dask.array<getitem, shape=(9, 1920, 1680), dtype=float64, chunksize=(1, 1920, 1680), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 13kB 11.6 11.6 11.6 11.6 ... 13.0 13.0 13.0 13.0
  * y            (y) float64 15kB 46.8 46.8 46.8 46.8 ... 45.2 45.2 45.2 45.2
  * variable     (variable) int64 72B 0 1 2 3 4 5 6 7 8
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:  Area
    _FillValue:     0

xarray.DataArray

'basins'

• variable: 9
• y: 1920
• x: 1680

• dask.array<chunksize=(1, 1920, 1680), meta=np.ndarray>

  Array Chunk Bytes 221.48 MiB 24.61 MiB Shape (9, 1920, 1680) (1, 1920, 1680) Dask graph 9 chunks in 47 graph layers Data type float64 numpy.ndarray

• Coordinates: (4)
  • x
    (x)
    float64
    11.6 11.6 11.6 ... 13.0 13.0 13.0

    array([11.600417, 11.60125 , 11.602083, ..., 12.997917, 12.99875 , 12.999583])

  • y
    (y)
    float64
    46.8 46.8 46.8 ... 45.2 45.2 45.2

    array([46.799583, 46.79875 , 46.797917, ..., 45.202083, 45.20125 , 45.200417])

  • variable
    (variable)
    int64
    0 1 2 3 4 5 6 7 8

    array([0, 1, 2, 3, 4, 5, 6, 7, 8])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    GeoTransform :

    11.600000000000023 0.0008333333333333203 0.0 46.8 0.0 -0.0008333333333333345

    x_dim :

    x

    y_dim :

    y

    array(0)

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Index([11.600416666666689, 11.601250000000022, 11.602083333333356,
           11.602916666666689, 11.603750000000023, 11.604583333333355,
            11.60541666666669, 11.606250000000022, 11.607083333333355,
            11.60791666666669,
           ...
           12.992083333333333, 12.992916666666666,           12.99375,
           12.994583333333333, 12.995416666666667,           12.99625,
           12.997083333333332, 12.997916666666667,           12.99875,
           12.999583333333334],
          dtype='float64', name='x', length=1680))

  • y
    PandasIndex

    PandasIndex(Index([ 46.79958333333333,           46.79875, 46.797916666666666,
            46.79708333333333,           46.79625,  46.79541666666666,
            46.79458333333333, 46.793749999999996,  46.79291666666666,
            46.79208333333333,
           ...
            45.20791666666666,  45.20708333333333,           45.20625,
           45.205416666666665,  45.20458333333333,           45.20375,
            45.20291666666666,  45.20208333333333, 45.201249999999995,
            45.20041666666666],
          dtype='float64', name='y', length=1920))

  • variable
    PandasIndex

    PandasIndex(Index([0, 1, 2, 3, 4, 5, 6, 7, 8], dtype='int64', name='variable'))

• Attributes: (2)

  AREA_OR_POINT :

  Area

  _FillValue :

  0

TIP: To write a dataset back to a stack of raster in a single folder use the .raster.to_mapstack() method.

Netcdf driver

Many gridded datasets with a third dimension (e.g. time) are saved in netcdf or zarr files, which can be read with the netcdf and zarr drivers respectively. This data is read using the xarray.open_mfdataset() method. These formats are flexible and therefore HydroMT is not always able to read the geospatial attributes such as the CRS from the data and it has to be set through the data catalog yaml file.

If the data is stored per year or month, the {year} and {month} keys can be used in the path argument of a data source in the data catalog yaml file to speed up the reading of a temporal subset of the data using the date_tuple argument of DataCatalog.get_rasterdataset() (not in this example).

As example we use the ERA5 dataset.

# Note the crs argument as this is missing in the original data
data_catalog["era5_hourly"]

crs: 4326
data_type: RasterDataset
driver: netcdf
meta:
  category: meteo
  history: Extracted from Copernicus Climate Data Store
  paper_doi: 10.1002/qj.3803
  paper_ref: Hersbach et al. (2019)
  source_license: https://cds.climate.copernicus.eu/cdsapp/#!/terms/licence-to-use-copernicus-products
  source_url: https://doi.org/10.24381/cds.bd0915c6
  source_version: ERA5 hourly data on pressure levels
path: /home/runner/.hydromt_data/artifact_data/v0.0.8/era5_hourly.nc
unit_add:
  temp: -273.15
unit_mult:
  kin: 0.000277778
  kout: 0.000277778
  precip: 1000
  press_msl: 0.01

# Note that the some units are converted
ds = data_catalog.get_rasterdataset("era5_hourly")
ds

2024-02-26 15:12:12,322 - read raster data - rasterdataset - INFO - Reading era5_hourly netcdf data from /home/runner/.hydromt_data/artifact_data/v0.0.8/era5_hourly.nc

2024-02-26 15:12:12,356 - read raster data - rasterdataset - DEBUG - Convert units for 5 variables.

<xarray.Dataset> Size: 285kB
Dimensions:      (longitude: 6, latitude: 7, time: 336)
Coordinates:
  * longitude    (longitude) float32 24B 11.75 12.0 12.25 12.5 12.75 13.0
  * latitude     (latitude) float32 28B 46.75 46.5 46.25 46.0 45.75 45.5 45.25
  * time         (time) datetime64[ns] 3kB 2010-02-01 ... 2010-02-14T23:00:00
    spatial_ref  int32 4B ...
Data variables:
    press_msl    (time, latitude, longitude) float32 56kB dask.array<chunksize=(336, 7, 6), meta=np.ndarray>
    kin          (time, latitude, longitude) float32 56kB dask.array<chunksize=(336, 7, 6), meta=np.ndarray>
    temp         (time, latitude, longitude) float32 56kB dask.array<chunksize=(336, 7, 6), meta=np.ndarray>
    kout         (time, latitude, longitude) float32 56kB dask.array<chunksize=(336, 7, 6), meta=np.ndarray>
    precip       (time, latitude, longitude) float32 56kB dask.array<chunksize=(336, 7, 6), meta=np.ndarray>
Attributes:
    Conventions:     CF-1.6
    history:         Extracted from Copernicus Climate Data Store
    category:        meteo
    paper_doi:       10.1002/qj.3803
    paper_ref:       Hersbach et al. (2019)
    source_license:  https://cds.climate.copernicus.eu/cdsapp/#!/terms/licenc...
    source_url:      https://doi.org/10.24381/cds.bd0915c6
    source_version:  ERA5 hourly data on pressure levels

xarray.Dataset

• Dimensions:
  • longitude: 6
  • latitude: 7
  • time: 336
• Coordinates: (4)
  • longitude
    (longitude)
    float32
    11.75 12.0 12.25 12.5 12.75 13.0

    array([11.75, 12.  , 12.25, 12.5 , 12.75, 13.  ], dtype=float32)

  • latitude
    (latitude)
    float32
    46.75 46.5 46.25 ... 45.5 45.25

    array([46.75, 46.5 , 46.25, 46.  , 45.75, 45.5 , 45.25], dtype=float32)

  • time
    (time)
    datetime64[ns]
    2010-02-01 ... 2010-02-14T23:00:00

    array(['2010-02-01T00:00:00.000000000', '2010-02-01T01:00:00.000000000',
           '2010-02-01T02:00:00.000000000', ..., '2010-02-14T21:00:00.000000000',
           '2010-02-14T22:00:00.000000000', '2010-02-14T23:00:00.000000000'],
          dtype='datetime64[ns]')

  • spatial_ref
    ()
    int32
    ...

    crs_wkt :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    x_dim :

    longitude

    y_dim :

    latitude

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]

    [1 values with dtype=int32]

• Data variables: (5)
  • press_msl
    (time, latitude, longitude)
    float32
    dask.array<chunksize=(336, 7, 6), meta=np.ndarray>

    units :

    Pa

    long_name :

    Mean sea level pressure

    standard_name :

    air_pressure_at_mean_sea_level

    _FillValue :

    nan

    Array Chunk Bytes 55.12 kiB 55.12 kiB Shape (336, 7, 6) (336, 7, 6) Dask graph 1 chunks in 7 graph layers Data type float32 numpy.ndarray

  • kin
    (time, latitude, longitude)
    float32
    dask.array<chunksize=(336, 7, 6), meta=np.ndarray>

    units :

    J m**-2

    long_name :

    Surface solar radiation downwards

    standard_name :

    surface_downwelling_shortwave_flux_in_air

    _FillValue :

    nan

    Array Chunk Bytes 55.12 kiB 55.12 kiB Shape (336, 7, 6) (336, 7, 6) Dask graph 1 chunks in 7 graph layers Data type float32 numpy.ndarray

  • temp
    (time, latitude, longitude)
    float32
    dask.array<chunksize=(336, 7, 6), meta=np.ndarray>

    units :

    K

    long_name :

    2 metre temperature

    _FillValue :

    nan

    Array Chunk Bytes 55.12 kiB 55.12 kiB Shape (336, 7, 6) (336, 7, 6) Dask graph 1 chunks in 7 graph layers Data type float32 numpy.ndarray

  • kout
    (time, latitude, longitude)
    float32
    dask.array<chunksize=(336, 7, 6), meta=np.ndarray>

    units :

    J m**-2

    long_name :

    TOA incident solar radiation

    _FillValue :

    nan

    Array Chunk Bytes 55.12 kiB 55.12 kiB Shape (336, 7, 6) (336, 7, 6) Dask graph 1 chunks in 7 graph layers Data type float32 numpy.ndarray

  • precip
    (time, latitude, longitude)
    float32
    dask.array<chunksize=(336, 7, 6), meta=np.ndarray>

    units :

    m

    long_name :

    Total precipitation

    _FillValue :

    nan

    Array Chunk Bytes 55.12 kiB 55.12 kiB Shape (336, 7, 6) (336, 7, 6) Dask graph 1 chunks in 7 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • longitude
    PandasIndex

    PandasIndex(Index([11.75, 12.0, 12.25, 12.5, 12.75, 13.0], dtype='float32', name='longitude'))

  • latitude
    PandasIndex

    PandasIndex(Index([46.75, 46.5, 46.25, 46.0, 45.75, 45.5, 45.25], dtype='float32', name='latitude'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2010-02-01 00:00:00', '2010-02-01 01:00:00',
                   '2010-02-01 02:00:00', '2010-02-01 03:00:00',
                   '2010-02-01 04:00:00', '2010-02-01 05:00:00',
                   '2010-02-01 06:00:00', '2010-02-01 07:00:00',
                   '2010-02-01 08:00:00', '2010-02-01 09:00:00',
                   ...
                   '2010-02-14 14:00:00', '2010-02-14 15:00:00',
                   '2010-02-14 16:00:00', '2010-02-14 17:00:00',
                   '2010-02-14 18:00:00', '2010-02-14 19:00:00',
                   '2010-02-14 20:00:00', '2010-02-14 21:00:00',
                   '2010-02-14 22:00:00', '2010-02-14 23:00:00'],
                  dtype='datetime64[ns]', name='time', length=336, freq=None))

• Attributes: (8)

  Conventions :

  CF-1.6

  history :

  Extracted from Copernicus Climate Data Store

  category :

  meteo

  paper_doi :

  10.1002/qj.3803

  paper_ref :

  Hersbach et al. (2019)

  source_license :

  https://cds.climate.copernicus.eu/cdsapp/#!/terms/licence-to-use-copernicus-products

  source_url :

  https://doi.org/10.24381/cds.bd0915c6

  source_version :

  ERA5 hourly data on pressure levels

Raster_tindex driver

If the raster data is tiled but for each tile a different CRS is used (for instance a different UTM projection for each UTM zone), this dataset cannot be described using a VRT file. In this case a vector file can be build to use a raster tile index using gdaltindex and read using open_raster_from_tindex(). To read the data into a single xarray.Dataset the data needs to be reprojected and mosaiced to a single CRS while reading. As this type of data cannot be loaded lazily the method is typically used with an area of interest for which the data is loaded and combined.

As example we use the GRWL mask raster tiles for which we have created a tileindex using the aforementioned gdaltindex command line tool. Note that the path points to the GeoPackage output of the gdaltindex tool.

data_catalog["grwl_mask"]

data_type: RasterDataset
driver: raster_tindex
driver_kwargs:
  chunks:
    x: 3000
    y: 3000
  mosaic_kwargs:
    method: nearest
  tileindex: location
meta:
  paper_doi: 10.1126/science.aat0636
  paper_ref: Allen and Pavelsky (2018)
  source_license: CC BY 4.0
  source_url: https://doi.org/10.5281/zenodo.1297434
  source_version: 1.01
nodata: 0
path: /home/runner/.hydromt_data/artifact_data/v0.0.8/grwl_tindex.gpkg

# the tileindex is a GeoPackage vector file
# with an attribute column 'location' (see also the tileindex argument under driver_kwargs) containing the (relative) paths to the raster file data
import geopandas as gpd

fn_tindex = data_catalog["grwl_mask"].path
print(fn_tindex)
gpd.read_file(fn_tindex, rows=5)

/home/runner/.hydromt_data/artifact_data/v0.0.8/grwl_tindex.gpkg

	location	EPSG	geometry
0	GRWL_mask_V01.01//NA17.tif	EPSG:32617	POLYGON ((-80.92667 4.00635, -77.99333 4.00082...
1	GRWL_mask_V01.01//NA18.tif	EPSG:32618	POLYGON ((-78.00811 4.00082, -71.99333 4.00082...
2	GRWL_mask_V01.01//NA19.tif	EPSG:32619	POLYGON ((-72.00811 4.00082, -65.99333 4.00082...
3	GRWL_mask_V01.01//NA20.tif	EPSG:32620	POLYGON ((-66.00811 4.00082, -59.99333 4.00082...
4	GRWL_mask_V01.01//NA21.tif	EPSG:32621	POLYGON ((-60.00811 4.00082, -53.99333 4.00082...

# this returns a DataArray (single variable) wit a mosaic of several files (see source_file attribute)
ds = data_catalog.get_rasterdataset("grwl_mask", bbox=bbox)
ds

2024-02-26 15:12:12,447 - read raster data - rasterdataset - INFO - Reading grwl_mask raster_tindex data from /home/runner/.hydromt_data/artifact_data/v0.0.8/grwl_tindex.gpkg

2024-02-26 15:12:13,545 - read raster data - rasterdataset - DEBUG - Clip to [241499.448, 5023876.339, 343279.303, 5177117.241] (epsg:32633))

<xarray.DataArray 'grwl_tindex' (y: 5108, x: 3393)> Size: 17MB
array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)
Coordinates:
  * y            (y) float64 41kB 5.177e+06 5.177e+06 ... 5.024e+06 5.024e+06
  * x            (x) float64 27kB 2.415e+05 2.415e+05 ... 3.432e+05 3.433e+05
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:       Area
    DataType:            Generic
    RepresentationType:  ATHEMATIC
    _FillValue:          0
    source_file:         NL33.tif; NL32.tif

xarray.DataArray

'grwl_tindex'

• y: 5108
• x: 3393

• 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

  array([[0, 0, 0, ..., 0, 0, 0],
         [0, 0, 0, ..., 0, 0, 0],
         [0, 0, 0, ..., 0, 0, 0],
         ...,
         [0, 0, 0, ..., 0, 0, 0],
         [0, 0, 0, ..., 0, 0, 0],
         [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)

• Coordinates: (3)
  • y
    (y)
    float64
    5.177e+06 5.177e+06 ... 5.024e+06

    axis :

    Y

    long_name :

    Northing

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([5177115., 5177085., 5177055., ..., 5023965., 5023935., 5023905.])

  • x
    (x)
    float64
    2.415e+05 2.415e+05 ... 3.433e+05

    axis :

    X

    long_name :

    Easting

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([241515., 241545., 241575., ..., 343215., 343245., 343275.])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCRS["UTM_Zone_33N",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 33N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",15,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1]],ID["EPSG",32633]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    UTM_Zone_33N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    15.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCRS["UTM_Zone_33N",BASEGEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 33N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",15,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["easting",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["northing",north,ORDER[2],LENGTHUNIT["metre",1]],ID["EPSG",32633]]

    GeoTransform :

    241440.0 30.0 0.0 5177190.0 0.0 -30.0

    x_dim :

    x

    y_dim :

    y

    array(0)

• Indexes: (2)
  • y
    PandasIndex

    PandasIndex(Index([5177115.0, 5177085.0, 5177055.0, 5177025.0, 5176995.0, 5176965.0,
           5176935.0, 5176905.0, 5176875.0, 5176845.0,
           ...
           5024175.0, 5024145.0, 5024115.0, 5024085.0, 5024055.0, 5024025.0,
           5023995.0, 5023965.0, 5023935.0, 5023905.0],
          dtype='float64', name='y', length=5108))

  • x
    PandasIndex

    PandasIndex(Index([241515.0, 241545.0, 241575.0, 241605.0, 241635.0, 241665.0, 241695.0,
           241725.0, 241755.0, 241785.0,
           ...
           343005.0, 343035.0, 343065.0, 343095.0, 343125.0, 343155.0, 343185.0,
           343215.0, 343245.0, 343275.0],
          dtype='float64', name='x', length=3393))

• Attributes: (5)

  AREA_OR_POINT :

  Area

  DataType :

  Generic

  RepresentationType :

  ATHEMATIC

  _FillValue :

  0

  source_file :

  NL33.tif; NL32.tif