diff --git a/src/geogenalg/application/generalize_contours.py b/src/geogenalg/application/generalize_contours.py new file mode 100644 index 0000000..4fd77ca --- /dev/null +++ b/src/geogenalg/application/generalize_contours.py @@ -0,0 +1,110 @@ +# Copyright (c) 2026 National Land Survey of Finland (Maanmittauslaitos) +# +# This file is part of geogen-algorithms. +# +# SPDX-License-Identifier: MIT + +from typing import ClassVar + +from cartagen.utils import gaussian_smoothing +from geopandas import GeoDataFrame +from pydantic import Field + +from geogenalg.application import ( + BaseAlgorithm, + ReferenceDataInformation, + supports_identity, +) +from geogenalg.continuity import ( + add_contiguous_lines_information, + smooth_linestring_connections, +) +from geogenalg.core.geometry import assign_z_from_attribute +from geogenalg.identity import hash_duplicate_indexes +from geogenalg.split import split_lines_by_points + +SNAP_DISTANCE = 1.0 + + +@supports_identity +class GeneralizeContours(BaseAlgorithm): + """Generalize contour line geometries. + + Input should contain LineString geometries representing contours. + + Reference data should contain Point geometries representing + positions of slope lines. + + Output is a GeoDataFrame containing generalized contour lines. + + The algorithm does the following steps: + 1. Filter contours based on elevation interval. + 2. Split contours at slope reference points to preserve fixed locations. + 3. Apply Gaussian smoothing to contour geometries. + 4. Smooth connections between adjacent contour segments. + 5. Remove short contour lines. + """ + + interval: float = Field(5, gt=0) + """Elevation interval used for contour filtering.""" + gaussian_filter_strength: float = Field(8, ge=0) + """Sigma value used for Gaussian smoothing.""" + length_threshold: float = Field(200, ge=0) + """Minimum length for contour line.""" + level_attribute: str = Field("elevation_value") + """Attribute containing contour elevation values.""" + reference_key: str = Field("slope") + """Reference Point or MultiPoint data key for slope lines.""" + + valid_input_geometry_types: ClassVar = {"LineString"} + + reference_data_schema: ClassVar = { + "reference_key": ReferenceDataInformation( + required=True, + valid_geometry_types={"Point"}, + ), + } + + def _execute( + self, + data: GeoDataFrame, + reference_data: dict[str, GeoDataFrame], + ) -> GeoDataFrame: + gdf = data.copy() + reference_gdf = reference_data[self.reference_key] + gdf.geometry = gdf.geometry.force_2d() + + # Filter contours by elevation interval + gdf = gdf[gdf[self.level_attribute] % self.interval == 0] + + if gdf.empty: + return gdf.copy() + + # Split contours at slope line positions + gdf = split_lines_by_points(gdf, reference_gdf, SNAP_DISTANCE) + + # Gaussian smoothing + gdf.geometry = gdf.geometry.apply( + lambda geom: gaussian_smoothing(geom, sigma=self.gaussian_filter_strength) + ) + + # Smooth line connections + gdf = smooth_linestring_connections(gdf) + + # TODO: Handle potentially intersecting contours after smoothing + + # Add information about the total length of the continuous contour + gdf = add_contiguous_lines_information( + gdf, GeoDataFrame(geometry=[], crs=gdf.crs) + ) + + # Remove short contours + gdf = gdf[gdf["contiguous_length"] >= self.length_threshold] + + gdf = hash_duplicate_indexes(gdf, "contour") + gdf = assign_z_from_attribute(gdf, self.level_attribute, overwrite_z=True) + return gdf.drop( + [column for column in gdf.columns if column not in data.columns], axis=1 + ) + + # TODO: reduce the number of vertices diff --git a/src/geogenalg/split.py b/src/geogenalg/split.py index 3186858..f5ef786 100644 --- a/src/geogenalg/split.py +++ b/src/geogenalg/split.py @@ -119,4 +119,4 @@ def split_lines_by_points( return copy_gdf_as_empty(lines_gdf) result = combine_gdfs(result_rows, ignore_index=False) - return GeoDataFrame(result, crs=lines_gdf.crs) + return GeoDataFrame(result, geometry=lines_gdf.geometry.name, crs=lines_gdf.crs) diff --git a/test/application/test_generalize_contours.py b/test/application/test_generalize_contours.py new file mode 100644 index 0000000..6708392 --- /dev/null +++ b/test/application/test_generalize_contours.py @@ -0,0 +1,53 @@ +# Copyright (c) 2026 National Land Survey of Finland (Maanmittauslaitos) +# +# This file is part of geogen-algorithms. +# +# SPDX-License-Identifier: MIT + + +from pathlib import Path + +from conftest import IntegrationTest + +from geogenalg.application.generalize_contours import GeneralizeContours +from geogenalg.testing import GeoPackagePath +from geogenalg.utility.dataframe_processing import read_gdf_from_file_and_set_index + +UNIQUE_ID_COLUMN = "kmtk_id" + + +def test_generalize_contours( + testdata_path: Path, + tmp_path: Path, +) -> None: + gpkg = GeoPackagePath(testdata_path / "contours.gpkg") + + slope_data = read_gdf_from_file_and_set_index( + testdata_path / "contours.gpkg", + UNIQUE_ID_COLUMN, + layer="slope_line", + ) + + slope_path = tmp_path / "slope.gpkg" + slope_data.to_file(slope_path, layer="slope") + + IntegrationTest( + input_uri=gpkg.to_input("contour"), + control_uri=gpkg.to_input("contour_control"), + algorithm=GeneralizeContours( + interval=5, + gaussian_filter_strength=8, + length_threshold=200, + level_attribute="n60_elevation_value", + reference_key="slope", + ), + unique_id_column=UNIQUE_ID_COLUMN, + check_missing_reference=True, + reference_uris={ + "slope": GeoPackagePath(slope_path).to_input("slope"), + }, + assert_function_arguments={ + "check_less_precise": True, + }, + dummy_data_mandatory_columns=["n60_elevation_value"], + ).run() diff --git a/test/testdata/contours.gpkg b/test/testdata/contours.gpkg new file mode 100644 index 0000000..0f48bf5 Binary files /dev/null and b/test/testdata/contours.gpkg differ