diff --git a/.gitignore b/.gitignore index 0a19790..b92f823 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ +# GeoPackage auxiliary files +*.gpkg-shm +*.gpkg-wal + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/pyproject.toml b/pyproject.toml index a44bf53..73f32cf 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,12 +29,12 @@ dynamic = [ ] dependencies = [ "cartagen>=1.0.2", - "geopandas==1.0.1", + "geopandas==1.1.1", "matplotlib==3.9.2", "networkx==3.4.2", "pygeoops==0.4.0", "scikit-learn==1.6.1", - "shapely==2.0.6", + "shapely==2.1.2", "shapelysmooth==0.2.1", "simplification==0.7.13", "typer>=0.19.2", @@ -190,11 +190,6 @@ select = ["SC"] per-file-ignores = ["test/*:SC200"] # TODO: remove exclude section when qgis dependencies removed/function refactored exclude = """ - src/geogenalg/application/generalize_lakes.py - src/geogenalg/application/generalize_seas.py - src/geogenalg/buffer/exaggerate_thin_parts.py - src/geogenalg/buffer/generalize_islands.py - src/geogenalg/core/extract_interior_rings.py src/geogenalg/application/remove_short_watercourses.py src/geogenalg/application/remove_dense_watercourses.py src/geogenalg/application/generalize_watercourses_50k.py @@ -228,11 +223,6 @@ source_pkgs = ["geogenalg"] # TODO: remove omit section when qgis dependencies removed/function refactored omit = [ "src/geogenalg/main.py", - "src/geogenalg/application/generalize_lakes.py", - "src/geogenalg/application/generalize_seas.py", - "src/geogenalg/buffer/exaggerate_thin_parts.py", - "src/geogenalg/buffer/generalize_islands.py", - "src/geogenalg/core/extract_interior_rings.py", "src/geogenalg/application/remove_short_watercourses.py", "src/geogenalg/application/remove_dense_watercourses.py", "src/geogenalg/application/generalize_watercourses_50k.py", diff --git a/src/geogenalg/application/generalize_lakes.py b/src/geogenalg/application/generalize_lakes.py deleted file mode 100644 index 2f8381a..0000000 --- a/src/geogenalg/application/generalize_lakes.py +++ /dev/null @@ -1,507 +0,0 @@ -# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) -# -# This file is part of geogen-algorithms. -# -# This source code is licensed under the MIT license found in the -# LICENSE file in the root directory of this source tree. - -from yleistys_qgis_plugin.core.utils import ( - VISVALINGAM_ENUM_VERSION_CUTOFF, - extract_interior_rings_from_layer, - get_subalgorithm_result_layer, -) - -try: - import processing -except ImportError: - from qgis import processing - -from typing import Any - -from qgis.core import ( - Qgis, - QgsExpression, - QgsFeature, - QgsFeatureIterator, - QgsFeatureRequest, - QgsFeatureSink, - QgsMapToPixelSimplifier, - QgsProcessing, - QgsProcessingAlgorithm, - QgsProcessingContext, - QgsProcessingFeatureSourceDefinition, - QgsProcessingFeedback, - QgsProcessingParameterBoolean, - QgsProcessingParameterFeatureSink, - QgsProcessingParameterFeatureSource, - QgsProcessingParameterNumber, - QgsVectorLayer, - QgsWkbTypes, - edit, -) -from qgis.PyQt.QtCore import QCoreApplication - - -class GeneralizeLakes(QgsProcessingAlgorithm): - """This is an example algorithm that takes a vector layer and - creates a new identical one. - - It is meant to be used as an example of how to create your own - algorithms and explain methods and variables used to do it. An - algorithm like this will be available in all elements, and there - is not need for additional work. - - All Processing algorithms should extend the QgsGeneralizeLakes - class. - """ - - # Constants used to refer to parameters and outputs. They will be - # used when calling the algorithm from another algorithm, or when - # calling from the QGIS console. - - INPUT = "INPUT" - OUTPUT = "OUTPUT" - MIN_LAKE_AREA = "MIN_LAKE_AREA" - MIN_ISLAND_AREA = "MIN_ISLAND_AREA" - MAX_SIMPLIFICATION_TOLERANCE = "MAX_SIMPLIFICATION_TOLERANCE" - EXAGGERATE_THIN_PARTS = "EXAGGERATE_THIN_PARTS" - - def __init__(self) -> None: - super().__init__() - - self._name = "generalizelakes" - self._display_name = "Generalize lakes" - self._group_id = "generalization" - self._group = "Generalization" - self._short_help_string = "" - - def tr(self, string) -> str: - """Returns a translatable string with the self.tr() function.""" - return QCoreApplication.translate("Processing", string) - - def createInstance(self): # noqa: N802 - return GeneralizeLakes() - - def name(self) -> str: - """Returns the algorithm name, used for identifying the algorithm. This - string should be fixed for the algorithm, and must not be localised. - The name should be unique within each provider. Names should contain - lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._name - - def displayName(self) -> str: # noqa: N802 - """Returns the translated algorithm name, which should be used for any - user-visible display of the algorithm name. - """ - return self.tr(self._display_name) - - def groupId(self) -> str: # noqa: N802 - """Returns the unique ID of the group this algorithm belongs to. This - string should be fixed for the algorithm, and must not be localised. - The group id should be unique within each provider. Group id should - contain lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._group_id - - def group(self) -> str: - """Returns the name of the group this algorithm belongs to. This string - should be localised. - """ - return self.tr(self._group) - - def shortHelpString(self) -> str: # noqa: N802 - """Returns a localised short helper string for the algorithm. This string - should provide a basic description about what the algorithm does and the - parameters and outputs associated with it.. - """ - return self.tr(self._short_help_string) - - def initAlgorithm(self, configuration=None) -> None: # noqa: N802 - """Here we define the inputs and output of the algorithm, along - with some other properties. - """ - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT, - self.tr("Input layer"), - [QgsProcessing.SourceType.TypeVectorPolygon], - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MIN_LAKE_AREA, - self.tr("Remove lakes with an area under"), - defaultValue=0, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MIN_ISLAND_AREA, - self.tr("Remove islands with an area under"), - defaultValue=100, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MAX_SIMPLIFICATION_TOLERANCE, - self.tr("Maximum simplification tolerance"), - defaultValue=10, - minValue=0, - ) - ) - - # used for comparing results - self.addParameter( - QgsProcessingParameterBoolean( - self.EXAGGERATE_THIN_PARTS, - self.tr("Exaggerate thin parts"), - defaultValue=True, - ) - ) - - self.addParameter( - QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr("Generalized")) - ) - - def processAlgorithm( # noqa: N802 - self, - parameters: dict[str, Any], - context: QgsProcessingContext, - feedback: QgsProcessingFeedback, - ) -> dict: - """Here is where the processing itself takes place.""" - # Initialize feedback if it is None - if feedback is None: - feedback = QgsProcessingFeedback() - - def check_feature_count_difference( - layer_1: QgsVectorLayer, layer_2: QgsVectorLayer, *, warn_if_different: bool - ) -> None: - difference = layer_1.featureCount() - layer_2.featureCount() - - log = ( - feedback.pushWarning - if warn_if_different and difference != 0 - else feedback.pushInfo - ) - - log( - f"Difference in feature count is: {difference} ({layer_1.featureCount()}/{layer_2.featureCount()})" - ) - - input_layer: QgsVectorLayer = self.parameterAsVectorLayer( - parameters, self.INPUT, context - ) - min_lake_area: int = self.parameterAsInt( - parameters, self.MIN_LAKE_AREA, context - ) - min_island_area: int = self.parameterAsInt( - parameters, self.MIN_ISLAND_AREA, context - ) - max_simplification_tolerance: float = self.parameterAsDouble( - parameters, self.MAX_SIMPLIFICATION_TOLERANCE, context - ) - exaggerate_thin_parts: bool = self.parameterAsBoolean( - parameters, self.EXAGGERATE_THIN_PARTS, context - ) - - # copy layer since we might have to delete some features and we don't want - # to affect the actual data source - copied_layer = input_layer.materialize(QgsFeatureRequest()) - - (sink, dest_id) = self.parameterAsSink( - parameters, - self.OUTPUT, - context, - copied_layer.fields(), - copied_layer.wkbType(), - copied_layer.sourceCrs(), - ) - - feedback.setProgress(0) - feedback.setProgressText("Selecting features to merge") - - # Select all features that touch another feature with processing script - processing.run( # pyright: ignore - "native:selectbylocation", - { - "INPUT": input_layer, - "PREDICATE": [4, 5], # touch or overlap - "INTERSECT": input_layer, # layer to compare features to - "METHOD": 0, # create a new selection - }, - context=context, - feedback=feedback, - is_child_algorithm=True, - ) - - feedback.pushInfo( - f"Selected {input_layer.selectedFeatureCount()}/{input_layer.featureCount()} features to be merged." - ) - - feedback.setProgressText("Merging features") - - # Looks like we have to use the input layer here because using the - # temporary copied layer causes problems, or would have to be - # added to the project layer tree which is not wanted - - if input_layer.selectedFeatureCount() > 0: - dissolve_layer = get_subalgorithm_result_layer( - "native:dissolve", - { - "INPUT": QgsProcessingFeatureSourceDefinition( - input_layer.source(), - selectedFeaturesOnly=True, - ), - "SEPARATE_DISJOINT": False, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - combined_features_layer = get_subalgorithm_result_layer( - "native:multiparttosingleparts", - { - "INPUT": dissolve_layer, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - with edit(copied_layer): - copied_layer.deleteFeatures(input_layer.selectedFeatureIds()) - - fid_field_idx: int = copied_layer.fields().indexOf("fid") - change_fid: bool = fid_field_idx != -1 - - running_fid = 0 - if change_fid: - running_fid = copied_layer.maximumValue(fid_field_idx) - - combined_feature: QgsFeature - for combined_feature in combined_features_layer.getFeatures(): - if change_fid: - combined_feature.setAttribute("fid", running_fid + 1) - running_fid += 1 - copied_layer.addFeature(QgsFeature(combined_feature)) - - input_layer.removeSelection() - - check_feature_count_difference( - input_layer, copied_layer, warn_if_different=False - ) - - feedback.setProgressText("Eliminating features") - - # Delete small lakes - eliminate_expression = QgsExpression(f"$area >= {min_lake_area}") - no_small_lakes = copied_layer.materialize( - QgsFeatureRequest(eliminate_expression) - ) - - # Delete small islands - eliminated_layer = get_subalgorithm_result_layer( - "native:deleteholes", - { - "INPUT": no_small_lakes, - "MIN_AREA": min_island_area, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - check_feature_count_difference( - no_small_lakes, eliminated_layer, warn_if_different=True - ) - - feedback.setProgressText("Extracting islands") - - island_layer = extract_interior_rings_from_layer(eliminated_layer) - # Extract all islands from the layer - # island_layer = get_subalgorithm_result_layer( - # "yleistys_qgis_plugin:extractinteriorrings", - # { - # "INPUT": eliminated_layer, - # "FEATURE_PER_RING": True, - # "OUTPUT": "TEMPORARY_OUTPUT", - # }, - # context, - # feedback, - # ) - - feedback.pushInfo(f"Extracted {island_layer.featureCount()} islands.") - - generalized_islands = get_subalgorithm_result_layer( - "yleistys_qgis_plugin:generalizeislands", - { - "INPUT": island_layer, - "MIN_WIDTH": 185, - "MIN_ELONGATION": 3.349, - "MAX_SIMPLIFICATION_TOLERANCE": max_simplification_tolerance, - "EXAGGERATE_BY": 3, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - check_feature_count_difference( - island_layer, generalized_islands, warn_if_different=True - ) - - # Delete all islands from the lakes. This is done because - # islands and lakes need to be processed differently. - # The islands will be added back in later. - exterior_ring_lakes = get_subalgorithm_result_layer( - "native:deleteholes", - { - "INPUT": eliminated_layer, - "MIN_AREA": 0.0, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - check_feature_count_difference( - exterior_ring_lakes, eliminated_layer, warn_if_different=True - ) - - if exaggerate_thin_parts: - feedback.setProgressText("Exaggerating thin parts") - - exterior_ring_lakes = get_subalgorithm_result_layer( - "yleistys_qgis_plugin:exaggeratethinparts", - { - "INPUT": exterior_ring_lakes, - "WIDTH_TOLERANCE": 20, - "EXAGGERATE_BY": 3, - "MAX_THINNESS": 0.5, - "MIN_AREA": 200, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - check_feature_count_difference( - exterior_ring_lakes, exterior_ring_lakes, warn_if_different=True - ) - - # After exaggerating thin parts simplify and smooth lakes - - feedback.setProgressText("Simplifying and smoothing lakes") - - if Qgis.QGIS_VERSION_INT < VISVALINGAM_ENUM_VERSION_CUTOFF: # pyright: ignore - simplifier = QgsMapToPixelSimplifier( - QgsMapToPixelSimplifier.SimplifyFlag.SimplifyGeometry, - 1, - QgsMapToPixelSimplifier.SimplifyAlgorithm.Visvalingam, - ) - else: - simplifier = QgsMapToPixelSimplifier( - QgsMapToPixelSimplifier.SimplifyFlag.SimplifyGeometry, - 1, - Qgis.VectorSimplificationAlgorithm.Visvalingam, # pyright: ignore - ) - - total = ( - 100.0 / exterior_ring_lakes.featureCount() - if exterior_ring_lakes.featureCount() - else 0 - ) - - with edit(exterior_ring_lakes): - features: QgsFeatureIterator = exterior_ring_lakes.getFeatures() - for current, feature in enumerate(features): - if feedback.isCanceled(): - break - - geom = feature.geometry() - - vertex_count: int = geom.constGet().nCoordinates() - add_to_tolerance: float = vertex_count / 200 - add_to_tolerance = min( - add_to_tolerance, max_simplification_tolerance - 1 - ) - - simplifier.setTolerance(1 + add_to_tolerance) - - new_geom = simplifier.simplify(geom) - new_geom.removeDuplicateNodes() - smoothed_geom = new_geom.smooth() - - exterior_ring_lakes.changeGeometry(feature.id(), smoothed_geom) - - feedback.setProgress(int(current * total)) - - feedback.pushInfo( - f"Lake features after simplifying and smoothing: {exterior_ring_lakes.featureCount()}" - ) - - # Both lakes and islands have been processed, - # add islands back in - - feedback.setProgressText("Adding generalized islands") - - difference_layer = get_subalgorithm_result_layer( - "native:difference", - { - "INPUT": exterior_ring_lakes, - "OVERLAY": generalized_islands, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - check_feature_count_difference( - difference_layer, exterior_ring_lakes, warn_if_different=True - ) - - if QgsWkbTypes.isSingleType(input_layer.wkbType()): - feedback.setProgressText("Converting to single polygons") - # the difference algorithm changes the geometries to - # multipolygons, convert back to single polygon - generalized_lakes = get_subalgorithm_result_layer( - "native:multiparttosingleparts", - { - "INPUT": difference_layer, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - check_feature_count_difference( - difference_layer, generalized_lakes, warn_if_different=True - ) - else: - generalized_lakes = difference_layer - - feedback.setProgress(0) - total = ( - 100.0 / generalized_lakes.featureCount() - if generalized_lakes.featureCount() - else 0 - ) - - for current, feature in enumerate(generalized_lakes.getFeatures()): - if feedback.isCanceled(): - break - - sink.addFeature(QgsFeature(feature), QgsFeatureSink.Flag.FastInsert) - - feedback.setProgress(int(current * total)) - - return {self.OUTPUT: dest_id} diff --git a/src/geogenalg/application/generalize_seas.py b/src/geogenalg/application/generalize_seas.py deleted file mode 100644 index 4638d97..0000000 --- a/src/geogenalg/application/generalize_seas.py +++ /dev/null @@ -1,440 +0,0 @@ -# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) -# -# This file is part of geogen-algorithms. -# -# This source code is licensed under the MIT license found in the -# LICENSE file in the root directory of this source tree. - -from typing import Any - -from qgis.core import ( - Qgis, - QgsExpression, - QgsFeatureIterator, - QgsFeatureRequest, - QgsFeatureSink, - QgsFeatureSource, - QgsGeometry, - QgsMapToPixelSimplifier, - QgsProcessing, - QgsProcessingAlgorithm, - QgsProcessingContext, - QgsProcessingFeedback, - QgsProcessingParameterBoolean, - QgsProcessingParameterFeatureSink, - QgsProcessingParameterFeatureSource, - QgsProcessingParameterNumber, - QgsVectorLayer, - QgsWkbTypes, - edit, -) -from qgis.PyQt.QtCore import QCoreApplication -from yleistys_qgis_plugin.core.utils import ( - VISVALINGAM_ENUM_VERSION_CUTOFF, - copy_feature_with_geometry, - get_subalgorithm_result_layer, -) - - -class GeneralizeSeas(QgsProcessingAlgorithm): - """This is an example algorithm that takes a vector layer and - creates a new identical one. - - It is meant to be used as an example of how to create your own - algorithms and explain methods and variables used to do it. An - algorithm like this will be available in all elements, and there - is not need for additional work. - - All Processing algorithms should extend the QgsProcessingAlgorithm - class. - """ - - # Constants used to refer to parameters and outputs. They will be - # used when calling the algorithm from another algorithm, or when - # calling from the QGIS console. - - INPUT = "INPUT" - OUTPUT = "OUTPUT" - MIN_ISLAND_AREA = "MIN_ISLAND_AREA" - MAX_SIMPLIFICATION_TOLERANCE = "MAX_SIMPLIFICATION_TOLERANCE" - SIMPLIFY_TOLERANCE = "SIMPLIFY_TOLERANCE" - EXAGGERATE_THIN_PARTS = "EXAGGERATE_THIN_PARTS" - - def __init__(self) -> None: - super().__init__() - - self._name = "generalizeseas" - self._display_name = "Generalize seas" - self._group_id = "generalization" - self._group = "Generalization" - self._short_help_string = "" - - def tr(self, string) -> str: - """Returns a translatable string with the self.tr() function.""" - return QCoreApplication.translate("Processing", string) - - def createInstance(self): # noqa: N802 - return GeneralizeSeas() - - def name(self) -> str: - """Returns the algorithm name, used for identifying the algorithm. This - string should be fixed for the algorithm, and must not be localised. - The name should be unique within each provider. Names should contain - lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._name - - def displayName(self) -> str: # noqa: N802 - """Returns the translated algorithm name, which should be used for any - user-visible display of the algorithm name. - """ - return self.tr(self._display_name) - - def groupId(self) -> str: # noqa: N802 - """Returns the unique ID of the group this algorithm belongs to. This - string should be fixed for the algorithm, and must not be localised. - The group id should be unique within each provider. Group id should - contain lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._group_id - - def group(self) -> str: - """Returns the name of the group this algorithm belongs to. This string - should be localised. - """ - return self.tr(self._group) - - def shortHelpString(self) -> str: # noqa: N802 - """Returns a localised short helper string for the algorithm. This string - should provide a basic description about what the algorithm does and the - parameters and outputs associated with it.. - """ - return self.tr(self._short_help_string) - - def initAlgorithm(self, configuration=None) -> None: # noqa: N802 - """Here we define the inputs and output of the algorithm, along - with some other properties. - """ - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT, - self.tr("Input layer"), - [QgsProcessing.SourceType.TypeVectorPolygon], - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.SIMPLIFY_TOLERANCE, - self.tr("Simplification tolerance"), - defaultValue=10, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MIN_ISLAND_AREA, - self.tr("Island area threshold"), - defaultValue=100, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MAX_SIMPLIFICATION_TOLERANCE, - self.tr("Maximum simplification tolerance"), - defaultValue=10, - minValue=0, - ) - ) - - # used for comparing results - self.addParameter( - QgsProcessingParameterBoolean( - self.EXAGGERATE_THIN_PARTS, - self.tr("Exaggerate thin parts"), - defaultValue=True, - ) - ) - - self.addParameter( - QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr("Generalized")) - ) - - def processAlgorithm( # noqa: N802 - self, - parameters: dict[str, Any], - context: QgsProcessingContext, - feedback: QgsProcessingFeedback, - ) -> dict: - """Here is where the processing itself takes place.""" - input_layer: QgsFeatureSource = self.parameterAsVectorLayer( - parameters, self.INPUT, context - ) - min_island_area: int = self.parameterAsInt( - parameters, self.MIN_ISLAND_AREA, context - ) - max_simplification_tolerance: float = self.parameterAsDouble( - parameters, self.MAX_SIMPLIFICATION_TOLERANCE, context - ) - simplify_tolerance: int = self.parameterAsInt( - parameters, self.SIMPLIFY_TOLERANCE, context - ) - exaggerate_thin_parts: bool = self.parameterAsBoolean( - parameters, self.EXAGGERATE_THIN_PARTS, context - ) - - sink, dest_id = self.parameterAsSink( - parameters, - self.OUTPUT, - context, - input_layer.fields(), - input_layer.wkbType(), - input_layer.sourceCrs(), - ) - - if feedback is None: - feedback = QgsProcessingFeedback() - - # There's multiple stages where dissolving is required, - # save for convenience and clarity - def _dissolve(layer: QgsVectorLayer) -> QgsVectorLayer: - params = { - "INPUT": layer, - "SEPARATE_DISJOINT": False, - "OUTPUT": "TEMPORARY_OUTPUT", - } - - return get_subalgorithm_result_layer( - "native:dissolve", params, context, feedback - ) - - # There's multiple stages where the multipart to single part algorithm - # is required, save for convenience and clarity - def _multiparttosingleparts(layer: QgsVectorLayer) -> QgsVectorLayer:# noqa: SC200 - params = { - "INPUT": layer, - "OUTPUT": "TEMPORARY_OUTPUT", - } - - return get_subalgorithm_result_layer( - "native:multiparttosingleparts", params, context, feedback - ) - - # Part of the process is to simplify and smooth the polygons. - # The sea features also include the territorial water borders. - # These features are large along with having a small number - # of vertices which causes the smoothing operation to distort - # the borders significantly. - # For now handle this issue by extracting the relevant features - # and combine them with the result at the end. - expression = QgsExpression("territorial_waters_category_id > 1") - categories_dissolved = _dissolve( - input_layer.materialize(QgsFeatureRequest(expression)) - ) - territorial_exterior_ring = get_subalgorithm_result_layer( - "native:deleteholes", - { - "INPUT": categories_dissolved, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - # Processing features one-by-one or in clusters introduces - # many issues, such as having to deal with islands which - # are not holes but instead between multiple features and the - # fact that exaggerating thin parts may modify the features' - # geometries and topology such that it'll create new holes - # between the features. - - # Dissolving all adjacent features to one large feature will - # make processing significantly slower, but it solves most of - # these issues - layer = _dissolve(input_layer) - layer = _multiparttosingleparts(layer) - - feedback.setProgressText("Eliminating features") - - # Remove small islands first - eliminated_layer = get_subalgorithm_result_layer( - "native:deleteholes", - { - "INPUT": layer, - "MIN_AREA": min_island_area, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - # Exaggerate thin parts. The islands are added back in - # as interior rings later and this operation may change - # the geometry of an island so this has to be done before - # extracting the islands, otherwise its effect would be - # nullified. - if exaggerate_thin_parts: - feedback.setProgressText("Exaggerating thin parts") - - eliminated_layer = get_subalgorithm_result_layer( - "yleistys_qgis_plugin:exaggeratethinparts", - { - "INPUT": eliminated_layer, - "WIDTH_TOLERANCE": 20, - "EXAGGERATE_BY": 3, - "MAX_THINNESS": 0.5, - "MIN_AREA": 200, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - # The exterior ring and islands need to be processed differently, - # f.e. small islands should't be simplified as much, so extract - # all islands at this point and process them separately - extracted_islands = get_subalgorithm_result_layer( - "yleistys_qgis_plugin:extractinteriorrings", - { - "INPUT": eliminated_layer, - "FEATURE_PER_RING": True, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - generalized_islands = get_subalgorithm_result_layer( - "yleistys_qgis_plugin:generalizeislands", - { - "INPUT": extracted_islands, - "MIN_WIDTH": 185, - "MIN_ELONGATION": 3.349, - "MAX_SIMPLIFICATION_TOLERANCE": max_simplification_tolerance, - "EXAGGERATE_BY": 3, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - # Delete the islands to process the exterior ring separately - exterior_ring = get_subalgorithm_result_layer( - "native:deleteholes", - { - "INPUT": eliminated_layer, - "MIN_AREA": 0.0, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - if Qgis.QGIS_VERSION_INT < VISVALINGAM_ENUM_VERSION_CUTOFF: # pyright: ignore - simplifier = QgsMapToPixelSimplifier( - QgsMapToPixelSimplifier.SimplifyFlag.SimplifyGeometry, - simplify_tolerance, - QgsMapToPixelSimplifier.SimplifyAlgorithm.Visvalingam, - ) - else: - simplifier = QgsMapToPixelSimplifier( - QgsMapToPixelSimplifier.SimplifyFlag.SimplifyGeometry, - simplify_tolerance, - Qgis.VectorSimplificationAlgorithm.Visvalingam, # pyright: ignore - ) - - feedback.setProgressText("Simplifying and smoothing") - - total = ( - 100.0 / exterior_ring.featureCount() if exterior_ring.featureCount() else 0 - ) - - # Perform simplification and smoothing feature-by-feature - with edit(exterior_ring): - features: QgsFeatureIterator = exterior_ring.getFeatures() - for current, feature in enumerate(features): - if feedback.isCanceled(): - break - - geom = feature.geometry() - - # simplify geometries with more vertices a greater amount, - # capping at max_simplification_tolerance - vertex_count: int = geom.constGet().nCoordinates() - add_to_tolerance: float = vertex_count / 200 - add_to_tolerance = min( - add_to_tolerance, max_simplification_tolerance - 1 - ) - - simplifier.setTolerance(1 + add_to_tolerance) - - new_geom = simplifier.simplify(geom) - - # simplification may introduce duplicate nodes, delete - new_geom.removeDuplicateNodes() - - new_geom = new_geom.smooth() - - exterior_ring.changeGeometry(feature.id(), new_geom) - - feedback.setProgress(int(current * total)) - - feedback.setProgressText("Adding generalized islands") - - generalized_sea = get_subalgorithm_result_layer( - "native:difference", - { - "INPUT": exterior_ring, - "OVERLAY": generalized_islands, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - diff = get_subalgorithm_result_layer( - "native:difference", - { - "INPUT": generalized_sea, - "OVERLAY": territorial_exterior_ring, - "OUTPUT": "TEMPORARY_OUTPUT", - }, - context, - feedback, - ) - - diff_singles = _multiparttosingleparts(diff) - - biggest_geom: QgsGeometry = QgsGeometry() - for feature in diff_singles.getFeatures(): - geom = feature.geometry() - if geom.area() > biggest_geom.area(): - biggest_geom = geom - - combined = territorial_exterior_ring.getGeometry(1).combine(biggest_geom) - - with edit(generalized_sea): - res = generalized_sea.changeGeometry(1, combined) - feedback.pushInfo(f"result: {res}") - - # the difference algorithm changes the geometry type to - # multipolygon, convert back to single polygon - if QgsWkbTypes.isSingleType(input_layer.wkbType()): - generalized_sea = _multiparttosingleparts(generalized_sea) - - for gen_feature in generalized_sea.getFeatures(): - if feedback.isCanceled(): - break - - geom = gen_feature.geometry() - sink.addFeature( - copy_feature_with_geometry(gen_feature, geom), - QgsFeatureSink.Flag.FastInsert, - ) - - return {self.OUTPUT: dest_id} diff --git a/src/geogenalg/application/generalize_water_areas.py b/src/geogenalg/application/generalize_water_areas.py new file mode 100644 index 0000000..c4f7022 --- /dev/null +++ b/src/geogenalg/application/generalize_water_areas.py @@ -0,0 +1,186 @@ +# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) +# +# This file is part of geogen-algorithms. +# +# This source code is licensed under the MIT license found in the +# LICENSE file in the root directory of this source tree. + +import operator +from dataclasses import dataclass +from typing import override + +from geopandas import GeoDataFrame, GeoSeries +from shapely import MultiPoint, MultiPolygon, Polygon, force_2d + +from geogenalg.application import BaseAlgorithm, supports_identity +from geogenalg.core.exceptions import GeometryTypeError +from geogenalg.core.geometry import ( + assign_nearest_z, + chaikin_smooth_keep_topology, + extract_interior_rings, + perforate_polygon_with_gdf_exteriors, +) +from geogenalg.exaggeration import ( + exaggerate_thin_polygons, + extract_narrow_polygon_parts, +) +from geogenalg.utility.validation import check_gdf_geometry_type + + +@supports_identity +@dataclass(frozen=True) +class GeneralizeWaterAreas(BaseAlgorithm): + """Generalizes polygonal water areas. + + The algorithm does the following: + - Removes interior rings (islands) under given size + - Removes water areas under given size + - Exaggerates (buffer) thin sections of areas + - Exaggerates thin islands + - Simplifies the areas + - Smooths the simplified areas, while retaining topology between areas + + A shoreline reference data may optionally be entered. This affects the + algorithm so that any vertices not present in the shoreline data will not + be modified while smoothing. This makes sense to use with sea part + features, so that territorial water borders are not modified. + """ + + min_area: float = 4000.0 + """Features under this area will be removed.""" + area_simplification_tolerance: float = 10.0 + """Simplification tolerance used for water areas.""" + thin_section_width: float = 20.0 + """Sections under this width will be exaggerated.""" + thin_section_min_size: float = 200.0 + """Don't exaggerate thin sections under this size.""" + thin_section_exaggerate_by: float = 3.0 + """By how many CRS units thin sections will be exaggerated.""" + island_min_area: float = 100.0 + """Islands under this area will be removed.""" + island_min_width: float = 185.0 + """Islands under this width will be considered for exaggeration.""" + island_min_elongation: float = 0.25 + """Islands under this elongation will be considered for exaggeration.""" + island_exaggerate_by: float = 3.0 + """By how many CRS units thin islands will be exaggerated.""" + island_simplification_tolerance: float = 10.0 + """Simplification tolerance used for islands.""" + smoothing_passes: int = 3 + """How many smoothing passes will be performed. Each smoothing passes + (nearly) doubles the vertex count.""" + reference_key: str = "shoreline" + """Reference data key to use as a shoreline layer for preventing smoothing of + non-shoreline vertices.""" + + @override + def _execute( + self, + data: GeoDataFrame, + reference_data: dict[str, GeoDataFrame], + ) -> GeoDataFrame: + if not check_gdf_geometry_type(data, ["Polygon"]): + msg = "Input data must only contain Polygons." + raise GeometryTypeError(msg) + + if not 0.0 <= self.island_min_elongation <= 1.0: + msg = "Minimum island elongation must be between 0.0 and 1.0." + raise ValueError(msg) + + if self.reference_key in reference_data: + shoreline_gdf = reference_data[self.reference_key] + if not check_gdf_geometry_type(shoreline_gdf, ["LineString"]): + msg = "Reference data must only contain LineStrings." + raise GeometryTypeError(msg) + + # Extract out any points which are found in the input (Polygon) + # data, but not the shoreline LineString data. This allows to + # determine non-shoreline vertices (e.g. territorial water borders + # in case of sea areas) and skip smoothing them later. + data_points = data.extract_unique_points().union_all() + shoreline_points = shoreline_gdf.extract_unique_points().union_all() + skip_coords = force_2d(data_points.difference(shoreline_points)) + else: + skip_coords = MultiPoint() + + gdf = data.copy() + + thin_sections = ( + extract_narrow_polygon_parts(gdf, self.thin_section_width) + .explode(index_parts=True) + .geometry + ) + thin_sections = thin_sections.loc[ + thin_sections.geometry.area > self.thin_section_min_size + ].buffer(self.thin_section_exaggerate_by) + + def add_exaggerated_parts(geom: Polygon) -> Polygon: + intersecting_geoms = thin_sections.loc[ + thin_sections.geometry.intersects(geom) + ].union_all() + return geom.union(intersecting_geoms) + + gdf.geometry = gdf.geometry.apply(add_exaggerated_parts) + + # Extract from unary union to catch islands which are between lake parts + islands_geom = extract_interior_rings(gdf.union_all()) + islands = GeoDataFrame( + geometry=GeoSeries(islands_geom).explode(index_parts=True) + ) + + gdf.geometry = gdf.geometry.apply(lambda geom: Polygon(geom.exterior)) + + # When islands are extracted, we get just an exterior ring. Because of + # recursive lakes and islands, it is required to add any lakes as + # interior rings to the extracted islands. + islands.geometry = islands.geometry.apply( + lambda geom: perforate_polygon_with_gdf_exteriors(geom, gdf) + ) + + islands = exaggerate_thin_polygons( + islands, + self.island_min_width, + self.island_min_elongation, + self.island_exaggerate_by, + ) + + gdf.geometry = gdf.geometry.simplify_coverage( + self.area_simplification_tolerance, + ) + islands.geometry = islands.geometry.simplify( + self.island_simplification_tolerance, + ) + islands = islands.loc[islands.geometry.area > self.island_min_area] + + # Delete small areas + gdf = gdf.loc[gdf.geometry.area > self.min_area] + + # Add islands back in + gdf.geometry = gdf.geometry.difference(islands.geometry.union_all()) + + # The difference between the islands and areas might result in tiny + # slivers and therefore change the geometry to a MultiPolygon. Resolve + # this by transforming any MultiPolygons to Polygons, retaining only + # the largest part. + def multipolygon_to_single_keep_largest( + geom: MultiPolygon | Polygon, + ) -> Polygon: + if isinstance(geom, Polygon): + return geom + + areas = [(geom.area, geom) for geom in geom.geoms] + areas.sort(key=operator.itemgetter(0), reverse=True) + + return areas[0][1] + + mask = gdf.geometry.geom_type == "MultiPolygon" + gdf.loc[mask, gdf.geometry.name] = gdf.loc[mask].geometry.apply( + multipolygon_to_single_keep_largest + ) + + gdf.geometry = chaikin_smooth_keep_topology( + gdf.geometry, + extra_skip_coords=skip_coords, + ) + + return assign_nearest_z(data, gdf) diff --git a/src/geogenalg/application/watercourse_areas.py b/src/geogenalg/application/watercourse_areas.py index e1c889d..c52857c 100644 --- a/src/geogenalg/application/watercourse_areas.py +++ b/src/geogenalg/application/watercourse_areas.py @@ -28,10 +28,10 @@ elongation, explode_line, extend_line_to_nearest, - extract_interior_rings, + extract_interior_rings_gdf, lines_to_segments, move_to_point, - rectangle_dimensions, + oriented_envelope_dimensions, scale_line_to_length, ) @@ -79,7 +79,7 @@ def generalize_watercourse_areas( # noqa: C901, PLR0913, PLR0915 def exaggerate_thin_island(island: Polygon) -> Polygon: try: - width = rectangle_dimensions(island.oriented_envelope).width + width = oriented_envelope_dimensions(island.oriented_envelope).width if elongation(island) >= min_island_elongation and width < min_island_width: return island.buffer(distance=exaggerate_island_by, quadsegs=2) # noqa: SC200 except GeometryOperationError: @@ -95,7 +95,7 @@ def exaggerate_thin_island(island: Polygon) -> Polygon: watercourses.geometry = watercourses.geometry.apply( lambda x: remove_inner_rings(x, 0, crs=watercourses.crs) ) - islands = extract_interior_rings(areas) + islands = extract_interior_rings_gdf(areas) islands = islands.loc[islands.geometry.area >= min_island_area] islands.geometry = islands.geometry.apply(exaggerate_thin_island) watercourses.geometry = watercourses.geometry.difference(islands.unary_union) diff --git a/src/geogenalg/buffer/__init__.py b/src/geogenalg/buffer/__init__.py deleted file mode 100644 index 161bae7..0000000 --- a/src/geogenalg/buffer/__init__.py +++ /dev/null @@ -1,6 +0,0 @@ -# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) -# -# This file is part of geogen-algorithms. -# -# This source code is licensed under the MIT license found in the -# LICENSE file in the root directory of this source tree. diff --git a/src/geogenalg/buffer/exaggerate_thin_parts.py b/src/geogenalg/buffer/exaggerate_thin_parts.py deleted file mode 100644 index 9bd9041..0000000 --- a/src/geogenalg/buffer/exaggerate_thin_parts.py +++ /dev/null @@ -1,213 +0,0 @@ -# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) -# -# This file is part of geogen-algorithms. -# -# This source code is licensed under the MIT license found in the -# LICENSE file in the root directory of this source tree. - - - -from typing import Any - -from qgis.core import ( - QgsFeature, - QgsFeatureSink, - QgsFeatureSource, - QgsGeometry, - QgsProcessing, - QgsProcessingAlgorithm, - QgsProcessingContext, - QgsProcessingFeedback, - QgsProcessingParameterBoolean, - QgsProcessingParameterFeatureSink, - QgsProcessingParameterFeatureSource, - QgsProcessingParameterNumber, -) -from qgis.PyQt.QtCore import QCoreApplication -from yleistys_qgis_plugin.core.exaggeration import buffer_elimination, filter_and_buffer -from yleistys_qgis_plugin.core.utils import copy_feature_with_geometry - - -class ExaggerateThinParts(QgsProcessingAlgorithm): - INPUT = "INPUT" - OUTPUT = "OUTPUT" - WIDTH_TOLERANCE = "WIDTH_TOLERANCE" - MAX_THINNESS = "MAX_THINNESS" - MIN_AREA = "MIN_AREA" - EXAGGERATE_BY = "EXAGGERATE_BY" - COMBINE = "COMBINE" - - def __init__(self) -> None: - super().__init__() - - self._name = "exaggeratethinparts" - self._display_name = "Exaggerate thin parts" - self._group_id = "utility" - self._group = "Utility" - self._short_help_string = "" - - def tr(self, string) -> str: - """Returns a translatable string with the self.tr() function.""" - return QCoreApplication.translate("Processing", string) - - def createInstance(self): # noqa: N802 - return ExaggerateThinParts() - - def name(self) -> str: - """Returns the algorithm name, used for identifying the algorithm. This - string should be fixed for the algorithm, and must not be localised. - The name should be unique within each provider. Names should contain - lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._name - - def displayName(self) -> str: # noqa: N802 - """Returns the translated algorithm name, which should be used for any - user-visible display of the algorithm name. - """ - return self.tr(self._display_name) - - def groupId(self) -> str: # noqa: N802 - """Returns the unique ID of the group this algorithm belongs to. This - string should be fixed for the algorithm, and must not be localised. - The group id should be unique within each provider. Group id should - contain lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._group_id - - def group(self) -> str: - """Returns the name of the group this algorithm belongs to. This string - should be localised. - """ - return self.tr(self._group) - - def shortHelpString(self) -> str: # noqa: N802 - """Returns a localised short helper string for the algorithm. This string - should provide a basic description about what the algorithm does and the - parameters and outputs associated with it.. - """ - return self.tr(self._short_help_string) - - def initAlgorithm(self, configuration=None) -> None: # noqa: N802 - """Here we define the inputs and output of the algorithm, along - with some other properties. - """ - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT, - self.tr("Input layer"), - [QgsProcessing.SourceType.TypeVectorPolygon], - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.WIDTH_TOLERANCE, - self.tr("Width tolerance"), - defaultValue=10, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.EXAGGERATE_BY, - self.tr("Exaggerate by"), - defaultValue=10, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MAX_THINNESS, - self.tr("Maximum thinness"), - defaultValue=0.5, - type=QgsProcessingParameterNumber.Type.Double, - minValue=0.0, - maxValue=1.0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MIN_AREA, - self.tr("Minimum area"), - defaultValue=200, - type=QgsProcessingParameterNumber.Type.Double, - ) - ) - - self.addParameter( - QgsProcessingParameterBoolean( - self.COMBINE, - self.tr("Combine into geometries"), - defaultValue=True, - ) - ) - - self.addParameter( - QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr("Exaggerated")) - ) - - def processAlgorithm( # noqa: N802 - self, - parameters: dict[str, Any], - context: QgsProcessingContext, - feedback: QgsProcessingFeedback, - ) -> dict: - """Here is where the processing itself takes place.""" - if feedback is None: - feedback = QgsProcessingFeedback() - - source: QgsFeatureSource = self.parameterAsSource( - parameters, self.INPUT, context - ) - width_tolerance = self.parameterAsInt(parameters, self.WIDTH_TOLERANCE, context) - exaggerate_by = self.parameterAsInt(parameters, self.EXAGGERATE_BY, context) - max_thinness = self.parameterAsInt(parameters, self.MAX_THINNESS, context) - min_area = self.parameterAsInt(parameters, self.MIN_AREA, context) - combine = self.parameterAsBool(parameters, self.COMBINE, context) - - sink, dest_id = self.parameterAsSink( - parameters, - self.OUTPUT, - context, - source.fields(), - source.wkbType(), - source.sourceCrs(), - ) - - total = 100.0 / source.featureCount() if source.featureCount() else 0 - - feature: QgsFeature - for current, feature in enumerate(source.getFeatures()): - geom = feature.geometry() - - eliminated_buffer = buffer_elimination(geom, width_tolerance / 2) - diff = geom.difference(eliminated_buffer) - - exaggerated_areas: QgsGeometry | None = filter_and_buffer( - diff, exaggerate_by, min_area, max_thinness - ) - - feedback.setProgress(int(current * total)) - - if combine: - final_geom = ( - geom.combine(exaggerated_areas) if exaggerated_areas else geom - ) - else: - final_geom = exaggerated_areas or None - - if not final_geom: - continue - - sink.addFeature( - copy_feature_with_geometry(feature, final_geom), - QgsFeatureSink.Flag.FastInsert, - ) - - return {self.OUTPUT: dest_id} diff --git a/src/geogenalg/buffer/generalize_islands.py b/src/geogenalg/buffer/generalize_islands.py deleted file mode 100644 index 72d2b1f..0000000 --- a/src/geogenalg/buffer/generalize_islands.py +++ /dev/null @@ -1,226 +0,0 @@ -# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) -# -# This file is part of geogen-algorithms. -# -# This source code is licensed under the MIT license found in the -# LICENSE file in the root directory of this source tree. - - - -from typing import Any - -from qgis.core import ( - Qgis, - QgsFeature, - QgsFeatureSink, - QgsFeatureSource, - QgsMapToPixelSimplifier, - QgsProcessing, - QgsProcessingAlgorithm, - QgsProcessingContext, - QgsProcessingFeedback, - QgsProcessingParameterFeatureSink, - QgsProcessingParameterFeatureSource, - QgsProcessingParameterNumber, -) -from qgis.PyQt.QtCore import QCoreApplication -from yleistys_qgis_plugin.core.utils import ( - VISVALINGAM_ENUM_VERSION_CUTOFF, - copy_feature_with_geometry, -) - - -class GeneralizeIslands(QgsProcessingAlgorithm): - INPUT = "INPUT" - OUTPUT = "OUTPUT" - MIN_WIDTH = "MIN_WIDTH" - MIN_ELONGATION = "MIN_ELONGATION" - MAX_SIMPLIFICATION_TOLERANCE = "MAX_SIMPLIFICATION_TOLERANCE" - EXAGGERATE_BY = "EXAGGERATE_BY" - - def __init__(self) -> None: - super().__init__() - - self._name = "generalizeislands" - self._display_name = "Generalize islands" - self._group_id = "utility" - self._group = "Utility" - self._short_help_string = "" - - def tr(self, string) -> str: - """Returns a translatable string with the self.tr() function.""" - return QCoreApplication.translate("Processing", string) - - def createInstance(self): # noqa: N802 - return GeneralizeIslands() - - def name(self) -> str: - """Returns the algorithm name, used for identifying the algorithm. This - string should be fixed for the algorithm, and must not be localised. - The name should be unique within each provider. Names should contain - lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._name - - def displayName(self) -> str: # noqa: N802 - """Returns the translated algorithm name, which should be used for any - user-visible display of the algorithm name. - """ - return self.tr(self._display_name) - - def groupId(self) -> str: # noqa: N802 - """Returns the unique ID of the group this algorithm belongs to. This - string should be fixed for the algorithm, and must not be localised. - The group id should be unique within each provider. Group id should - contain lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._group_id - - def group(self) -> str: - """Returns the name of the group this algorithm belongs to. This string - should be localised. - """ - return self.tr(self._group) - - def shortHelpString(self) -> str: # noqa: N802 - """Returns a localised short helper string for the algorithm. This string - should provide a basic description about what the algorithm does and the - parameters and outputs associated with it.. - """ - return self.tr(self._short_help_string) - - def initAlgorithm(self, configuration=None) -> None: # noqa: N802 - """Here we define the inputs and output of the algorithm, along - with some other properties. - """ - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT, - self.tr("Input layer"), - [QgsProcessing.SourceType.TypeVectorPolygon], - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MIN_WIDTH, - self.tr("Exaggerate only if max width under"), - defaultValue=185, - type=QgsProcessingParameterNumber.Type.Double, - minValue=0.0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MIN_ELONGATION, - self.tr("Exaggerate only if elongation over"), - defaultValue=3.349, - type=QgsProcessingParameterNumber.Type.Double, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.EXAGGERATE_BY, - self.tr("Exaggerate by"), - defaultValue=3, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterNumber( - self.MAX_SIMPLIFICATION_TOLERANCE, - self.tr("Maximum simplification tolerance"), - defaultValue=10, - minValue=0, - ) - ) - - self.addParameter( - QgsProcessingParameterFeatureSink( - self.OUTPUT, self.tr("Generalized islands") - ) - ) - - def processAlgorithm( # noqa: N802 - self, - parameters: dict[str, Any], - context: QgsProcessingContext, - feedback: QgsProcessingFeedback, - ) -> dict: - """Here is where the processing itself takes place.""" - if feedback is None: - feedback = QgsProcessingFeedback() - - source: QgsFeatureSource = self.parameterAsSource( - parameters, self.INPUT, context - ) - min_width: float = self.parameterAsDouble(parameters, self.MIN_WIDTH, context) - min_elongation: float = self.parameterAsDouble( - parameters, self.MIN_ELONGATION, context - ) - max_simplification_tolerance: float = self.parameterAsDouble( - parameters, self.MAX_SIMPLIFICATION_TOLERANCE, context - ) - exaggerate_by: float = self.parameterAsDouble( - parameters, self.EXAGGERATE_BY, context - ) - - sink, dest_id = self.parameterAsSink( - parameters, - self.OUTPUT, - context, - source.fields(), - source.wkbType(), - source.sourceCrs(), - ) - - total = 100.0 / source.featureCount() if source.featureCount() else 0 - - if Qgis.QGIS_VERSION_INT < VISVALINGAM_ENUM_VERSION_CUTOFF: # pyright: ignore - simplifier = QgsMapToPixelSimplifier( - QgsMapToPixelSimplifier.SimplifyFlag.SimplifyGeometry, - 1, - QgsMapToPixelSimplifier.SimplifyAlgorithm.Visvalingam, - ) - else: - simplifier = QgsMapToPixelSimplifier( - QgsMapToPixelSimplifier.SimplifyFlag.SimplifyGeometry, - 1, - Qgis.VectorSimplificationAlgorithm.Visvalingam, # pyright: ignore - ) - - feature: QgsFeature - for current, feature in enumerate(source.getFeatures()): - if feedback.isCanceled(): - break - - geom = feature.geometry() - - _, _, _, bbox_width, bbox_height = geom.orientedMinimumBoundingBox() - elongation: float = bbox_height / bbox_width - - # check whether the island should be exaggerated - if elongation >= min_elongation and bbox_width < min_width: - geom = geom.buffer(distance=exaggerate_by, segments=2) - - vertex_count: int = geom.constGet().nCoordinates() - add_to_tolerance: float = vertex_count / 200 - add_to_tolerance = min(add_to_tolerance, max_simplification_tolerance - 1) - - simplifier.setTolerance(1 + add_to_tolerance) - - simplified_geom = simplifier.simplify(geom) - smoothed_geom = simplified_geom.smooth() - sink.addFeature( - copy_feature_with_geometry(feature, smoothed_geom), - QgsFeatureSink.Flag.FastInsert, - ) - - feedback.setProgress(current * total) - - return {self.OUTPUT: dest_id} diff --git a/src/geogenalg/core/extract_interior_rings.py b/src/geogenalg/core/extract_interior_rings.py deleted file mode 100644 index 6c7d328..0000000 --- a/src/geogenalg/core/extract_interior_rings.py +++ /dev/null @@ -1,175 +0,0 @@ -# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) -# -# This file is part of geogen-algorithms. -# -# This source code is licensed under the MIT license found in the -# LICENSE file in the root directory of this source tree. - - - -from typing import Any - -from qgis.core import ( - Qgis, - QgsFeature, - QgsFeatureRequest, - QgsFeatureSink, - QgsFeatureSource, - QgsGeometry, - QgsPointXY, - QgsProcessing, - QgsProcessingAlgorithm, - QgsProcessingContext, - QgsProcessingFeedback, - QgsProcessingParameterFeatureSink, - QgsProcessingParameterFeatureSource, - QgsRectangle, - QgsWkbTypes, -) -from qgis.PyQt.QtCore import QCoreApplication -from geogenalg.core.utils import ( - extract_interior_rings, -) - - -class ExtractInteriorRings(QgsProcessingAlgorithm): - INPUT = "INPUT" - OUTPUT = "OUTPUT" - - def __init__(self) -> None: - super().__init__() - - self._name = "extractinteriorrings" - self._display_name = "Extract interior rings (holes) from layer" - self._group_id = "utility" - self._group = "Utility" - self._short_help_string = "" - - def tr(self, string) -> str: - """Returns a translatable string with the self.tr() function.""" - return QCoreApplication.translate("Processing", string) - - def createInstance(self): # noqa: N802 - return ExtractInteriorRings() - - def name(self) -> str: - """Returns the algorithm name, used for identifying the algorithm. This - string should be fixed for the algorithm, and must not be localised. - The name should be unique within each provider. Names should contain - lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._name - - def displayName(self) -> str: # noqa: N802 - """Returns the translated algorithm name, which should be used for any - user-visible display of the algorithm name. - """ - return self.tr(self._display_name) - - def groupId(self) -> str: # noqa: N802 - """Returns the unique ID of the group this algorithm belongs to. This - string should be fixed for the algorithm, and must not be localised. - The group id should be unique within each provider. Group id should - contain lowercase alphanumeric characters only and no spaces or other - formatting characters. - """ - return self._group_id - - def group(self) -> str: - """Returns the name of the group this algorithm belongs to. This string - should be localised. - """ - return self.tr(self._group) - - def shortHelpString(self) -> str: # noqa: N802 - """Returns a localised short helper string for the algorithm. This string - should provide a basic description about what the algorithm does and the - parameters and outputs associated with it.. - """ - return self.tr(self._short_help_string) - - def initAlgorithm(self, configuration=None) -> None: # noqa: N802 - """Here we define the inputs and output of the algorithm, along - with some other properties. - """ - self.addParameter( - QgsProcessingParameterFeatureSource( - self.INPUT, - self.tr("Input layer"), - [QgsProcessing.SourceType.TypeVectorPolygon], - ) - ) - - self.addParameter( - QgsProcessingParameterFeatureSink(self.OUTPUT, self.tr("Extracted")) - ) - - def processAlgorithm( # noqa: N802 - self, - parameters: dict[str, Any], - context: QgsProcessingContext, - feedback: QgsProcessingFeedback, - ) -> dict: - """Here is where the processing itself takes place.""" - if feedback is None: - feedback = QgsProcessingFeedback() - - source: QgsFeatureSource = self.parameterAsSource( - parameters, self.INPUT, context - ) - - sink, dest_id = self.parameterAsSink( - parameters, - self.OUTPUT, - context, - source.fields(), - source.wkbType(), - source.sourceCrs(), - ) - - geom_rings: list[QgsGeometry] = [] - - for feature in source.getFeatures(): - if feedback.isCanceled(): - return {} - geom_rings.append(extract_interior_rings(feature.geometry())) - - for geom_ring in geom_rings: - if feedback.isCanceled(): - return {} - for primitive_ring in geom_ring.constParts(): - if feedback.isCanceled(): - return {} - - interior_ring: QgsGeometry = QgsGeometry(primitive_ring.clone()) - - # let's check if there are any other exterior rings inside the interior - # ring and add those as rings to handle recursive islands - rect: QgsRectangle = interior_ring.boundingBox() - - request = QgsFeatureRequest().setFilterRect(rect) - - for feat in source.getFeatures(request): - geom = feat.geometry() - if geom.within(interior_ring): - rings: list[list[QgsPointXY]] - - if QgsWkbTypes.isSingleType(geom.wkbType()): - rings = [geom.asPolygon()[0]] - else: - rings = [ring[0] for ring in geom.asMultiPolygon()] - - for ring in rings: - res = interior_ring.addRing(ring) - if res != Qgis.GeometryOperationResult.Success: - feedback.reportError( - "Could not add ring to extracted interior ring!", - fatalError=True, - ) - - feature = QgsFeature(source.fields()) - feature.setGeometry(interior_ring) - sink.addFeature(feature, QgsFeatureSink.Flag.FastInsert) - - return {self.OUTPUT: dest_id} diff --git a/src/geogenalg/core/geometry.py b/src/geogenalg/core/geometry.py index ae2057b..cdb2d90 100644 --- a/src/geogenalg/core/geometry.py +++ b/src/geogenalg/core/geometry.py @@ -6,7 +6,6 @@ # LICENSE file in the root directory of this source tree. -import math from enum import Enum from statistics import mean from typing import NamedTuple @@ -22,7 +21,7 @@ Point, Polygon, affinity, - count_coordinates, + force_2d, get_coordinates, polygonize, shortest_line, @@ -39,17 +38,255 @@ ) -class LineExtendFrom(Enum): # noqa: D101 +class LineExtendFrom(Enum): + """An enumeration describing from which end a linestring should be extended.""" + END = -1 START = 0 BOTH = 1 -class Dimensions(NamedTuple): # noqa: D101 +class Dimensions(NamedTuple): + """Simple named tuple describing the dimensions of a rectangle.""" + width: float height: float +def chaikin_smooth_skip_coords( + geom: LineString | Polygon, + skip_coords: list[Point] | MultiPoint, + *, + iterations: int = 1, +) -> LineString | Polygon: + """Smooth input linestring or polygon. + + This is an implementation of Chaikin's corner cutting line smoothing + algorithm, with the addition of being able to skip the smoothing of + specific points. + + Args: + ---- + geom: Geometry to smooth. + iterations: Number of iterations i.e. smoothing passes. Each pass + doubles the number of vertices (minus the skipped coordinates) i.e. + increase with caution, growth is exponential. + skip_coords: List of points to skip. If a corresponding coordinate exists + in the input geometry, it is guaranteed to remain unchanged. + + Returns: + ------- + Smoothed geometry. + + Note: + ---- + If skipping of coordinates is not required, prefer shapelysmooth's + function as it is likely to be more efficient. + + + """ + if isinstance(skip_coords, MultiPoint): + skip_coords_ = [point.coords[0] for point in skip_coords.geoms] + else: + skip_coords_ = [point.coords[0] for point in skip_coords] + + # TODO: allow processing 2.5D geometries and multigeometries? + + def _process_coord_sequence( + seq: CoordinateSequence, + output: list[tuple[float, ...]], + ) -> list[tuple[float, ...]]: + for i in range(len(seq) - 1): + current_coord = seq[i] + next_coord = seq[i + 1] + + current_x = current_coord[0] + current_y = current_coord[1] + next_x = next_coord[0] + next_y = next_coord[1] + + q = ((0.75 * current_x + 0.25 * next_x), (0.75 * current_y + 0.25 * next_y)) + r = ((0.25 * current_x + 0.75 * next_x), (0.25 * current_y + 0.75 * next_y)) + + if current_coord in skip_coords_: + output.append(current_coord) + if next_coord not in skip_coords_: + output.append(r) + continue + + output.append(q) + + if next_coord not in skip_coords_: + output.append(r) + + return output + + def _process_linestring(geom: LineString) -> LineString: + coords = [geom.coords[0]] + + processed_coords = _process_coord_sequence(geom.coords, coords) + processed_coords.append(geom.coords[-1]) + + return LineString(processed_coords) + + def _process_polygon(geom: Polygon) -> Polygon: + exterior = _process_coord_sequence(geom.exterior.coords, []) + interiors = [ + _process_coord_sequence(interior.coords, []) for interior in geom.interiors + ] + + return Polygon(exterior, interiors) + + process_function = ( + _process_linestring if isinstance(geom, LineString) else _process_polygon + ) + + result = geom + for _ in range(iterations): + result = process_function(result) + + return result + + +def get_topological_points(geoseries: GeoSeries) -> list[Point]: + """Find all topological points in a GeoSeries. + + Topological point referring to a point which is shared by two or more + geometries in the GeoSeries. + + Args: + ---- + geoseries: The GeoSeries to find topological points in. + + Returns: + ------- + List of all the topological points (if any). + + Raises: + ------ + GeometryOperationError: If union could not be performed on unique points + in the GeoSeries. + + """ + if geoseries.empty: + return [] + + unique_points = geoseries.force_2d().extract_unique_points().union_all() + + if isinstance(unique_points, Point): + unique_points = MultiPoint([unique_points]) + + if not isinstance(unique_points, MultiPoint): + msg = "Unique points not a Point or a MultiPoint" + raise GeometryOperationError(msg) + + # TODO: maybe there's a faster way to do this + topo_points: list[Point] = [] + for point in unique_points.geoms: + intersections = geoseries.intersects(point) + + # If there are more than 1 intersections + if len(intersections.loc[intersections]) > 1: + topo_points.append(point) + + return topo_points + + +def chaikin_smooth_keep_topology( + geoseries: GeoSeries, + iterations: int = 3, + *, + extra_skip_coords: list[Point] | MultiPoint | None, +) -> GeoSeries: + """Apply smoothing algorithm while keeping topological points unchanged. + + Args: + ---- + geoseries: GeoSeries to be smoothed. + iterations: Number of smoothing passes, increase for a smoother result. + Note that each pass roughly doubles the number of vertices, so + growth will be exponential. Anything above 6-7 is unlikely to + produce cartographically meaningful differences and therefore + unnecessarily increase vertex count. + extra_skip_coords: Any additional coordinates in addition to + topological points which should not be smoothed. + + Returns: + ------- + GeoSeries with smoothed geometries, with shared topological points + unchanged. + + Note: + ---- + If the input has 3D geometries, they will be changed to 2D geometries. + + """ + copy = geoseries.copy() + + skipped_coords = get_topological_points(copy) + + if isinstance(extra_skip_coords, MultiPoint): + extra_skip_coords = list(extra_skip_coords.geoms) + + if extra_skip_coords is not None: + skipped_coords += extra_skip_coords + + return copy.apply( + lambda geom: chaikin_smooth_skip_coords( + force_2d(geom), + iterations=iterations, + skip_coords=skipped_coords, + ), + ) + + +def perforate_polygon_with_gdf_exteriors( + geometry: Polygon, + gdf: GeoDataFrame, +) -> Polygon: + """Add the exteriors of a polygon GeoDataFrame as holes to a Polygon. + + The purpose of this function over a simple difference operation is to + handle cases where there are polygon(s) inside interior rings of the input + geometry and you need to add holes only for the exteriors. This is useful + f.e. for recursive islands and lakes. With difference, the result would be + a multipolygon with all the contained geometries as parts. + + Args: + ---- + geometry: A Polygon to which holes will be added + gdf: A GeoDataFrame whose geometry exteriors will be added as holes + + Returns: + ------- + A new Polygon with the holes added. + + """ + features_within_geom = gdf.geometry.loc[gdf.geometry.within(geometry)] + exteriors = features_within_geom.apply(lambda geom: Polygon(geom.exterior)) + + return geometry.difference(exteriors.union_all()) + + +def extract_interior_rings(geometry: Polygon | MultiPolygon) -> MultiPolygon: + """Extract interior rings (holes) of a geometry as a new geometry. + + Returns + ------- + Interior rings as a MultiPolygon. + + """ + if isinstance(geometry, Polygon): + polygons = [Polygon(interior) for interior in geometry.interiors] + else: + polygons = [] + for geom in geometry.geoms: + for interior in geom.interiors: + polygons.append(Polygon(interior)) + + return MultiPolygon(polygons) + + def mean_z(geom: BaseGeometry, nodata_value: float = -999.999) -> float: """Calculate the mean of z values in a geometry's vertices. @@ -93,31 +330,25 @@ def mean_z(geom: BaseGeometry, nodata_value: float = -999.999) -> float: return mean(z_values) -def rectangle_dimensions(rect: Polygon) -> Dimensions: - """Calculate the dimensions of a rectangular polygon. +def oriented_envelope_dimensions(geom: Polygon) -> Dimensions: + """Calculate the dimensions of a geometry from its oriented envelope. - Returns: - The width and height of a rectangular polygon. + Returns + ------- + The width and height of the oriented envelope. - Raises: - GeometryOperationError: If input polygon is invalid. - InvalidGeometryError: If input polygon is not a rectangle. + Raises + ------ + InvalidGeometryError: If input polygon is invalid. """ - required_coords = 5 - if count_coordinates(rect) != required_coords or not math.isclose( # noqa: SC200 - rect.area, - rect.oriented_envelope.area, - rel_tol=1e-06, # noqa: SC200 - ): - msg = f"Not a rectangle: {rect}" - raise GeometryOperationError(msg) - - if not rect.is_valid: + if not geom.is_valid: msg = "Rectangle not valid" raise InvalidGeometryError(msg) - x, y = rect.exterior.coords.xy + envelope = geom.oriented_envelope + + x, y = envelope.exterior.coords.xy point_1 = Point(x[0], y[0]) point_2 = Point(x[1], y[1]) @@ -125,36 +356,31 @@ def rectangle_dimensions(rect: Polygon) -> Dimensions: lengths = (point_1.distance(point_2), point_2.distance(point_3)) - return Dimensions(min(lengths), max(lengths)) + return Dimensions(width=min(lengths), height=max(lengths)) def elongation(polygon: Polygon) -> float: """Calculate the elongation of a polygon. - Returns: + Returns + ------- The ratio of the height (longest side) to the width (shorter side) of its minimum bounding oriented envelope. - Raises: - TypeError: if created envelope is not of type polygon. - """ - envelope: BaseGeometry = polygon.oriented_envelope - if not isinstance(envelope, Polygon): - msg = "Envelope is not a Polygon" - raise TypeError(msg) - - dimensions = rectangle_dimensions(envelope) - return dimensions.height / dimensions.width + dimensions = oriented_envelope_dimensions(polygon) + return dimensions.width / dimensions.height -def extract_interior_rings(areas: GeoDataFrame) -> GeoDataFrame: +def extract_interior_rings_gdf(areas: GeoDataFrame) -> GeoDataFrame: """Extract the interior rings of a polygon geodataframe. - Returns: + Returns + ------- A new geodataframe containing the interior rings. - Raises: + Raises + ------ GeometryTypeError: If geometries are not polygons or multipolygons. """ @@ -181,7 +407,8 @@ def extract_interior_rings(areas: GeoDataFrame) -> GeoDataFrame: def explode_line(line: LineString) -> list[LineString]: """Explode a line geometry to all its segments. - Returns: + Returns + ------- All segments of a line geometry. """ @@ -193,11 +420,13 @@ def explode_line(line: LineString) -> list[LineString]: def lines_to_segments(lines: GeoSeries) -> GeoSeries: """Convert lines to segments. - Returns: + Returns + ------- A GeoSeries of line geometries containing all line segments from the input GeoSeries. - Raises: + Raises + ------ GeometryTypeError: if param `lines` contains other geometry types than lines. @@ -218,11 +447,13 @@ def scale_line_to_length( ) -> LineString: """Scale line to given length. - Returns: + Returns + ------- A version of the input geometry which has been scaled to the given length. - Raises: + Raises + ------ ValueError: If param `geom` has no length. ValueError: If param `length` is less than or equal to 0. @@ -246,7 +477,8 @@ def scale_line_to_length( def move_to_point(geom: BaseGeometry, point: Point) -> BaseGeometry: """Move geometry to given point position. - Returns: + Returns + ------- A version of the input geometry which has been moved to the given point position. @@ -265,12 +497,14 @@ def extend_line_to_nearest( ) -> LineString: """Extend line to nearest point. - Returns: + Returns + ------- A version of the input line which has been extended to the nearest point of another geometry. The extension can be done from the end or start of the line, or alternatively both ends. - Raises: + Raises + ------ GeometryOperationError: If invalid result got from line union or merge operations. @@ -320,12 +554,14 @@ def extend_line_to_nearest( def point_on_line(line: LineString, distance: float) -> Point: """Get point along line. - Returns: + Returns + ------- A point which is the given distance away and collinear to either the first or last segment of the input line. If distance is > 0 the point is calculated from the last segment and if < 0 the first segment. - Raises: + Raises + ------ ValueError: If line does not have exactly 2 vertices. """ @@ -371,11 +607,13 @@ def extend_line_by( ) -> LineString: """Extend line by given distance. - Returns: + Returns + ------- Version of the input line which has been extended by the given distance either from its start or end segment, or both. - Raises: + Raises + ------ ValueError: if extend_from is less or equal than 0 """ diff --git a/src/geogenalg/exaggeration.py b/src/geogenalg/exaggeration.py index 3adce4a..4dc9207 100644 --- a/src/geogenalg/exaggeration.py +++ b/src/geogenalg/exaggeration.py @@ -6,15 +6,54 @@ # LICENSE file in the root directory of this source tree. -import geopandas as gpd +from dataclasses import dataclass +from typing import Literal + +from geopandas import GeoDataFrame +from shapely import BufferCapStyle, BufferJoinStyle, Polygon +from shapely.geometry.base import BaseGeometry from geogenalg.core.exceptions import GeometryTypeError +from geogenalg.core.geometry import elongation, oriented_envelope_dimensions from geogenalg.utility.validation import check_gdf_geometry_type +@dataclass +class BufferOptions: + """Data class for passing options to buffer function.""" + + quad_segments: int = 16 + cap_style: BufferCapStyle | Literal["round", "square", "flat"] = "round" + join_style: BufferJoinStyle | Literal["round", "mitre", "bevel"] = "round" + mitre_limit: float = 5 + single_sided: bool = False + + +def buffer_with_options( + geom: BaseGeometry, + distance: float, + options: BufferOptions, +) -> Polygon: + """Buffer geometry with options from BufferOptions object. + + Returns + ------- + Buffered geometry. + + """ + return geom.buffer( + distance=distance, + quad_segs=options.quad_segments, # noqa: SC200 + cap_style=options.cap_style, + join_style=options.join_style, + mitre_limit=options.mitre_limit, + single_sided=options.single_sided, + ) + + def extract_narrow_polygon_parts( - input_gdf: gpd.GeoDataFrame, threshold: float -) -> gpd.GeoDataFrame: + input_gdf: GeoDataFrame, threshold: float +) -> GeoDataFrame: """Extract polygon parts narrower than the threshold. Args: @@ -61,3 +100,47 @@ def extract_narrow_polygon_parts( ).buffer(0.1, cap_style="flat", join_style="mitre") return narrow_parts_gdf + + +def exaggerate_thin_polygons( + input_gdf: GeoDataFrame, + width_threshold: float, + elongation_threshold: float, + exaggerate_by: float, + buffer_options: BufferOptions | None = None, +) -> GeoDataFrame: + """Exaggerate polygon if it's considered thin. + + If the polygon's oriented bounding box is under the width and elongation + threshold, it will be considered thin. + + Args: + ---- + input_gdf: GeoDataFrame containing polygon geometries. + width_threshold: maximum value for width for polygon to be considered thin + elongation_threshold: maximum value for elongation for polygon to be + considered thin + exaggerate_by: buffer distance for exaggeration + buffer_options: options to pass to buffer function + + Returns: + ------- + GeoDataFrame containing the input features with the thin polygons exaggerated. + + """ + if buffer_options is None: + buffer_options = BufferOptions() + + modified = input_gdf.copy() + + def _process_polygon(geom: Polygon) -> Polygon: + width = oriented_envelope_dimensions(geom).width + + if elongation(geom) <= elongation_threshold and width <= width_threshold: + return buffer_with_options(geom, exaggerate_by, buffer_options) + + return geom + + modified.geometry = modified.geometry.apply(_process_polygon) + + return modified diff --git a/src/geogenalg/main.py b/src/geogenalg/main.py index ded6fed..b3c651b 100644 --- a/src/geogenalg/main.py +++ b/src/geogenalg/main.py @@ -25,6 +25,7 @@ ) from geogenalg.application.generalize_fences import GeneralizeFences from geogenalg.application.generalize_landcover import GeneralizeLandcover +from geogenalg.application.generalize_water_areas import GeneralizeWaterAreas from geogenalg.utility.dataframe_processing import read_gdf_from_file_and_set_index GEOPACKAGE_URI_HELP = ( @@ -274,6 +275,7 @@ def build_app() -> None: "clusters_to_centroids": GeneralizePointClustersAndPolygonsToCentroids, "fences": GeneralizeFences, "landcover": GeneralizeLandcover, + "water_areas": GeneralizeWaterAreas, } for cli_command_name, alg in commands_and_algs.items(): diff --git a/src/geogenalg/testing.py b/src/geogenalg/testing.py new file mode 100644 index 0000000..264e3cf --- /dev/null +++ b/src/geogenalg/testing.py @@ -0,0 +1,72 @@ +# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) +# +# This file is part of geogen-algorithms. +# +# This source code is licensed under the MIT license found in the +# LICENSE file in the root directory of this source tree. + + +from logging import getLogger +from pathlib import Path +from tempfile import gettempdir + +from geopandas import GeoDataFrame, overlay +from geopandas.testing import assert_geodataframe_equal + +logger = getLogger(__name__) + + +def assert_gdf_equal_save_diff( + result: GeoDataFrame, + control: GeoDataFrame, + *, + directory: Path | None = None, + save_control: bool = True, +) -> None: + """Assert GeoDataFrames are equal and if not, save results. + + Intended mainly for convenience in debugging. In tests, the + geopandas.testing and pandas.testing modules should be used. + + Args: + ---- + result: Result GeoDataFrame. + control: Control GeoDataFrame. + directory: Specify directory in which results will be saved. If None, a + temporary directory will be created. + save_control: Whether to save control dataset to directory as well. + + Raises: + ------ + AssertionError: If GeoDataFrames are different. + + """ + if directory is None: + directory = Path(gettempdir()) / "geogenalg_tests" + + if not directory.exists(): + directory.mkdir() + + try: + assert_geodataframe_equal(result, control) + except AssertionError: + result_path = directory / "result.gpkg" + geometry_difference_path = directory / "geomdiff.gpkg" + directory / "compare.gpkg" + + result.to_file(result_path) + + geom_diff = overlay(result, control, how="difference") + if not geom_diff.union_all().is_empty: + geom_diff.to_file(geometry_difference_path) + + # TODO: support saving dataframe comparison to file as well + # something like (comparison = result.compare(control)) + + if save_control: + control_path = directory / "control.gpkg" + control.to_file(control_path) + + logger.info("Test results saved to %s.", directory) + + raise diff --git a/test/application/test_generalize_water_areas.py b/test/application/test_generalize_water_areas.py new file mode 100644 index 0000000..fcb77a8 --- /dev/null +++ b/test/application/test_generalize_water_areas.py @@ -0,0 +1,127 @@ +# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) +# +# This file is part of geogen-algorithms. +# +# This source code is licensed under the MIT license found in the +# LICENSE file in the root directory of this source tree. + +import re +from pathlib import Path +from tempfile import TemporaryDirectory + +import pytest +from geopandas import read_file +from pandas.testing import assert_frame_equal +from pygeoops import GeoDataFrame +from shapely import Point, box + +from geogenalg.application.generalize_water_areas import GeneralizeWaterAreas +from geogenalg.core.exceptions import GeometryTypeError +from geogenalg.utility.dataframe_processing import read_gdf_from_file_and_set_index + +UNIQUE_ID_COLUMN = "kmtk_id" + + +def test_generalize_lakes(testdata_path: Path): + input_path = testdata_path / "lakes.gpkg" + input_data = read_gdf_from_file_and_set_index( + input_path, + UNIQUE_ID_COLUMN, + layer="lake_part", + ) + + temp_dir = TemporaryDirectory() + output_path = Path(temp_dir.name) / "output.gpkg" + + algorithm = GeneralizeWaterAreas( + min_area=4000.0, + island_min_area=100.0, + area_simplification_tolerance=10.0, + thin_section_width=20.0, + thin_section_min_size=200.0, + thin_section_exaggerate_by=3.0, + island_min_width=185.0, + island_min_elongation=0.25, + island_exaggerate_by=3.0, + island_simplification_tolerance=10.0, + smoothing_passes=3, + ) + + control = read_gdf_from_file_and_set_index( + input_path, + UNIQUE_ID_COLUMN, + layer="control", + ) + algorithm.execute(input_data, {}).to_file(output_path, layer="result") + + result = read_gdf_from_file_and_set_index( + output_path, UNIQUE_ID_COLUMN, layer="result" + ) + + assert_frame_equal(result, control) + + +def test_generalize_sea(testdata_path: Path): + input_path = testdata_path / "sea.gpkg" + input_data = read_gdf_from_file_and_set_index( + input_path, + UNIQUE_ID_COLUMN, + layer="sea_part", + ) + + shoreline = read_file(input_path, layer="shoreline") + + temp_dir = TemporaryDirectory() + output_path = Path(temp_dir.name) / "output.gpkg" + + algorithm = GeneralizeWaterAreas( + min_area=4000.0, + area_simplification_tolerance=10.0, + thin_section_width=20.0, + thin_section_min_size=200.0, + thin_section_exaggerate_by=3.0, + island_min_area=100.0, + island_min_width=185.0, + island_min_elongation=0.25, + island_exaggerate_by=3.0, + island_simplification_tolerance=10.0, + smoothing_passes=3, + ) + + control = read_gdf_from_file_and_set_index( + input_path, + UNIQUE_ID_COLUMN, + layer="control", + ) + algorithm.execute(input_data, {"shoreline": shoreline}).to_file( + output_path, layer="result" + ) + + result = read_gdf_from_file_and_set_index( + output_path, UNIQUE_ID_COLUMN, layer="result" + ) + + assert_frame_equal(control, result) + + +def test_invalid_geom_type() -> None: + gdf = GeoDataFrame({"id": [1]}, geometry=[Point(0, 0)]) + + with pytest.raises( + GeometryTypeError, + match=re.escape("Input data must only contain Polygons."), + ): + GeneralizeWaterAreas().execute(data=gdf, reference_data={}) + + +def test_invalid_shoreline_geom_type() -> None: + gdf = GeoDataFrame({"id": [1]}, geometry=[box(0, 0, 1, 1)]) + reference_gdf = GeoDataFrame({"id": [1]}, geometry=[Point(0, 0)]) + + with pytest.raises( + GeometryTypeError, + match=re.escape("Reference data must only contain LineStrings."), + ): + GeneralizeWaterAreas().execute( + data=gdf, reference_data={"shoreline": reference_gdf} + ) diff --git a/test/core/test_geometry.py b/test/core/test_geometry.py index 635d103..8a90e47 100644 --- a/test/core/test_geometry.py +++ b/test/core/test_geometry.py @@ -5,33 +5,530 @@ # This source code is licensed under the MIT license found in the # LICENSE file in the root directory of this source tree. +import math + import pytest -from geopandas import GeoDataFrame, gpd +from geopandas import GeoDataFrame, GeoSeries +from geopandas.testing import assert_geoseries_equal from pandas import DataFrame from pandas.testing import assert_frame_equal -from shapely import GeometryCollection, MultiPolygon, from_wkt, to_wkt +from shapely import GeometryCollection, MultiPolygon, box, from_wkt, to_wkt from shapely.geometry import LineString, MultiLineString, MultiPoint, Point, Polygon from shapely.geometry.base import BaseGeometry from geogenalg.core.exceptions import ( - GeometryOperationError, GeometryTypeError, ) from geogenalg.core.geometry import ( + Dimensions, LineExtendFrom, assign_nearest_z, + chaikin_smooth_keep_topology, + chaikin_smooth_skip_coords, elongation, explode_line, extend_line_to_nearest, extract_interior_rings, + extract_interior_rings_gdf, + get_topological_points, lines_to_segments, mean_z, move_to_point, - rectangle_dimensions, + oriented_envelope_dimensions, + perforate_polygon_with_gdf_exteriors, scale_line_to_length, ) +@pytest.mark.parametrize( + ("input_geoseries", "extra_skip_coords", "expected_geoseries"), + [ + ( + GeoSeries(), + [], + GeoSeries(), + ), + ( + GeoSeries( + [ + box(0, 0, 2, 2), + box(0, 2, 2, 4), + box(0, 4, 2, 5), + ] + ), + None, + GeoSeries( + [ + Polygon( + [ + [2, 0.5], + [2, 2], + [0, 2], + [0, 0.5], + [0.5, 0], + [1.5, 0], + [2, 0.5], + ] + ), + box(0, 2, 2, 4), + Polygon( + [ + [2, 4], + [2, 4.75], + [1.5, 5], + [0.5, 5], + [0, 4.75], + [0, 4], + [2, 4], + ] + ), + ] + ), + ), + ], + ids=[ + "empty series", + "middle box not smoothed", + ], +) +def test_chaikin_smooth_keep_topology( + input_geoseries: GeoSeries, + extra_skip_coords: list[Point] | MultiPoint | None, + expected_geoseries: GeoSeries, +): + assert_geoseries_equal( + chaikin_smooth_keep_topology( + input_geoseries, + 1, + extra_skip_coords=extra_skip_coords, + ), + expected_geoseries, + ) + + +@pytest.mark.parametrize( + ("input_geometry", "skip_coords", "iterations", "expected_geometry"), + [ + ( + LineString( + [ + [0, 0], + [1, 1], + [1, 2], + ] + ), + [], + 1, + LineString( + [ + [0, 0], + [0.25, 0.25], + [0.75, 0.75], + [1, 1.25], + [1, 1.75], + [1, 2], + ] + ), + ), + ( + LineString( + [ + [0, 0], + [1, 1], + [1, 2], + ] + ), + [Point(1, 1)], + 1, + LineString( + [ + [0, 0], + [0.25, 0.25], + [1, 1], + [1, 1.75], + [1, 2], + ] + ), + ), + ( + LineString( + [ + [0, 0], + [1, 1], + [1, 2], + ] + ), + [Point(1, 1)], + 2, + LineString( + [ + [0, 0], + [0.0625, 0.0625], + [0.1875, 0.1875], + [0.4375, 0.4375], + [1, 1], + [1, 1.5625], + [1, 1.8125], + [1, 1.9375], + [1, 2], + ] + ), + ), + ( + box(0, 0, 1, 1), + [], + 1, + Polygon( + [ + [1, 0.25], + [1, 0.75], + [0.75, 1], + [0.25, 1], + [0, 0.75], + [0, 0.25], + [0.25, 0], + [0.75, 0], + [1, 0.25], + ] + ), + ), + ( + box(0, 0, 1, 1), + [Point(0, 0)], + 1, + Polygon( + [ + [1, 0.25], + [1, 0.75], + [0.75, 1], + [0.25, 1], + [0, 0.75], + [0, 0], + [0.75, 0], + [1, 0.25], + ] + ), + ), + ( + box(0, 0, 1, 1), + MultiPoint([Point(0, 0), Point(1, 1)]), + 1, + Polygon( + [ + [1, 0.25], + [1, 1], + [0.25, 1], + [0, 0.75], + [0, 0], + [0.75, 0], + [1, 0.25], + ] + ), + ), + ], + ids=[ + "three point linestring, no skip", + "three point linestring, skip one", + "three point linestring, skip one, two iterations", + "square polygon, no skip", + "square polygon, skip one", + "square polygon, skip two", + ], +) +def test_chaikin_smooth_skip_coords( + input_geometry: LineString | Polygon, + skip_coords: list[Point] | MultiPoint, + iterations: int, + expected_geometry: LineString | Polygon, +): + assert ( + chaikin_smooth_skip_coords( + input_geometry, + skip_coords, + iterations=iterations, + ) + == expected_geometry + ) + + +@pytest.mark.parametrize( + ("input_geoseries", "expected_points"), + [ + ( + GeoSeries(), + [], + ), + ( + GeoSeries(Point(0, 0)), + [], + ), + ( + GeoSeries( + [ + LineString([[0, 0], [0, 1]]), + LineString([[0, 0], [0, -1]]), + ] + ), + [Point(0, 0)], + ), + ( + GeoSeries( + [ + LineString([[0, 0], [0, 1]]), + LineString([[0, 0], [0, -1]]), + LineString([[0, 0], [-1, 0]]), + ] + ), + [Point(0, 0)], + ), + ( + GeoSeries( + [ + box(0, 0, 2, 2), + box(0, 2, 2, 4), + box(0, 4, 2, 5), + ] + ), + [ + Point(0, 2), + Point(0, 4), + Point(2, 2), + Point(2, 4), + ], + ), + ], + ids=[ + "empty series", + "no points", + "two geoms, one shared point", + "three geoms, one shared point", + "polygons, many shared points", + ], +) +def test_get_topological_points( + input_geoseries: GeoSeries, + expected_points: list[Point], +): + assert get_topological_points(input_geoseries) == expected_points + + +@pytest.mark.parametrize( + ("input_geometry", "gdf", "expected_geometry"), + [ + ( + box(0, 0, 20, 20), + GeoDataFrame( + geometry=[ + box(5, 5, 6, 6), + box(7, 7, 8, 8), + ], + ), + Polygon( + shell=[[0, 0], [0, 20], [20, 20], [20, 0], [0, 0]], + holes=[ + [[6, 5], [6, 6], [5, 6], [5, 5], [6, 5]], + [[8, 7], [8, 8], [7, 8], [7, 7], [8, 7]], + ], + ), + ), + ( + box(0, 0, 20, 20), + GeoDataFrame( + geometry=[ + # Feature within ring + # Feature with ring + Polygon( + shell=[[1, 1], [1, 10], [10, 10], [10, 1], [1, 1]], + holes=[[[2, 2], [2, 9], [9, 9], [9, 2], [2, 2]]], + ), + Polygon( + shell=[[3, 3], [3, 8], [8, 8], [8, 3], [3, 3]], + holes=[[[4, 4], [4, 7], [7, 7], [7, 4], [4, 4]]], + ), + ## + # Another feature within ring + # Feature with ring + Polygon( + shell=[[11, 11], [11, 19], [19, 19], [19, 11], [11, 11]], + holes=[[[12, 12], [12, 18], [18, 18], [18, 12], [12, 12]]], + ), + # Feature within ring + Polygon( + shell=[[13, 13], [13, 17], [17, 17], [17, 13], [13, 13]], + holes=[[[14, 14], [14, 16], [16, 16], [16, 14], [14, 14]]], + ), + ## + # Regular standalone feature + box(1, 11, 9, 19), + ], + ), + Polygon( + shell=[[0, 0], [0, 20], [20, 20], [20, 0]], + holes=[ + [[1, 10], [1, 1], [10, 1], [10, 10], [1, 10]], + [[1, 11], [9, 11], [9, 19], [1, 19], [1, 11]], + [[19, 19], [11, 19], [11, 11], [19, 11], [19, 19]], + ], + ), + ), + ], + ids=[ + "two exteriors", + "recursive features", + ], +) +def test_perforate_polygon_with_gdf_exteriors( + input_geometry: Polygon, + gdf: GeoDataFrame, + expected_geometry: Polygon, +): + assert ( + perforate_polygon_with_gdf_exteriors(input_geometry, gdf) == expected_geometry + ) + + +@pytest.mark.parametrize( + ("input_geometry", "expected_geometry"), + [ + ( + Polygon( + shell=[[0, 0, 4.0], [0, 2, 5.5], [2, 2, 6.5], [2, 0, 1.5]], + holes=[ + [ + [0.5, 0.5, 4.0], + [0.5, 1, 1.0], + [1, 1, 2.3], + [1, 0.5, 0.5], + ] + ], + ), + MultiPolygon( + [ + Polygon( + shell=[ + [0.5, 0.5, 4.0], + [0.5, 1, 1.0], + [1, 1, 2.3], + [1, 0.5, 0.5], + ], + ) + ] + ), + ), + ( + Polygon( + shell=[[0, 0], [0, 20], [20, 20], [20, 0]], + holes=[ + [ + [0.5, 0.5], + [0.5, 1], + [1, 1], + [1, 0.5], + ], + [ + [6, 6], + [6, 8], + [8, 8], + [8, 6], + ], + ], + ), + MultiPolygon( + [ + Polygon( + shell=[ + [0.5, 0.5], + [0.5, 1], + [1, 1], + [1, 0.5], + ], + ), + Polygon( + shell=[ + [6, 6], + [6, 8], + [8, 8], + [8, 6], + ], + ), + ] + ), + ), + ( + MultiPolygon( + [ + Polygon( + shell=[[0, 0], [0, 20], [20, 20], [20, 0]], + holes=[ + [[0.5, 0.5], [0.5, 1], [1, 1], [1, 0.5]], + [[6.5, 6.5], [6.5, 7.5], [7.5, 7.5], [7.5, 6.5]], + ], + ) + ] + ), + MultiPolygon( + [ + Polygon( + shell=[[0.5, 0.5], [0.5, 1], [1, 1], [1, 0.5]], + ), + Polygon( + shell=[[6.5, 6.5], [6.5, 7.5], [7.5, 7.5], [7.5, 6.5]], + ), + ] + ), + ), + ( + MultiPolygon( + [ + Polygon( + shell=[[0, 0], [0, 20], [20, 20], [20, 0]], + holes=[ + [[0.5, 0.5], [0.5, 1], [1, 1], [1, 0.5]], + ], + ), + Polygon( + shell=[[100, 100], [100, 110], [110, 110], [110, 100]], + holes=[[[105, 105], [105, 106], [106, 106], [106, 105]]], + ), + box(-10, -10, -5, -5), + ] + ), + MultiPolygon( + [ + Polygon( + shell=[[0.5, 0.5], [0.5, 1], [1, 1], [1, 0.5]], + ), + Polygon( + shell=[[105, 105], [105, 106], [106, 106], [106, 105]], + ), + ] + ), + ), + ( + Polygon(), + MultiPolygon(), + ), + ( + box(10, 10, 20, 20), + MultiPolygon(), + ), + ], + ids=[ + "polygon, single ring", + "polygon, two rings", + "multipolygon, two rings", + "multipolygon, three parts", + "empty polygon", + "no rings", + ], +) +def test_extract_interior_rings( + input_geometry: Polygon | MultiPolygon, + expected_geometry: MultiPolygon, +): + assert extract_interior_rings(input_geometry) == expected_geometry + + @pytest.mark.parametrize( ("geometry", "expected_mean"), [ @@ -169,24 +666,73 @@ def test_mean_z_no_geometrycollection(): mean_z(GeometryCollection([Point(0, 0, 0), Point(0, 0, 1)])) -def test_rectangle_dimensions(): - rect = from_wkt("POLYGON ((0 0, 0 4, 1 4, 1 0, 0 0))") - - assert rectangle_dimensions(rect).width == 1.0 - assert rectangle_dimensions(rect).height == 4.0 - - rect2 = from_wkt("POLYGON ((0 0, 0 4, 1 4, 1 5, 1 0, 0 0))") - with pytest.raises(GeometryOperationError): - rectangle_dimensions(rect2) - +@pytest.mark.parametrize( + ("geom", "expected_dimensions"), + [ + (box(0, 0, 1, 1), Dimensions(width=1, height=1)), + (box(0, 1, 2, 2), Dimensions(width=1.0, height=2.0)), + ( + Polygon( + [ + [0, 0], + [2, 0], + [2, 1], + [4, 1], + [4, -2], + [2, -2], + [2, -1], + [0, -1], + [0, 0], + ], + ), + Dimensions(width=3.0, height=4.0), + ), + ], + ids=[ + "square", + "rectangle", + "irregular polygon", + ], +) +def test_oriented_envelope_dimensions(geom: Polygon, expected_dimensions: Dimensions): + assert oriented_envelope_dimensions(geom) == expected_dimensions -def test_elongation(): - rect = from_wkt("POLYGON ((0 0, 0 4, 1 4, 1 0, 0 0))") - assert elongation(rect) == 4.0 +@pytest.mark.parametrize( + ("polygon", "expected_elongation"), + [ + (Polygon([[0, 0], [0, 4], [1, 4], [1, 0], [0, 0]]), 0.25), + (Polygon([[0, 0], [0, 2], [2, 2], [2, 0], [0, 0]]), 1.0), + ( + Polygon( + [ + [0, 0], + [2, 0], + [2, 1], + [4, 1], + [4, -2], + [2, -2], + [2, -1], + [0, -1], + [0, 0], + ] + ), + 0.75, + ), + (Polygon([[0, 0], [0, 19], [1, 19], [1, 0], [0, 0]]), 0.052631579), + ], + ids=[ + "elongated", + "square", + "irregular shape", + "very elongated", + ], +) +def test_elongation(polygon: Polygon, expected_elongation: float): + assert math.isclose(elongation(polygon), expected_elongation, rel_tol=1e-08) -def test_extract_interior_rings(): +def test_extract_interior_rings_gdf(): polygon_with_holes = from_wkt( """POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0), (2 3, 3 3, 3 2, 2 2, 2 3), @@ -194,12 +740,12 @@ def test_extract_interior_rings(): (7 8, 8 8, 8 7, 7 7, 7 8))""" ) - gdf = gpd.GeoDataFrame( + gdf = GeoDataFrame( geometry=[polygon_with_holes], crs="EPSG:3857", ) - extracted_holes = extract_interior_rings(gdf) + extracted_holes = extract_interior_rings_gdf(gdf) assert len(extracted_holes.index) == 3 assert extracted_holes.crs == "EPSG:3857" @@ -224,7 +770,7 @@ def test_explode_line_to_segments(): def test_lines_to_segments(): - lines = gpd.GeoSeries( + lines = GeoSeries( [ from_wkt("LINESTRING (0 0, 0 1)"), from_wkt("LINESTRING (10 0, 10 15)"), @@ -242,7 +788,7 @@ def test_lines_to_segments(): assert segments[3].wkt == "LINESTRING (21 21, 22 22)" assert segments[4].wkt == "LINESTRING (22 22, 23 23)" - invalid_series = gpd.GeoSeries( + invalid_series = GeoSeries( [ from_wkt("LINESTRING (0 0, 0 1)"), from_wkt("LINESTRING (10 0, 10 15)"), diff --git a/test/test_exaggeration.py b/test/test_exaggeration.py index a1a3801..bbbe015 100644 --- a/test/test_exaggeration.py +++ b/test/test_exaggeration.py @@ -6,8 +6,9 @@ # LICENSE file in the root directory of this source tree. import re -import geopandas as gpd import pytest +from geopandas import GeoDataFrame +from geopandas.testing import assert_geodataframe_equal from pandas.testing import assert_frame_equal from shapely import equals_exact from shapely.geometry import LineString, MultiPolygon, Point, Polygon @@ -16,33 +17,115 @@ from geogenalg.core.exceptions import GeometryTypeError +@pytest.mark.parametrize( + ( + "input_gdf", + "width_threshold", + "elongation_threshold", + "expected_gdf", + ), + [ + ( + GeoDataFrame( + geometry=[ + Polygon( + [[0, 0], [0, 4], [1, 4], [1, 0], [0, 0]] + ), # elongated enough, too wide + Polygon( + [[0, 0], [0, 0.5], [0.5, 0.5], [0.5, 0], [0, 0]] + ), # thin enough, not elongated enough + Polygon( + [[0, 0], [0, 4], [4, 4], [4, 0], [0, 0]] + ), # not elongated enough, too wide + ], + ), + 0.55, # width_threshold + 0.5, # elongation_threshold + GeoDataFrame( + geometry=[ + Polygon([[0, 0], [0, 4], [1, 4], [1, 0], [0, 0]]), + Polygon([[0, 0], [0, 0.5], [0.5, 0.5], [0.5, 0], [0, 0]]), + Polygon([[0, 0], [0, 4], [4, 4], [4, 0], [0, 0]]), + ], + ), + ), + ( + GeoDataFrame( + geometry=[ + Polygon( + [[0, 0], [0, 4], [1, 4], [1, 0], [0, 0]] + ), # should exaggerate + Polygon( + [[0, 0], [0, 0.5], [0.5, 0.5], [0.5, 0], [0, 0]] + ), # not elongated enough + Polygon( + [[0, 0], [0, 4], [4, 4], [4, 0], [0, 0]] + ), # too wide and not elongated enough + ], + ), + 2, # width_threshold + 0.5, # elongation_threshold + GeoDataFrame( + geometry=[ + Polygon([[-1, -1], [-1, 5], [2, 5], [2, -1], [-1, -1]]), + Polygon([[0, 0], [0, 0.5], [0.5, 0.5], [0.5, 0], [0, 0]]), + Polygon([[0, 0], [0, 4], [4, 4], [4, 0], [0, 0]]), + ], + ), + ), + ], + ids=[ + "no exaggerations", + "one exaggerated", + ], +) +def test_exaggerate_thin_polygons( + input_gdf: GeoDataFrame, + width_threshold: float, + elongation_threshold: float, + expected_gdf: GeoDataFrame, +): + result = exaggeration.exaggerate_thin_polygons( + input_gdf, + width_threshold, + elongation_threshold, + 1.0, + buffer_options=exaggeration.BufferOptions( + cap_style="square", + join_style="mitre", + ), + ) + + assert_geodataframe_equal(result, expected_gdf) + + @pytest.mark.parametrize( ("input_gdf", "threshold", "expected_gdf"), [ ( - gpd.GeoDataFrame( + GeoDataFrame( {"id": [1]}, geometry=[Polygon([(0, 0, 5), (0, 0.4, 5), (2, 0.4, 5), (2, 0, 5)])], ), 1.0, - gpd.GeoDataFrame( + GeoDataFrame( {"id": [1]}, geometry=[Polygon([(0, 0, 5), (0, 0.4, 5), (2, 0.4, 5), (2, 0, 5)])], ), ), ( - gpd.GeoDataFrame( + GeoDataFrame( {"id": [2]}, geometry=[Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])], ), 1.0, - gpd.GeoDataFrame( + GeoDataFrame( {"id": [2]}, geometry=[Polygon()], ), ), ( - gpd.GeoDataFrame( + GeoDataFrame( {"id": [3]}, geometry=[ MultiPolygon( @@ -54,36 +137,36 @@ ], ), 1.0, - gpd.GeoDataFrame( + GeoDataFrame( {"id": [3]}, geometry=[Polygon([(0, 0.4, 2), (2, 0.4, 3), (2, 0, 4), (0, 0, 1)])], ), ), ( - gpd.GeoDataFrame(columns=["id", "geometry"], geometry=[]), + GeoDataFrame(columns=["id", "geometry"], geometry=[]), 1.0, - gpd.GeoDataFrame(columns=["id", "geometry"], geometry=[]), + GeoDataFrame(columns=["id", "geometry"], geometry=[]), ), ( - gpd.GeoDataFrame( + GeoDataFrame( {"id": [4], "name": ["test-area"]}, geometry=[Polygon([(0, 0), (0, 0.4), (3, 0.4), (3, 0)])], ), 1.0, - gpd.GeoDataFrame( + GeoDataFrame( {"id": [4], "name": ["test-area"]}, geometry=[Polygon([(0, 0), (0, 0.4), (3, 0.4), (3, 0)])], ), ), ( - gpd.GeoDataFrame( + GeoDataFrame( {"id": [5]}, geometry=[ Polygon([(0, 0), (4, 0), (4, 2), (2.2, 2), (2.2, 0.4), (0, 0.4)]) ], ), 0.5, - gpd.GeoDataFrame( + GeoDataFrame( {"id": [5]}, geometry=[ Polygon( @@ -98,7 +181,7 @@ ), ), ( - gpd.GeoDataFrame( + GeoDataFrame( {"id": [6]}, geometry=[ Polygon( @@ -116,7 +199,7 @@ ], ), 1, - gpd.GeoDataFrame( + GeoDataFrame( {"id": [6]}, geometry=[ MultiPolygon( @@ -140,7 +223,7 @@ ], ) def test_extract_narrow_polygon_parts( - input_gdf: gpd.GeoDataFrame, threshold: float, expected_gdf: gpd.GeoDataFrame + input_gdf: GeoDataFrame, threshold: float, expected_gdf: GeoDataFrame ): result_gdf = exaggeration.extract_narrow_polygon_parts(input_gdf, threshold) @@ -155,7 +238,9 @@ def test_extract_narrow_polygon_parts( strict=False, ) ): - assert equals_exact(result_geometry, expected_geometry, tolerance=1e-3) + assert equals_exact( + result_geometry, expected_geometry, tolerance=1e-3, normalize=True + ) # Check attributes result_attrs = result_sorted.drop(columns="geometry") @@ -164,7 +249,7 @@ def test_extract_narrow_polygon_parts( def test_extract_narrow_polygon_parts_raises_on_non_polygon_geometry(): - invalid_gdf = gpd.GeoDataFrame( + invalid_gdf = GeoDataFrame( {"id": [1, 2, 3]}, geometry=[ Point(0, 0), diff --git a/test/test_testing.py b/test/test_testing.py new file mode 100644 index 0000000..bc902b7 --- /dev/null +++ b/test/test_testing.py @@ -0,0 +1,81 @@ +# Copyright (c) 2025 National Land Survey of Finland (Maanmittauslaitos) +# +# This file is part of geogen-algorithms. +# +# This source code is licensed under the MIT license found in the +# LICENSE file in the root directory of this source tree. + + +from pathlib import Path +from tempfile import TemporaryDirectory + +import pytest +from geopandas import GeoDataFrame, read_file +from geopandas.testing import assert_geodataframe_equal +from shapely import Point, Polygon, box + +from geogenalg.testing import assert_gdf_equal_save_diff + + +def test_assert_gdf_equal_save_diff_geom(): + temp_dir = TemporaryDirectory() + temp_dir_path = Path(temp_dir.name) + + gdf_result = GeoDataFrame( + geometry=[box(0.5, 0.5, 1.5, 1.5)], + ) + gdf_control = GeoDataFrame( + geometry=[box(0, 0, 1, 1)], + ) + + with pytest.raises(AssertionError): + assert_gdf_equal_save_diff(gdf_result, gdf_control, directory=temp_dir_path) + + result_path = temp_dir_path / "result.gpkg" + geomdiff_path = temp_dir_path / "geomdiff.gpkg" + control_path = temp_dir_path / "control.gpkg" + + assert result_path.exists() + assert geomdiff_path.exists() + assert control_path.exists() + + control = read_file(control_path) + result = read_file(result_path) + assert_geodataframe_equal(control, gdf_control) + assert_geodataframe_equal(result, gdf_result) + + geomdiff = read_file(geomdiff_path) + geomdiff_expected = GeoDataFrame( + geometry=[ + Polygon( + [ + [0.5, 1.5], + [1.5, 1.5], + [1.5, 0.5], + [1.0, 0.5], + [1.0, 1.0], + [0.5, 1.0], + [0.5, 1.5], + ] + ) + ], + ) + assert_geodataframe_equal(geomdiff, geomdiff_expected) + + +def test_assert_gdf_equal_save_diff_is_equal(): + temp_dir = TemporaryDirectory() + temp_dir_path = Path(temp_dir.name) + + gdf1 = GeoDataFrame( + {id: [1, 2]}, + geometry=[Point(0, 0), Point(1, 0)], + ) + + gdf2 = gdf1.copy() + + assert_gdf_equal_save_diff(gdf1, gdf2, directory=temp_dir_path) + + assert not (temp_dir_path / "result.gpkg").exists() + assert not (temp_dir_path / "geomdiff.gpkg").exists() + assert not (temp_dir_path / "control.gpkg").exists() diff --git a/test/testdata/buildings_generalized_100k_helsinki.gpkg b/test/testdata/buildings_generalized_100k_helsinki.gpkg index 74fb049..47337bd 100644 Binary files a/test/testdata/buildings_generalized_100k_helsinki.gpkg and b/test/testdata/buildings_generalized_100k_helsinki.gpkg differ diff --git a/test/testdata/buildings_generalized_50k_helsinki.gpkg b/test/testdata/buildings_generalized_50k_helsinki.gpkg index a2792ea..4626f5d 100644 Binary files a/test/testdata/buildings_generalized_50k_helsinki.gpkg and b/test/testdata/buildings_generalized_50k_helsinki.gpkg differ diff --git a/test/testdata/lakes.gpkg b/test/testdata/lakes.gpkg new file mode 100644 index 0000000..05aae88 Binary files /dev/null and b/test/testdata/lakes.gpkg differ diff --git a/test/testdata/marshes_small_area.gpkg b/test/testdata/marshes_small_area.gpkg index 586a642..8f6a59a 100644 Binary files a/test/testdata/marshes_small_area.gpkg and b/test/testdata/marshes_small_area.gpkg differ diff --git a/test/testdata/sea.gpkg b/test/testdata/sea.gpkg new file mode 100644 index 0000000..3db2d49 Binary files /dev/null and b/test/testdata/sea.gpkg differ diff --git a/uv.lock b/uv.lock index b332ade..e6d2600 100644 --- a/uv.lock +++ b/uv.lock @@ -319,12 +319,12 @@ lint = [ [package.metadata] requires-dist = [ { name = "cartagen", specifier = ">=1.0.2" }, - { name = "geopandas", specifier = "==1.0.1" }, + { name = "geopandas", specifier = "==1.1.1" }, { name = "matplotlib", specifier = "==3.9.2" }, { name = "networkx", specifier = "==3.4.2" }, { name = "pygeoops", specifier = "==0.4.0" }, { name = "scikit-learn", specifier = "==1.6.1" }, - { name = "shapely", specifier = "==2.0.6" }, + { name = "shapely", specifier = "==2.1.2" }, { name = "shapelysmooth", specifier = "==0.2.1" }, { name = "simplification", specifier = "==0.7.13" }, { name = "typer", specifier = ">=0.19.2" }, @@ -352,7 +352,7 @@ lint = [ [[package]] name = "geopandas" -version = "1.0.1" +version = "1.1.1" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "numpy" }, @@ -362,9 +362,9 @@ dependencies = [ { name = "pyproj" }, { name = "shapely" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/39/08/2cf5d85356e45b10b8d066cf4c3ba1e9e3185423c48104eed87e8afd0455/geopandas-1.0.1.tar.gz", hash = "sha256:b8bf70a5534588205b7a56646e2082fb1de9a03599651b3d80c99ea4c2ca08ab", size = 317736, upload-time = "2024-07-02T12:26:52.928Z" } +sdist = { url = "https://files.pythonhosted.org/packages/8c/76/e1960ba846f153ab109575242abf89dc98f8e057faa32f3decf4cce9247a/geopandas-1.1.1.tar.gz", hash = "sha256:1745713f64d095c43e72e08e753dbd271678254b24f2e01db8cdb8debe1d293d", size = 332655, upload-time = "2025-06-26T21:04:56.57Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/c4/64/7d344cfcef5efddf9cf32f59af7f855828e9d74b5f862eddf5bfd9f25323/geopandas-1.0.1-py3-none-any.whl", hash = "sha256:01e147d9420cc374d26f51fc23716ac307f32b49406e4bd8462c07e82ed1d3d6", size = 323587, upload-time = "2024-07-02T12:26:50.876Z" }, + { url = "https://files.pythonhosted.org/packages/0b/70/d5cd0696eff08e62fdbdebe5b46527facb4e7220eabe0ac6225efab50168/geopandas-1.1.1-py3-none-any.whl", hash = "sha256:589e61aaf39b19828843df16cb90234e72897e2579be236f10eee0d052ad98e8", size = 338365, upload-time = "2025-06-26T21:04:55.139Z" }, ] [[package]] @@ -1208,31 +1208,61 @@ wheels = [ [[package]] name = "shapely" -version = "2.0.6" +version = "2.1.2" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "numpy" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/4a/89/0d20bac88016be35ff7d3c0c2ae64b477908f1b1dfa540c5d69ac7af07fe/shapely-2.0.6.tar.gz", hash = "sha256:997f6159b1484059ec239cacaa53467fd8b5564dabe186cd84ac2944663b0bf6", size = 282361, upload-time = "2024-08-19T21:57:22.303Z" } -wheels = [ - { url = "https://files.pythonhosted.org/packages/37/15/269d8e1f7f658a37e61f7028683c546f520e4e7cedba1e32c77ff9d3a3c7/shapely-2.0.6-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:5aeb0f51a9db176da9a30cb2f4329b6fbd1e26d359012bb0ac3d3c7781667a9e", size = 1449578, upload-time = "2024-08-19T21:56:24.058Z" }, - { url = "https://files.pythonhosted.org/packages/37/63/e182e43081fffa0a2d970c480f2ef91647a6ab94098f61748c23c2a485f2/shapely-2.0.6-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:9a7a78b0d51257a367ee115f4d41ca4d46edbd0dd280f697a8092dd3989867b2", size = 1296792, upload-time = "2024-08-19T21:56:26.044Z" }, - { url = "https://files.pythonhosted.org/packages/6e/5a/d019f69449329dcd517355444fdb9ddd58bec5e080b8bdba007e8e4c546d/shapely-2.0.6-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f32c23d2f43d54029f986479f7c1f6e09c6b3a19353a3833c2ffb226fb63a855", size = 2443997, upload-time = "2024-08-19T22:00:54.836Z" }, - { url = "https://files.pythonhosted.org/packages/25/aa/53f145e5a610a49af9ac49f2f1be1ec8659ebd5c393d66ac94e57c83b00e/shapely-2.0.6-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b3dc9fb0eb56498912025f5eb352b5126f04801ed0e8bdbd867d21bdbfd7cbd0", size = 2528334, upload-time = "2024-08-19T21:56:27.53Z" }, - { url = "https://files.pythonhosted.org/packages/64/64/0c7b0a22b416d36f6296b92bb4219d82b53d0a7c47e16fd0a4c85f2f117c/shapely-2.0.6-cp311-cp311-win32.whl", hash = "sha256:d93b7e0e71c9f095e09454bf18dad5ea716fb6ced5df3cb044564a00723f339d", size = 1294669, upload-time = "2024-08-19T21:56:29.509Z" }, - { url = "https://files.pythonhosted.org/packages/b1/5a/6a67d929c467a1973b6bb9f0b00159cc343b02bf9a8d26db1abd2f87aa23/shapely-2.0.6-cp311-cp311-win_amd64.whl", hash = "sha256:c02eb6bf4cfb9fe6568502e85bb2647921ee49171bcd2d4116c7b3109724ef9b", size = 1442032, upload-time = "2024-08-19T21:56:31.158Z" }, - { url = "https://files.pythonhosted.org/packages/46/77/efd9f9d4b6a762f976f8b082f54c9be16f63050389500fb52e4f6cc07c1a/shapely-2.0.6-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:cec9193519940e9d1b86a3b4f5af9eb6910197d24af02f247afbfb47bcb3fab0", size = 1450326, upload-time = "2024-08-19T21:56:33.166Z" }, - { url = "https://files.pythonhosted.org/packages/68/53/5efa6e7a4036a94fe6276cf7bbb298afded51ca3396b03981ad680c8cc7d/shapely-2.0.6-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:83b94a44ab04a90e88be69e7ddcc6f332da7c0a0ebb1156e1c4f568bbec983c3", size = 1298480, upload-time = "2024-08-19T21:56:35.317Z" }, - { url = "https://files.pythonhosted.org/packages/88/a2/1be1db4fc262e536465a52d4f19d85834724fedf2299a1b9836bc82fe8fa/shapely-2.0.6-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:537c4b2716d22c92036d00b34aac9d3775e3691f80c7aa517c2c290351f42cd8", size = 2439311, upload-time = "2024-08-19T22:01:00.611Z" }, - { url = "https://files.pythonhosted.org/packages/d5/7d/9a57e187cbf2fbbbdfd4044a4f9ce141c8d221f9963750d3b001f0ec080d/shapely-2.0.6-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:98fea108334be345c283ce74bf064fa00cfdd718048a8af7343c59eb40f59726", size = 2524835, upload-time = "2024-08-19T21:56:36.87Z" }, - { url = "https://files.pythonhosted.org/packages/6d/0a/f407509ab56825f39bf8cfce1fb410238da96cf096809c3e404e5bc71ea1/shapely-2.0.6-cp312-cp312-win32.whl", hash = "sha256:42fd4cd4834747e4990227e4cbafb02242c0cffe9ce7ef9971f53ac52d80d55f", size = 1295613, upload-time = "2024-08-19T21:56:38.962Z" }, - { url = "https://files.pythonhosted.org/packages/7b/b3/857afd9dfbfc554f10d683ac412eac6fa260d1f4cd2967ecb655c57e831a/shapely-2.0.6-cp312-cp312-win_amd64.whl", hash = "sha256:665990c84aece05efb68a21b3523a6b2057e84a1afbef426ad287f0796ef8a48", size = 1442539, upload-time = "2024-08-19T21:56:40.686Z" }, - { url = "https://files.pythonhosted.org/packages/34/e8/d164ef5b0eab86088cde06dee8415519ffd5bb0dd1bd9d021e640e64237c/shapely-2.0.6-cp313-cp313-macosx_10_13_x86_64.whl", hash = "sha256:42805ef90783ce689a4dde2b6b2f261e2c52609226a0438d882e3ced40bb3013", size = 1445344, upload-time = "2024-08-19T21:56:42.714Z" }, - { url = "https://files.pythonhosted.org/packages/ce/e2/9fba7ac142f7831757a10852bfa465683724eadbc93d2d46f74a16f9af04/shapely-2.0.6-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:6d2cb146191a47bd0cee8ff5f90b47547b82b6345c0d02dd8b25b88b68af62d7", size = 1296182, upload-time = "2024-08-19T21:56:44.64Z" }, - { url = "https://files.pythonhosted.org/packages/cf/dc/790d4bda27d196cd56ec66975eaae3351c65614cafd0e16ddde39ec9fb92/shapely-2.0.6-cp313-cp313-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:e3fdef0a1794a8fe70dc1f514440aa34426cc0ae98d9a1027fb299d45741c381", size = 2423426, upload-time = "2024-08-19T22:01:05.167Z" }, - { url = "https://files.pythonhosted.org/packages/af/b0/f8169f77eac7392d41e231911e0095eb1148b4d40c50ea9e34d999c89a7e/shapely-2.0.6-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:2c665a0301c645615a107ff7f52adafa2153beab51daf34587170d85e8ba6805", size = 2513249, upload-time = "2024-08-19T21:56:47.1Z" }, - { url = "https://files.pythonhosted.org/packages/f6/1d/a8c0e9ab49ff2f8e4dedd71b0122eafb22a18ad7e9d256025e1f10c84704/shapely-2.0.6-cp313-cp313-win32.whl", hash = "sha256:0334bd51828f68cd54b87d80b3e7cee93f249d82ae55a0faf3ea21c9be7b323a", size = 1294848, upload-time = "2024-08-19T21:56:48.914Z" }, - { url = "https://files.pythonhosted.org/packages/23/38/2bc32dd1e7e67a471d4c60971e66df0bdace88656c47a9a728ace0091075/shapely-2.0.6-cp313-cp313-win_amd64.whl", hash = "sha256:d37d070da9e0e0f0a530a621e17c0b8c3c9d04105655132a87cfff8bd77cc4c2", size = 1441371, upload-time = "2024-08-19T21:56:51.108Z" }, +sdist = { url = "https://files.pythonhosted.org/packages/4d/bc/0989043118a27cccb4e906a46b7565ce36ca7b57f5a18b78f4f1b0f72d9d/shapely-2.1.2.tar.gz", hash = "sha256:2ed4ecb28320a433db18a5bf029986aa8afcfd740745e78847e330d5d94922a9", size = 315489, upload-time = "2025-09-24T13:51:41.432Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/8f/8d/1ff672dea9ec6a7b5d422eb6d095ed886e2e523733329f75fdcb14ee1149/shapely-2.1.2-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:91121757b0a36c9aac3427a651a7e6567110a4a67c97edf04f8d55d4765f6618", size = 1820038, upload-time = "2025-09-24T13:50:15.628Z" }, + { url = "https://files.pythonhosted.org/packages/4f/ce/28fab8c772ce5db23a0d86bf0adaee0c4c79d5ad1db766055fa3dab442e2/shapely-2.1.2-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:16a9c722ba774cf50b5d4541242b4cce05aafd44a015290c82ba8a16931ff63d", size = 1626039, upload-time = "2025-09-24T13:50:16.881Z" }, + { url = "https://files.pythonhosted.org/packages/70/8b/868b7e3f4982f5006e9395c1e12343c66a8155c0374fdc07c0e6a1ab547d/shapely-2.1.2-cp311-cp311-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:cc4f7397459b12c0b196c9efe1f9d7e92463cbba142632b4cc6d8bbbbd3e2b09", size = 3001519, upload-time = "2025-09-24T13:50:18.606Z" }, + { url = "https://files.pythonhosted.org/packages/13/02/58b0b8d9c17c93ab6340edd8b7308c0c5a5b81f94ce65705819b7416dba5/shapely-2.1.2-cp311-cp311-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:136ab87b17e733e22f0961504d05e77e7be8c9b5a8184f685b4a91a84efe3c26", size = 3110842, upload-time = "2025-09-24T13:50:21.77Z" }, + { url = "https://files.pythonhosted.org/packages/af/61/8e389c97994d5f331dcffb25e2fa761aeedfb52b3ad9bcdd7b8671f4810a/shapely-2.1.2-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:16c5d0fc45d3aa0a69074979f4f1928ca2734fb2e0dde8af9611e134e46774e7", size = 4021316, upload-time = "2025-09-24T13:50:23.626Z" }, + { url = "https://files.pythonhosted.org/packages/d3/d4/9b2a9fe6039f9e42ccf2cb3e84f219fd8364b0c3b8e7bbc857b5fbe9c14c/shapely-2.1.2-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:6ddc759f72b5b2b0f54a7e7cde44acef680a55019eb52ac63a7af2cf17cb9cd2", size = 4178586, upload-time = "2025-09-24T13:50:25.443Z" }, + { url = "https://files.pythonhosted.org/packages/16/f6/9840f6963ed4decf76b08fd6d7fed14f8779fb7a62cb45c5617fa8ac6eab/shapely-2.1.2-cp311-cp311-win32.whl", hash = "sha256:2fa78b49485391224755a856ed3b3bd91c8455f6121fee0db0e71cefb07d0ef6", size = 1543961, upload-time = "2025-09-24T13:50:26.968Z" }, + { url = "https://files.pythonhosted.org/packages/38/1e/3f8ea46353c2a33c1669eb7327f9665103aa3a8dfe7f2e4ef714c210b2c2/shapely-2.1.2-cp311-cp311-win_amd64.whl", hash = "sha256:c64d5c97b2f47e3cd9b712eaced3b061f2b71234b3fc263e0fcf7d889c6559dc", size = 1722856, upload-time = "2025-09-24T13:50:28.497Z" }, + { url = "https://files.pythonhosted.org/packages/24/c0/f3b6453cf2dfa99adc0ba6675f9aaff9e526d2224cbd7ff9c1a879238693/shapely-2.1.2-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:fe2533caae6a91a543dec62e8360fe86ffcdc42a7c55f9dfd0128a977a896b94", size = 1833550, upload-time = "2025-09-24T13:50:30.019Z" }, + { url = "https://files.pythonhosted.org/packages/86/07/59dee0bc4b913b7ab59ab1086225baca5b8f19865e6101db9ebb7243e132/shapely-2.1.2-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:ba4d1333cc0bc94381d6d4308d2e4e008e0bd128bdcff5573199742ee3634359", size = 1643556, upload-time = "2025-09-24T13:50:32.291Z" }, + { url = "https://files.pythonhosted.org/packages/26/29/a5397e75b435b9895cd53e165083faed5d12fd9626eadec15a83a2411f0f/shapely-2.1.2-cp312-cp312-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:0bd308103340030feef6c111d3eb98d50dc13feea33affc8a6f9fa549e9458a3", size = 2988308, upload-time = "2025-09-24T13:50:33.862Z" }, + { url = "https://files.pythonhosted.org/packages/b9/37/e781683abac55dde9771e086b790e554811a71ed0b2b8a1e789b7430dd44/shapely-2.1.2-cp312-cp312-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:1e7d4d7ad262a48bb44277ca12c7c78cb1b0f56b32c10734ec9a1d30c0b0c54b", size = 3099844, upload-time = "2025-09-24T13:50:35.459Z" }, + { url = "https://files.pythonhosted.org/packages/d8/f3/9876b64d4a5a321b9dc482c92bb6f061f2fa42131cba643c699f39317cb9/shapely-2.1.2-cp312-cp312-musllinux_1_2_aarch64.whl", hash = "sha256:e9eddfe513096a71896441a7c37db72da0687b34752c4e193577a145c71736fc", size = 3988842, upload-time = "2025-09-24T13:50:37.478Z" }, + { url = "https://files.pythonhosted.org/packages/d1/a0/704c7292f7014c7e74ec84eddb7b109e1fbae74a16deae9c1504b1d15565/shapely-2.1.2-cp312-cp312-musllinux_1_2_x86_64.whl", hash = "sha256:980c777c612514c0cf99bc8a9de6d286f5e186dcaf9091252fcd444e5638193d", size = 4152714, upload-time = "2025-09-24T13:50:39.9Z" }, + { url = "https://files.pythonhosted.org/packages/53/46/319c9dc788884ad0785242543cdffac0e6530e4d0deb6c4862bc4143dcf3/shapely-2.1.2-cp312-cp312-win32.whl", hash = "sha256:9111274b88e4d7b54a95218e243282709b330ef52b7b86bc6aaf4f805306f454", size = 1542745, upload-time = "2025-09-24T13:50:41.414Z" }, + { url = "https://files.pythonhosted.org/packages/ec/bf/cb6c1c505cb31e818e900b9312d514f381fbfa5c4363edfce0fcc4f8c1a4/shapely-2.1.2-cp312-cp312-win_amd64.whl", hash = "sha256:743044b4cfb34f9a67205cee9279feaf60ba7d02e69febc2afc609047cb49179", size = 1722861, upload-time = "2025-09-24T13:50:43.35Z" }, + { url = "https://files.pythonhosted.org/packages/c3/90/98ef257c23c46425dc4d1d31005ad7c8d649fe423a38b917db02c30f1f5a/shapely-2.1.2-cp313-cp313-macosx_10_13_x86_64.whl", hash = "sha256:b510dda1a3672d6879beb319bc7c5fd302c6c354584690973c838f46ec3e0fa8", size = 1832644, upload-time = "2025-09-24T13:50:44.886Z" }, + { url = "https://files.pythonhosted.org/packages/6d/ab/0bee5a830d209adcd3a01f2d4b70e587cdd9fd7380d5198c064091005af8/shapely-2.1.2-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:8cff473e81017594d20ec55d86b54bc635544897e13a7cfc12e36909c5309a2a", size = 1642887, upload-time = "2025-09-24T13:50:46.735Z" }, + { url = "https://files.pythonhosted.org/packages/2d/5e/7d7f54ba960c13302584c73704d8c4d15404a51024631adb60b126a4ae88/shapely-2.1.2-cp313-cp313-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:fe7b77dc63d707c09726b7908f575fc04ff1d1ad0f3fb92aec212396bc6cfe5e", size = 2970931, upload-time = "2025-09-24T13:50:48.374Z" }, + { url = "https://files.pythonhosted.org/packages/f2/a2/83fc37e2a58090e3d2ff79175a95493c664bcd0b653dd75cb9134645a4e5/shapely-2.1.2-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:7ed1a5bbfb386ee8332713bf7508bc24e32d24b74fc9a7b9f8529a55db9f4ee6", size = 3082855, upload-time = "2025-09-24T13:50:50.037Z" }, + { url = "https://files.pythonhosted.org/packages/44/2b/578faf235a5b09f16b5f02833c53822294d7f21b242f8e2d0cf03fb64321/shapely-2.1.2-cp313-cp313-musllinux_1_2_aarch64.whl", hash = "sha256:a84e0582858d841d54355246ddfcbd1fce3179f185da7470f41ce39d001ee1af", size = 3979960, upload-time = "2025-09-24T13:50:51.74Z" }, + { url = "https://files.pythonhosted.org/packages/4d/04/167f096386120f692cc4ca02f75a17b961858997a95e67a3cb6a7bbd6b53/shapely-2.1.2-cp313-cp313-musllinux_1_2_x86_64.whl", hash = "sha256:dc3487447a43d42adcdf52d7ac73804f2312cbfa5d433a7d2c506dcab0033dfd", size = 4142851, upload-time = "2025-09-24T13:50:53.49Z" }, + { url = "https://files.pythonhosted.org/packages/48/74/fb402c5a6235d1c65a97348b48cdedb75fb19eca2b1d66d04969fc1c6091/shapely-2.1.2-cp313-cp313-win32.whl", hash = "sha256:9c3a3c648aedc9f99c09263b39f2d8252f199cb3ac154fadc173283d7d111350", size = 1541890, upload-time = "2025-09-24T13:50:55.337Z" }, + { url = "https://files.pythonhosted.org/packages/41/47/3647fe7ad990af60ad98b889657a976042c9988c2807cf322a9d6685f462/shapely-2.1.2-cp313-cp313-win_amd64.whl", hash = "sha256:ca2591bff6645c216695bdf1614fca9c82ea1144d4a7591a466fef64f28f0715", size = 1722151, upload-time = "2025-09-24T13:50:57.153Z" }, + { url = "https://files.pythonhosted.org/packages/3c/49/63953754faa51ffe7d8189bfbe9ca34def29f8c0e34c67cbe2a2795f269d/shapely-2.1.2-cp313-cp313t-macosx_10_13_x86_64.whl", hash = "sha256:2d93d23bdd2ed9dc157b46bc2f19b7da143ca8714464249bef6771c679d5ff40", size = 1834130, upload-time = "2025-09-24T13:50:58.49Z" }, + { url = "https://files.pythonhosted.org/packages/7f/ee/dce001c1984052970ff60eb4727164892fb2d08052c575042a47f5a9e88f/shapely-2.1.2-cp313-cp313t-macosx_11_0_arm64.whl", hash = "sha256:01d0d304b25634d60bd7cf291828119ab55a3bab87dc4af1e44b07fb225f188b", size = 1642802, upload-time = "2025-09-24T13:50:59.871Z" }, + { url = "https://files.pythonhosted.org/packages/da/e7/fc4e9a19929522877fa602f705706b96e78376afb7fad09cad5b9af1553c/shapely-2.1.2-cp313-cp313t-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:8d8382dd120d64b03698b7298b89611a6ea6f55ada9d39942838b79c9bc89801", size = 3018460, upload-time = "2025-09-24T13:51:02.08Z" }, + { url = "https://files.pythonhosted.org/packages/a1/18/7519a25db21847b525696883ddc8e6a0ecaa36159ea88e0fef11466384d0/shapely-2.1.2-cp313-cp313t-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:19efa3611eef966e776183e338b2d7ea43569ae99ab34f8d17c2c054d3205cc0", size = 3095223, upload-time = "2025-09-24T13:51:04.472Z" }, + { url = "https://files.pythonhosted.org/packages/48/de/b59a620b1f3a129c3fecc2737104a0a7e04e79335bd3b0a1f1609744cf17/shapely-2.1.2-cp313-cp313t-musllinux_1_2_aarch64.whl", hash = "sha256:346ec0c1a0fcd32f57f00e4134d1200e14bf3f5ae12af87ba83ca275c502498c", size = 4030760, upload-time = "2025-09-24T13:51:06.455Z" }, + { url = "https://files.pythonhosted.org/packages/96/b3/c6655ee7232b417562bae192ae0d3ceaadb1cc0ffc2088a2ddf415456cc2/shapely-2.1.2-cp313-cp313t-musllinux_1_2_x86_64.whl", hash = "sha256:6305993a35989391bd3476ee538a5c9a845861462327efe00dd11a5c8c709a99", size = 4170078, upload-time = "2025-09-24T13:51:08.584Z" }, + { url = "https://files.pythonhosted.org/packages/a0/8e/605c76808d73503c9333af8f6cbe7e1354d2d238bda5f88eea36bfe0f42a/shapely-2.1.2-cp313-cp313t-win32.whl", hash = "sha256:c8876673449f3401f278c86eb33224c5764582f72b653a415d0e6672fde887bf", size = 1559178, upload-time = "2025-09-24T13:51:10.73Z" }, + { url = "https://files.pythonhosted.org/packages/36/f7/d317eb232352a1f1444d11002d477e54514a4a6045536d49d0c59783c0da/shapely-2.1.2-cp313-cp313t-win_amd64.whl", hash = "sha256:4a44bc62a10d84c11a7a3d7c1c4fe857f7477c3506e24c9062da0db0ae0c449c", size = 1739756, upload-time = "2025-09-24T13:51:12.105Z" }, + { url = "https://files.pythonhosted.org/packages/fc/c4/3ce4c2d9b6aabd27d26ec988f08cb877ba9e6e96086eff81bfea93e688c7/shapely-2.1.2-cp314-cp314-macosx_10_13_x86_64.whl", hash = "sha256:9a522f460d28e2bf4e12396240a5fc1518788b2fcd73535166d748399ef0c223", size = 1831290, upload-time = "2025-09-24T13:51:13.56Z" }, + { url = "https://files.pythonhosted.org/packages/17/b9/f6ab8918fc15429f79cb04afa9f9913546212d7fb5e5196132a2af46676b/shapely-2.1.2-cp314-cp314-macosx_11_0_arm64.whl", hash = "sha256:1ff629e00818033b8d71139565527ced7d776c269a49bd78c9df84e8f852190c", size = 1641463, upload-time = "2025-09-24T13:51:14.972Z" }, + { url = "https://files.pythonhosted.org/packages/a5/57/91d59ae525ca641e7ac5551c04c9503aee6f29b92b392f31790fcb1a4358/shapely-2.1.2-cp314-cp314-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:f67b34271dedc3c653eba4e3d7111aa421d5be9b4c4c7d38d30907f796cb30df", size = 2970145, upload-time = "2025-09-24T13:51:16.961Z" }, + { url = "https://files.pythonhosted.org/packages/8a/cb/4948be52ee1da6927831ab59e10d4c29baa2a714f599f1f0d1bc747f5777/shapely-2.1.2-cp314-cp314-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:21952dc00df38a2c28375659b07a3979d22641aeb104751e769c3ee825aadecf", size = 3073806, upload-time = "2025-09-24T13:51:18.712Z" }, + { url = "https://files.pythonhosted.org/packages/03/83/f768a54af775eb41ef2e7bec8a0a0dbe7d2431c3e78c0a8bdba7ab17e446/shapely-2.1.2-cp314-cp314-musllinux_1_2_aarch64.whl", hash = "sha256:1f2f33f486777456586948e333a56ae21f35ae273be99255a191f5c1fa302eb4", size = 3980803, upload-time = "2025-09-24T13:51:20.37Z" }, + { url = "https://files.pythonhosted.org/packages/9f/cb/559c7c195807c91c79d38a1f6901384a2878a76fbdf3f1048893a9b7534d/shapely-2.1.2-cp314-cp314-musllinux_1_2_x86_64.whl", hash = "sha256:cf831a13e0d5a7eb519e96f58ec26e049b1fad411fc6fc23b162a7ce04d9cffc", size = 4133301, upload-time = "2025-09-24T13:51:21.887Z" }, + { url = "https://files.pythonhosted.org/packages/80/cd/60d5ae203241c53ef3abd2ef27c6800e21afd6c94e39db5315ea0cbafb4a/shapely-2.1.2-cp314-cp314-win32.whl", hash = "sha256:61edcd8d0d17dd99075d320a1dd39c0cb9616f7572f10ef91b4b5b00c4aeb566", size = 1583247, upload-time = "2025-09-24T13:51:23.401Z" }, + { url = "https://files.pythonhosted.org/packages/74/d4/135684f342e909330e50d31d441ace06bf83c7dc0777e11043f99167b123/shapely-2.1.2-cp314-cp314-win_amd64.whl", hash = "sha256:a444e7afccdb0999e203b976adb37ea633725333e5b119ad40b1ca291ecf311c", size = 1773019, upload-time = "2025-09-24T13:51:24.873Z" }, + { url = "https://files.pythonhosted.org/packages/a3/05/a44f3f9f695fa3ada22786dc9da33c933da1cbc4bfe876fe3a100bafe263/shapely-2.1.2-cp314-cp314t-macosx_10_13_x86_64.whl", hash = "sha256:5ebe3f84c6112ad3d4632b1fd2290665aa75d4cef5f6c5d77c4c95b324527c6a", size = 1834137, upload-time = "2025-09-24T13:51:26.665Z" }, + { url = "https://files.pythonhosted.org/packages/52/7e/4d57db45bf314573427b0a70dfca15d912d108e6023f623947fa69f39b72/shapely-2.1.2-cp314-cp314t-macosx_11_0_arm64.whl", hash = "sha256:5860eb9f00a1d49ebb14e881f5caf6c2cf472c7fd38bd7f253bbd34f934eb076", size = 1642884, upload-time = "2025-09-24T13:51:28.029Z" }, + { url = "https://files.pythonhosted.org/packages/5a/27/4e29c0a55d6d14ad7422bf86995d7ff3f54af0eba59617eb95caf84b9680/shapely-2.1.2-cp314-cp314t-manylinux2014_aarch64.manylinux_2_17_aarch64.whl", hash = "sha256:b705c99c76695702656327b819c9660768ec33f5ce01fa32b2af62b56ba400a1", size = 3018320, upload-time = "2025-09-24T13:51:29.903Z" }, + { url = "https://files.pythonhosted.org/packages/9f/bb/992e6a3c463f4d29d4cd6ab8963b75b1b1040199edbd72beada4af46bde5/shapely-2.1.2-cp314-cp314t-manylinux2014_x86_64.manylinux_2_17_x86_64.whl", hash = "sha256:a1fd0ea855b2cf7c9cddaf25543e914dd75af9de08785f20ca3085f2c9ca60b0", size = 3094931, upload-time = "2025-09-24T13:51:32.699Z" }, + { url = "https://files.pythonhosted.org/packages/9c/16/82e65e21070e473f0ed6451224ed9fa0be85033d17e0c6e7213a12f59d12/shapely-2.1.2-cp314-cp314t-musllinux_1_2_aarch64.whl", hash = "sha256:df90e2db118c3671a0754f38e36802db75fe0920d211a27481daf50a711fdf26", size = 4030406, upload-time = "2025-09-24T13:51:34.189Z" }, + { url = "https://files.pythonhosted.org/packages/7c/75/c24ed871c576d7e2b64b04b1fe3d075157f6eb54e59670d3f5ffb36e25c7/shapely-2.1.2-cp314-cp314t-musllinux_1_2_x86_64.whl", hash = "sha256:361b6d45030b4ac64ddd0a26046906c8202eb60d0f9f53085f5179f1d23021a0", size = 4169511, upload-time = "2025-09-24T13:51:36.297Z" }, + { url = "https://files.pythonhosted.org/packages/b1/f7/b3d1d6d18ebf55236eec1c681ce5e665742aab3c0b7b232720a7d43df7b6/shapely-2.1.2-cp314-cp314t-win32.whl", hash = "sha256:b54df60f1fbdecc8ebc2c5b11870461a6417b3d617f555e5033f1505d36e5735", size = 1602607, upload-time = "2025-09-24T13:51:37.757Z" }, + { url = "https://files.pythonhosted.org/packages/9a/f6/f09272a71976dfc138129b8faf435d064a811ae2f708cb147dccdf7aacdb/shapely-2.1.2-cp314-cp314t-win_amd64.whl", hash = "sha256:0036ac886e0923417932c2e6369b6c52e38e0ff5d9120b90eef5cd9a5fc5cae9", size = 1796682, upload-time = "2025-09-24T13:51:39.233Z" }, ] [[package]] diff --git a/whitelist.txt b/whitelist.txt index 3e02ca0..a12d616 100644 --- a/whitelist.txt +++ b/whitelist.txt @@ -1,6 +1,5 @@ 2d alg -alg algs arctan2 argmax @@ -14,6 +13,7 @@ changelog cleandoc concat const +coord coords crs dataframe @@ -38,6 +38,7 @@ geopandas geoseries getfullargspec getsource +gettempdir gpd hypot idx @@ -45,6 +46,7 @@ inplace landcover linalg linemerge +linestring loc multipolygon nodata @@ -59,6 +61,7 @@ sindex sklearn starmap tolist +topo typer unary utils