439 add proximity to anomaly tool (#440)

Author: nmaarnio

docs/raster_processing/proximity_to_anomaly.md

@@ -0,0 +1,3 @@
+# Proximity to anomaly
+
+::: eis_toolkit.raster_processing.proximity_to_anomaly

eis_toolkit/cli.py

@@ -1315,6 +1315,65 @@ def distance_to_anomaly_cli(
     typer.echo(f"Computing distance to anomaly completed, writing raster to {output_raster}.")
 
 
+# PROXIMITY TO ANOMALY
+@app.command()
+def proximity_to_anomaly_cli(
+    input_raster: INPUT_FILE_OPTION,
+    output_raster: OUTPUT_FILE_OPTION,
+    threshold_criteria: Annotated[ThresholdCriteria, typer.Option(case_sensitive=False)],
+    first_threshold_criteria_value: float = typer.Option(),
+    second_threshold_criteria_value: float = None,
+    max_distance: float = typer.Option(),
+    max_distance_value: float = 0.0,
+    anomaly_value: float = 1.0,
+):
+    """
+    Calculate proximity from each raster cell to nearest anomaly cell.
+
+    Uses only the first band of the raster.
+    """
+    # from sys import platform
+
+    from eis_toolkit.raster_processing.proximity_to_anomaly import proximity_to_anomaly
+
+    typer.echo("Progress: 10%")
+
+    if second_threshold_criteria_value is not None:
+        threshold_criteria_value = (first_threshold_criteria_value, second_threshold_criteria_value)
+    else:
+        threshold_criteria_value = first_threshold_criteria_value
+
+    with rasterio.open(input_raster) as raster:
+        typer.echo("Progress: 25%")
+        # Use optimized version if Windows
+        # if platform == "win32":
+        #     out_image, out_meta = proximity_to_anomaly_gdal(
+        #         anomaly_raster_profile=raster.profile,
+        #         anomaly_raster_data=raster.read(1),
+        #         threshold_criteria_value=threshold_criteria_value,
+        #         threshold_criteria=get_enum_values(threshold_criteria),
+        #         max_distance=max_distance,
+        #         scaling_range=(anomaly_value, max_distance_value),
+        #     )
+        # else:
+        out_image, out_meta = proximity_to_anomaly(
+            anomaly_raster_profile=raster.profile,
+            anomaly_raster_data=raster.read(1),
+            threshold_criteria_value=threshold_criteria_value,
+            threshold_criteria=get_enum_values(threshold_criteria),
+            max_distance=max_distance,
+            scaling_range=(anomaly_value, max_distance_value),
+        )
+
+    typer.echo("Progress: 75%")
+
+    with rasterio.open(output_raster, "w", **out_meta) as dest:
+        dest.write(out_image, 1)
+    typer.echo("Progress: 100%")
+
+    typer.echo(f"Computing proximity to anomaly completed, writing raster to {output_raster}.")
+
+
 # EXTRACT VALUES FROM RASTER
 @app.command()
 def extract_values_from_raster_cli(

eis_toolkit/raster_processing/proximity_to_anomaly.py

@@ -0,0 +1,106 @@
+from numbers import Number
+
+import numpy as np
+from beartype import beartype
+from beartype.typing import Literal, Tuple, Union
+from rasterio import profiles
+
+from eis_toolkit.raster_processing.distance_to_anomaly import distance_to_anomaly, distance_to_anomaly_gdal
+from eis_toolkit.transformations.linear import _min_max_scaling
+
+
+@beartype
+def proximity_to_anomaly(
+    anomaly_raster_profile: Union[profiles.Profile, dict],
+    anomaly_raster_data: np.ndarray,
+    threshold_criteria_value: Union[Tuple[Number, Number], Number],
+    threshold_criteria: Literal["lower", "higher", "in_between", "outside"],
+    max_distance: Number,
+    scaling_range: Tuple[Number, Number] = (1, 0),
+) -> Tuple[np.ndarray, Union[profiles.Profile, dict]]:
+    """Calculate proximity from raster cell to nearest anomaly.
+
+    The criteria for what is anomalous can be defined as a single number and
+    criteria text of "higher" or "lower". Alternatively, the definition can be
+    a range where values inside (criteria text of "within") or outside are
+    marked as anomalous (criteria text of "outside"). If anomaly_raster_profile does
+    contain "nodata" key, np.nan is assumed to correspond to nodata values.
+
+    Scales proximity values linearly in the given range. The first number in scale_range
+    denotes the value at the anomaly cells, the second at given maximum_distance.
+
+    Args:
+        anomaly_raster_profile: The raster profile in which the distances
+            to the nearest anomalous value are determined.
+        anomaly_raster_data: The raster data in which the distances
+            to the nearest anomalous value are determined.
+        threshold_criteria_value: Value(s) used to define anomalous.
+            If the threshold criteria requires a tuple of values,
+            the first value should be the minimum and the second
+            the maximum value.
+        threshold_criteria: Method to define anomalous.
+        max_distance: The maximum distance in the output array beyond which
+            proximity is considered 0.
+        scaling_range: Min and max values used for scaling the proximity values.
+            Defaults to (1, 0).
+
+    Returns:
+        A 2D numpy array with the distances to anomalies computed
+        and the original anomaly raster profile.
+    """
+    out_image, anomaly_raster_profile = distance_to_anomaly(
+        anomaly_raster_profile, anomaly_raster_data, threshold_criteria_value, threshold_criteria, max_distance
+    )
+    out_image = _min_max_scaling(out_image, scaling_range)
+
+    return out_image, anomaly_raster_profile
+
+
+@beartype
+def proximity_to_anomaly_gdal(
+    anomaly_raster_profile: Union[profiles.Profile, dict],
+    anomaly_raster_data: np.ndarray,
+    threshold_criteria_value: Union[Tuple[Number, Number], Number],
+    threshold_criteria: Literal["lower", "higher", "in_between", "outside"],
+    max_distance: Number,
+    scaling_range: Tuple[Number, Number] = (1, 0),
+) -> Tuple[np.ndarray, Union[profiles.Profile, dict]]:
+    """Calculate proximity from raster cell to nearest anomaly.
+
+    Uses an optimized, faster version of `distance_to_anomaly` in the background.
+    Available only on Windows.
+
+    The criteria for what is anomalous can be defined as a single number and
+    criteria text of "higher" or "lower". Alternatively, the definition can be
+    a range where values inside (criteria text of "within") or outside are
+    marked as anomalous (criteria text of "outside"). If anomaly_raster_profile does
+    contain "nodata" key, np.nan is assumed to correspond to nodata values.
+
+    Scales proximity values linearly in the given range. The first number in scale_range
+    denotes the value at the anomaly cells, the second at given maximum_distance.
+
+    Args:
+        anomaly_raster_profile: The raster profile in which the distances
+            to the nearest anomalous value are determined.
+        anomaly_raster_data: The raster data in which the distances
+            to the nearest anomalous value are determined.
+        threshold_criteria_value: Value(s) used to define anomalous.
+            If the threshold criteria requires a tuple of values,
+            the first value should be the minimum and the second
+            the maximum value.
+        threshold_criteria: Method to define anomalous.
+        max_distance: The maximum distance in the output array beyond which
+            proximity is considered 0.
+        scaling_range: Min and max values used for scaling the proximity values.
+            Defaults to (1, 0).
+
+    Returns:
+        A 2D numpy array with the distances to anomalies computed
+        and the original anomaly raster profile.
+    """
+    out_image, anomaly_raster_profile = distance_to_anomaly_gdal(
+        anomaly_raster_profile, anomaly_raster_data, threshold_criteria_value, threshold_criteria, max_distance
+    )
+    out_image = _min_max_scaling(out_image, scaling_range)
+
+    return out_image, anomaly_raster_profile