Auxiliary Data and Datacubes

DARTS uses several auxiliary data - data which does not change between different scenes and / or time steps. Raster auxiliary data is stored in Zarr Datacubes.

Currently, the following auxiliary data is used:

• ArcticDEM
• Tasseled Cap indices (Brightness, Greenness, Wetness)

with more to come.

ArcticDEM

The ArcticDEM is downloaded via their STAC server using these extend files.

The user can specify the download directory, where the ArcticDEM will be procedurally stored in a Zarr Datacube. The user can also specify the resolution of the ArcticDEM, which is either 2m, 10m or 32m. Each resolution is stored in their own Zarr Datacube.

darts_acquisition.load_arcticdem

load_arcticdem(
    geobox: odc.geo.geobox.GeoBox,
    data_dir: pathlib.Path | str,
    resolution: darts_acquisition.arcticdem.RESOLUTIONS,
    buffer: int = 0,
    persist: bool = True,
) -> xarray.Dataset

Load the ArcticDEM for the given geobox, fetch new data from the STAC server if necessary.

Parameters:

• geobox (odc.geo.geobox.GeoBox) –

  The geobox for which the tile should be loaded.

• data_dir (pathlib.Path | str) –

  The directory where the ArcticDEM data is stored.

• resolution (typing.Literal[2, 10, 32]) –

  The resolution of the ArcticDEM data in m.

• buffer (int, default: 0 ) –

  The buffer around the projected (epsg:3413) geobox in pixels. Defaults to 0.

• persist (bool, default: True ) –

  If the data should be persisted in memory. If not, this will return a Dask backed Dataset. Defaults to True.

Returns:

• xarray.Dataset –

  xr.Dataset: The ArcticDEM tile, with a buffer applied. Note: The buffer is applied in the arcticdem dataset's CRS, hence the orientation might be different. Final dataset is NOT matched to the reference CRS and resolution.

Warning

Geobox must be in a meter based CRS.

Usage

Since the API of the load_arcticdem is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

import xarray as xr
import odc.geo.xr

from darts_aquisition import load_arcticdem

# Assume "optical" is an already loaded s2 based dataarray

arcticdem = load_arcticdem(
    optical.odc.geobox,
    "/path/to/arcticdem-parent-directory",
    resolution=2,
    buffer=ceil(self.tpi_outer_radius / 2 * sqrt(2))
)

# Now we can for example match the resolution and extent of the optical data:
arcticdem = arcticdem.odc.reproject(optical.odc.geobox, resampling="cubic")

The buffer parameter is used to extend the region of interest by a certain amount of pixels. This comes handy when calculating e.g. the Topographic Position Index (TPI), which requires a buffer around the region of interest to remove edge effects.

Raises:

• ValueError –

  If the resolution is not supported.

Source code in darts-acquisition/src/darts_acquisition/arcticdem.py

@stopuhr.funkuhr("Loading ArcticDEM", printer=logger.debug, print_kwargs=True)
def load_arcticdem(
    geobox: GeoBox, data_dir: Path | str, resolution: RESOLUTIONS, buffer: int = 0, persist: bool = True
) -> xr.Dataset:
    """Load the ArcticDEM for the given geobox, fetch new data from the STAC server if necessary.

    Args:
        geobox (GeoBox): The geobox for which the tile should be loaded.
        data_dir (Path | str): The directory where the ArcticDEM data is stored.
        resolution (Literal[2, 10, 32]): The resolution of the ArcticDEM data in m.
        buffer (int, optional): The buffer around the projected (epsg:3413) geobox in pixels. Defaults to 0.
        persist (bool, optional): If the data should be persisted in memory.
            If not, this will return a Dask backed Dataset. Defaults to True.

    Returns:
        xr.Dataset: The ArcticDEM tile, with a buffer applied.
            Note: The buffer is applied in the arcticdem dataset's CRS, hence the orientation might be different.
            Final dataset is NOT matched to the reference CRS and resolution.

    Warning:
        Geobox must be in a meter based CRS.

    Usage:
        Since the API of the `load_arcticdem` is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

        ```python
        import xarray as xr
        import odc.geo.xr

        from darts_aquisition import load_arcticdem

        # Assume "optical" is an already loaded s2 based dataarray

        arcticdem = load_arcticdem(
            optical.odc.geobox,
            "/path/to/arcticdem-parent-directory",
            resolution=2,
            buffer=ceil(self.tpi_outer_radius / 2 * sqrt(2))
        )

        # Now we can for example match the resolution and extent of the optical data:
        arcticdem = arcticdem.odc.reproject(optical.odc.geobox, resampling="cubic")
        ```

        The `buffer` parameter is used to extend the region of interest by a certain amount of pixels.
        This comes handy when calculating e.g. the Topographic Position Index (TPI), which requires a buffer around the region of interest to remove edge effects.

    Raises:
        ValueError: If the resolution is not supported.

    """  # noqa: E501
    odc.stac.configure_rio(cloud_defaults=True, aws={"aws_unsigned": True})

    match resolution:
        case 2:
            accessor = smart_geocubes.ArcticDEM2m(data_dir)
        case 10:
            accessor = smart_geocubes.ArcticDEM10m(data_dir)
        case 32:
            accessor = smart_geocubes.ArcticDEM32m(data_dir)
        case _:
            raise ValueError(f"Resolution {resolution} not supported, only 2m, 10m and 32m are supported")

    accessor.assert_created()

    arcticdem = accessor.load(geobox, buffer=buffer, persist=persist)

    # Change dtype of the datamask to uint8 for later reproject_match
    arcticdem["datamask"] = arcticdem.datamask.astype("uint8")

    return arcticdem

Tasseled Cap indices (TCVIS)

The TCVIS data is downloaded from Google Earth-Engine (GEE) using the TCVIS collection from Ingmar Nitze: "users/ingmarnitze/TCTrend_SR_2000-2019_TCVIS".

darts_acquisition.load_tcvis

load_tcvis(
    geobox: odc.geo.geobox.GeoBox,
    data_dir: pathlib.Path | str,
    buffer: int = 0,
    persist: bool = True,
) -> xarray.Dataset

Load the TCVIS for the given geobox, fetch new data from GEE if necessary.

Parameters:

• geobox (odc.geo.geobox.GeoBox) –

  The geobox to load the data for.

• data_dir (pathlib.Path | str) –

  The directory to store the downloaded data for faster access for consecutive calls.

• buffer (int, default: 0 ) –

  The buffer around the geobox in pixels. Defaults to 0.

• persist (bool, default: True ) –

  If the data should be persisted in memory. If not, this will return a Dask backed Dataset. Defaults to True.

Returns:

• xarray.Dataset –

  xr.Dataset: The TCVIS dataset.

Usage

Since the API of the load_tcvis is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

import xarray as xr
import odc.geo.xr

from darts_aquisition import load_tcvis

# Assume "optical" is an already loaded s2 based dataarray

tcvis = load_tcvis(
    optical.odc.geobox,
    "/path/to/tcvis-parent-directory",
)

# Now we can for example match the resolution and extent of the optical data:
tcvis = tcvis.odc.reproject(optical.odc.geobox, resampling="cubic")

Source code in darts-acquisition/src/darts_acquisition/tcvis.py

@stopuhr.funkuhr("Loading TCVIS", printer=logger.debug, print_kwargs=True)
def load_tcvis(
    geobox: GeoBox,
    data_dir: Path | str,
    buffer: int = 0,
    persist: bool = True,
) -> xr.Dataset:
    """Load the TCVIS for the given geobox, fetch new data from GEE if necessary.

    Args:
        geobox (GeoBox): The geobox to load the data for.
        data_dir (Path | str): The directory to store the downloaded data for faster access for consecutive calls.
        buffer (int, optional): The buffer around the geobox in pixels. Defaults to 0.
        persist (bool, optional): If the data should be persisted in memory.
            If not, this will return a Dask backed Dataset. Defaults to True.

    Returns:
        xr.Dataset: The TCVIS dataset.

    Usage:
        Since the API of the `load_tcvis` is based on GeoBox, one can load a specific ROI based on an existing Xarray DataArray:

        ```python
        import xarray as xr
        import odc.geo.xr

        from darts_aquisition import load_tcvis

        # Assume "optical" is an already loaded s2 based dataarray

        tcvis = load_tcvis(
            optical.odc.geobox,
            "/path/to/tcvis-parent-directory",
        )

        # Now we can for example match the resolution and extent of the optical data:
        tcvis = tcvis.odc.reproject(optical.odc.geobox, resampling="cubic")
        ```

    """  # noqa: E501
    accessor = smart_geocubes.TCTrend(data_dir, create_icechunk_storage=False)

    # We want to assume that the datacube is already created to be save in a multi-process environment
    accessor.assert_created()

    tcvis = accessor.load(geobox, buffer=buffer, persist=persist)

    # Rename to follow our conventions
    tcvis = tcvis.rename_vars(
        {
            "TCB_slope": "tc_brightness",
            "TCG_slope": "tc_greenness",
            "TCW_slope": "tc_wetness",
        }
    )

    return tcvis

Why Zarr Datacubes?

Zarr is a file format for storing chunked, compressed, N-dimensional arrays. It is designed to store large arrays of data, and to facilitate fast and efficient IO. Zarr works well integrated with Dask and Xarray.

By storing the auxiliary data in Zarr Datacubes, it is much easier and faster to access the data of interest. If we would use GeoTiffs, we would have to first create a Cloud-Optimized GeoTiff (COG), which is basically an ensemble (mosaic) of multiple GeoTiffs. Then we would have to read from the COG, which behind the scenes would open multiple GeoTiffs and crops them to fit the region of interest. E.g. Opening a specific region of interest 10km x 10km from a 2m resolution COG would take up to 2 minutes, if the COGs extend is panarctic. Opening the same region from a Zarr Datacube takes less than 1 second.

Inspiration

This implementation and concept is heavily inspired by EarthMovers implementation of serverless datacube generation.

Procedural download

Info

The currently used auxiliary data is downloaded on demand, only data actually used is downloaded and stored on your local machine. Hence, the stored datacubes can be thought of as a cache, which is filled with data as needed.

There are currently two implementations of the procedural download used: a cloud based STAC download and a download via Google Earth-Engine.

Because the single tiles of the STAC mosaic can be overlapping and intersect with multiple Zarr chunks, the STAC download is slightly more complicated. Since Google Earth-Engine allows for exact geoboxes, download of the exact chunks is possible. This reduces the complexity of the download.

	STAC	GEE
1. ROI
2. ROI

The above graphics shows the difference between loading data from STAC (left) and Google Earth-Engine (right). With the STAC download, the data is downloaded from a mosaic of tiles, which can be overlapping with each other and cover multiple Zarr chunks. It may occur that a chunk is not fully covered by the STAC mosaic, which results in only partial loaded chunks. In such cases, the missing data in these chunks will be updated if the other intersecting tile is downloaded, which may occur to a later time if a connected ROI is requested. The download process is much easier for GEE, since one can request the exact geoboxes of the Zarr chunks and GEE will handle the rest. Hence, chunks will always be fully covered by the downloaded data.

Regarding the open ROI process, both implementations follow the same principle:

1. Check which Tiles / Chunks intersect with the region of interest
2. Dowload all new Tiles / Chunks
3. Store the new Tiles / Chunks in their specific Zarr chunks
4. Return the region of interest of the Zarr Datacube

STAC download

Google Earth-Engine download