From 4f9128c416860c8b1aa58c312c25e3d18f546df6 Mon Sep 17 00:00:00 2001 From: Niko Aarnio Date: Tue, 25 Mar 2025 12:13:24 +0200 Subject: [PATCH 1/4] feat(utils): add utility function for generating a full base raster profile --- eis_toolkit/cli.py | 83 ++++++--------------------------- eis_toolkit/utilities/raster.py | 51 ++++++++++++++++++++ 2 files changed, 64 insertions(+), 70 deletions(-) diff --git a/eis_toolkit/cli.py b/eis_toolkit/cli.py index c34be4c1..069b9488 100644 --- a/eis_toolkit/cli.py +++ b/eis_toolkit/cli.py @@ -15,6 +15,7 @@ import numpy as np import pandas as pd import rasterio +import rasterio.profiles import typer from beartype.typing import List, Optional, Sequence, Tuple, Union from typing_extensions import Annotated @@ -1963,23 +1964,14 @@ def idw_interpolation_cli( search_radius: Optional[float] = None, ): """Apply inverse distance weighting (IDW) interpolation to input vector file.""" - from eis_toolkit.exceptions import InvalidParameterValueException - from eis_toolkit.utilities.raster import profile_from_extent_and_pixel_size + from eis_toolkit.utilities.raster import base_profile from eis_toolkit.vector_processing.idw_interpolation import idw with ProgressLog.reading_input_files(): geodataframe = gpd.read_file(input_vector) if base_raster is None or base_raster == "": - if any(bound is None for bound in extent) or pixel_size is None or pixel_size <= 0: - raise InvalidParameterValueException( - "Expected positive pixel size and defined extent in absence of base raster. " - + f"Pixel size: {pixel_size}, extent: {extent}." - ) - profile = profile_from_extent_and_pixel_size(extent, (pixel_size, pixel_size)) - profile["crs"] = geodataframe.crs - profile["driver"] = "GTiff" - profile["dtype"] = "float32" + profile = base_profile(extent, pixel_size, geodataframe.crs) else: with rasterio.open(base_raster) as raster: profile = raster.profile.copy() @@ -1993,7 +1985,6 @@ def idw_interpolation_cli( search_radius=search_radius, ) - profile["count"] = 1 with ProgressLog.saving_output_files(output_raster): with rasterio.open(output_raster, "w", **profile) as dst: dst.write(out_image, 1) @@ -2015,23 +2006,14 @@ def kriging_interpolation_cli( method: Annotated[KrigingMethod, typer.Option(case_sensitive=False)] = KrigingMethod.ordinary, ): """Apply kriging interpolation to input vector file.""" - from eis_toolkit.exceptions import InvalidParameterValueException - from eis_toolkit.utilities.raster import profile_from_extent_and_pixel_size + from eis_toolkit.utilities.raster import base_profile from eis_toolkit.vector_processing.kriging_interpolation import kriging with ProgressLog.reading_input_files(): geodataframe = gpd.read_file(input_vector) if base_raster is None or base_raster == "": - if any(bound is None for bound in extent) or pixel_size is None or pixel_size <= 0: - raise InvalidParameterValueException( - "Expected positive pixel size and defined extent in absence of base raster. " - + f"Pixel size: {pixel_size}, extent: {extent}." - ) - profile = profile_from_extent_and_pixel_size(extent, (pixel_size, pixel_size)) - profile["crs"] = geodataframe.crs - profile["driver"] = "GTiff" - profile["dtype"] = "float32" + profile = base_profile(extent, pixel_size, geodataframe.crs) else: with rasterio.open(base_raster) as raster: profile = raster.profile.copy() @@ -2046,7 +2028,6 @@ def kriging_interpolation_cli( method=get_enum_values(method), ) - profile["count"] = 1 with ProgressLog.saving_output_files(output_raster): with rasterio.open(output_raster, "w", **profile) as dst: dst.write(out_image, 1) @@ -2073,23 +2054,14 @@ def rasterize_cli( Either base raster or pixel size + extent must be provided. """ - from eis_toolkit.exceptions import InvalidParameterValueException - from eis_toolkit.utilities.raster import profile_from_extent_and_pixel_size + from eis_toolkit.utilities.raster import base_profile from eis_toolkit.vector_processing.rasterize_vector import rasterize_vector with ProgressLog.reading_input_files(): geodataframe = gpd.read_file(input_vector) if base_raster is None or base_raster == "": - if any(bound is None for bound in extent) or pixel_size is None or pixel_size <= 0: - raise InvalidParameterValueException( - "Expected positive pixel size and defined extent in absence of base raster. " - + f"Pixel size: {pixel_size}, extent: {extent}." - ) - profile = profile_from_extent_and_pixel_size(extent, (pixel_size, pixel_size)) - profile["crs"] = geodataframe.crs - profile["driver"] = "GTiff" - profile["dtype"] = "float32" + profile = base_profile(extent, pixel_size, geodataframe.crs) else: with rasterio.open(base_raster) as raster: profile = raster.profile.copy() @@ -2105,7 +2077,6 @@ def rasterize_cli( get_enum_values(merge_strategy), ) - profile["count"] = 1 with ProgressLog.saving_output_files(output_raster): with rasterio.open(output_raster, "w", **profile) as dst: dst.write(out_image, 1) @@ -2151,23 +2122,14 @@ def vector_density_cli( Either base raster or pixel size + extent must be provided. """ - from eis_toolkit.exceptions import InvalidParameterValueException - from eis_toolkit.utilities.raster import profile_from_extent_and_pixel_size + from eis_toolkit.utilities.raster import base_profile from eis_toolkit.vector_processing.vector_density import vector_density with ProgressLog.reading_input_files(): geodataframe = gpd.read_file(input_vector) if base_raster is None or base_raster == "": - if any(bound is None for bound in extent) or pixel_size is None or pixel_size <= 0: - raise InvalidParameterValueException( - "Expected positive pixel size and defined extent in absence of base raster. " - + f"Pixel size: {pixel_size}, extent: {extent}." - ) - profile = profile_from_extent_and_pixel_size(extent, (pixel_size, pixel_size)) - profile["crs"] = geodataframe.crs - profile["driver"] = "GTiff" - profile["dtype"] = "float32" + profile = base_profile(extent, pixel_size, geodataframe.crs) else: with rasterio.open(base_raster) as raster: profile = raster.profile.copy() @@ -2180,7 +2142,6 @@ def vector_density_cli( statistic=get_enum_values(statistic), ) - profile["count"] = 1 with ProgressLog.saving_output_files(output_raster): with rasterio.open(output_raster, "w", **profile) as dst: dst.write(out_image, 1) @@ -2199,23 +2160,14 @@ def distance_computation_cli( max_distance: float = None, ): """Calculate distance from raster cell to nearest geometry.""" - from eis_toolkit.exceptions import InvalidParameterValueException - from eis_toolkit.utilities.raster import profile_from_extent_and_pixel_size + from eis_toolkit.utilities.raster import base_profile from eis_toolkit.vector_processing.distance_computation import distance_computation with ProgressLog.reading_input_files(): geodataframe = gpd.read_file(input_vector) if base_raster is None or base_raster == "": - if any(bound is None for bound in extent) or pixel_size is None or pixel_size <= 0: - raise InvalidParameterValueException( - "Expected positive pixel size and defined extent in absence of base raster. " - + f"Pixel size: {pixel_size}, extent: {extent}." - ) - profile = profile_from_extent_and_pixel_size(extent, (pixel_size, pixel_size)) - profile["crs"] = geodataframe.crs - profile["driver"] = "GTiff" - profile["dtype"] = "float32" + profile = base_profile(extent, pixel_size, geodataframe.crs) mask = None else: with rasterio.open(base_raster) as raster: @@ -2289,23 +2241,14 @@ def proximity_computation_cli( extent: Tuple[float, float, float, float] = (None, None, None, None), ): """Calculate proximity from raster cell to nearest geometry.""" - from eis_toolkit.exceptions import InvalidParameterValueException - from eis_toolkit.utilities.raster import profile_from_extent_and_pixel_size + from eis_toolkit.utilities.raster import base_profile from eis_toolkit.vector_processing.proximity_computation import proximity_computation with ProgressLog.reading_input_files(): geodataframe = gpd.read_file(input_vector) if base_raster is None or base_raster == "": - if any(bound is None for bound in extent) or pixel_size is None or pixel_size <= 0: - raise InvalidParameterValueException( - "Expected positive pixel size and defined extent in absence of base raster. " - + f"Pixel size: {pixel_size}, extent: {extent}." - ) - profile = profile_from_extent_and_pixel_size(extent, (pixel_size, pixel_size)) - profile["crs"] = geodataframe.crs - profile["driver"] = "GTiff" - profile["dtype"] = "float32" + profile = base_profile(extent, pixel_size, geodataframe.crs) mask = None else: with rasterio.open(base_raster) as raster: diff --git a/eis_toolkit/utilities/raster.py b/eis_toolkit/utilities/raster.py index 2228002f..648e36da 100644 --- a/eis_toolkit/utilities/raster.py +++ b/eis_toolkit/utilities/raster.py @@ -5,6 +5,7 @@ import rasterio from beartype import beartype from beartype.typing import Literal, Sequence, Tuple, Union +from pyproj import CRS from rasterio import profiles, transform from eis_toolkit.exceptions import ( @@ -170,3 +171,53 @@ def profile_from_extent_and_pixel_size( } raster_profile = profiles.Profile(raster_meta) return raster_profile + + +@beartype +def base_profile( + extent: Tuple[Number, Number, Number, Number], + pixel_size: Number, + crs: CRS, + dtype=np.float32, + count: int = 1, + nodata: float = -9999.0, + driver: str = "GTiff", + compress: str = "lzw", + round_strategy: Literal["nearest", "up", "down"] = "up", +) -> profiles.Profile: + """ + Create a raster profile from given extent, pixel size, CRS and other optional parameters. + + If extent and pixel size do not match exactly, i.e. raster width and height + calcalated from bounds and pixel size are not integers, rounding for the width and + height is performed. + + Args: + extent: Raster extent in the form (coord_west, coord_east, coord_south, coord_north). + pixel_size: Desired pixel size. If two values are provided, first is used for x and second for y. + If one value is provided, the value is used for both directions. + crs: Coordinate reference system. + dtype: Data type. + count: Number of raster bands. + nodata: Nodata value. + driver: Driver, defaults to GeoTIFF file format. + compress: Compression mode. + round_strategy: The rounding strategy if extent and pixel size do not match exactly. + Defaults to "up". + + Returns: + Rasterio profile. + """ + if any(bound is None for bound in extent) or pixel_size is None or pixel_size <= 0: + raise InvalidParameterValueException( + "Expected positive pixel size and defined extent in absence of base raster. " + + f"Pixel size: {pixel_size}, extent: {extent}." + ) + profile = profile_from_extent_and_pixel_size(extent, (pixel_size, pixel_size), round_strategy) + profile["crs"] = crs + profile["dtype"] = dtype + profile["count"] = count + profile["nodata"] = nodata + profile["driver"] = driver + profile["compress"] = compress + return profile From bd6ad780b87f47df973d1a890483af7ba82a617d Mon Sep 17 00:00:00 2001 From: Niko Aarnio Date: Tue, 25 Mar 2025 12:41:34 +0200 Subject: [PATCH 2/4] style: small docstring and param name modifications --- .../training_data_tools/points_to_raster.py | 1 + eis_toolkit/training_data_tools/random_sampling.py | 14 +++++++------- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/eis_toolkit/training_data_tools/points_to_raster.py b/eis_toolkit/training_data_tools/points_to_raster.py index 1e5689ca..bf5d8b61 100644 --- a/eis_toolkit/training_data_tools/points_to_raster.py +++ b/eis_toolkit/training_data_tools/points_to_raster.py @@ -124,6 +124,7 @@ def points_to_raster( Args: geodataframe: The geodataframe points set to be converted into raster. + raster_profile: The raster profile determining the output raster grid properties. attribute: Values to be be assigned to the geodataframe. radius: Radius to be applied around the geodataframe in [m]. buffer: Buffers the matrix value when two or more radii with different attribute value overlap. diff --git a/eis_toolkit/training_data_tools/random_sampling.py b/eis_toolkit/training_data_tools/random_sampling.py index 82114d56..8dc979b3 100644 --- a/eis_toolkit/training_data_tools/random_sampling.py +++ b/eis_toolkit/training_data_tools/random_sampling.py @@ -33,7 +33,7 @@ def _random_sampling( @beartype def generate_negatives( raster_array: np.ndarray, - raster_meta: Union[profiles.Profile, dict], + raster_profile: Union[profiles.Profile, dict], sample_number: Number, random_seed: int = 48, ) -> Tuple[gpd.GeoDataFrame, np.ndarray, Union[profiles.Profile, dict]]: @@ -41,15 +41,15 @@ def generate_negatives( Args: raster_array: Raster array with marked positives. - raster_meta: Raster metadata. - sample_number: maximum number of negatives to be generated. + raster_profile: The raster profile determining the output raster grid properties. + sample_number: Maximum number of negatives to be generated. random_seed: Seed for generating random negatives. Returns: A tuple containing the shapely points, output raster as a NumPy array and updated metadata. Raises: - EmptyDataException: The raster array is empty. + EmptyDataException: The raster array is empty. """ if raster_array.size == 0: @@ -80,11 +80,11 @@ def generate_negatives( out_array[row, col] = -1 - x, y = rasterio.transform.xy(raster_meta["transform"], row, col) + x, y = rasterio.transform.xy(raster_profile["transform"], row, col) points = [Point(x[i], y[i]) for i in range(len(x))] sample_negative = gpd.GeoDataFrame(geometry=points) - sample_negative.set_crs(raster_meta["crs"], allow_override=True, inplace=True) + sample_negative.set_crs(raster_profile["crs"], allow_override=True, inplace=True) - return sample_negative, out_array, raster_meta + return sample_negative, out_array, raster_profile From 59ee53b8f44e783d33f868888ab2aaac711e06b1 Mon Sep 17 00:00:00 2001 From: Niko Aarnio Date: Tue, 25 Mar 2025 12:42:23 +0200 Subject: [PATCH 3/4] feat(CLI): add Points to raster CLI function --- eis_toolkit/cli.py | 54 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/eis_toolkit/cli.py b/eis_toolkit/cli.py index 069b9488..b68f6b3b 100644 --- a/eis_toolkit/cli.py +++ b/eis_toolkit/cli.py @@ -356,6 +356,14 @@ class ReplaceCondition(str, Enum): greater_than_or_equal = "greater_than_or_equal" +class BufferOption(str, Enum): + """Buffer options.""" + + avg = "avg" + min = "min" + max = "max" + + INPUT_FILE_OPTION = Annotated[ Path, typer.Option( @@ -2275,6 +2283,52 @@ def proximity_computation_cli( ProgressLog.finish() +# --- TRAINING DATA TOOLS --- + + +@app.command() +def points_to_raster_cli( + input_vector: INPUT_FILE_OPTION, + output_raster: OUTPUT_FILE_OPTION, + base_raster: INPUT_FILE_OPTION = None, + pixel_size: float = None, + extent: Tuple[float, float, float, float] = (None, None, None, None), + attribute: Optional[str] = None, + radius: Optional[int] = None, + buffer: Annotated[BufferOption, typer.Option(case_sensitive=False)] = None, +): + """Convert a point data set into a binary raster.""" + from eis_toolkit.training_data_tools.points_to_raster import points_to_raster + from eis_toolkit.utilities.raster import base_profile + + with ProgressLog.reading_input_files(): + geodataframe = gpd.read_file(input_vector) + + if base_raster is None or base_raster == "": + profile = base_profile(extent, pixel_size, geodataframe.crs) + mask = None + else: + with rasterio.open(base_raster) as raster: + profile = raster.profile.copy() + raster_array = raster.read(1) + mask = (raster_array == profile["nodata"]) | np.isnan(raster_array) + + with ProgressLog.running_algorithm(): + out_image, out_profile = points_to_raster( + geodataframe=geodataframe, raster_profile=profile, attribute=attribute, radius=radius, buffer=buffer + ) + + # Apply nodata mask + if mask is not None: + out_image[mask] = out_profile["nodata"] + + with ProgressLog.saving_output_files(output_raster): + with rasterio.open(output_raster, "w", **out_profile) as dst: + dst.write(out_image, 1) + + ProgressLog.finish() + + # --- PREDICTION --- From a4a9c3b108dd1efdf9705d5a21b0d0b477f5e776 Mon Sep 17 00:00:00 2001 From: Niko Aarnio Date: Tue, 25 Mar 2025 13:11:41 +0200 Subject: [PATCH 4/4] feat(CLI): add Generate negatives CLI function --- eis_toolkit/cli.py | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/eis_toolkit/cli.py b/eis_toolkit/cli.py index b68f6b3b..8c35164a 100644 --- a/eis_toolkit/cli.py +++ b/eis_toolkit/cli.py @@ -447,7 +447,8 @@ def saving_output_files(savepath: Union[str, Sequence[str]]): # noqa: D102 yield if isinstance(savepath, Sequence): for file in savepath: - typer.echo(f"✅ Output file(s) saved to {file}\n") + typer.echo(f"✅ Output file(s) saved to {file}") + typer.echo(" ") else: typer.echo(f"✅ Output file(s) saved to {savepath}\n") @@ -2285,7 +2286,7 @@ def proximity_computation_cli( # --- TRAINING DATA TOOLS --- - +# POINTS TO RASTER @app.command() def points_to_raster_cli( input_vector: INPUT_FILE_OPTION, @@ -2329,6 +2330,38 @@ def points_to_raster_cli( ProgressLog.finish() +# GENERATE NEGATIVES +@app.command() +def generate_negatives_cli( + input_raster: INPUT_FILE_OPTION, + output_raster: OUTPUT_FILE_OPTION, + output_vector: OUTPUT_FILE_OPTION, + sample_number: Annotated[int, typer.Option()], + random_seed: int = 48, +): + """Generate probable negatives from raster array with marked positives.""" + from eis_toolkit.training_data_tools.random_sampling import generate_negatives + + with ProgressLog.reading_input_files(): + with rasterio.open(input_raster) as raster: + raster_array = raster.read(1) + profile = raster.profile.copy() + mask = (raster_array == profile["nodata"]) | np.isnan(raster_array) + + with ProgressLog.running_algorithm(): + out_points, out_image, out_profile = generate_negatives( + raster_array=raster_array, raster_profile=profile, sample_number=sample_number, random_seed=random_seed + ) + out_image[mask] = out_profile["nodata"] + + with ProgressLog.saving_output_files([output_raster, output_vector]): + with rasterio.open(output_raster, "w", **out_profile) as dst: + dst.write(out_image, 1) + out_points.to_file(output_vector) + + ProgressLog.finish() + + # --- PREDICTION ---