Skip to content

v2

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 the aspect of the terrain surface from an ArcticDEM Dataset.

Parameters:

  • arcticdem_ds (xarray.Dataset) –

    The ArcticDEM Dataset containing the 'dem' variable.

Returns:

  • xarray.Dataset –

    xr.Dataset: The input Dataset with the calculated aspect added as a new variable 'aspect'.

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py 
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
@stopwatch("Calculating aspect", printer=logger.debug)
def calculate_aspect(arcticdem_ds: xr.Dataset) -> xr.Dataset:
    """Calculate the aspect of the terrain surface from an ArcticDEM Dataset.

    Args:
        arcticdem_ds (xr.Dataset): The ArcticDEM Dataset containing the 'dem' variable.

    Returns:
        xr.Dataset: The input Dataset with the calculated aspect added as a new variable 'aspect'.

    """
    aspect_deg = aspect(arcticdem_ds.dem)
    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 the curvature of the terrain surface from an ArcticDEM Dataset.

Parameters:

  • arcticdem_ds (xarray.Dataset) –

    The ArcticDEM Dataset containing the 'dem' variable.

Returns:

  • xarray.Dataset –

    xr.Dataset: The input Dataset with the calculated curvature added as a new variable 'curvature'.

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py 
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
@stopwatch("Calculating curvature", printer=logger.debug)
def calculate_curvature(arcticdem_ds: xr.Dataset) -> xr.Dataset:
    """Calculate the curvature of the terrain surface from an ArcticDEM Dataset.

    Args:
        arcticdem_ds (xr.Dataset): The ArcticDEM Dataset containing the 'dem' variable.

    Returns:
        xr.Dataset: The input Dataset with the calculated curvature added as a new variable 'curvature'.

    """
    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 the hillshade of the terrain surface from an ArcticDEM Dataset.

Parameters:

  • arcticdem_ds (xarray.Dataset) –

    The ArcticDEM Dataset containing the 'dem' variable.

  • azimuth (int, default: 225 ) –

    The azimuth angle of the light source in degrees. Defaults to 225.

  • angle_altitude (int, default: 25 ) –

    The altitude angle of the light source in degrees. Defaults to 25.

Returns:

  • xarray.Dataset –

    xr.Dataset: The input Dataset with the calculated hillshade added as a new variable 'hillshade'.

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
@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 the hillshade of the terrain surface from an ArcticDEM Dataset.

    Args:
        arcticdem_ds (xr.Dataset): The ArcticDEM Dataset containing the 'dem' variable.
        azimuth (int, optional): The azimuth angle of the light source in degrees. Defaults to 225.
        angle_altitude (int, optional): The altitude angle of the light source in degrees. Defaults to 25.

    Returns:
        xr.Dataset: The input Dataset with the calculated hillshade added as a new variable 'hillshade'.

    """
    hillshade_da = hillshade(arcticdem_ds.dem, azimuth=azimuth, 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 from an xarray Dataset containing spectral bands.

This function will clip the NIR and Red bands to the range [0, 1] before calculating NDVI to avoid potential numerical instabilities from negative reflections.

Parameters:

  • optical (xarray.Dataset) –

    The xarray Dataset containing the spectral bands.

Returns:

  • xarray.Dataset –

    xr.DataArray: A new DataArray containing the calculated NDVI values.

Notes

NDVI (Normalized Difference Vegetation Index) is calculated using the formula: NDVI = (NIR - Red) / (NIR + Red)

Source code in darts-preprocessing/src/darts_preprocessing/engineering/indices.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
@stopwatch("Calculating NDVI", printer=logger.debug)
def calculate_ndvi(optical: xr.Dataset) -> xr.Dataset:
    """Calculate NDVI from an xarray Dataset containing spectral bands.

    This function will clip the NIR and Red bands to the range [0, 1] before calculating NDVI to avoid
    potential numerical instabilities from negative reflections.

    Args:
        optical (xr.Dataset): The xarray Dataset containing the spectral bands.

    Returns:
        xr.DataArray: A new DataArray containing the calculated NDVI values.

    Notes:
        NDVI (Normalized Difference Vegetation Index) is calculated using the formula:
            NDVI = (NIR - Red) / (NIR + Red)

    """
    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 the slope of the terrain surface from an ArcticDEM Dataset.

Parameters:

  • arcticdem_ds (xarray.Dataset) –

    The ArcticDEM Dataset containing the 'dem' variable.

Returns:

  • xarray.Dataset –

    xr.Dataset: The input Dataset with the calculated slope added as a new variable 'slope'.

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
@stopwatch("Calculating slope", printer=logger.debug)
def calculate_slope(arcticdem_ds: xr.Dataset) -> xr.Dataset:
    """Calculate the slope of the terrain surface from an ArcticDEM Dataset.

    Args:
        arcticdem_ds (xr.Dataset): The ArcticDEM Dataset containing the 'dem' variable.

    Returns:
        xr.Dataset: The input Dataset with the calculated slope added as a new variable 'slope'.

    """
    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.

Parameters:

  • arcticdem_ds (xarray.Dataset) –

    The ArcticDEM Dataset containing the 'dem' variable.

  • outer_radius (int) –

    The outer radius of the annulus kernel in m.

  • inner_radius (int) –

    The inner radius of the annulus kernel in m.

Returns:

  • xarray.Dataset –

    xr.Dataset: The input Dataset with the calculated TPI added as a new variable 'tpi'.

Source code in darts-preprocessing/src/darts_preprocessing/engineering/arcticdem.py 
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
@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.

    Args:
        arcticdem_ds (xr.Dataset): The ArcticDEM Dataset containing the 'dem' variable.
        outer_radius (int, optional): The outer radius of the annulus kernel in m.
        inner_radius (int, optional): The inner radius of the annulus kernel in m.

    Returns:
        xr.Dataset: The input Dataset with the calculated TPI added as a new variable 'tpi'.

    """
    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:

Returns:

Source code in darts-preprocessing/src/darts_preprocessing/v2.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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 
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
@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,
    ds_tcvis: xarray.Dataset,
    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.

The processing steps are: - Calculate NDVI - Calculate slope, hillshade, aspect, curvature and relative elevation from ArcticDEM - Merge everything into a single ds.

Parameters:

  • ds_optical (xarray.Dataset) –

    The Planet scene optical data or Sentinel 2 scene optical dataset including data_masks.

  • ds_arcticdem (xarray.Dataset) –

    The ArcticDEM dataset.

  • ds_tcvis (xarray.Dataset) –

    The TCVIS dataset.

  • tpi_outer_radius (int, default: 100 ) –

    The outer radius of the annulus kernel for the tpi calculation in m. Defaults to 100m.

  • tpi_inner_radius (int, default: 0 ) –

    The inner radius of the annulus kernel for the tpi calculation in m. Defaults to 0.

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

    The device to run the tpi and slope calculations on. If "cuda" take the first device (0), if int take the specified device. Defaults to "cuda" if cuda is available, else "cpu".

Returns:

Source code in darts-preprocessing/src/darts_preprocessing/v2.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
@stopwatch("Preprocessing", printer=logger.debug)
def preprocess_v2(
    ds_optical: xr.Dataset,
    ds_arcticdem: xr.Dataset,
    ds_tcvis: xr.Dataset,
    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.

    The processing steps are:
    - Calculate NDVI
    - Calculate slope, hillshade, aspect, curvature and relative elevation from ArcticDEM
    - Merge everything into a single ds.

    Args:
        ds_optical (xr.Dataset): The Planet scene optical data or Sentinel 2 scene optical dataset including data_masks.
        ds_arcticdem (xr.Dataset): The ArcticDEM dataset.
        ds_tcvis (xr.Dataset): The TCVIS dataset.
        tpi_outer_radius (int, optional): The outer radius of the annulus kernel for the tpi calculation
            in m. Defaults to 100m.
        tpi_inner_radius (int, optional): The inner radius of the annulus kernel for the tpi calculation
            in m. Defaults to 0.
        device (Literal["cuda", "cpu"] | int, optional): The device to run the tpi and slope calculations on.
            If "cuda" take the first device (0), if int take the specified device.
            Defaults to "cuda" if cuda is available, else "cpu".

    Returns:
        xr.Dataset: The preprocessed dataset.

    """
    # 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)

    # 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 delte 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

    # 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)

    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

    return ds_optical