diff --git a/docs/raster_processing/proximity_to_anomaly.md b/docs/raster_processing/proximity_to_anomaly.md new file mode 100644 index 00000000..e4f9a2e1 --- /dev/null +++ b/docs/raster_processing/proximity_to_anomaly.md @@ -0,0 +1,3 @@ +# Proximity to anomaly + +::: eis_toolkit.raster_processing.proximity_to_anomaly diff --git a/eis_toolkit/cli.py b/eis_toolkit/cli.py index 3e3f0c21..9d77bcb1 100644 --- a/eis_toolkit/cli.py +++ b/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( diff --git a/eis_toolkit/raster_processing/proximity_to_anomaly.py b/eis_toolkit/raster_processing/proximity_to_anomaly.py new file mode 100644 index 00000000..50a1a0a4 --- /dev/null +++ b/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