darts_preprocessing.v2

PLANET scene based preprocessing.

logger module-attribute

logger = logging.getLogger(
    __name__.replace("darts_", "darts.")
)

calculate_aspect

calculate_aspect(
    arcticdem_ds: xarray.Dataset,
) -> xarray.Dataset

Calculate aspect (compass direction) of the terrain surface from an ArcticDEM Dataset.

Aspect indicates the downslope direction of the maximum rate of change in elevation.

Parameters:

• arcticdem_ds (xarray.Dataset) –

  Dataset containing: - dem (float32): Digital Elevation Model

Returns:

• xarray.Dataset –

  xr.Dataset: Input Dataset with new data variable added:

• xarray.Dataset –

  • aspect (float32): Aspect in degrees clockwise from north [0-360], or -1 for flat areas.

  • long_name: "Aspect"

  • units: "degrees"
  • description: Compass direction of slope
  • source: "ArcticDEM"

Note

Aspect values:

• 0° or 360°: North-facing
• 90°: East-facing
• 180°: South-facing
• 270°: West-facing
• -1: Flat (no dominant direction)

Example

from darts_preprocessing import calculate_aspect

arcticdem_with_aspect = calculate_aspect(arcticdem_ds)

# Identify south-facing slopes (135-225 degrees)
south_facing = (arcticdem_with_aspect.aspect > 135) & (arcticdem_with_aspect.aspect < 225)

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py

@stopwatch("Calculating aspect", printer=logger.debug)
def calculate_aspect(arcticdem_ds: xr.Dataset) -> xr.Dataset:
    """Calculate aspect (compass direction) of the terrain surface from an ArcticDEM Dataset.

    Aspect indicates the downslope direction of the maximum rate of change in elevation.

    Args:
        arcticdem_ds (xr.Dataset): Dataset containing:
            - dem (float32): Digital Elevation Model

    Returns:
        xr.Dataset: Input Dataset with new data variable added:

        - aspect (float32): Aspect in degrees clockwise from north [0-360], or -1 for flat areas.

            - long_name: "Aspect"
            - units: "degrees"
            - description: Compass direction of slope
            - source: "ArcticDEM"

    Note:
        Aspect values:

        - 0° or 360°: North-facing
        - 90°: East-facing
        - 180°: South-facing
        - 270°: West-facing
        - -1: Flat (no dominant direction)

    Example:
        ```python
        from darts_preprocessing import calculate_aspect

        arcticdem_with_aspect = calculate_aspect(arcticdem_ds)

        # Identify south-facing slopes (135-225 degrees)
        south_facing = (arcticdem_with_aspect.aspect > 135) & (arcticdem_with_aspect.aspect < 225)
        ```

    """
    aspect_deg = aspect(arcticdem_ds.dem)

    # Aspect is always calculated in the projection - thus "north" is rather an "up"
    # To get the true north, we need to correct the aspect based on the coordinates
    x = arcticdem_ds.x.expand_dims({"y": arcticdem_ds.y})
    y = arcticdem_ds.y.expand_dims({"x": arcticdem_ds.x})
    if arcticdem_ds.cupy.is_cupy:
        x = cp.asarray(x)
        y = cp.asarray(y)
        correction_offset = cp.arctan2(x, y) * (180 / np.pi) + 90
    else:
        correction_offset = np.arctan2(x, y) * (180 / np.pi) + 90
    aspect_deg = (aspect_deg + correction_offset) % 360

    aspect_deg.attrs = {
        "long_name": "Aspect",
        "units": "degrees",
        "description": "The compass direction that the slope faces, in degrees clockwise from north.",
        "source": "ArcticDEM",
    }
    arcticdem_ds["aspect"] = aspect_deg.compute()
    return arcticdem_ds

calculate_curvature

calculate_curvature(
    arcticdem_ds: xarray.Dataset,
) -> xarray.Dataset

Calculate curvature of the terrain surface from an ArcticDEM Dataset.

Curvature measures the rate of change of slope, indicating terrain convexity or concavity.

Parameters:

• arcticdem_ds (xarray.Dataset) –

  Dataset containing: - dem (float32): Digital Elevation Model

Returns:

• xarray.Dataset –

  xr.Dataset: Input Dataset with new data variable added:

• xarray.Dataset –

  • curvature (float32): Curvature values.

  • long_name: "Curvature"

  • description: Rate of change of slope
  • source: "ArcticDEM"

Note

Curvature interpretation:

• Positive values: Convex terrain (hills, ridges)
• Negative values: Concave terrain (valleys, depressions)
• Near zero: Planar terrain

Example

from darts_preprocessing import calculate_curvature

arcticdem_with_curv = calculate_curvature(arcticdem_ds)

# Identify ridges (convex areas)
ridges = arcticdem_with_curv.curvature > 0.1

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py

@stopwatch("Calculating curvature", printer=logger.debug)
def calculate_curvature(arcticdem_ds: xr.Dataset) -> xr.Dataset:
    """Calculate curvature of the terrain surface from an ArcticDEM Dataset.

    Curvature measures the rate of change of slope, indicating terrain convexity or concavity.

    Args:
        arcticdem_ds (xr.Dataset): Dataset containing:
            - dem (float32): Digital Elevation Model

    Returns:
        xr.Dataset: Input Dataset with new data variable added:

        - curvature (float32): Curvature values.

            - long_name: "Curvature"
            - description: Rate of change of slope
            - source: "ArcticDEM"

    Note:
        Curvature interpretation:

        - Positive values: Convex terrain (hills, ridges)
        - Negative values: Concave terrain (valleys, depressions)
        - Near zero: Planar terrain

    Example:
        ```python
        from darts_preprocessing import calculate_curvature

        arcticdem_with_curv = calculate_curvature(arcticdem_ds)

        # Identify ridges (convex areas)
        ridges = arcticdem_with_curv.curvature > 0.1
        ```

    """
    curvature_da = curvature(arcticdem_ds.dem)
    curvature_da.attrs = {
        "long_name": "Curvature",
        "units": "",
        "description": "The curvature of the terrain surface.",
        "source": "ArcticDEM",
    }
    arcticdem_ds["curvature"] = curvature_da.compute()
    return arcticdem_ds

calculate_hillshade

calculate_hillshade(
    arcticdem_ds: xarray.Dataset,
    azimuth: int = 225,
    angle_altitude: int = 25,
) -> xarray.Dataset

Calculate hillshade of the terrain surface from an ArcticDEM Dataset.

Hillshade simulates illumination of terrain from a specified sun position, useful for visualization and terrain analysis.

Parameters:

• arcticdem_ds (xarray.Dataset) –

  Dataset containing: - dem (float32): Digital Elevation Model

• azimuth (int, default: 225 ) –

  Light source azimuth in degrees clockwise from north [0-360]. Defaults to 225 (southwest).

• angle_altitude (int, default: 25 ) –

  Light source altitude angle in degrees above horizon [0-90]. Defaults to 25.

Returns:

• xarray.Dataset –

  xr.Dataset: Input Dataset with new data variable added:

• xarray.Dataset –

  • hillshade (float32): Illumination values [0-255], where 0 is shadow and 255 is fully lit.

  • long_name: "Hillshade"

  • description: Documents azimuth and angle_altitude used
  • source: "ArcticDEM"

Note

Common azimuth/altitude combinations:

• 315°/45°: Classic northwest illumination (default for many GIS applications)
• 225°/25°: Southwest with low sun (better for visualizing subtle features)

The hillshade calculation accounts for both slope and aspect of the terrain.

Example

from darts_preprocessing import calculate_hillshade

# Default southwest illumination
arcticdem_with_hs = calculate_hillshade(arcticdem_ds)

# Custom sun position
arcticdem_custom = calculate_hillshade(
    arcticdem_ds,
    azimuth=315,
    angle_altitude=45
)

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py

@stopwatch.f("Calculating hillshade", printer=logger.debug, print_kwargs=["azimuth", "angle_altitude"])
def calculate_hillshade(arcticdem_ds: xr.Dataset, azimuth: int = 225, angle_altitude: int = 25) -> xr.Dataset:
    """Calculate hillshade of the terrain surface from an ArcticDEM Dataset.

    Hillshade simulates illumination of terrain from a specified sun position, useful
    for visualization and terrain analysis.

    Args:
        arcticdem_ds (xr.Dataset): Dataset containing:
            - dem (float32): Digital Elevation Model
        azimuth (int, optional): Light source azimuth in degrees clockwise from north [0-360].
            Defaults to 225 (southwest).
        angle_altitude (int, optional): Light source altitude angle in degrees above horizon [0-90].
            Defaults to 25.

    Returns:
        xr.Dataset: Input Dataset with new data variable added:

        - hillshade (float32): Illumination values [0-255], where 0 is shadow and 255 is fully lit.

            - long_name: "Hillshade"
            - description: Documents azimuth and angle_altitude used
            - source: "ArcticDEM"

    Note:
        Common azimuth/altitude combinations:

        - 315°/45°: Classic northwest illumination (default for many GIS applications)
        - 225°/25°: Southwest with low sun (better for visualizing subtle features)

        The hillshade calculation accounts for both slope and aspect of the terrain.

    Example:
        ```python
        from darts_preprocessing import calculate_hillshade

        # Default southwest illumination
        arcticdem_with_hs = calculate_hillshade(arcticdem_ds)

        # Custom sun position
        arcticdem_custom = calculate_hillshade(
            arcticdem_ds,
            azimuth=315,
            angle_altitude=45
        )
        ```

    """
    x, y = arcticdem_ds.x.mean().item(), arcticdem_ds.y.mean().item()
    correction_offset = np.arctan2(x, y) * (180 / np.pi) + 90
    azimuth_corrected = (azimuth - correction_offset + 360) % 360

    hillshade_da = hillshade(arcticdem_ds.dem, azimuth=azimuth_corrected, angle_altitude=angle_altitude)
    hillshade_da.attrs = {
        "long_name": "Hillshade",
        "units": "",
        "description": f"The hillshade based on azimuth {azimuth} and angle_altitude {angle_altitude}.",
        "source": "ArcticDEM",
    }
    arcticdem_ds["hillshade"] = hillshade_da.compute()
    return arcticdem_ds

calculate_ndvi

calculate_ndvi(optical: xarray.Dataset) -> xarray.Dataset

Calculate NDVI (Normalized Difference Vegetation Index) from spectral bands.

NDVI is a widely-used vegetation index that indicates photosynthetic activity and vegetation health. Values range from -1 to 1, with higher values indicating denser, healthier vegetation.

Parameters:

• optical (xarray.Dataset) –

  Dataset containing spectral bands: - nir (float32): Near-infrared reflectance [0-1] - red (float32): Red reflectance [0-1]

Returns:

• xarray.Dataset –

  xr.DataArray: NDVI values with attributes: - long_name: "NDVI"

Note

Formula: NDVI = (NIR - Red) / (NIR + Red)

Input bands are clipped to [0, 1] before calculation to avoid numerical instabilities from negative reflectance values or sensor artifacts. The final result is also clipped to ensure values remain in the valid [-1, 1] range.

Example

Calculate NDVI from optical data:

from darts_preprocessing import calculate_ndvi

# optical contains 'nir' and 'red' bands
ndvi = calculate_ndvi(optical)

# Mask vegetation
vegetation_mask = ndvi > 0.3

Source code in darts-preprocessing/src/darts_preprocessing/engineering/indices.py

@stopwatch("Calculating NDVI", printer=logger.debug)
def calculate_ndvi(optical: xr.Dataset) -> xr.Dataset:
    """Calculate NDVI (Normalized Difference Vegetation Index) from spectral bands.

    NDVI is a widely-used vegetation index that indicates photosynthetic activity and
    vegetation health. Values range from -1 to 1, with higher values indicating denser,
    healthier vegetation.

    Args:
        optical (xr.Dataset): Dataset containing spectral bands:
            - nir (float32): Near-infrared reflectance [0-1]
            - red (float32): Red reflectance [0-1]

    Returns:
        xr.DataArray: NDVI values with attributes:
            - long_name: "NDVI"

    Note:
        Formula: NDVI = (NIR - Red) / (NIR + Red)

        Input bands are clipped to [0, 1] before calculation to avoid numerical instabilities
        from negative reflectance values or sensor artifacts. The final result is also clipped
        to ensure values remain in the valid [-1, 1] range.

    Example:
        Calculate NDVI from optical data:

        ```python
        from darts_preprocessing import calculate_ndvi

        # optical contains 'nir' and 'red' bands
        ndvi = calculate_ndvi(optical)

        # Mask vegetation
        vegetation_mask = ndvi > 0.3
        ```

    """
    nir = optical["nir"].clip(0, 1)
    r = optical["red"].clip(0, 1)
    ndvi = (nir - r) / (nir + r)
    ndvi = ndvi.clip(-1, 1)
    ndvi = ndvi.assign_attrs({"long_name": "NDVI"}).rename("ndvi")
    return ndvi

calculate_slope

calculate_slope(
    arcticdem_ds: xarray.Dataset,
) -> xarray.Dataset

Calculate slope of the terrain surface from an ArcticDEM Dataset.

Slope represents the rate of change of elevation, indicating terrain steepness.

Parameters:

• arcticdem_ds (xarray.Dataset) –

  Dataset containing: - dem (float32): Digital Elevation Model

Returns:

• xarray.Dataset –

  xr.Dataset: Input Dataset with new data variable added:

• xarray.Dataset –

  • slope (float32): Slope in degrees [0-90].

  • long_name: "Slope"

  • units: "degrees"
  • source: "ArcticDEM"

Note

Slope is calculated using finite difference methods on the DEM. Values approaching 90° indicate near-vertical terrain.

Example

from darts_preprocessing import calculate_slope

arcticdem_with_slope = calculate_slope(arcticdem_ds)

# Mask steep terrain
steep_areas = arcticdem_with_slope.slope > 30

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py

@stopwatch("Calculating slope", printer=logger.debug)
def calculate_slope(arcticdem_ds: xr.Dataset) -> xr.Dataset:
    """Calculate slope of the terrain surface from an ArcticDEM Dataset.

    Slope represents the rate of change of elevation, indicating terrain steepness.

    Args:
        arcticdem_ds (xr.Dataset): Dataset containing:
            - dem (float32): Digital Elevation Model

    Returns:
        xr.Dataset: Input Dataset with new data variable added:

        - slope (float32): Slope in degrees [0-90].

            - long_name: "Slope"
            - units: "degrees"
            - source: "ArcticDEM"

    Note:
        Slope is calculated using finite difference methods on the DEM.
        Values approaching 90° indicate near-vertical terrain.

    Example:
        ```python
        from darts_preprocessing import calculate_slope

        arcticdem_with_slope = calculate_slope(arcticdem_ds)

        # Mask steep terrain
        steep_areas = arcticdem_with_slope.slope > 30
        ```

    """
    slope_deg = slope(arcticdem_ds.dem)
    slope_deg.attrs = {
        "long_name": "Slope",
        "units": "degrees",
        "description": "The slope of the terrain surface in degrees.",
        "source": "ArcticDEM",
    }
    arcticdem_ds["slope"] = slope_deg.compute()
    return arcticdem_ds

calculate_topographic_position_index

calculate_topographic_position_index(
    arcticdem_ds: xarray.Dataset,
    outer_radius: int,
    inner_radius: int,
) -> xarray.Dataset

Calculate the Topographic Position Index (TPI) from an ArcticDEM Dataset.

TPI measures the relative topographic position of a point by comparing its elevation to the mean elevation of the surrounding neighborhood. Positive values indicate higher positions (ridges), negative values indicate lower positions (valleys).

Parameters:

• arcticdem_ds (xarray.Dataset) –

  The ArcticDEM Dataset containing the 'dem' variable (float32).

• outer_radius (int) –

  The outer radius of the neighborhood in meters. Can also be specified as string with units (e.g., "100m" or "10px").

• inner_radius (int) –

  The inner radius of the annulus kernel in meters. If > 0, creates an annulus (ring) instead of a circle. Set to 0 for a circular kernel. Can also be specified as string with units (e.g., "50m" or "5px").

Returns:

• xarray.Dataset –

  xr.Dataset: The input Dataset with a new data variable added:

• xarray.Dataset –

  • tpi (float32): Topographic Position Index values.

  • long_name: "Topographic Position Index (TPI)"

  • description: Details about the kernel used

Note

Kernel shape combinations:

• inner_radius=0: Circular kernel comparing each cell to all neighbors within outer_radius
• inner_radius>0: Annulus kernel comparing each cell to neighbors in a ring between inner_radius and outer_radius. Useful for multi-scale terrain analysis.

The actual radii used are rounded to the nearest pixel based on the DEM resolution.

Example

Calculate TPI with circular and annulus kernels:

from darts_preprocessing import calculate_topographic_position_index

# Circular kernel (100m radius)
arcticdem_with_tpi = calculate_topographic_position_index(
    arcticdem_ds=arcticdem,
    outer_radius=100,
    inner_radius=0
)

# Annulus kernel (50-100m ring)
arcticdem_multi_scale = calculate_topographic_position_index(
    arcticdem_ds=arcticdem,
    outer_radius=100,
    inner_radius=50
)

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py

@stopwatch.f("Calculating TPI", printer=logger.debug, print_kwargs=["outer_radius", "inner_radius"])
def calculate_topographic_position_index(arcticdem_ds: xr.Dataset, outer_radius: int, inner_radius: int) -> xr.Dataset:
    """Calculate the Topographic Position Index (TPI) from an ArcticDEM Dataset.

    TPI measures the relative topographic position of a point by comparing its elevation to
    the mean elevation of the surrounding neighborhood. Positive values indicate higher
    positions (ridges), negative values indicate lower positions (valleys).

    Args:
        arcticdem_ds (xr.Dataset): The ArcticDEM Dataset containing the 'dem' variable (float32).
        outer_radius (int): The outer radius of the neighborhood in meters.
            Can also be specified as string with units (e.g., "100m" or "10px").
        inner_radius (int): The inner radius of the annulus kernel in meters.
            If > 0, creates an annulus (ring) instead of a circle. Set to 0 for a circular kernel.
            Can also be specified as string with units (e.g., "50m" or "5px").

    Returns:
        xr.Dataset: The input Dataset with a new data variable added:

        - tpi (float32): Topographic Position Index values.

            - long_name: "Topographic Position Index (TPI)"
            - description: Details about the kernel used

    Note:
        Kernel shape combinations:

        - inner_radius=0: Circular kernel comparing each cell to all neighbors within outer_radius
        - inner_radius>0: Annulus kernel comparing each cell to neighbors in a ring between
          inner_radius and outer_radius. Useful for multi-scale terrain analysis.

        The actual radii used are rounded to the nearest pixel based on the DEM resolution.

    Example:
        Calculate TPI with circular and annulus kernels:

        ```python
        from darts_preprocessing import calculate_topographic_position_index

        # Circular kernel (100m radius)
        arcticdem_with_tpi = calculate_topographic_position_index(
            arcticdem_ds=arcticdem,
            outer_radius=100,
            inner_radius=0
        )

        # Annulus kernel (50-100m ring)
        arcticdem_multi_scale = calculate_topographic_position_index(
            arcticdem_ds=arcticdem,
            outer_radius=100,
            inner_radius=50
        )
        ```

    """
    cellsize_x, cellsize_y = convolution.calc_cellsize(arcticdem_ds.dem)  # Should be equal to the resolution of the DEM
    # Use an annulus kernel if inner_radius is greater than 0
    outer_radius: Distance = Distance.parse(outer_radius, cellsize_x)
    if inner_radius > 0:
        inner_radius: Distance = Distance.parse(inner_radius, cellsize_x)
        kernel = convolution.annulus_kernel(cellsize_x, cellsize_y, outer_radius.meter, inner_radius.meter)
        attr_cell_description = (
            f"within a ring at a distance of {inner_radius}-{outer_radius} cells away from the focal cell."
        )
        logger.debug(
            f"Calculating Topographic Position Index with annulus kernel of {inner_radius}-{outer_radius} cells."
        )
    else:
        kernel = convolution.circle_kernel(cellsize_x, cellsize_y, outer_radius.meter)
        attr_cell_description = f"within a circle at a distance of {outer_radius} cells away from the focal cell."
        logger.debug(f"Calculating Topographic Position Index with circle kernel of {outer_radius} cells.")

    # Change dtype of kernel to float32 since we don't need the precision and f32 is faster
    kernel = kernel.astype("float32")

    if has_cuda_and_cupy() and arcticdem_ds.cupy.is_cupy:
        kernel = cp.asarray(kernel)

    tpi = arcticdem_ds.dem - convolution.convolution_2d(arcticdem_ds.dem, kernel) / kernel.sum()
    tpi.attrs = {
        "long_name": "Topographic Position Index",
        "units": "m",
        "description": "The difference between the elevation of a cell and the mean elevation of the surrounding"
        f"cells {attr_cell_description}",
        "source": "ArcticDEM",
    }

    arcticdem_ds["tpi"] = tpi.compute()

    return arcticdem_ds

get_azimuth_and_elevation

get_azimuth_and_elevation(
    ds_optical: xarray.Dataset,
) -> tuple[float, float]

Get the azimuth and elevation from the optical dataset attributes.

Parameters:

• ds_optical (xarray.Dataset) –

  The optical dataset.

Returns:

• tuple[float, float] –

  tuple[float, float]: The azimuth and elevation.

Source code in darts-preprocessing/src/darts_preprocessing/v2.py

def get_azimuth_and_elevation(ds_optical: xr.Dataset) -> tuple[float, float]:
    """Get the azimuth and elevation from the optical dataset attributes.

    Args:
        ds_optical (xr.Dataset): The optical dataset.

    Returns:
        tuple[float, float]: The azimuth and elevation.

    """
    azimuth = ds_optical.attrs.get("azimuth", float("nan"))
    elevation = ds_optical.attrs.get("elevation", float("nan"))
    if isnan(azimuth):
        azimuth = 225
        logger.warning("No azimuth found in optical dataset attributes. Using default value of 225 degrees.")
    if isnan(elevation):
        elevation = 25
        logger.warning("No sun elevation found in optical dataset attributes. Using default value of 25 degrees.")
    if not isinstance(azimuth, (int, float)):
        azimuth = 225
        logger.warning(
            f"Azimuth found in optical dataset is {azimuth}, which is not a number. Using default value of 225 degrees."
        )
    if not isinstance(elevation, (int, float)):
        elevation = 25
        logger.warning(
            f"Sun elevation found in optical dataset is {elevation}, which is not a number."
            " Using default value of 25 degrees."
        )

    azimuth = round(azimuth)
    elevation = round(elevation)
    return azimuth, elevation

preprocess_arcticdem

preprocess_arcticdem(
    ds_arcticdem: xarray.Dataset,
    tpi_outer_radius: int,
    tpi_inner_radius: int,
    azimuth: int,
    angle_altitude: int,
) -> xarray.Dataset

Preprocess the ArcticDEM data with mdoern (DARTS v2) preprocessing steps.

Parameters:

• ds_arcticdem (xarray.Dataset) –

  The ArcticDEM dataset.

• tpi_outer_radius (int) –

  The outer radius of the annulus kernel for the tpi calculation in number of cells.

• tpi_inner_radius (int) –

  The inner radius of the annulus kernel for the tpi calculation in number of cells.

• azimuth (int) –

  The azimuth angle of the light source in degrees for hillshade calculation.

• angle_altitude (int) –

  The altitude angle of the light source in degrees for hillshade

Returns:

• xarray.Dataset –

  xr.Dataset: The preprocessed ArcticDEM dataset.

Source code in darts-preprocessing/src/darts_preprocessing/v2.py

@stopwatch.f(
    "Preprocessing arcticdem",
    printer=logger.debug,
    print_kwargs=["tpi_outer_radius", "tpi_inner_radius", "azimuth", "angle_altitude"],
)
def preprocess_arcticdem(
    ds_arcticdem: xr.Dataset,
    tpi_outer_radius: int,
    tpi_inner_radius: int,
    azimuth: int,
    angle_altitude: int,
) -> xr.Dataset:
    """Preprocess the ArcticDEM data with mdoern (DARTS v2) preprocessing steps.

    Args:
        ds_arcticdem (xr.Dataset): The ArcticDEM dataset.
        tpi_outer_radius (int): The outer radius of the annulus kernel for the tpi calculation in number of cells.
        tpi_inner_radius (int): The inner radius of the annulus kernel for the tpi calculation in number of cells.
        azimuth (int): The azimuth angle of the light source in degrees for hillshade calculation.
        angle_altitude (int): The altitude angle of the light source in degrees for hillshade

    Returns:
        xr.Dataset: The preprocessed ArcticDEM dataset.

    """
    ds_arcticdem = calculate_topographic_position_index(ds_arcticdem, tpi_outer_radius, tpi_inner_radius)
    ds_arcticdem = calculate_slope(ds_arcticdem)
    ds_arcticdem = calculate_hillshade(ds_arcticdem, azimuth=azimuth, angle_altitude=angle_altitude)
    ds_arcticdem = calculate_aspect(ds_arcticdem)
    ds_arcticdem = calculate_curvature(ds_arcticdem)

    return ds_arcticdem

preprocess_v2

preprocess_v2(
    ds_optical: xarray.Dataset,
    ds_arcticdem: xarray.Dataset | None,
    ds_tcvis: xarray.Dataset | None,
    tpi_outer_radius: int = 100,
    tpi_inner_radius: int = 0,
    device: typing.Literal["cuda", "cpu"]
    | int = darts_utils.cuda.DEFAULT_DEVICE,
) -> xarray.Dataset

Preprocess optical data with modern (DARTS v2) preprocessing steps.

This function combines optical imagery with terrain (ArcticDEM) and temporal vegetation indices (TCVIS) to create a multi-source feature dataset for segmentation. All auxiliary data sources are reprojected and cropped to match the optical data's extent and resolution.

Processing steps

1. Calculate NDVI from optical bands
2. If TCVIS provided: Reproject and merge Tasseled Cap trends
3. If ArcticDEM provided: Calculate terrain features (TPI, slope, hillshade, aspect, curvature) using solar geometry from optical data attributes

Parameters:

• ds_optical (xarray.Dataset) –

  Optical imagery dataset (PlanetScope or Sentinel-2) containing: - Required variables: blue, green, red, nir (float32, reflectance values) - Required variables: quality_data_mask, valid_data_mask (uint8) - Required attributes: azimuth (float), elevation (float) for hillshade calculation

• ds_arcticdem (xarray.Dataset | None) –

  ArcticDEM dataset containing 'dem' (float32) and 'arcticdem_data_mask' (uint8). If None, terrain features are skipped.

• ds_tcvis (xarray.Dataset | None) –

  TCVIS dataset containing tc_brightness, tc_greenness, tc_wetness (float). If None, TCVIS features are skipped.

• tpi_outer_radius (int, default: 100 ) –

  Outer radius for TPI calculation in meters. Defaults to 100m.

• tpi_inner_radius (int, default: 0 ) –

  Inner radius for TPI annulus kernel in meters. Set to 0 for circular kernel. Defaults to 0.

• device (typing.Literal['cuda', 'cpu'] | int, default: darts_utils.cuda.DEFAULT_DEVICE ) –

  Device for GPU-accelerated computations (NDVI, TPI, slope). Use "cuda" for first GPU, int for specific GPU, or "cpu". Defaults to "cuda" if available, else "cpu".

Returns:

• xarray.Dataset –

  xr.Dataset: Preprocessed dataset with all input optical variables plus:

• xarray.Dataset –

  Added from optical processing: - ndvi (float32): Normalized Difference Vegetation Index Attributes: long_name="NDVI"

• xarray.Dataset –

  Added from TCVIS (if ds_tcvis provided): - tc_brightness (float): Tasseled Cap brightness trend - tc_greenness (float): Tasseled Cap greenness trend - tc_wetness (float): Tasseled Cap wetness trend

• xarray.Dataset –

  Added from ArcticDEM (if ds_arcticdem provided): - dem (float32): Elevation in meters - relative_elevation (float32): Topographic Position Index (TPI) Attributes: long_name="Topographic Position Index (TPI)" - slope (float32): Slope in degrees [0-90] Attributes: long_name="Slope" - hillshade (uint8): Hillshade values [0-255] Attributes: long_name="Hillshade" - aspect (float32): Aspect in degrees [0-360] Attributes: long_name="Aspect" - curvature (float32): Surface curvature Attributes: long_name="Curvature" - arcticdem_data_mask (uint8): DEM validity mask

Note

Attribute usage: - azimuth attribute from ds_optical: Used for hillshade calculation (solar azimuth angle). Falls back to 225° if missing or invalid. - elevation attribute from ds_optical: Used for hillshade calculation (solar elevation angle). Falls back to 25° if missing or invalid.

Processing behavior: - If both ds_tcvis and ds_arcticdem are None, only NDVI is calculated. - ArcticDEM is buffered by tpi_outer_radius before reprojection to avoid edge effects, then cropped back to optical extent after terrain feature calculation. - Reprojection uses cubic resampling for smooth terrain features. - GPU acceleration (if device="cuda") significantly speeds up TPI and slope calculations.

Example

Complete preprocessing with all data sources:

from darts_preprocessing import preprocess_v2
from darts_acquisition import load_cdse_s2_sr_scene, load_arcticdem, load_tcvis

# Load optical data
optical = load_cdse_s2_sr_scene(s2_scene_id, ...)

# Load auxiliary data
arcticdem = load_arcticdem(optical.odc.geobox, ...)
tcvis = load_tcvis(optical.odc.geobox, ...)

# Preprocess
preprocessed = preprocess_v2(
    ds_optical=optical,
    ds_arcticdem=arcticdem,
    ds_tcvis=tcvis,
    tpi_outer_radius=100,
    tpi_inner_radius=0,
    device="cuda"
)

# Result contains: blue, green, red, nir, ndvi, tc_brightness, tc_greenness,
# tc_wetness, dem, relative_elevation, slope, hillshade, aspect, curvature

Source code in darts-preprocessing/src/darts_preprocessing/v2.py

@stopwatch("Preprocessing", printer=logger.debug)
def preprocess_v2(
    ds_optical: xr.Dataset,
    ds_arcticdem: xr.Dataset | None,
    ds_tcvis: xr.Dataset | None,
    tpi_outer_radius: int = 100,
    tpi_inner_radius: int = 0,
    device: Literal["cuda", "cpu"] | int = DEFAULT_DEVICE,
) -> xr.Dataset:
    """Preprocess optical data with modern (DARTS v2) preprocessing steps.

    This function combines optical imagery with terrain (ArcticDEM) and temporal vegetation
    indices (TCVIS) to create a multi-source feature dataset for segmentation. All auxiliary
    data sources are reprojected and cropped to match the optical data's extent and resolution.

    Processing steps:
        1. Calculate NDVI from optical bands
        2. If TCVIS provided: Reproject and merge Tasseled Cap trends
        3. If ArcticDEM provided: Calculate terrain features (TPI, slope, hillshade, aspect, curvature)
           using solar geometry from optical data attributes

    Args:
        ds_optical (xr.Dataset): Optical imagery dataset (PlanetScope or Sentinel-2) containing:
            - Required variables: blue, green, red, nir (float32, reflectance values)
            - Required variables: quality_data_mask, valid_data_mask (uint8)
            - Required attributes: azimuth (float), elevation (float) for hillshade calculation
        ds_arcticdem (xr.Dataset | None): ArcticDEM dataset containing 'dem' (float32) and
            'arcticdem_data_mask' (uint8). If None, terrain features are skipped.
        ds_tcvis (xr.Dataset | None): TCVIS dataset containing tc_brightness, tc_greenness,
            tc_wetness (float). If None, TCVIS features are skipped.
        tpi_outer_radius (int, optional): Outer radius for TPI calculation in meters.
            Defaults to 100m.
        tpi_inner_radius (int, optional): Inner radius for TPI annulus kernel in meters.
            Set to 0 for circular kernel. Defaults to 0.
        device (Literal["cuda", "cpu"] | int, optional): Device for GPU-accelerated computations
            (NDVI, TPI, slope). Use "cuda" for first GPU, int for specific GPU, or "cpu".
            Defaults to "cuda" if available, else "cpu".

    Returns:
        xr.Dataset: Preprocessed dataset with all input optical variables plus:

        Added from optical processing:
            - ndvi (float32): Normalized Difference Vegetation Index
              Attributes: long_name="NDVI"

        Added from TCVIS (if ds_tcvis provided):
            - tc_brightness (float): Tasseled Cap brightness trend
            - tc_greenness (float): Tasseled Cap greenness trend
            - tc_wetness (float): Tasseled Cap wetness trend

        Added from ArcticDEM (if ds_arcticdem provided):
            - dem (float32): Elevation in meters
            - relative_elevation (float32): Topographic Position Index (TPI)
              Attributes: long_name="Topographic Position Index (TPI)"
            - slope (float32): Slope in degrees [0-90]
              Attributes: long_name="Slope"
            - hillshade (uint8): Hillshade values [0-255]
              Attributes: long_name="Hillshade"
            - aspect (float32): Aspect in degrees [0-360]
              Attributes: long_name="Aspect"
            - curvature (float32): Surface curvature
              Attributes: long_name="Curvature"
            - arcticdem_data_mask (uint8): DEM validity mask

    Note:
        Attribute usage:
        - `azimuth` attribute from ds_optical: Used for hillshade calculation (solar azimuth angle).
          Falls back to 225° if missing or invalid.
        - `elevation` attribute from ds_optical: Used for hillshade calculation (solar elevation angle).
          Falls back to 25° if missing or invalid.

        Processing behavior:
        - If both ds_tcvis and ds_arcticdem are None, only NDVI is calculated.
        - ArcticDEM is buffered by tpi_outer_radius before reprojection to avoid edge effects,
          then cropped back to optical extent after terrain feature calculation.
        - Reprojection uses cubic resampling for smooth terrain features.
        - GPU acceleration (if device="cuda") significantly speeds up TPI and slope calculations.

    Example:
        Complete preprocessing with all data sources:

        ```python
        from darts_preprocessing import preprocess_v2
        from darts_acquisition import load_cdse_s2_sr_scene, load_arcticdem, load_tcvis

        # Load optical data
        optical = load_cdse_s2_sr_scene(s2_scene_id, ...)

        # Load auxiliary data
        arcticdem = load_arcticdem(optical.odc.geobox, ...)
        tcvis = load_tcvis(optical.odc.geobox, ...)

        # Preprocess
        preprocessed = preprocess_v2(
            ds_optical=optical,
            ds_arcticdem=arcticdem,
            ds_tcvis=tcvis,
            tpi_outer_radius=100,
            tpi_inner_radius=0,
            device="cuda"
        )

        # Result contains: blue, green, red, nir, ndvi, tc_brightness, tc_greenness,
        # tc_wetness, dem, relative_elevation, slope, hillshade, aspect, curvature
        ```

    """
    # Move to GPU for faster calculations
    ds_optical = move_to_device(ds_optical, device)
    # Calculate NDVI
    ds_optical["ndvi"] = calculate_ndvi(ds_optical)
    ds_optical = move_to_host(ds_optical)

    if ds_tcvis is None and ds_arcticdem is None:
        logger.debug("No auxiliary data provided. Only NDVI was calculated.")
        return ds_optical

    if ds_tcvis is not None:
        # Reproject TCVIS to optical data
        with stopwatch("Reprojecting TCVIS", printer=logger.debug):
            # *: Reprojecting this way will not alter the datatype of the data!
            # Should be uint8 before and after reprojection.
            ds_tcvis = ds_tcvis.odc.reproject(ds_optical.odc.geobox, resampling="cubic")

        # !: Reprojecting with f64 coordinates and values behind the decimal point can result in a coordinate missmatch:
        # E.g. ds_optical has x coordinates [2.123, ...] then is can happen that the
        # reprojected ds_tcvis coordinates are [2.12300001, ...]
        # This results is all-nan assigments later when adding the variables of the reprojected dataset to the original
        assert (ds_optical.x == ds_tcvis.x).all(), "x coordinates do not match! See code comment above"
        assert (ds_optical.y == ds_tcvis.y).all(), "y coordinates do not match! See code comment above"

        # ?: Do ds_tcvis and ds_optical now share the same memory on the GPU?
        #  or do I need to delete ds_tcvis to free memory?
        # Same question for ArcticDEM
        ds_optical["tc_brightness"] = ds_tcvis.tc_brightness
        ds_optical["tc_greenness"] = ds_tcvis.tc_greenness
        ds_optical["tc_wetness"] = ds_tcvis.tc_wetness

    if ds_arcticdem is not None:
        # Calculate TPI and slope from ArcticDEM
        with stopwatch("Reprojecting ArcticDEM", printer=logger.debug):
            ds_arcticdem = ds_arcticdem.odc.reproject(
                ds_optical.odc.geobox.buffered(tpi_outer_radius), resampling="cubic"
            )
        # Move to same device as optical
        ds_arcticdem = move_to_device(ds_arcticdem, device)

        assert (ds_optical.x == ds_arcticdem.x).all(), "x coordinates do not match! See code comment above"
        assert (ds_optical.y == ds_arcticdem.y).all(), "y coordinates do not match! See code comment above"

        azimuth, angle_altitude = get_azimuth_and_elevation(ds_optical)
        ds_arcticdem = preprocess_arcticdem(
            ds_arcticdem,
            tpi_outer_radius,
            tpi_inner_radius,
            azimuth,
            angle_altitude,
        )
        ds_arcticdem = move_to_host(ds_arcticdem)

        # TODO: Check if crop can be done with apply_mask = False
        # -> Then the type conversion of the arcticdem data mask would not be necessary anymore
        # -> And this would also allow to keep the data on the GPU
        ds_arcticdem = ds_arcticdem.odc.crop(ds_optical.odc.geobox.extent)
        # For some reason, we need to reindex, because the reproject + crop of the arcticdem sometimes results
        # in floating point errors. These error are at the order of 1e-10, hence, way below millimeter precision.
        ds_arcticdem["x"] = ds_optical.x
        ds_arcticdem["y"] = ds_optical.y

        ds_optical["dem"] = ds_arcticdem.dem
        ds_optical["relative_elevation"] = ds_arcticdem.tpi
        ds_optical["slope"] = ds_arcticdem.slope
        ds_optical["hillshade"] = ds_arcticdem.hillshade
        ds_optical["aspect"] = ds_arcticdem.aspect
        ds_optical["curvature"] = ds_arcticdem.curvature
        ds_optical["arcticdem_data_mask"] = ds_arcticdem.arcticdem_data_mask.astype("uint8")

    return ds_optical