From 0d85e85b70babcaabcf4edbd72824c89192ad4fa Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 12 Apr 2024 15:42:36 +0200 Subject: [PATCH 01/33] add fucntion merge mask --- lidro/vectors/merge_vector.py | 42 +++++++++++++++++++++++++++++++ test/vectors/test_merge_vector.py | 21 ++++++++++++++++ 2 files changed, 63 insertions(+) create mode 100644 lidro/vectors/merge_vector.py create mode 100644 test/vectors/test_merge_vector.py diff --git a/lidro/vectors/merge_vector.py b/lidro/vectors/merge_vector.py new file mode 100644 index 00000000..64f735f8 --- /dev/null +++ b/lidro/vectors/merge_vector.py @@ -0,0 +1,42 @@ +# -*- coding: utf-8 -*- +""" Merge +""" +import os + +import geopandas as gpd +from shapely.ops import unary_union + + +def merge_geom(input_folder: str, output_folder: str, crs: str): + """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, + filter mask (keep only water's area > 1000 m²) and save it as a GeoJSON file. + + Args: + input_folder (str): folder who contains severals Masks Hydro + output_folder (str): output folder + crs (str): a pyproj CRS object used to create the output GeoJSON file + """ + # List for stocking all GeoDataFrame + polys = [] + + # Parcourir tous les fichiers dans le dossier + for fichier in os.listdir(input_folder): + if fichier.endswith(".GeoJSON"): + # Charger chaque fichier GeoJSON en tant que GeoDataFrame + geojson = gpd.read_file(os.path.join(input_folder, fichier)) + merge_geom = geojson["geometry"].unary_union + # Ajouter le GeoDataFrame à la liste + polys.append(merge_geom) + + # Union geometry + mergedPolys = unary_union(polys) + + geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) + + # keep only water's area > 1000 m² + filter_geometry = [geom for geom in geometry if geom.area > 1000] + + gdf = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) + + # save the result + gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py new file mode 100644 index 00000000..d8175a46 --- /dev/null +++ b/test/vectors/test_merge_vector.py @@ -0,0 +1,21 @@ +import os +import shutil +from pathlib import Path + +from lidro.vectors.merge_vector import merge_geom + +INPUT_FOLDER = Path("./tmp/main/main_lidro_default") +TMP_PATH = Path("./tmp/vectors/merge_vector") + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_merge_mask_default(): + output = "tmp/vectors/merge_vector/MaskHydro_merge.geojson" + merge_geom(INPUT_FOLDER, TMP_PATH, "epsg:2154") + + assert Path(output).exists() From 35b00c68386c7a19b275fd147b74442198c2bf98 Mon Sep 17 00:00:00 2001 From: mdupays Date: Tue, 23 Apr 2024 11:54:30 +0200 Subject: [PATCH 02/33] add function merge all mask hydro --- configs/configs_lidro.yaml | 2 +- lidro/main.py | 18 ++++-- lidro/vectors/check_rectify_geometry.py | 84 +++++++++++++++++++++++++ lidro/vectors/convert_to_vector.py | 15 +---- lidro/vectors/merge_vector.py | 33 ++++++++-- test/vectors/test_convert_to_vector.py | 2 +- 6 files changed, 127 insertions(+), 27 deletions(-) create mode 100644 lidro/vectors/check_rectify_geometry.py diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index ac84fc6d..4660c541 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -19,7 +19,7 @@ io: filter: # Classes to be considered as "non-water" - keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 65, 66, 67] + keep_classes: [0, 1, 2, 3, 6, 17, 65, 66, 67] hydra: output_subdir: null diff --git a/lidro/main.py b/lidro/main.py index f645152b..ccd11012 100644 --- a/lidro/main.py +++ b/lidro/main.py @@ -9,15 +9,18 @@ from pyproj import CRS from lidro.vectors.convert_to_vector import create_hydro_vector_mask +from lidro.vectors.merge_vector import merge_geom @hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") def main(config: DictConfig): - """Create a vector mask of hydro surfaces from the points classification of the input LAS/LAZ file, + """Step 1 : Create a vector mask of hydro surfaces from the points classification of the input LAS/LAZ file, and save it as a GeoJSON file. It can run either on a single file, or on each file of a folder + Step 2 : Merge all vector masks of hydro surfaces + Args: config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default) It contains the algorithm parameters and the input/output parameters @@ -53,9 +56,9 @@ def main_on_one_tile(filename): Args: filename (str): filename to the LAS file """ - tilename = os.path.splitext(filename)[0] # filename to the LAS file - input_file = os.path.join(input_dir, filename) # path to the LAS file - output_file = os.path.join(output_dir, f"MaskHydro_{tilename}.GeoJSON") # path to the Mask Hydro file + tilename = os.path.splitext(filename)[0] # filename to the LAS file + input_file = os.path.join(input_dir, filename) # path to the LAS file + output_file = os.path.join(output_dir, f"MaskHydro_{tilename}.GeoJSON") # path to the Mask Hydro file logging.info(f"\nCreate Mask Hydro 1 for tile : {tilename}") create_hydro_vector_mask(input_file, output_file, pixel_size, tile_size, classe, crs) @@ -64,9 +67,14 @@ def main_on_one_tile(filename): main_on_one_tile(initial_las_filename) else: - # Lauch creating mask hydro tile by tile + # Lauch creating Mask Hydro tile by tile for file in os.listdir(input_dir): main_on_one_tile(file) + if os.path.isdir(output_dir): + output_dir_merge = os.path.join(output_dir, "merge") # folder with Mask Hydro file + os.makedirs(output_dir_merge, exist_ok=True) # Create folder "merge" + # Merge all Mash Hydro + merge_geom(output_dir, output_dir_merge, crs) if __name__ == "__main__": diff --git a/lidro/vectors/check_rectify_geometry.py b/lidro/vectors/check_rectify_geometry.py new file mode 100644 index 00000000..9cc9343c --- /dev/null +++ b/lidro/vectors/check_rectify_geometry.py @@ -0,0 +1,84 @@ +# -*- coding: utf-8 -*- +""" Rectify and check geometry +""" +from shapely.geometry import CAP_STYLE, MultiPolygon, Polygon +from shapely.validation import make_valid + + +def remove_hole(multipoly): + """Remove small holes (surface < 50 m²) + + Args: + - multipoly (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry without holes (< 50 m²) + """ + list_parts = [] + eps = 50 + + for polygon in multipoly.geoms: + list_interiors = [] + + for interior in polygon.interiors: + p = Polygon(interior) + + if p.area > eps: + list_interiors.append(interior) + + temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) + list_parts.append(temp_pol) + + new_multipolygon = MultiPolygon(list_parts) + + return new_multipolygon + + +def simplify_geometry(polygon): + """Simplify certain shapes to facilitate calculations + + Args: + polygon (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry simplify + """ + # Positive buffer half the size of the HYDRO mask resolution (1m resolution) + _geom = polygon.buffer(0.5, cap_style=CAP_STYLE.square) + geom = _geom.buffer(-1, cap_style=CAP_STYLE.square) + return geom + + +def fix_invalid_geometry(geometry): + """Set invalid geoemtries + + Args: + geometry (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry valid + """ + + if not geometry.is_valid: + return make_valid(geometry) + else: + return geometry + + +def check_geometry(initial_gdf): + """Check topology + + Args: + initial_gdf (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry valid + """ + # Obtain simple geometries + gdf_simple = initial_gdf.explode(ignore_index=True) + # Delete duplicate geoemtries if any + gdf_without_duplicates = gdf_simple.drop_duplicates(ignore_index=True) + # Identify invalid geometries and keep only valid ones + gdf_valid = gdf_without_duplicates.copy() + gdf_valid.geometry = gdf_valid.geometry.apply(lambda geom: fix_invalid_geometry(geom)) + return gdf_valid diff --git a/lidro/vectors/convert_to_vector.py b/lidro/vectors/convert_to_vector.py index b3a7b065..993982ac 100644 --- a/lidro/vectors/convert_to_vector.py +++ b/lidro/vectors/convert_to_vector.py @@ -40,19 +40,6 @@ def create_hydro_vector_mask(filename: str, output: str, pixel_size: float, tile ) if value != 0 ] - # keep only water's area > 1000 m² - filter_geometry = [geom for geom in geometry if geom.area > 1000] - nb_filter_geometry = len(filter_geometry) - - gdf = gpd.GeoDataFrame( - { - "layer": [ii for ii, geom in enumerate(filter_geometry)] * np.ones(nb_filter_geometry), - "dalles": filename, - "geometry": filter_geometry, - }, - geometry="geometry", - crs=crs, - ) - # save the result + gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs) gdf.to_file(output, driver="GeoJSON", crs=crs) diff --git a/lidro/vectors/merge_vector.py b/lidro/vectors/merge_vector.py index 64f735f8..067e42f0 100644 --- a/lidro/vectors/merge_vector.py +++ b/lidro/vectors/merge_vector.py @@ -4,8 +4,15 @@ import os import geopandas as gpd +from shapely.geometry import MultiPolygon from shapely.ops import unary_union +from lidro.vectors.check_rectify_geometry import ( + check_geometry, + remove_hole, + simplify_geometry, +) + def merge_geom(input_folder: str, output_folder: str, crs: str): """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, @@ -19,13 +26,13 @@ def merge_geom(input_folder: str, output_folder: str, crs: str): # List for stocking all GeoDataFrame polys = [] - # Parcourir tous les fichiers dans le dossier + # Browse all files in folder for fichier in os.listdir(input_folder): if fichier.endswith(".GeoJSON"): - # Charger chaque fichier GeoJSON en tant que GeoDataFrame + # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(os.path.join(input_folder, fichier)) merge_geom = geojson["geometry"].unary_union - # Ajouter le GeoDataFrame à la liste + # Add the GeoDataFrame to the list polys.append(merge_geom) # Union geometry @@ -33,10 +40,24 @@ def merge_geom(input_folder: str, output_folder: str, crs: str): geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) - # keep only water's area > 1000 m² - filter_geometry = [geom for geom in geometry if geom.area > 1000] + # keep only water's area > 100 m² + filter_geometry = [geom for geom in geometry if geom.area > 100] + gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) + + # Check and rectify the invalid geometry + gdf = check_geometry(gdf_filter) + + # Correcting geometric errors: simplifying certain shapes to make calculations easier + buffer_geom = gdf["geometry"].apply(simplify_geometry) + simplify_geom = buffer_geom.simplify(tolerance=0.5, preserve_topology=True) - gdf = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) + # Correction of holes (< 50m²) in Hydrological Masks + list_parts = [] + for poly in simplify_geom: + list_parts.append(poly) + not_hole_geom = remove_hole(MultiPolygon(list_parts)) + geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False) + gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs) # save the result gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/test/vectors/test_convert_to_vector.py b/test/vectors/test_convert_to_vector.py index b481db03..882f264e 100644 --- a/test/vectors/test_convert_to_vector.py +++ b/test/vectors/test_convert_to_vector.py @@ -29,7 +29,7 @@ def test_create_hydro_vector_mask_default(): AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2154"]]' Path(output).parent.mkdir(exist_ok=True) - create_hydro_vector_mask(las_file, output, 1, 1000, [0, 1, 2, 3, 4, 5, 6, 17, 66], crs) + create_hydro_vector_mask(las_file, output, 1, 1000, [0, 1, 2, 3, 6, 17, 66], crs) assert Path(output).exists() From 7a512775bccdf43494d9f017b64e167f1f73007b Mon Sep 17 00:00:00 2001 From: mdupays Date: Tue, 23 Apr 2024 17:51:16 +0200 Subject: [PATCH 03/33] update merge mask with other parameters v2 --- configs/configs_lidro.yaml | 2 +- lidro/vectors/check_rectify_geometry.py | 7 ++++--- lidro/vectors/merge_vector.py | 8 ++++---- test/vectors/test_convert_to_vector.py | 2 +- 4 files changed, 10 insertions(+), 9 deletions(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 4660c541..ac84fc6d 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -19,7 +19,7 @@ io: filter: # Classes to be considered as "non-water" - keep_classes: [0, 1, 2, 3, 6, 17, 65, 66, 67] + keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 65, 66, 67] hydra: output_subdir: null diff --git a/lidro/vectors/check_rectify_geometry.py b/lidro/vectors/check_rectify_geometry.py index 9cc9343c..48134298 100644 --- a/lidro/vectors/check_rectify_geometry.py +++ b/lidro/vectors/check_rectify_geometry.py @@ -12,10 +12,10 @@ def remove_hole(multipoly): - multipoly (GeoJSON): Hydro Mask geometry Returns: - GeoJSON: Hydro Mask geometry without holes (< 50 m²) + GeoJSON: Hydro Mask geometry without holes (< 100 m²) """ list_parts = [] - eps = 50 + eps = 100 for polygon in multipoly.geoms: list_interiors = [] @@ -43,7 +43,8 @@ def simplify_geometry(polygon): Returns: GeoJSON: Hydro Mask geometry simplify """ - # Positive buffer half the size of the HYDRO mask resolution (1m resolution) + # Positive buffer half the size of the Hydro Mask resolution (0.5 m resolution) + # + Negative buffer the size of the Hydro Mask Resolution (1m) _geom = polygon.buffer(0.5, cap_style=CAP_STYLE.square) geom = _geom.buffer(-1, cap_style=CAP_STYLE.square) return geom diff --git a/lidro/vectors/merge_vector.py b/lidro/vectors/merge_vector.py index 067e42f0..24538bcf 100644 --- a/lidro/vectors/merge_vector.py +++ b/lidro/vectors/merge_vector.py @@ -16,7 +16,7 @@ def merge_geom(input_folder: str, output_folder: str, crs: str): """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, - filter mask (keep only water's area > 1000 m²) and save it as a GeoJSON file. + filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file. Args: input_folder (str): folder who contains severals Masks Hydro @@ -40,8 +40,8 @@ def merge_geom(input_folder: str, output_folder: str, crs: str): geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) - # keep only water's area > 100 m² - filter_geometry = [geom for geom in geometry if geom.area > 100] + # keep only water's area > 150 m² + filter_geometry = [geom for geom in geometry if geom.area > 150] gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) # Check and rectify the invalid geometry @@ -51,7 +51,7 @@ def merge_geom(input_folder: str, output_folder: str, crs: str): buffer_geom = gdf["geometry"].apply(simplify_geometry) simplify_geom = buffer_geom.simplify(tolerance=0.5, preserve_topology=True) - # Correction of holes (< 50m²) in Hydrological Masks + # Correction of holes (< 100m²) in Hydrological Masks list_parts = [] for poly in simplify_geom: list_parts.append(poly) diff --git a/test/vectors/test_convert_to_vector.py b/test/vectors/test_convert_to_vector.py index 882f264e..9ef8c849 100644 --- a/test/vectors/test_convert_to_vector.py +++ b/test/vectors/test_convert_to_vector.py @@ -36,7 +36,7 @@ def test_create_hydro_vector_mask_default(): with open(output, "r") as f: geojson_data = json.load(f) - # Vérifie la structure globale du fichier GeoJSON + # CHECK STRUCTURE assert "type" in geojson_data assert geojson_data["type"] == "FeatureCollection" assert "features" in geojson_data From 1919c4b65f36b32975f242245ce2ec70299a27da Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 26 Apr 2024 17:06:41 +0200 Subject: [PATCH 04/33] MAJ from branch fix-refacto-mask --- Dockerfile | 16 ++++ Makefile | 36 +++++++- README.md | 50 +++++++---- ci/test.sh | 1 + configs/configs_lidro.yaml | 4 + ..._tile.sh => example_create_mask_by_tile.sh | 2 +- ...fault.sh => example_create_mask_default.sh | 2 +- .../pointcloud/filter_las.py | 0 .../{ => create_mask_hydro}/pointcloud/io.py | 0 .../pointcloud/read_las.py | 0 .../rasters/create_mask_raster.py | 15 ++-- .../vectors/convert_to_vector.py | 7 +- lidro/{main.py => main_create_mask.py} | 16 +--- lidro/vectors/check_rectify_geometry.py | 85 ------------------- lidro/vectors/merge_vector.py | 63 -------------- test/pointcloud/test_filter_las.py | 4 +- test/pointcloud/test_io.py | 4 +- test/pointcloud/test_read_las.py | 4 +- test/rasters/test_create_mask_raster.py | 6 +- ...{test_main.py => test_main_create_mask.py} | 10 +-- test/vectors/test_convert_to_vector.py | 37 ++++---- test/vectors/test_merge_vector.py | 21 ----- 22 files changed, 145 insertions(+), 238 deletions(-) create mode 100644 Dockerfile create mode 100644 ci/test.sh rename example_lidro_by_tile.sh => example_create_mask_by_tile.sh (81%) rename example_lidro_default.sh => example_create_mask_default.sh (71%) rename lidro/{ => create_mask_hydro}/pointcloud/filter_las.py (100%) rename lidro/{ => create_mask_hydro}/pointcloud/io.py (100%) rename lidro/{ => create_mask_hydro}/pointcloud/read_las.py (100%) rename lidro/{ => create_mask_hydro}/rasters/create_mask_raster.py (80%) rename lidro/{ => create_mask_hydro}/vectors/convert_to_vector.py (86%) rename lidro/{main.py => main_create_mask.py} (79%) delete mode 100644 lidro/vectors/check_rectify_geometry.py delete mode 100644 lidro/vectors/merge_vector.py rename test/{test_main.py => test_main_create_mask.py} (93%) delete mode 100644 test/vectors/test_merge_vector.py diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 00000000..e69d4f0a --- /dev/null +++ b/Dockerfile @@ -0,0 +1,16 @@ +FROM mambaorg/micromamba:latest as mamba_pdal +COPY environment.yml /environment.yml +USER root +RUN micromamba env create -n lidro -f /environment.yml + +FROM debian:bullseye-slim + + +WORKDIR /lidro +RUN mkdir tmp +COPY lidro lidro +COPY test test +COPY configs configs + +# Copy test data that are stored directly in the lidro-data repository ("http://gitlab.forge-idi.ign.fr/Lidar/lidro-data.git") +COPY data/pointcloud data/pointcloud \ No newline at end of file diff --git a/Makefile b/Makefile index 6e6d0c54..001ed431 100644 --- a/Makefile +++ b/Makefile @@ -1,6 +1,40 @@ +# Makefile to manage main tasks +# cf. https://blog.ianpreston.ca/conda/python/bash/2020/05/13/conda_envs.html#makefile +# Oneshell means I can run multiple lines in a recipe in the same shell, so I don't have to +# chain commands together with semicolon +.ONESHELL: +SHELL = /bin/bash install: - mamba env update -n lidro -f environment.yml + pip install -e . install-precommit: pre-commit install + +testing: + ./ci/test.sh + +mamba-env-create: + mamba env create -n lidro -f environment.yml + +mamba-env-update: + mamba env update -n lidro -f environment.yml + +############################## +# Docker +############################## + +PROJECT_NAME=ignimagelidar/lidro +VERSION=`python -m lidro.version` + +docker-build: + docker build -t ${PROJECT_NAME}:${VERSION} -f Dockerfile . + +docker-test: + docker run --rm -it ${PROJECT_NAME}:${VERSION} python -m pytest -s + +docker-remove: + docker rmi -f `docker images | grep ${PROJECT_NAME} | tr -s ' ' | cut -d ' ' -f 3` + +docker-deploy: + docker push ${PROJECT_NAME}:${VERSION} \ No newline at end of file diff --git a/README.md b/README.md index 9ea13653..b9724a1d 100644 --- a/README.md +++ b/README.md @@ -1,30 +1,47 @@ -### LIDRO ### -Lidro (Aplanissement des surfaces d'eaux) est un outil permettant de créer automatiquement des points virtuels le long des surfaces d'eaux afin de créer des modèles numériques cohérents avec les modèles hydrologiques. Le jeu de données en entrée correspond à un nuage de point lidar classifiés. +### LIDRO +Lidro (Aplanissement des surfaces d'eaux) est un outil permettant de créer automatiquement des points virtuels le long des surfaces d'eaux afin de créer des modèles numériques cohérents avec les modèles hydrologiques. Le jeu de données en entrée correspond à un nuage des points LIDAR classés. ## Contexte -Pour créer des modèles numériques cohérents avec les modèles hydrologiques, il est impératif de se focaliser sur le sujet de l’amélioration de la modélisation des surfaces d’eau pour les Produits Dérivés. ​ +Pour créer des modèles numériques cohérents avec les modèles hydrologiques, il est impératif de se focaliser sur l’amélioration de la modélisation des surfaces d’eau. ​ Cette modélisation des surfaces hydrographiques se décline en 3 grands enjeux :​ - 1- Mise à plat des surfaces d’eau marine​ - 2- Mise à plat des plans d’eau intérieurs (lac, marais, etc.)​ - 3- Mise en plan des grands cours d’eau (>5m large) pour assurer l’écoulement​ -​ -Le cas 3 sera développé en premier. +* Mise à plat des surfaces d’eau marine​ +* Mise à plat des plans d’eau intérieurs (lac, marais, etc.)​ +* Mise en plan des grands cours d’eau (>5m large) pour assurer l’écoulement​. A noter que cette étape sera développée en premier. ## Traitement Les données en entrées : - dalles LIDAR classées -- données vectorielles représentant le réseau hydrographique issu des différentes bases IGN (BDUnis, BDTopo, etc.) +- données vectorielles représentant le réseau hydrographique issu des différentes bases de données IGN (BDUnis, BDTopo, etc.) Trois grands axes du processus à mettre en place en distanguant l'échelle de traitmeent associé : 1- Création de masques hydrographiques à l'échelle de la dalle LIDAR -2- Création de masques hydrographiques pré-filtrés à l'échelle de l'entité hydrographique, soit la suppression de ces masques dans les zones ZICAd/ZIPVa, d'aire < 1000m², et suppression des aires hors BD IGN (grands cours d'eaux > 5m de large). -3.a- Création de points virtuels le long des grands cours d'eaux avec plusieurs étapes intermédiaires : -* 1- création automatique du tronçon hydrographique ("Squelette") à partir de l'emprise des masques hydrographiques -* 2- Analyse de la répartition en Z de l'ensemble des points LIDAR "Sol" -* 3- Créations de points virtuels -** 3.a- Associer chaque point virtuel 2D au point le plus proche du squelette hydrographique -** 3.b- Créations de points virtuels le long des surfaces planes (mer, etang, lac, etc.): analyse statisque de l'ensmeble des points LIDAR "Sol" le long des côtes/berges afin d'obtenir une surface plane +2- Création de masques hydrographiques pré-filtrés à l'échelle du chantier, soit : +* la suppression de ces masques dans les zones ZICAD/ZIPVA +* la suppression des aires < 150 m² +* la suppression des aires < 1000 m² hors BD IGN (grands cours d'eaux > 5m de large) + +3- Création de points virtuels le long de deux entités hydrographiques : +* ##### Grands cours d'eaux (> 5 m de large dans la BD Unis). +Il existe plusieurs étapes intermédiaires : +1- création automatique du tronçon hydrographique ("Squelette hydrographique", soit les tronçons hydrographiques dans la BD Unid) à partir de l'emprise du masque hydrographique "écoulement" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle) +2- Analyse de la répartition en Z de l'ensemble des points LIDAR "Sol" +3- Créations de points virtuels nécéssitant plusieurs étapes intermédiaires : +* Associer chaque point virtuel 2D au point le plus proche du squelette hydrographique +* Traitement "Z" du squelette : + * Analyser la répartition en Z de l’ensemble des points LIDAR extrait à l’étape précédente afin d’extraire une seule valeur d’altitude au point fictif Z le plus proche du squelette hydrographique. A noter que l’altitude correspond donc à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes. Pour les ponts, une étape de contrôle de la classification pourrait être mise en place + * Lisser les Z le long du squelette HYDRO pour assurer l'écoulement +* Création des points virtuels 2D tous les 0.5 mètres le long des bords du masque hydrographique "écoulement" en cohérence en "Z" avec le squelette hydrographique + +* ##### Surfaces planes (mer, lac, étang, etc.) +Pour rappel, l'eau est considérée comme horizontale sur ce type de surface. +Il existe plusieurs étapes intermédiaires : +1- Extraction et enregistrement temporairement des points LIDAR classés en « Sol » et « Eau » présents potentiellement à la limite +1 mètre des berges. Pour cela, on s'appuie sur 'emprise du masque hydrographique "surface plane" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle). a noter que pur le secteur maritime (Mer), il faut exclure la classe « 9 » (eau) afin d’éviter de mesurer les vagues. +2- Analyse statisque de l'ensemble des points LIDAR "Sol / Eau" le long des côtes/berges afin d'obtenir une surface plane. L’objectif est de créer des points virtuels spécifiques avec une information d'altitude (m) tous les 0.5 m sur les bords des surfaces maritimes et des plans d’eaux à partir du masque hydrographique "surface plane". Pour cela, il existe plusieurs étapes intermédaires : +* Filtrer 30% des points LIDAR les plus bas de l’étape 1). afin de supprimer les altitudes trop élevées +* Analyser la répartition en Z de ces points afin d’extraire une seule valeur d’altitude selon l’objet hydrographique : + * Pour les Plans d’eau : l’altitude correspond à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes, + * Pour la Mer : l’altitude correspond à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes ## Installation des dépendances (conda) @@ -47,6 +64,7 @@ mamba env update -n lidro -f environment.yml conda activate lidro ``` +## Contribuer Installer pre-commit ``` pre-commit install diff --git a/ci/test.sh b/ci/test.sh new file mode 100644 index 00000000..823d3d84 --- /dev/null +++ b/ci/test.sh @@ -0,0 +1 @@ +python -m pytest -s ./test -v \ No newline at end of file diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index ac84fc6d..c03cc87d 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -17,6 +17,10 @@ io: no_data_value: -9999 tile_size: 1000 +raster: + # size for dilatation + dilatation_size: 3 + filter: # Classes to be considered as "non-water" keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 65, 66, 67] diff --git a/example_lidro_by_tile.sh b/example_create_mask_by_tile.sh similarity index 81% rename from example_lidro_by_tile.sh rename to example_create_mask_by_tile.sh index f9c4de70..fd2dd27b 100644 --- a/example_lidro_by_tile.sh +++ b/example_create_mask_by_tile.sh @@ -1,5 +1,5 @@ # For lauching Mask Hydro -python -m lidro.main \ +python -m lidro.main_create_mask \ io.input_filename=LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.las \ io.input_dir=./data/pointcloud/ \ io.output_dir=./tmp/ \ diff --git a/example_lidro_default.sh b/example_create_mask_default.sh similarity index 71% rename from example_lidro_default.sh rename to example_create_mask_default.sh index 3b639bee..8818c5b4 100644 --- a/example_lidro_default.sh +++ b/example_create_mask_default.sh @@ -1,5 +1,5 @@ # For lauching Mask Hydro -python -m lidro.main \ +python -m lidro.main_create_mask \ io.input_dir=./data/pointcloud/ \ io.output_dir=./tmp/ \ diff --git a/lidro/pointcloud/filter_las.py b/lidro/create_mask_hydro/pointcloud/filter_las.py similarity index 100% rename from lidro/pointcloud/filter_las.py rename to lidro/create_mask_hydro/pointcloud/filter_las.py diff --git a/lidro/pointcloud/io.py b/lidro/create_mask_hydro/pointcloud/io.py similarity index 100% rename from lidro/pointcloud/io.py rename to lidro/create_mask_hydro/pointcloud/io.py diff --git a/lidro/pointcloud/read_las.py b/lidro/create_mask_hydro/pointcloud/read_las.py similarity index 100% rename from lidro/pointcloud/read_las.py rename to lidro/create_mask_hydro/pointcloud/read_las.py diff --git a/lidro/rasters/create_mask_raster.py b/lidro/create_mask_hydro/rasters/create_mask_raster.py similarity index 80% rename from lidro/rasters/create_mask_raster.py rename to lidro/create_mask_hydro/rasters/create_mask_raster.py index 3174ee1a..797119e4 100644 --- a/lidro/rasters/create_mask_raster.py +++ b/lidro/create_mask_hydro/rasters/create_mask_raster.py @@ -6,9 +6,9 @@ import numpy as np import scipy.ndimage -from lidro.pointcloud.filter_las import filter_pointcloud -from lidro.pointcloud.io import get_pointcloud_origin -from lidro.pointcloud.read_las import read_pointcloud +from lidro.create_mask_hydro.pointcloud.filter_las import filter_pointcloud +from lidro.create_mask_hydro.pointcloud.io import get_pointcloud_origin +from lidro.create_mask_hydro.pointcloud.read_las import read_pointcloud def create_occupancy_map(points: np.array, tile_size: int, pixel_size: float, origin: Tuple[int, int]): @@ -33,7 +33,7 @@ def create_occupancy_map(points: np.array, tile_size: int, pixel_size: float, or return bins -def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, classes: List[int]): +def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, classes: List[int], dilatation_size: int): """ "Detect hydrographic surfaces in a tile from the classified points of the input pointcloud An hydrographic surface is define as a surface where there is no points from any class different from water @@ -43,6 +43,7 @@ def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, class pixel_size (float): distance between each node of the raster grid (in meters) classes (List[int]): List of classes to use for the binarisation (points with other classification values are ignored) + dilatation_size (int): size for dilatation raster Returns: smoothed_water (np.array): 2D binary array (x, y) of the water presence from the point cloud @@ -64,7 +65,9 @@ def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, class # Revert occupancy map to keep pixels where there is no point of the selected classes detected_water = np.logical_not(occupancy) - # Apply a mathematical morphology operations: dilation (/ ! \ closing : reduct size of image) - morphology_bins = scipy.ndimage.binary_dilation(detected_water, structure=np.ones((3, 3))).astype(np.uint8) + # Apply a mathematical morphology operations: DILATATION + # / ! \ NOT "CLOSING", due to the reduction in the size of hydro masks, particularly at the tile borders. + # / ! \ WITH "CLOSING" => Masks Hydro are no longer continuous, when they are merged + morphology_bins = scipy.ndimage.binary_dilation(detected_water, structure=np.ones((dilatation_size, dilatation_size))).astype(np.uint8) return morphology_bins, pcd_origin diff --git a/lidro/vectors/convert_to_vector.py b/lidro/create_mask_hydro/vectors/convert_to_vector.py similarity index 86% rename from lidro/vectors/convert_to_vector.py rename to lidro/create_mask_hydro/vectors/convert_to_vector.py index 993982ac..3b0657d5 100644 --- a/lidro/vectors/convert_to_vector.py +++ b/lidro/create_mask_hydro/vectors/convert_to_vector.py @@ -7,10 +7,10 @@ from rasterio.transform import from_origin from shapely.geometry import shape as shapely_shape -from lidro.rasters.create_mask_raster import detect_hydro_by_tile +from lidro.create_mask_hydro.rasters.create_mask_raster import detect_hydro_by_tile -def create_hydro_vector_mask(filename: str, output: str, pixel_size: float, tile_size: int, classes: list, crs: str): +def create_hydro_vector_mask(filename: str, output: str, pixel_size: float, tile_size: int, classes: list, crs: str, dilatation_size: int): """Create a vector mask of hydro surfaces in a tile from the points classification of the input LAS/LAZ file, and save it as a GeoJSON file. @@ -22,9 +22,10 @@ def create_hydro_vector_mask(filename: str, output: str, pixel_size: float, tile tile_size (int): size of the intermediate raster grid (in meters) classes (list): List of classes to consider as water (points with other classification values are ignored) crs (str): a pyproj CRS object used to create the output GeoJSON file + dilatation_size (int): size for dilatation raster """ # Read a binary image representing hydrographic surface(s) - binary_image, pcd_origin = detect_hydro_by_tile(filename, tile_size, pixel_size, classes) + binary_image, pcd_origin = detect_hydro_by_tile(filename, tile_size, pixel_size, classes, dilatation_size) # Extract origin origin_x = pcd_origin[0] diff --git a/lidro/main.py b/lidro/main_create_mask.py similarity index 79% rename from lidro/main.py rename to lidro/main_create_mask.py index ccd11012..84313aa8 100644 --- a/lidro/main.py +++ b/lidro/main_create_mask.py @@ -8,19 +8,16 @@ from omegaconf import DictConfig from pyproj import CRS -from lidro.vectors.convert_to_vector import create_hydro_vector_mask -from lidro.vectors.merge_vector import merge_geom +from lidro.create_mask_hydro.vectors.convert_to_vector import create_hydro_vector_mask @hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") def main(config: DictConfig): - """Step 1 : Create a vector mask of hydro surfaces from the points classification of the input LAS/LAZ file, + """Create a vector mask of hydro surfaces from the points classification of the input LAS/LAZ file, and save it as a GeoJSON file. It can run either on a single file, or on each file of a folder - Step 2 : Merge all vector masks of hydro surfaces - Args: config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default) It contains the algorithm parameters and the input/output parameters @@ -49,6 +46,7 @@ def main(config: DictConfig): tile_size = config.io.tile_size crs = CRS.from_user_input(config.io.srid) classe = config.filter.keep_classes + dilatation_size = config.raster.dilatation_size def main_on_one_tile(filename): """Lauch main.py on one tile @@ -60,7 +58,7 @@ def main_on_one_tile(filename): input_file = os.path.join(input_dir, filename) # path to the LAS file output_file = os.path.join(output_dir, f"MaskHydro_{tilename}.GeoJSON") # path to the Mask Hydro file logging.info(f"\nCreate Mask Hydro 1 for tile : {tilename}") - create_hydro_vector_mask(input_file, output_file, pixel_size, tile_size, classe, crs) + create_hydro_vector_mask(input_file, output_file, pixel_size, tile_size, classe, crs, dilatation_size) if initial_las_filename: # Lauch creating mask by one tile: @@ -70,12 +68,6 @@ def main_on_one_tile(filename): # Lauch creating Mask Hydro tile by tile for file in os.listdir(input_dir): main_on_one_tile(file) - if os.path.isdir(output_dir): - output_dir_merge = os.path.join(output_dir, "merge") # folder with Mask Hydro file - os.makedirs(output_dir_merge, exist_ok=True) # Create folder "merge" - # Merge all Mash Hydro - merge_geom(output_dir, output_dir_merge, crs) - if __name__ == "__main__": main() diff --git a/lidro/vectors/check_rectify_geometry.py b/lidro/vectors/check_rectify_geometry.py deleted file mode 100644 index 48134298..00000000 --- a/lidro/vectors/check_rectify_geometry.py +++ /dev/null @@ -1,85 +0,0 @@ -# -*- coding: utf-8 -*- -""" Rectify and check geometry -""" -from shapely.geometry import CAP_STYLE, MultiPolygon, Polygon -from shapely.validation import make_valid - - -def remove_hole(multipoly): - """Remove small holes (surface < 50 m²) - - Args: - - multipoly (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry without holes (< 100 m²) - """ - list_parts = [] - eps = 100 - - for polygon in multipoly.geoms: - list_interiors = [] - - for interior in polygon.interiors: - p = Polygon(interior) - - if p.area > eps: - list_interiors.append(interior) - - temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) - list_parts.append(temp_pol) - - new_multipolygon = MultiPolygon(list_parts) - - return new_multipolygon - - -def simplify_geometry(polygon): - """Simplify certain shapes to facilitate calculations - - Args: - polygon (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry simplify - """ - # Positive buffer half the size of the Hydro Mask resolution (0.5 m resolution) - # + Negative buffer the size of the Hydro Mask Resolution (1m) - _geom = polygon.buffer(0.5, cap_style=CAP_STYLE.square) - geom = _geom.buffer(-1, cap_style=CAP_STYLE.square) - return geom - - -def fix_invalid_geometry(geometry): - """Set invalid geoemtries - - Args: - geometry (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry valid - """ - - if not geometry.is_valid: - return make_valid(geometry) - else: - return geometry - - -def check_geometry(initial_gdf): - """Check topology - - Args: - initial_gdf (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry valid - """ - # Obtain simple geometries - gdf_simple = initial_gdf.explode(ignore_index=True) - # Delete duplicate geoemtries if any - gdf_without_duplicates = gdf_simple.drop_duplicates(ignore_index=True) - # Identify invalid geometries and keep only valid ones - gdf_valid = gdf_without_duplicates.copy() - gdf_valid.geometry = gdf_valid.geometry.apply(lambda geom: fix_invalid_geometry(geom)) - return gdf_valid diff --git a/lidro/vectors/merge_vector.py b/lidro/vectors/merge_vector.py deleted file mode 100644 index 24538bcf..00000000 --- a/lidro/vectors/merge_vector.py +++ /dev/null @@ -1,63 +0,0 @@ -# -*- coding: utf-8 -*- -""" Merge -""" -import os - -import geopandas as gpd -from shapely.geometry import MultiPolygon -from shapely.ops import unary_union - -from lidro.vectors.check_rectify_geometry import ( - check_geometry, - remove_hole, - simplify_geometry, -) - - -def merge_geom(input_folder: str, output_folder: str, crs: str): - """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, - filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file. - - Args: - input_folder (str): folder who contains severals Masks Hydro - output_folder (str): output folder - crs (str): a pyproj CRS object used to create the output GeoJSON file - """ - # List for stocking all GeoDataFrame - polys = [] - - # Browse all files in folder - for fichier in os.listdir(input_folder): - if fichier.endswith(".GeoJSON"): - # Load each GeoJSON file as GeoDataFrame - geojson = gpd.read_file(os.path.join(input_folder, fichier)) - merge_geom = geojson["geometry"].unary_union - # Add the GeoDataFrame to the list - polys.append(merge_geom) - - # Union geometry - mergedPolys = unary_union(polys) - - geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) - - # keep only water's area > 150 m² - filter_geometry = [geom for geom in geometry if geom.area > 150] - gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) - - # Check and rectify the invalid geometry - gdf = check_geometry(gdf_filter) - - # Correcting geometric errors: simplifying certain shapes to make calculations easier - buffer_geom = gdf["geometry"].apply(simplify_geometry) - simplify_geom = buffer_geom.simplify(tolerance=0.5, preserve_topology=True) - - # Correction of holes (< 100m²) in Hydrological Masks - list_parts = [] - for poly in simplify_geom: - list_parts.append(poly) - not_hole_geom = remove_hole(MultiPolygon(list_parts)) - geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False) - gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs) - - # save the result - gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/test/pointcloud/test_filter_las.py b/test/pointcloud/test_filter_las.py index 12e50d2f..15ce3bc7 100644 --- a/test/pointcloud/test_filter_las.py +++ b/test/pointcloud/test_filter_las.py @@ -4,9 +4,9 @@ import numpy as np -from lidro.pointcloud.filter_las import filter_pointcloud +from lidro.create_mask_hydro.pointcloud.filter_las import filter_pointcloud -TMP_PATH = Path("./tmp/pointcloud/filter_las") +TMP_PATH = Path("./tmp/create_mask_hydro/pointcloud/filter_las") def setup_module(module): diff --git a/test/pointcloud/test_io.py b/test/pointcloud/test_io.py index 727527f5..c8e79ef0 100644 --- a/test/pointcloud/test_io.py +++ b/test/pointcloud/test_io.py @@ -5,9 +5,9 @@ import laspy import numpy as np -from lidro.pointcloud.io import get_pointcloud_origin +from lidro.create_mask_hydro.pointcloud.io import get_pointcloud_origin -TMP_PATH = Path("./tmp/pointcloud/io") +TMP_PATH = Path("./tmp/create_mask_hydro/pointcloud/io") LAS_FILE = "./data/pointcloud/LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.las" diff --git a/test/pointcloud/test_read_las.py b/test/pointcloud/test_read_las.py index 92e0d2dc..8f3748b4 100644 --- a/test/pointcloud/test_read_las.py +++ b/test/pointcloud/test_read_las.py @@ -5,9 +5,9 @@ import numpy as np from pyproj import CRS -from lidro.pointcloud.read_las import read_pointcloud +from lidro.create_mask_hydro.pointcloud.read_las import read_pointcloud -TMP_PATH = Path("./tmp/pointcloud/io") +TMP_PATH = Path("./tmp/create_mask_hydro/pointcloud/io") LAS_FILE = "./data/pointcloud/LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.las" diff --git a/test/rasters/test_create_mask_raster.py b/test/rasters/test_create_mask_raster.py index ee5271b6..33865095 100644 --- a/test/rasters/test_create_mask_raster.py +++ b/test/rasters/test_create_mask_raster.py @@ -7,9 +7,9 @@ import rasterio from rasterio.transform import from_origin -from lidro.rasters.create_mask_raster import create_occupancy_map, detect_hydro_by_tile +from lidro.create_mask_hydro.rasters.create_mask_raster import create_occupancy_map, detect_hydro_by_tile -TMP_PATH = Path("./tmp/rasters/create_mask_raster") +TMP_PATH = Path("./tmp/create_mask_hydro/rasters/create_mask_raster") LAS_FILE = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" @@ -43,7 +43,7 @@ def test_create_occupancy_map_default(): @pytest.mark.returnfile def test_detect_hydro_by_tile_return_file(): output_tif = TMP_PATH / "Semis_2021_0830_6291_LA93_IGN69_size.tif" - array, origin = detect_hydro_by_tile(LAS_FILE, tile_size, pixel_size, classes=[0, 1, 2, 3, 4, 5, 6, 17, 66]) + array, origin = detect_hydro_by_tile(LAS_FILE, tile_size, pixel_size, classes=[0, 1, 2, 3, 4, 5, 6, 17, 66], dilatation_size=1) assert isinstance(array, np.ndarray) is True assert list(array.shape) == [tile_size / pixel_size] * 2 diff --git a/test/test_main.py b/test/test_main_create_mask.py similarity index 93% rename from test/test_main.py rename to test/test_main_create_mask.py index 8556970a..37e07e4f 100644 --- a/test/test_main.py +++ b/test/test_main_create_mask.py @@ -5,21 +5,21 @@ import pytest from hydra import compose, initialize -from lidro.main import main +from lidro.main_create_mask import main INPUT_DIR = Path("data") / "pointcloud" -OUTPUT_DIR = Path("tmp") / "main" +OUTPUT_DIR = Path("tmp") / "create_mask_hydro/main" def setup_module(module): - os.makedirs("tmp/main", exist_ok=True) + os.makedirs("tmp/create_mask_hydro/main", exist_ok=True) def test_main_run_okay(): repo_dir = Path.cwd().parent - cmd = f"""python -m lidro.main \ + cmd = f"""python -m lidro.main_create_mask \ io.input_dir="{repo_dir}/lidro/data/pointcloud/"\ - io.output_dir="{repo_dir}/lidro/tmp/main/" + io.output_dir="{repo_dir}/lidro/tmp/create_mask_hydro/main/" """ sp.run(cmd, shell=True, check=True) diff --git a/test/vectors/test_convert_to_vector.py b/test/vectors/test_convert_to_vector.py index 9ef8c849..4514674a 100644 --- a/test/vectors/test_convert_to_vector.py +++ b/test/vectors/test_convert_to_vector.py @@ -3,9 +3,14 @@ import shutil from pathlib import Path -from lidro.vectors.convert_to_vector import create_hydro_vector_mask +from lidro.create_mask_hydro.vectors.convert_to_vector import create_hydro_vector_mask + +TMP_PATH = Path("./tmp/create_mask_hydro/vectors/convert_to_vector") + +las_file = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" +output = "./tmp/create_mask_hydro/vectors/convert_to_vector/MaskHydro_Semis_2021_0830_6291_LA93_IGN69.GeoJSON" +output_main = "./tmp/create_mask_hydro/main/main_lidro_default/MaskHydro_Semis_2021_0830_6291_LA93_IGN69.GeoJSON" -TMP_PATH = Path("./tmp/vectors/convert_to_vector") def setup_module(module): @@ -15,9 +20,6 @@ def setup_module(module): def test_create_hydro_vector_mask_default(): - las_file = "./data/pointcloud/LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.las" - output = "./tmp/vectorize_bins/Test_MaskHydro_LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST.geojson" - crs = 'PROJCS["RGF93 v1 / Lambert-93",GEOGCS["RGF93 v1",DATUM["Reseau_Geodesique_Francais_1993_v1",\ SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],\ AUTHORITY["EPSG","6171"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],\ @@ -28,13 +30,17 @@ def test_create_hydro_vector_mask_default(): UNIT["metre",1,AUTHORITY["EPSG","9001"]],\ AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2154"]]' - Path(output).parent.mkdir(exist_ok=True) - create_hydro_vector_mask(las_file, output, 1, 1000, [0, 1, 2, 3, 6, 17, 66], crs) + create_hydro_vector_mask(las_file, output, 1, 1000, [0, 1, 2, 3, 4, 5, 6, 17, 66], crs, 3) assert Path(output).exists() +def test_check_structure_default(): + # Output with open(output, "r") as f: geojson_data = json.load(f) + + with open(output_main, "r") as f: + geojson_data_main = json.load(f) # CHECK STRUCTURE assert "type" in geojson_data @@ -42,14 +48,15 @@ def test_create_hydro_vector_mask_default(): assert "features" in geojson_data assert isinstance(geojson_data["features"], list) - # Vérifie le quatrième polygone correspondent bien + # CHECK POLYGON for feature in geojson_data["features"]: geometry = feature["geometry"] coordinates = geometry["coordinates"] - assert coordinates[3] == [ - [706231.0, 6626179.0], - [706232.0, 6626179.0], - [706232.0, 6626178.0], - [706231.0, 6626178.0], - [706231.0, 6626179.0], - ] + + for feature in geojson_data_main["features"]: + geometry_main = feature["geometry"] + coordinates_main = geometry_main["coordinates"] + + assert coordinates[0] == coordinates_main[0] + + diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py deleted file mode 100644 index d8175a46..00000000 --- a/test/vectors/test_merge_vector.py +++ /dev/null @@ -1,21 +0,0 @@ -import os -import shutil -from pathlib import Path - -from lidro.vectors.merge_vector import merge_geom - -INPUT_FOLDER = Path("./tmp/main/main_lidro_default") -TMP_PATH = Path("./tmp/vectors/merge_vector") - - -def setup_module(module): - if TMP_PATH.is_dir(): - shutil.rmtree(TMP_PATH) - os.makedirs(TMP_PATH) - - -def test_merge_mask_default(): - output = "tmp/vectors/merge_vector/MaskHydro_merge.geojson" - merge_geom(INPUT_FOLDER, TMP_PATH, "epsg:2154") - - assert Path(output).exists() From 83626499f9f84f615069796f387145f9083fd5e0 Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 26 Apr 2024 17:19:03 +0200 Subject: [PATCH 05/33] create function merge_mask --- configs/configs_lidro.yaml | 11 ++- lidro/main_merge_mask.py | 54 ++++++++++++ .../vectors/check_rectify_geometry.py | 87 +++++++++++++++++++ .../merge_mask_hydro/vectors/merge_vector.py | 67 ++++++++++++++ main_create_mask.log | 8 ++ 5 files changed, 226 insertions(+), 1 deletion(-) create mode 100644 lidro/main_merge_mask.py create mode 100644 lidro/merge_mask_hydro/vectors/check_rectify_geometry.py create mode 100644 lidro/merge_mask_hydro/vectors/merge_vector.py create mode 100644 main_create_mask.log diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index c03cc87d..26cdc1b6 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -20,7 +20,16 @@ io: raster: # size for dilatation dilatation_size: 3 - + +vector: + # Filter water's area (m²) + water_area: 150 + # Parameters for buffer + buffer_positive: 1 + buffer_negative: -0.5 + # Tolerance from DOuglas-Peucker + tolerance: 1 + filter: # Classes to be considered as "non-water" keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 65, 66, 67] diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py new file mode 100644 index 00000000..29cc435f --- /dev/null +++ b/lidro/main_merge_mask.py @@ -0,0 +1,54 @@ +""" Main script for calculate Mask HYDRO 1 +""" + +import logging +import os + +import hydra +from omegaconf import DictConfig +from pyproj import CRS +from lidro.merge_mask_hydro.vectors.merge_vector import merge_geom + + +@hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") +def main(config: DictConfig): + """ Merge all vector masks of hydro surfacesfrom the points classification of the input LAS/LAZ file, + and save it as a GeoJSON file. + + It can run either on each file of a folder + + Args: + config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default) + It contains the algorithm parameters and the input/output parameters + """ + logging.basicConfig(level=logging.INFO) + + # Check input/output files and folders + input_dir = config.io.input_dir + if input_dir is None: + raise ValueError("""config.io.input_dir is empty, please provide an input directory in the configuration""") + + if not os.path.isdir(input_dir): + raise FileNotFoundError(f"""The input directory ({input_dir}) doesn't exist.""") + + output_dir = config.io.output_dir + if output_dir is None: + raise ValueError("""config.io.output_dir is empty, please provide an input directory in the configuration""") + + os.makedirs(output_dir, exist_ok=True) + + # Parameters for cmerging Mask Hydro + water_area = config.vector.water_area # keep only water's area (> 150 m² by default) + buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro + buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro + tolerance = config.vector.tolerance # Tolerance from Douglas-Peucker + crs = CRS.from_user_input(config.io.srid) + + if os.path.isdir(input_dir): + os.makedirs(output_dir, exist_ok=True) # Create folder "merge" + # Merge all Mash Hydro + merge_geom(output_dir, output_dir, crs, water_area, buffer_positive, buffer_negative, tolerance) + + +if __name__ == "__main__": + main() diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py new file mode 100644 index 00000000..3bfa853d --- /dev/null +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -0,0 +1,87 @@ +# -*- coding: utf-8 -*- +""" Rectify and check geometry +""" +from shapely.geometry import CAP_STYLE, MultiPolygon, Polygon +from shapely.validation import make_valid + + +def remove_hole(multipoly): + """Remove small holes (surface < 100 m²) + + Args: + - multipoly (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry without holes (< 100 m²) + """ + list_parts = [] + eps = 100 + + for polygon in multipoly.geoms: + list_interiors = [] + + for interior in polygon.interiors: + p = Polygon(interior) + + if p.area > eps: + list_interiors.append(interior) + + temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) + list_parts.append(temp_pol) + + new_multipolygon = MultiPolygon(list_parts) + + return new_multipolygon + + +def simplify_geometry(polygon, buffer_positive: float, buffer_negative: float): + """Buffer geometry + Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks + + Args: + polygon (GeoJSON): Hydro Mask geometry + buffer_positive (float): positive buffer from Mask Hydro + buffer_negative (float): negative buffer from Mask Hydro + + Returns: + GeoJSON: Hydro Mask geometry simplify + """ + # Buffer + _geom = polygon.buffer(buffer_positive, cap_style=CAP_STYLE.square) + geom = _geom.buffer(buffer_negative, cap_style=CAP_STYLE.square) + return geom + + +def fix_invalid_geometry(geometry): + """Set invalid geoemtries + + Args: + geometry (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry valid + """ + + if not geometry.is_valid: + return make_valid(geometry) + else: + return geometry + + +def check_geometry(initial_gdf): + """Check topology + + Args: + initial_gdf (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry valid + """ + # Obtain simple geometries + gdf_simple = initial_gdf.explode(ignore_index=True) + # Delete duplicate geoemtries if any + gdf_without_duplicates = gdf_simple.drop_duplicates(ignore_index=True) + # Identify invalid geometries and keep only valid ones + gdf_valid = gdf_without_duplicates.copy() + gdf_valid.geometry = gdf_valid.geometry.apply(lambda geom: fix_invalid_geometry(geom)) + return gdf_valid diff --git a/lidro/merge_mask_hydro/vectors/merge_vector.py b/lidro/merge_mask_hydro/vectors/merge_vector.py new file mode 100644 index 00000000..40e0bc0a --- /dev/null +++ b/lidro/merge_mask_hydro/vectors/merge_vector.py @@ -0,0 +1,67 @@ +# -*- coding: utf-8 -*- +""" Merge +""" +import os + +import geopandas as gpd +from shapely.geometry import MultiPolygon +from shapely.ops import unary_union + +from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( + check_geometry, + remove_hole, + simplify_geometry, +) + + +def merge_geom(input_folder: str, output_folder: str, crs: str, water_area: int, buffer_positive: float, buffer_negative: float, tolerance: float): + """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, + filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file. + + Args: + input_folder (str): folder who contains severals Masks Hydro + output_folder (str): output folder + crs (str): a pyproj CRS object used to create the output GeoJSON file + water_area (int): filter Mask Hydro : keep only water's area (> 150 m² by default) + buffer_positive (int): positive buffer from Mask Hydro + buffer_negative (int): negative buffer from Mask Hydro + tolerance (float): All parts of a simplified geometry will be no more than tolerance distance from the original. + """ + # List for stocking all GeoDataFrame + polys = [] + + # Browse all files in folder + for fichier in os.listdir(input_folder): + if fichier.endswith(".GeoJSON"): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(os.path.join(input_folder, fichier)) + merge_geom = geojson["geometry"].unary_union + # Add the GeoDataFrame to the list + polys.append(merge_geom) + + # Union geometry + mergedPolys = unary_union(polys) + + geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) + + # keep only water's area ()> 150 m² by default) + filter_geometry = [geom for geom in geometry if geom.area > water_area] + gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) + + # Check and rectify the invalid geometry + gdf = check_geometry(gdf_filter) + + # Correcting geometric errors: simplifying certain shapes to make calculations easier + buffer_geom = gdf["geometry"].apply(simplify_geometry(buffer_positive, buffer_negative)) + simplify_geom = buffer_geom.simplify(tolerance=tolerance, preserve_topology=True) + + # Correction of holes (< 100m²) in Hydrological Masks + list_parts = [] + for poly in simplify_geom: + list_parts.append(poly) + not_hole_geom = remove_hole(MultiPolygon(list_parts)) + geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False) + gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs) + + # save the result + gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/main_create_mask.log b/main_create_mask.log new file mode 100644 index 00000000..7651dc75 --- /dev/null +++ b/main_create_mask.log @@ -0,0 +1,8 @@ +[2024-04-26 17:17:35,080][root][INFO] - +Create Mask Hydro 1 for tile : LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST +[2024-04-26 17:17:35,212][root][INFO] - +Create Mask Hydro 1 for tile : Semis_2021_0830_6291_LA93_IGN69 +[2024-04-26 17:18:25,871][root][INFO] - +Create Mask Hydro 1 for tile : LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST +[2024-04-26 17:18:25,978][root][INFO] - +Create Mask Hydro 1 for tile : Semis_2021_0830_6291_LA93_IGN69 From 3c46b8cba78aa1aebe9ab88e127cc558cc610109 Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 26 Apr 2024 17:22:24 +0200 Subject: [PATCH 06/33] create function merge_mask --- example_merge_mask_default.sh | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 example_merge_mask_default.sh diff --git a/example_merge_mask_default.sh b/example_merge_mask_default.sh new file mode 100644 index 00000000..c829f219 --- /dev/null +++ b/example_merge_mask_default.sh @@ -0,0 +1,7 @@ +# For lauching Mask Hydro +python -m lidro.main_merge_mask \ +io.input_dir=./tmp/main_lidro_default/ \ +io.output_dir=./tmp/mask_merge/ \ + + + From e3836b3176357bf19a8dd76351e9d3151ac60f2a Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 26 Apr 2024 17:23:59 +0200 Subject: [PATCH 07/33] refacto gitignore --- .gitignore | 3 ++- main_create_mask.log | 8 -------- 2 files changed, 2 insertions(+), 9 deletions(-) delete mode 100644 main_create_mask.log diff --git a/.gitignore b/.gitignore index 08b71877..15b8857d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ **/__pycache__ tmp -main.log \ No newline at end of file +main.log +*.log diff --git a/main_create_mask.log b/main_create_mask.log deleted file mode 100644 index 7651dc75..00000000 --- a/main_create_mask.log +++ /dev/null @@ -1,8 +0,0 @@ -[2024-04-26 17:17:35,080][root][INFO] - -Create Mask Hydro 1 for tile : LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST -[2024-04-26 17:17:35,212][root][INFO] - -Create Mask Hydro 1 for tile : Semis_2021_0830_6291_LA93_IGN69 -[2024-04-26 17:18:25,871][root][INFO] - -Create Mask Hydro 1 for tile : LHD_FXX_0706_6627_PTS_C_LAMB93_IGN69_TEST -[2024-04-26 17:18:25,978][root][INFO] - -Create Mask Hydro 1 for tile : Semis_2021_0830_6291_LA93_IGN69 From 1e7abfb0aac43b41cb750545dc8bdaa4f1ba3ec9 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 29 Apr 2024 16:50:30 +0200 Subject: [PATCH 08/33] refacto merge_mask from dev branche --- .gitignore | 3 +- Dockerfile | 16 ---- Makefile | 40 --------- README.md | 59 +++++++------ ci/test.sh | 1 - configs/configs_lidro.yaml | 5 +- example_merge_mask_default.sh | 7 -- .../rasters/create_mask_raster.py | 41 +++++---- .../vectors/convert_to_vector.py | 4 +- lidro/main_create_mask.py | 5 +- lidro/main_merge_mask.py | 54 ------------ .../vectors/check_rectify_geometry.py | 87 ------------------- .../merge_mask_hydro/vectors/merge_vector.py | 67 -------------- test/pointcloud/test_read_las.py | 2 +- test/rasters/test_create_mask_raster.py | 9 +- test/vectors/test_convert_to_vector.py | 10 +-- 16 files changed, 77 insertions(+), 333 deletions(-) delete mode 100644 Dockerfile delete mode 100644 Makefile delete mode 100644 ci/test.sh delete mode 100644 example_merge_mask_default.sh delete mode 100644 lidro/main_merge_mask.py delete mode 100644 lidro/merge_mask_hydro/vectors/check_rectify_geometry.py delete mode 100644 lidro/merge_mask_hydro/vectors/merge_vector.py diff --git a/.gitignore b/.gitignore index 15b8857d..08b71877 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,3 @@ **/__pycache__ tmp -main.log -*.log +main.log \ No newline at end of file diff --git a/Dockerfile b/Dockerfile deleted file mode 100644 index e69d4f0a..00000000 --- a/Dockerfile +++ /dev/null @@ -1,16 +0,0 @@ -FROM mambaorg/micromamba:latest as mamba_pdal -COPY environment.yml /environment.yml -USER root -RUN micromamba env create -n lidro -f /environment.yml - -FROM debian:bullseye-slim - - -WORKDIR /lidro -RUN mkdir tmp -COPY lidro lidro -COPY test test -COPY configs configs - -# Copy test data that are stored directly in the lidro-data repository ("http://gitlab.forge-idi.ign.fr/Lidar/lidro-data.git") -COPY data/pointcloud data/pointcloud \ No newline at end of file diff --git a/Makefile b/Makefile deleted file mode 100644 index 001ed431..00000000 --- a/Makefile +++ /dev/null @@ -1,40 +0,0 @@ -# Makefile to manage main tasks -# cf. https://blog.ianpreston.ca/conda/python/bash/2020/05/13/conda_envs.html#makefile - -# Oneshell means I can run multiple lines in a recipe in the same shell, so I don't have to -# chain commands together with semicolon -.ONESHELL: -SHELL = /bin/bash -install: - pip install -e . - -install-precommit: - pre-commit install - -testing: - ./ci/test.sh - -mamba-env-create: - mamba env create -n lidro -f environment.yml - -mamba-env-update: - mamba env update -n lidro -f environment.yml - -############################## -# Docker -############################## - -PROJECT_NAME=ignimagelidar/lidro -VERSION=`python -m lidro.version` - -docker-build: - docker build -t ${PROJECT_NAME}:${VERSION} -f Dockerfile . - -docker-test: - docker run --rm -it ${PROJECT_NAME}:${VERSION} python -m pytest -s - -docker-remove: - docker rmi -f `docker images | grep ${PROJECT_NAME} | tr -s ' ' | cut -d ' ' -f 3` - -docker-deploy: - docker push ${PROJECT_NAME}:${VERSION} \ No newline at end of file diff --git a/README.md b/README.md index b9724a1d..8446a3ab 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,5 @@ -### LIDRO +# LIDRO + Lidro (Aplanissement des surfaces d'eaux) est un outil permettant de créer automatiquement des points virtuels le long des surfaces d'eaux afin de créer des modèles numériques cohérents avec les modèles hydrologiques. Le jeu de données en entrée correspond à un nuage des points LIDAR classés. ## Contexte @@ -15,38 +16,43 @@ Les données en entrées : - données vectorielles représentant le réseau hydrographique issu des différentes bases de données IGN (BDUnis, BDTopo, etc.) Trois grands axes du processus à mettre en place en distanguant l'échelle de traitmeent associé : -1- Création de masques hydrographiques à l'échelle de la dalle LIDAR -2- Création de masques hydrographiques pré-filtrés à l'échelle du chantier, soit : -* la suppression de ces masques dans les zones ZICAD/ZIPVA -* la suppression des aires < 150 m² -* la suppression des aires < 1000 m² hors BD IGN (grands cours d'eaux > 5m de large) - -3- Création de points virtuels le long de deux entités hydrographiques : -* ##### Grands cours d'eaux (> 5 m de large dans la BD Unis). -Il existe plusieurs étapes intermédiaires : -1- création automatique du tronçon hydrographique ("Squelette hydrographique", soit les tronçons hydrographiques dans la BD Unid) à partir de l'emprise du masque hydrographique "écoulement" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle) -2- Analyse de la répartition en Z de l'ensemble des points LIDAR "Sol" -3- Créations de points virtuels nécéssitant plusieurs étapes intermédiaires : -* Associer chaque point virtuel 2D au point le plus proche du squelette hydrographique -* Traitement "Z" du squelette : +* 1- Création de masques hydrographiques à l'échelle de la dalle LIDAR +* 2- Création de masques hydrographiques pré-filtrés à l'échelle du chantier, soit : + * la suppression de ces masques dans les zones ZICAD/ZIPVA + * la suppression des aires < 150 m² + * la suppression des aires < 1000 m² hors BD IGN (grands cours d'eau < 5m de large) +* 3- Création de points virtuels le long de deux entités hydrographiques : + * Grands cours d'eau (> 5 m de large dans la BD Unis). + * Surfaces planes (mer, lac, étang, etc.) + +### Traitement des grands cours d'eau (> 5 m de large dans la BD Uns). + +Il existe plusieurs étapes intermédiaires : +* 1- création automatique du tronçon hydrographique ("Squelette hydrographique", soit les tronçons hydrographiques dans la BD Unid) à partir de l'emprise du masque hydrographique "écoulement" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle) +* 2- Analyse de la répartition en Z de l'ensemble des points LIDAR "Sol" +* 3- Création de points virtuels nécéssitant plusieurs étapes intermédiaires : + * Associer chaque point virtuel 2D au point le plus proche du squelette hydrographique + * Traitement "Z" du squelette : * Analyser la répartition en Z de l’ensemble des points LIDAR extrait à l’étape précédente afin d’extraire une seule valeur d’altitude au point fictif Z le plus proche du squelette hydrographique. A noter que l’altitude correspond donc à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes. Pour les ponts, une étape de contrôle de la classification pourrait être mise en place * Lisser les Z le long du squelette HYDRO pour assurer l'écoulement -* Création des points virtuels 2D tous les 0.5 mètres le long des bords du masque hydrographique "écoulement" en cohérence en "Z" avec le squelette hydrographique + * Création des points virtuels 2D tous les 0.5 mètres le long des bords du masque hydrographique "écoulement" en cohérence en "Z" avec le squelette hydrographique -* ##### Surfaces planes (mer, lac, étang, etc.) +### Traitement des surfaces planes (mer, lac, étang, etc.) Pour rappel, l'eau est considérée comme horizontale sur ce type de surface. -Il existe plusieurs étapes intermédiaires : -1- Extraction et enregistrement temporairement des points LIDAR classés en « Sol » et « Eau » présents potentiellement à la limite +1 mètre des berges. Pour cela, on s'appuie sur 'emprise du masque hydrographique "surface plane" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle). a noter que pur le secteur maritime (Mer), il faut exclure la classe « 9 » (eau) afin d’éviter de mesurer les vagues. -2- Analyse statisque de l'ensemble des points LIDAR "Sol / Eau" le long des côtes/berges afin d'obtenir une surface plane. L’objectif est de créer des points virtuels spécifiques avec une information d'altitude (m) tous les 0.5 m sur les bords des surfaces maritimes et des plans d’eaux à partir du masque hydrographique "surface plane". Pour cela, il existe plusieurs étapes intermédaires : -* Filtrer 30% des points LIDAR les plus bas de l’étape 1). afin de supprimer les altitudes trop élevées -* Analyser la répartition en Z de ces points afin d’extraire une seule valeur d’altitude selon l’objet hydrographique : + +Il existe plusieurs étapes intermédiaires : +* 1- Extraction et enregistrement temporairement des points LIDAR classés en « Sol » et « Eau » présents potentiellement à la limite +1 mètre des berges. Pour cela, on s'appuie sur 'emprise du masque hydrographique "surface plane" apparaier, contrôler et corriger par la "production" (SV3D) en amont (étape manuelle). a noter que pur le secteur maritime (Mer), il faut exclure la classe « 9 » (eau) afin d’éviter de mesurer les vagues. +* 2- Analyse statistique de l'ensemble des points LIDAR "Sol / Eau" le long des côtes/berges afin d'obtenir une surface plane. + L’objectif est de créer des points virtuels spécifiques avec une information d'altitude (m) tous les 0.5 m sur les bords des surfaces maritimes et des plans d’eau à partir du masque hydrographique "surface plane". Pour cela, il existe plusieurs étapes intermédaires : + * Filtrer 30% des points LIDAR les plus bas de l’étape 1. afin de supprimer les altitudes trop élevées + * Analyser la répartition en Z de ces points afin d’extraire une seule valeur d’altitude selon l’objet hydrographique : * Pour les Plans d’eau : l’altitude correspond à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes, * Pour la Mer : l’altitude correspond à la ligne basse de la boxplot, soit la valeur minimale en excluant les valeurs aberrantes ## Installation des dépendances (conda) pré-requis: installer Mamba -Cloner le dépôt +Cloner le dépôt ``` git clone https://github.com/IGNF/lidro.git ``` @@ -85,7 +91,8 @@ Lidro se lance sur un seul fichier LAS/LAZ ou sur un Dossier Voir les tests fonctionnels en bas du README. -## Tests unitaires +## Tests +### Tests fonctionnels Tester sur un seul fichier LAS/LAZ ``` example_lidro_by_tile.sh @@ -96,8 +103,8 @@ Tester sur un dossier example_lidro_default.sh ``` -## Tests -Pour lancer les tests : +### Tests unitaires +Pour lancer les tests : ``` python -m pytest -s ``` diff --git a/ci/test.sh b/ci/test.sh deleted file mode 100644 index 823d3d84..00000000 --- a/ci/test.sh +++ /dev/null @@ -1 +0,0 @@ -python -m pytest -s ./test -v \ No newline at end of file diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 26cdc1b6..522c3482 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -19,7 +19,8 @@ io: raster: # size for dilatation - dilatation_size: 3 + dilation_size: 3 + vector: # Filter water's area (m²) @@ -27,7 +28,7 @@ vector: # Parameters for buffer buffer_positive: 1 buffer_negative: -0.5 - # Tolerance from DOuglas-Peucker + # Tolerance from Douglas-Peucker tolerance: 1 filter: diff --git a/example_merge_mask_default.sh b/example_merge_mask_default.sh deleted file mode 100644 index c829f219..00000000 --- a/example_merge_mask_default.sh +++ /dev/null @@ -1,7 +0,0 @@ -# For lauching Mask Hydro -python -m lidro.main_merge_mask \ -io.input_dir=./tmp/main_lidro_default/ \ -io.output_dir=./tmp/mask_merge/ \ - - - diff --git a/lidro/create_mask_hydro/rasters/create_mask_raster.py b/lidro/create_mask_hydro/rasters/create_mask_raster.py index 797119e4..18fb41e0 100644 --- a/lidro/create_mask_hydro/rasters/create_mask_raster.py +++ b/lidro/create_mask_hydro/rasters/create_mask_raster.py @@ -33,22 +33,23 @@ def create_occupancy_map(points: np.array, tile_size: int, pixel_size: float, or return bins -def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, classes: List[int], dilatation_size: int): +def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, classes: List[int], dilation_size: int): """ "Detect hydrographic surfaces in a tile from the classified points of the input pointcloud - An hydrographic surface is define as a surface where there is no points from any class different from water - - Args: - filename (str): input pointcloud - tile_size (int): size of the raster grid (in meters) - pixel_size (float): distance between each node of the raster grid (in meters) - classes (List[int]): List of classes to use for the binarisation (points with other - classification values are ignored) - dilatation_size (int): size for dilatation raster - - Returns: - smoothed_water (np.array): 2D binary array (x, y) of the water presence from the point cloud - pcd_origin (list): top left corner of the tile containing the point cloud - (infered from pointcloud bounding box and input tile size) + An hydrographic surface is defined as a surface where there is no points from any class different from water + The output hydrographic surface mask is dilated to make sure that the masks are continuous when merged with their + neighbours + Args: + filename (str): input pointcloud + tile_size (int): size of the raster grid (in meters) + pixel_size (float): distance between each node of the raster grid (in meters) + classes (List[int]): List of classes to use for the binarisation (points with other + classification values are ignored) + dilation_size (int): size of the structuring element for dilation + + Returns: + smoothed_water (np.array): 2D binary array (x, y) of the water presence from the point cloud + pcd_origin (list): top left corner of the tile containing the point cloud + (infered from pointcloud bounding box and input tile size) """ # Read pointcloud, and extract coordinates (X, Y, Z, and classification) of all points array, crs = read_pointcloud(filename) @@ -65,9 +66,11 @@ def detect_hydro_by_tile(filename: str, tile_size: int, pixel_size: float, class # Revert occupancy map to keep pixels where there is no point of the selected classes detected_water = np.logical_not(occupancy) - # Apply a mathematical morphology operations: DILATATION - # / ! \ NOT "CLOSING", due to the reduction in the size of hydro masks, particularly at the tile borders. + # Apply a mathematical morphology operations: DILATION + # / ! \ NOT "CLOSING", due to the reduction in the size of hydro masks (tile borders) # / ! \ WITH "CLOSING" => Masks Hydro are no longer continuous, when they are merged - morphology_bins = scipy.ndimage.binary_dilation(detected_water, structure=np.ones((dilatation_size, dilatation_size))).astype(np.uint8) + water_mask = scipy.ndimage.binary_dilation( + detected_water, structure=np.ones((dilation_size, dilation_size)) + ).astype(np.uint8) - return morphology_bins, pcd_origin + return water_mask, pcd_origin diff --git a/lidro/create_mask_hydro/vectors/convert_to_vector.py b/lidro/create_mask_hydro/vectors/convert_to_vector.py index 3b0657d5..867e2f36 100644 --- a/lidro/create_mask_hydro/vectors/convert_to_vector.py +++ b/lidro/create_mask_hydro/vectors/convert_to_vector.py @@ -10,7 +10,9 @@ from lidro.create_mask_hydro.rasters.create_mask_raster import detect_hydro_by_tile -def create_hydro_vector_mask(filename: str, output: str, pixel_size: float, tile_size: int, classes: list, crs: str, dilatation_size: int): +def create_hydro_vector_mask( + filename: str, output: str, pixel_size: float, tile_size: int, classes: list, crs: str, dilatation_size: int +): """Create a vector mask of hydro surfaces in a tile from the points classification of the input LAS/LAZ file, and save it as a GeoJSON file. diff --git a/lidro/main_create_mask.py b/lidro/main_create_mask.py index 84313aa8..3e01fb1c 100644 --- a/lidro/main_create_mask.py +++ b/lidro/main_create_mask.py @@ -46,7 +46,7 @@ def main(config: DictConfig): tile_size = config.io.tile_size crs = CRS.from_user_input(config.io.srid) classe = config.filter.keep_classes - dilatation_size = config.raster.dilatation_size + dilation_size = config.raster.dilation_size def main_on_one_tile(filename): """Lauch main.py on one tile @@ -58,7 +58,7 @@ def main_on_one_tile(filename): input_file = os.path.join(input_dir, filename) # path to the LAS file output_file = os.path.join(output_dir, f"MaskHydro_{tilename}.GeoJSON") # path to the Mask Hydro file logging.info(f"\nCreate Mask Hydro 1 for tile : {tilename}") - create_hydro_vector_mask(input_file, output_file, pixel_size, tile_size, classe, crs, dilatation_size) + create_hydro_vector_mask(input_file, output_file, pixel_size, tile_size, classe, crs, dilation_size) if initial_las_filename: # Lauch creating mask by one tile: @@ -69,5 +69,6 @@ def main_on_one_tile(filename): for file in os.listdir(input_dir): main_on_one_tile(file) + if __name__ == "__main__": main() diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py deleted file mode 100644 index 29cc435f..00000000 --- a/lidro/main_merge_mask.py +++ /dev/null @@ -1,54 +0,0 @@ -""" Main script for calculate Mask HYDRO 1 -""" - -import logging -import os - -import hydra -from omegaconf import DictConfig -from pyproj import CRS -from lidro.merge_mask_hydro.vectors.merge_vector import merge_geom - - -@hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") -def main(config: DictConfig): - """ Merge all vector masks of hydro surfacesfrom the points classification of the input LAS/LAZ file, - and save it as a GeoJSON file. - - It can run either on each file of a folder - - Args: - config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default) - It contains the algorithm parameters and the input/output parameters - """ - logging.basicConfig(level=logging.INFO) - - # Check input/output files and folders - input_dir = config.io.input_dir - if input_dir is None: - raise ValueError("""config.io.input_dir is empty, please provide an input directory in the configuration""") - - if not os.path.isdir(input_dir): - raise FileNotFoundError(f"""The input directory ({input_dir}) doesn't exist.""") - - output_dir = config.io.output_dir - if output_dir is None: - raise ValueError("""config.io.output_dir is empty, please provide an input directory in the configuration""") - - os.makedirs(output_dir, exist_ok=True) - - # Parameters for cmerging Mask Hydro - water_area = config.vector.water_area # keep only water's area (> 150 m² by default) - buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro - buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro - tolerance = config.vector.tolerance # Tolerance from Douglas-Peucker - crs = CRS.from_user_input(config.io.srid) - - if os.path.isdir(input_dir): - os.makedirs(output_dir, exist_ok=True) # Create folder "merge" - # Merge all Mash Hydro - merge_geom(output_dir, output_dir, crs, water_area, buffer_positive, buffer_negative, tolerance) - - -if __name__ == "__main__": - main() diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py deleted file mode 100644 index 3bfa853d..00000000 --- a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py +++ /dev/null @@ -1,87 +0,0 @@ -# -*- coding: utf-8 -*- -""" Rectify and check geometry -""" -from shapely.geometry import CAP_STYLE, MultiPolygon, Polygon -from shapely.validation import make_valid - - -def remove_hole(multipoly): - """Remove small holes (surface < 100 m²) - - Args: - - multipoly (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry without holes (< 100 m²) - """ - list_parts = [] - eps = 100 - - for polygon in multipoly.geoms: - list_interiors = [] - - for interior in polygon.interiors: - p = Polygon(interior) - - if p.area > eps: - list_interiors.append(interior) - - temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) - list_parts.append(temp_pol) - - new_multipolygon = MultiPolygon(list_parts) - - return new_multipolygon - - -def simplify_geometry(polygon, buffer_positive: float, buffer_negative: float): - """Buffer geometry - Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks - - Args: - polygon (GeoJSON): Hydro Mask geometry - buffer_positive (float): positive buffer from Mask Hydro - buffer_negative (float): negative buffer from Mask Hydro - - Returns: - GeoJSON: Hydro Mask geometry simplify - """ - # Buffer - _geom = polygon.buffer(buffer_positive, cap_style=CAP_STYLE.square) - geom = _geom.buffer(buffer_negative, cap_style=CAP_STYLE.square) - return geom - - -def fix_invalid_geometry(geometry): - """Set invalid geoemtries - - Args: - geometry (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry valid - """ - - if not geometry.is_valid: - return make_valid(geometry) - else: - return geometry - - -def check_geometry(initial_gdf): - """Check topology - - Args: - initial_gdf (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry valid - """ - # Obtain simple geometries - gdf_simple = initial_gdf.explode(ignore_index=True) - # Delete duplicate geoemtries if any - gdf_without_duplicates = gdf_simple.drop_duplicates(ignore_index=True) - # Identify invalid geometries and keep only valid ones - gdf_valid = gdf_without_duplicates.copy() - gdf_valid.geometry = gdf_valid.geometry.apply(lambda geom: fix_invalid_geometry(geom)) - return gdf_valid diff --git a/lidro/merge_mask_hydro/vectors/merge_vector.py b/lidro/merge_mask_hydro/vectors/merge_vector.py deleted file mode 100644 index 40e0bc0a..00000000 --- a/lidro/merge_mask_hydro/vectors/merge_vector.py +++ /dev/null @@ -1,67 +0,0 @@ -# -*- coding: utf-8 -*- -""" Merge -""" -import os - -import geopandas as gpd -from shapely.geometry import MultiPolygon -from shapely.ops import unary_union - -from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( - check_geometry, - remove_hole, - simplify_geometry, -) - - -def merge_geom(input_folder: str, output_folder: str, crs: str, water_area: int, buffer_positive: float, buffer_negative: float, tolerance: float): - """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, - filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file. - - Args: - input_folder (str): folder who contains severals Masks Hydro - output_folder (str): output folder - crs (str): a pyproj CRS object used to create the output GeoJSON file - water_area (int): filter Mask Hydro : keep only water's area (> 150 m² by default) - buffer_positive (int): positive buffer from Mask Hydro - buffer_negative (int): negative buffer from Mask Hydro - tolerance (float): All parts of a simplified geometry will be no more than tolerance distance from the original. - """ - # List for stocking all GeoDataFrame - polys = [] - - # Browse all files in folder - for fichier in os.listdir(input_folder): - if fichier.endswith(".GeoJSON"): - # Load each GeoJSON file as GeoDataFrame - geojson = gpd.read_file(os.path.join(input_folder, fichier)) - merge_geom = geojson["geometry"].unary_union - # Add the GeoDataFrame to the list - polys.append(merge_geom) - - # Union geometry - mergedPolys = unary_union(polys) - - geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) - - # keep only water's area ()> 150 m² by default) - filter_geometry = [geom for geom in geometry if geom.area > water_area] - gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) - - # Check and rectify the invalid geometry - gdf = check_geometry(gdf_filter) - - # Correcting geometric errors: simplifying certain shapes to make calculations easier - buffer_geom = gdf["geometry"].apply(simplify_geometry(buffer_positive, buffer_negative)) - simplify_geom = buffer_geom.simplify(tolerance=tolerance, preserve_topology=True) - - # Correction of holes (< 100m²) in Hydrological Masks - list_parts = [] - for poly in simplify_geom: - list_parts.append(poly) - not_hole_geom = remove_hole(MultiPolygon(list_parts)) - geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False) - gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs) - - # save the result - gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/test/pointcloud/test_read_las.py b/test/pointcloud/test_read_las.py index 8f3748b4..38778af8 100644 --- a/test/pointcloud/test_read_las.py +++ b/test/pointcloud/test_read_las.py @@ -23,4 +23,4 @@ def test_read_pointcloud_return_format_okay(): assert isinstance(output, np.ndarray) is True - assert CRS.from_user_input("EPSG:2154") + assert crs == CRS.from_user_input("EPSG:2154") diff --git a/test/rasters/test_create_mask_raster.py b/test/rasters/test_create_mask_raster.py index 33865095..078dac2d 100644 --- a/test/rasters/test_create_mask_raster.py +++ b/test/rasters/test_create_mask_raster.py @@ -7,7 +7,10 @@ import rasterio from rasterio.transform import from_origin -from lidro.create_mask_hydro.rasters.create_mask_raster import create_occupancy_map, detect_hydro_by_tile +from lidro.create_mask_hydro.rasters.create_mask_raster import ( + create_occupancy_map, + detect_hydro_by_tile, +) TMP_PATH = Path("./tmp/create_mask_hydro/rasters/create_mask_raster") @@ -43,7 +46,9 @@ def test_create_occupancy_map_default(): @pytest.mark.returnfile def test_detect_hydro_by_tile_return_file(): output_tif = TMP_PATH / "Semis_2021_0830_6291_LA93_IGN69_size.tif" - array, origin = detect_hydro_by_tile(LAS_FILE, tile_size, pixel_size, classes=[0, 1, 2, 3, 4, 5, 6, 17, 66], dilatation_size=1) + array, origin = detect_hydro_by_tile( + LAS_FILE, tile_size, pixel_size, classes=[0, 1, 2, 3, 4, 5, 6, 17, 66], dilation_size=1 + ) assert isinstance(array, np.ndarray) is True assert list(array.shape) == [tile_size / pixel_size] * 2 diff --git a/test/vectors/test_convert_to_vector.py b/test/vectors/test_convert_to_vector.py index 4514674a..8678c381 100644 --- a/test/vectors/test_convert_to_vector.py +++ b/test/vectors/test_convert_to_vector.py @@ -12,7 +12,6 @@ output_main = "./tmp/create_mask_hydro/main/main_lidro_default/MaskHydro_Semis_2021_0830_6291_LA93_IGN69.GeoJSON" - def setup_module(module): if TMP_PATH.is_dir(): shutil.rmtree(TMP_PATH) @@ -34,11 +33,12 @@ def test_create_hydro_vector_mask_default(): assert Path(output).exists() + def test_check_structure_default(): # Output with open(output, "r") as f: geojson_data = json.load(f) - + with open(output_main, "r") as f: geojson_data_main = json.load(f) @@ -48,7 +48,7 @@ def test_check_structure_default(): assert "features" in geojson_data assert isinstance(geojson_data["features"], list) - # CHECK POLYGON + # CHECK POLYGON for feature in geojson_data["features"]: geometry = feature["geometry"] coordinates = geometry["coordinates"] @@ -56,7 +56,5 @@ def test_check_structure_default(): for feature in geojson_data_main["features"]: geometry_main = feature["geometry"] coordinates_main = geometry_main["coordinates"] - - assert coordinates[0] == coordinates_main[0] - + assert coordinates[0] == coordinates_main[0] From 78a067f32ac4b16ae88deeb96a0f8770d467048f Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 29 Apr 2024 16:52:10 +0200 Subject: [PATCH 09/33] add function merge_mask and test --- lidro/main_merge_mask.py | 54 +++++++++++ .../vectors/check_rectify_geometry.py | 89 +++++++++++++++++++ .../merge_mask_hydro/vectors/merge_vector.py | 67 ++++++++++++++ test/vectors/test_merge_vector.py | 53 +++++++++++ 4 files changed, 263 insertions(+) create mode 100644 lidro/main_merge_mask.py create mode 100644 lidro/merge_mask_hydro/vectors/check_rectify_geometry.py create mode 100644 lidro/merge_mask_hydro/vectors/merge_vector.py create mode 100644 test/vectors/test_merge_vector.py diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py new file mode 100644 index 00000000..29cc435f --- /dev/null +++ b/lidro/main_merge_mask.py @@ -0,0 +1,54 @@ +""" Main script for calculate Mask HYDRO 1 +""" + +import logging +import os + +import hydra +from omegaconf import DictConfig +from pyproj import CRS +from lidro.merge_mask_hydro.vectors.merge_vector import merge_geom + + +@hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") +def main(config: DictConfig): + """ Merge all vector masks of hydro surfacesfrom the points classification of the input LAS/LAZ file, + and save it as a GeoJSON file. + + It can run either on each file of a folder + + Args: + config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default) + It contains the algorithm parameters and the input/output parameters + """ + logging.basicConfig(level=logging.INFO) + + # Check input/output files and folders + input_dir = config.io.input_dir + if input_dir is None: + raise ValueError("""config.io.input_dir is empty, please provide an input directory in the configuration""") + + if not os.path.isdir(input_dir): + raise FileNotFoundError(f"""The input directory ({input_dir}) doesn't exist.""") + + output_dir = config.io.output_dir + if output_dir is None: + raise ValueError("""config.io.output_dir is empty, please provide an input directory in the configuration""") + + os.makedirs(output_dir, exist_ok=True) + + # Parameters for cmerging Mask Hydro + water_area = config.vector.water_area # keep only water's area (> 150 m² by default) + buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro + buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro + tolerance = config.vector.tolerance # Tolerance from Douglas-Peucker + crs = CRS.from_user_input(config.io.srid) + + if os.path.isdir(input_dir): + os.makedirs(output_dir, exist_ok=True) # Create folder "merge" + # Merge all Mash Hydro + merge_geom(output_dir, output_dir, crs, water_area, buffer_positive, buffer_negative, tolerance) + + +if __name__ == "__main__": + main() diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py new file mode 100644 index 00000000..ce62e21d --- /dev/null +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -0,0 +1,89 @@ +# -*- coding: utf-8 -*- +""" Rectify and check geometry +""" +import geopandas as gpd +from shapely.geometry import CAP_STYLE, MultiPolygon, Polygon +from shapely.validation import make_valid + + +def remove_hole(multipoly): + """Remove small holes (surface < 100 m²) + + Args: + - multipoly (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry without holes (< 100 m²) + """ + list_parts = [] + eps = 100 + + for polygon in multipoly.geoms: + list_interiors = [] + + for interior in polygon.interiors: + p = Polygon(interior) + + if p.area > eps: + list_interiors.append(interior) + + temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) + list_parts.append(temp_pol) + + new_multipolygon = MultiPolygon(list_parts) + + return new_multipolygon + + +def simplify_geometry(s, buffer_positive: float, buffer_negative: float): + """Buffer geometry + Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks + + Args: + s (GeoJSON): Hydro Mask geometry + buffer_positive (float): positive buffer from Mask Hydro + buffer_negative (float): negative buffer from Mask Hydro + + Returns: + GeoJSON: Hydro Mask geometry simplify + """ + # Buffer + _geom = s.buffer(buffer_positive, cap_style=CAP_STYLE.square) + #_geom = polygon.buffer(buffer_positive, cap_style=CAP_STYLE.square) + geom = _geom.buffer(buffer_negative, cap_style=CAP_STYLE.square) + return geom + + +def fix_invalid_geometry(geometry): + """Set invalid geoemtries + + Args: + geometry (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry valid + """ + + if not geometry.is_valid: + return make_valid(geometry) + else: + return geometry + + +def check_geometry(initial_gdf): + """Check topology + + Args: + initial_gdf (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry valid + """ + # Obtain simple geometries + gdf_simple = initial_gdf.explode(ignore_index=True) + # Delete duplicate geoemtries if any + gdf_without_duplicates = gdf_simple.drop_duplicates(ignore_index=True) + # Identify invalid geometries and keep only valid ones + gdf_valid = gdf_without_duplicates.copy() + gdf_valid.geometry = gdf_valid.geometry.apply(lambda geom: fix_invalid_geometry(geom)) + return gdf_valid diff --git a/lidro/merge_mask_hydro/vectors/merge_vector.py b/lidro/merge_mask_hydro/vectors/merge_vector.py new file mode 100644 index 00000000..b8778baa --- /dev/null +++ b/lidro/merge_mask_hydro/vectors/merge_vector.py @@ -0,0 +1,67 @@ +# -*- coding: utf-8 -*- +""" Merge +""" +import os + +import geopandas as gpd +from shapely.geometry import MultiPolygon +from shapely.ops import unary_union + +from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( + check_geometry, + remove_hole, + simplify_geometry, +) + + +def merge_geom(input_folder: str, output_folder: str, crs: str, water_area: int, buffer_positive: float, buffer_negative: float, tolerance: float): + """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, + filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file. + + Args: + input_folder (str): folder who contains severals Masks Hydro + output_folder (str): output folder + crs (str): a pyproj CRS object used to create the output GeoJSON file + water_area (int): filter Mask Hydro : keep only water's area (> 150 m² by default) + buffer_positive (int): positive buffer from Mask Hydro + buffer_negative (int): negative buffer from Mask Hydro + tolerance (float): All parts of a simplified geometry will be no more than tolerance distance from the original. + """ + # List for stocking all GeoDataFrame + polys = [] + + # Browse all files in folder + for fichier in os.listdir(input_folder): + if fichier.endswith(".GeoJSON"): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(os.path.join(input_folder, fichier)) + merge_geom = geojson["geometry"].unary_union + # Add the GeoDataFrame to the list + polys.append(merge_geom) + + # Union geometry + mergedPolys = unary_union(polys) + + geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) + + # keep only water's area (> 150 m² by default) + filter_geometry = [geom for geom in geometry if geom.area > water_area] + gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) + + # Check and rectify the invalid geometry + gdf = check_geometry(gdf_filter) + + # Correcting geometric errors: simplifying certain shapes to make calculations easier + buffer_geom = simplify_geometry(gdf["geometry"], buffer_positive, buffer_negative).explode(index_parts=False) + simplify_geom = buffer_geom.simplify(tolerance=tolerance, preserve_topology=True) + + # Correction of holes (< 100m²) in Hydrological Masks + list_parts = [] + for poly in simplify_geom: + list_parts.append(poly) + not_hole_geom = remove_hole(MultiPolygon(list_parts)) + geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False) + gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs) + + # save the result + gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py new file mode 100644 index 00000000..8dce952d --- /dev/null +++ b/test/vectors/test_merge_vector.py @@ -0,0 +1,53 @@ +import json +import os +import shutil +from pathlib import Path + +from lidro.merge_mask_hydro.vectors.merge_vector import merge_geom + +TMP_PATH = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro") + +input_folder = "./data/mask_hydro/" +output = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro/MaskHydro_merge.geojson") + + + +def setup_module(module): + if TMP_PATH.is_dir(): + shutil.rmtree(TMP_PATH) + os.makedirs(TMP_PATH) + + +def test_create_hydro_vector_mask_default(): + crs = 'PROJCS["RGF93 v1 / Lambert-93",GEOGCS["RGF93 v1",DATUM["Reseau_Geodesique_Francais_1993_v1",\ + SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],\ + AUTHORITY["EPSG","6171"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],\ + UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4171"]],\ + PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",46.5],\ + PARAMETER["central_meridian",3],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",44],\ + PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],\ + UNIT["metre",1,AUTHORITY["EPSG","9001"]],\ + AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2154"]]' + + merge_geom(input_folder, TMP_PATH, crs, 150, 0.5, -1.5, 1) + + assert Path(output).exists() + +# def test_check_structure_default(): +# # Output +# with open(output, "r") as f: +# geojson_data = json.load(f) + +# # CHECK STRUCTURE +# assert "type" in geojson_data +# assert geojson_data["type"] == "FeatureCollection" +# assert "features" in geojson_data +# assert isinstance(geojson_data["features"], list) + +# # CHECK POLYGON +# for feature in geojson_data["features"]: +# geometry = feature["geometry"] +# coordinates = geometry["coordinates"] + + #assert coordinates[0] == coordinates_main[0] + \ No newline at end of file From 5befd2856dcc83e465f4c48beb3542cd1d0f1599 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 29 Apr 2024 16:53:10 +0200 Subject: [PATCH 10/33] add script for lauching merge_mask --- example_merge_mask_default.sh | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 example_merge_mask_default.sh diff --git a/example_merge_mask_default.sh b/example_merge_mask_default.sh new file mode 100644 index 00000000..c829f219 --- /dev/null +++ b/example_merge_mask_default.sh @@ -0,0 +1,7 @@ +# For lauching Mask Hydro +python -m lidro.main_merge_mask \ +io.input_dir=./tmp/main_lidro_default/ \ +io.output_dir=./tmp/mask_merge/ \ + + + From 23757cc9ab0788141c03471c001a11a916bc20c9 Mon Sep 17 00:00:00 2001 From: mdupays Date: Fri, 3 May 2024 16:30:19 +0200 Subject: [PATCH 11/33] add severals tests for merge_mask_vector --- example_merge_mask_default.sh | 4 +- lidro/main_merge_mask.py | 13 ++- .../vectors/check_rectify_geometry.py | 33 +----- .../merge_mask_hydro/vectors/merge_vector.py | 20 +++- lidro/merge_mask_hydro/vectors/remove_hole.py | 32 ++++++ test/test_main_merge_mask.py | 103 ++++++++++++++++++ test/vectors/test_check_rectify_geometry.py | 33 ++++++ test/vectors/test_merge_vector.py | 45 ++++---- test/vectors/test_remove_hole.py | 15 +++ 9 files changed, 235 insertions(+), 63 deletions(-) create mode 100644 lidro/merge_mask_hydro/vectors/remove_hole.py create mode 100644 test/test_main_merge_mask.py create mode 100644 test/vectors/test_check_rectify_geometry.py create mode 100644 test/vectors/test_remove_hole.py diff --git a/example_merge_mask_default.sh b/example_merge_mask_default.sh index c829f219..a2d6e65f 100644 --- a/example_merge_mask_default.sh +++ b/example_merge_mask_default.sh @@ -1,7 +1,7 @@ # For lauching Mask Hydro python -m lidro.main_merge_mask \ -io.input_dir=./tmp/main_lidro_default/ \ -io.output_dir=./tmp/mask_merge/ \ +io.input_dir=./data/mask_hydro/ \ +io.output_dir=./tmp/merge_mask_hydro/ \ diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py index 29cc435f..faabb136 100644 --- a/lidro/main_merge_mask.py +++ b/lidro/main_merge_mask.py @@ -7,12 +7,13 @@ import hydra from omegaconf import DictConfig from pyproj import CRS + from lidro.merge_mask_hydro.vectors.merge_vector import merge_geom @hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") def main(config: DictConfig): - """ Merge all vector masks of hydro surfacesfrom the points classification of the input LAS/LAZ file, + """Merge all vector masks of hydro surfacesfrom the points classification of the input LAS/LAZ file, and save it as a GeoJSON file. It can run either on each file of a folder @@ -38,16 +39,16 @@ def main(config: DictConfig): os.makedirs(output_dir, exist_ok=True) # Parameters for cmerging Mask Hydro - water_area = config.vector.water_area # keep only water's area (> 150 m² by default) - buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro - buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro - tolerance = config.vector.tolerance # Tolerance from Douglas-Peucker + water_area = config.vector.water_area # keep only water's area (> 150 m² by default) + buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro + buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro + tolerance = config.vector.tolerance # Tolerance from Douglas-Peucker crs = CRS.from_user_input(config.io.srid) if os.path.isdir(input_dir): os.makedirs(output_dir, exist_ok=True) # Create folder "merge" # Merge all Mash Hydro - merge_geom(output_dir, output_dir, crs, water_area, buffer_positive, buffer_negative, tolerance) + merge_geom(input_dir, output_dir, crs, water_area, buffer_positive, buffer_negative, tolerance) if __name__ == "__main__": diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py index ce62e21d..e42f020e 100644 --- a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -1,40 +1,10 @@ # -*- coding: utf-8 -*- """ Rectify and check geometry """ -import geopandas as gpd -from shapely.geometry import CAP_STYLE, MultiPolygon, Polygon +from shapely.geometry import CAP_STYLE from shapely.validation import make_valid -def remove_hole(multipoly): - """Remove small holes (surface < 100 m²) - - Args: - - multipoly (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry without holes (< 100 m²) - """ - list_parts = [] - eps = 100 - - for polygon in multipoly.geoms: - list_interiors = [] - - for interior in polygon.interiors: - p = Polygon(interior) - - if p.area > eps: - list_interiors.append(interior) - - temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) - list_parts.append(temp_pol) - - new_multipolygon = MultiPolygon(list_parts) - - return new_multipolygon - - def simplify_geometry(s, buffer_positive: float, buffer_negative: float): """Buffer geometry Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks @@ -49,7 +19,6 @@ def simplify_geometry(s, buffer_positive: float, buffer_negative: float): """ # Buffer _geom = s.buffer(buffer_positive, cap_style=CAP_STYLE.square) - #_geom = polygon.buffer(buffer_positive, cap_style=CAP_STYLE.square) geom = _geom.buffer(buffer_negative, cap_style=CAP_STYLE.square) return geom diff --git a/lidro/merge_mask_hydro/vectors/merge_vector.py b/lidro/merge_mask_hydro/vectors/merge_vector.py index b8778baa..6bf4409e 100644 --- a/lidro/merge_mask_hydro/vectors/merge_vector.py +++ b/lidro/merge_mask_hydro/vectors/merge_vector.py @@ -9,12 +9,20 @@ from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( check_geometry, - remove_hole, simplify_geometry, ) +from lidro.merge_mask_hydro.vectors.remove_hole import remove_hole -def merge_geom(input_folder: str, output_folder: str, crs: str, water_area: int, buffer_positive: float, buffer_negative: float, tolerance: float): +def merge_geom( + input_folder: str, + output_folder: str, + crs: str, + water_area: int, + buffer_positive: float, + buffer_negative: float, + tolerance: float, +): """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file. @@ -25,7 +33,8 @@ def merge_geom(input_folder: str, output_folder: str, crs: str, water_area: int, water_area (int): filter Mask Hydro : keep only water's area (> 150 m² by default) buffer_positive (int): positive buffer from Mask Hydro buffer_negative (int): negative buffer from Mask Hydro - tolerance (float): All parts of a simplified geometry will be no more than tolerance distance from the original. + tolerance (float): All parts of a simplified geometry will be no more + than tolerance distance from the original. """ # List for stocking all GeoDataFrame polys = [] @@ -59,9 +68,12 @@ def merge_geom(input_folder: str, output_folder: str, crs: str, water_area: int, list_parts = [] for poly in simplify_geom: list_parts.append(poly) + not_hole_geom = remove_hole(MultiPolygon(list_parts)) geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False) - gdf = gpd.GeoDataFrame(geometry=geometry, crs=crs) + # keep only water's area (> 150 m² by default) : + filter_geometry_second = [geom for geom in geometry if geom.area > water_area] + gdf = gpd.GeoDataFrame(geometry=filter_geometry_second, crs=crs) # save the result gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/lidro/merge_mask_hydro/vectors/remove_hole.py b/lidro/merge_mask_hydro/vectors/remove_hole.py new file mode 100644 index 00000000..420acbbd --- /dev/null +++ b/lidro/merge_mask_hydro/vectors/remove_hole.py @@ -0,0 +1,32 @@ +# -*- coding: utf-8 -*- +""" Remove small holes """ +from shapely.geometry import MultiPolygon, Polygon + + +def remove_hole(multipoly): + """Remove small holes (surface < 100 m²) + + Args: + - multipoly (GeoJSON): Hydro Mask geometry + + Returns: + GeoJSON: Hydro Mask geometry without holes (< 100 m²) + """ + list_parts = [] + eps = 100 + + for polygon in multipoly.geoms: + list_interiors = [] + + for interior in polygon.interiors: + p = Polygon(interior) + + if p.area > eps: + list_interiors.append(interior) + + temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) + list_parts.append(temp_pol) + + new_multipolygon = MultiPolygon(list_parts) + + return new_multipolygon diff --git a/test/test_main_merge_mask.py b/test/test_main_merge_mask.py new file mode 100644 index 00000000..3f627249 --- /dev/null +++ b/test/test_main_merge_mask.py @@ -0,0 +1,103 @@ +import os +import subprocess as sp +from pathlib import Path + +import pytest +from hydra import compose, initialize + +from lidro.main_merge_mask import main + +INPUT_DIR = Path("data") / "mask_hydro" +OUTPUT_DIR = Path("tmp") / "merge_mask_hydro/main" + + +def setup_module(module): + os.makedirs("tmp/merge_mask_hydro/main", exist_ok=True) + + +def test_main_run_okay(): + repo_dir = Path.cwd().parent + cmd = f"""python -m lidro.main_merge_mask \ + io.input_dir="{repo_dir}/lidro/data/mask_hydro/"\ + io.output_dir="{repo_dir}/lidro/tmp/merge_mask_hydro/main/" + """ + sp.run(cmd, shell=True, check=True) + + +def test_main_lidro_fail_no_input_dir(): + output_dir = OUTPUT_DIR / "main_no_input_dir" + pixel_size = 1 + water_area = 150 + buffer_positive = 0.5 + buffer_negative = -1.5 + tolerance = 1 + srid = 2154 + with initialize(version_base="1.2", config_path="../configs"): + # config is relative to a module + cfg = compose( + config_name="configs_lidro", + overrides=[ + f"io.output_dir={output_dir}", + f"io.pixel_size={pixel_size}", + f"vector.water_area={water_area}", + f"vector.buffer_positive={buffer_positive}", + f"vector.buffer_negative={buffer_negative}", + f"vector.tolerance={tolerance}", + f"io.srid={srid}", + ], + ) + with pytest.raises(ValueError): + main(cfg) + + +def test_main_lidro_fail_wrong_input_dir(): + output_dir = OUTPUT_DIR / "main_wrong_input_dir" + pixel_size = 1 + water_area = 150 + buffer_positive = 0.5 + buffer_negative = -1.5 + tolerance = 1 + srid = 2154 + with initialize(version_base="1.2", config_path="../configs"): + # config is relative to a module + cfg = compose( + config_name="configs_lidro", + overrides=[ + "io.input_dir=some_random_input_dir", + f"io.output_dir={output_dir}", + f"io.pixel_size={pixel_size}", + f"vector.water_area={water_area}", + f"vector.buffer_positive={buffer_positive}", + f"vector.buffer_negative={buffer_negative}", + f"vector.tolerance={tolerance}", + f"io.srid={srid}", + ], + ) + with pytest.raises(FileNotFoundError): + main(cfg) + + +def test_main_lidro_fail_no_output_dir(): + input_dir = INPUT_DIR + pixel_size = 1 + water_area = 150 + buffer_positive = 0.5 + buffer_negative = -1.5 + tolerance = 1 + srid = 2154 + with initialize(version_base="1.2", config_path="../configs"): + # config is relative to a module + cfg = compose( + config_name="configs_lidro", + overrides=[ + f"io.input_dir={input_dir}", + f"io.pixel_size={pixel_size}", + f"vector.water_area={water_area}", + f"vector.buffer_positive={buffer_positive}", + f"vector.buffer_negative={buffer_negative}", + f"vector.tolerance={tolerance}", + f"io.srid={srid}", + ], + ) + with pytest.raises(ValueError): + main(cfg) diff --git a/test/vectors/test_check_rectify_geometry.py b/test/vectors/test_check_rectify_geometry.py new file mode 100644 index 00000000..e010c657 --- /dev/null +++ b/test/vectors/test_check_rectify_geometry.py @@ -0,0 +1,33 @@ +import geopandas as gpd +from shapely.geometry import MultiPolygon + +from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( + check_geometry, + fix_invalid_geometry, + simplify_geometry, +) + +input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" + + +def test_simplify_geometry_default(): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(input) + geometry = geojson["geometry"].unary_union + buffer = simplify_geometry(geometry, 1, -1) + assert isinstance(buffer, MultiPolygon) + + +def test_check_geometry_default(): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(input) + check_geom = check_geometry(geojson) + assert isinstance(check_geom.dtypes, object) + + +def test_fix_invalid_geometry_default(): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(input) + valid_geom = geojson.geometry.apply(lambda geom: fix_invalid_geometry(geom)) + assert isinstance(valid_geom.dtypes, object) + assert valid_geom[2] == geojson.geometry[2] diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py index 8dce952d..eee73225 100644 --- a/test/vectors/test_merge_vector.py +++ b/test/vectors/test_merge_vector.py @@ -9,7 +9,7 @@ input_folder = "./data/mask_hydro/" output = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro/MaskHydro_merge.geojson") - +output_main = "./data/merge_mask_hydro/MaskHydro_merge.geojson" def setup_module(module): @@ -33,21 +33,28 @@ def test_create_hydro_vector_mask_default(): assert Path(output).exists() -# def test_check_structure_default(): -# # Output -# with open(output, "r") as f: -# geojson_data = json.load(f) - -# # CHECK STRUCTURE -# assert "type" in geojson_data -# assert geojson_data["type"] == "FeatureCollection" -# assert "features" in geojson_data -# assert isinstance(geojson_data["features"], list) - -# # CHECK POLYGON -# for feature in geojson_data["features"]: -# geometry = feature["geometry"] -# coordinates = geometry["coordinates"] - - #assert coordinates[0] == coordinates_main[0] - \ No newline at end of file + +def test_check_structure_default(): + # Output + with open(output, "r") as f: + geojson_data = json.load(f) + + with open(output_main, "r") as f: + geojson_data_main = json.load(f) + + # CHECK STRUCTURE + assert "type" in geojson_data + assert geojson_data["type"] == "FeatureCollection" + assert "features" in geojson_data + assert isinstance(geojson_data["features"], list) + + # CHECK POLYGON + for feature in geojson_data["features"]: + geometry = feature["geometry"] + coordinates = geometry["coordinates"] + + for feature in geojson_data_main["features"]: + geometry_main = feature["geometry"] + coordinates_main = geometry_main["coordinates"] + + assert coordinates[0] == coordinates_main[0] diff --git a/test/vectors/test_remove_hole.py b/test/vectors/test_remove_hole.py new file mode 100644 index 00000000..2c066927 --- /dev/null +++ b/test/vectors/test_remove_hole.py @@ -0,0 +1,15 @@ +import geopandas as gpd +from shapely.geometry import MultiPolygon + +from lidro.merge_mask_hydro.vectors.remove_hole import remove_hole + +input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" + + +def test_remove_hole_default(): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(input) + geometry = geojson["geometry"].unary_union + mask_without_hole = remove_hole(geometry) + + assert isinstance(mask_without_hole, MultiPolygon) From ab04c229f5063cf319f367264a7827c7853d1fd7 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 6 May 2024 10:55:51 +0200 Subject: [PATCH 12/33] update folder data --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 48ae4faf..36d9aab8 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 48ae4faf6a8a2c4233a22e927443e78a98c73d31 +Subproject commit 36d9aab8ae1e8eeb43449c89c2cbb136a0b1caa2 From 1a9d825d3ebc07e121115071e70f892ed1d08f99 Mon Sep 17 00:00:00 2001 From: mdupaysign <162304206+mdupaysign@users.noreply.github.com> Date: Mon, 6 May 2024 14:08:11 +0200 Subject: [PATCH 13/33] Update lidro/main_merge_mask.py Co-authored-by: leavauchier <120112647+leavauchier@users.noreply.github.com> --- lidro/main_merge_mask.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py index faabb136..26309b64 100644 --- a/lidro/main_merge_mask.py +++ b/lidro/main_merge_mask.py @@ -48,7 +48,14 @@ def main(config: DictConfig): if os.path.isdir(input_dir): os.makedirs(output_dir, exist_ok=True) # Create folder "merge" # Merge all Mash Hydro - merge_geom(input_dir, output_dir, crs, water_area, buffer_positive, buffer_negative, tolerance) + merge_geom( + input_dir, + output_dir, + CRS.from_user_input(config.io.srid), + config.vector.water_area, + config.vector.buffer_positive, + config.vector.buffer_negative, + config.vector.tolerance) if __name__ == "__main__": From 91ea20cc7dacdee86b2d5b8175c9a75361893f55 Mon Sep 17 00:00:00 2001 From: mdupaysign <162304206+mdupaysign@users.noreply.github.com> Date: Mon, 6 May 2024 14:08:30 +0200 Subject: [PATCH 14/33] Update lidro/main_merge_mask.py Co-authored-by: leavauchier <120112647+leavauchier@users.noreply.github.com> --- lidro/main_merge_mask.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py index 26309b64..7066990e 100644 --- a/lidro/main_merge_mask.py +++ b/lidro/main_merge_mask.py @@ -38,7 +38,7 @@ def main(config: DictConfig): os.makedirs(output_dir, exist_ok=True) - # Parameters for cmerging Mask Hydro + # Parameters for merging Mask Hydro water_area = config.vector.water_area # keep only water's area (> 150 m² by default) buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro From ed15b2f199c4cf17f82855bed4cbe5f20b15e1ea Mon Sep 17 00:00:00 2001 From: mdupaysign <162304206+mdupaysign@users.noreply.github.com> Date: Mon, 6 May 2024 14:08:40 +0200 Subject: [PATCH 15/33] Update lidro/main_merge_mask.py Co-authored-by: leavauchier <120112647+leavauchier@users.noreply.github.com> --- lidro/main_merge_mask.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py index 7066990e..af56b21b 100644 --- a/lidro/main_merge_mask.py +++ b/lidro/main_merge_mask.py @@ -13,7 +13,7 @@ @hydra.main(config_path="../configs/", config_name="configs_lidro.yaml", version_base="1.2") def main(config: DictConfig): - """Merge all vector masks of hydro surfacesfrom the points classification of the input LAS/LAZ file, + """Merge all vector masks of hydro surfaces from the points classification of the input LAS/LAZ file, and save it as a GeoJSON file. It can run either on each file of a folder From ef3c84d840ef386a79140d40afe18fdbc05ffe89 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 6 May 2024 14:15:53 +0200 Subject: [PATCH 16/33] update main_merge_mask.py --- configs/configs_lidro.yaml | 2 +- lidro/main_merge_mask.py | 6 ++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index a7caaa4d..9963df7b 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -23,7 +23,7 @@ raster: vector: # Filter water's area (m²) - water_area: 150 + min_water_area: 150 # Parameters for buffer buffer_positive: 1 buffer_negative: -0.5 diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py index af56b21b..be6dfcae 100644 --- a/lidro/main_merge_mask.py +++ b/lidro/main_merge_mask.py @@ -16,8 +16,6 @@ def main(config: DictConfig): """Merge all vector masks of hydro surfaces from the points classification of the input LAS/LAZ file, and save it as a GeoJSON file. - It can run either on each file of a folder - Args: config (DictConfig): hydra configuration (configs/configs_lidro.yaml by default) It contains the algorithm parameters and the input/output parameters @@ -39,7 +37,7 @@ def main(config: DictConfig): os.makedirs(output_dir, exist_ok=True) # Parameters for merging Mask Hydro - water_area = config.vector.water_area # keep only water's area (> 150 m² by default) + min_water_area = config.vector.min_water_area # keep only water's area (> 150 m² by default) buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro tolerance = config.vector.tolerance # Tolerance from Douglas-Peucker @@ -52,7 +50,7 @@ def main(config: DictConfig): input_dir, output_dir, CRS.from_user_input(config.io.srid), - config.vector.water_area, + config.vector.min_water_area, config.vector.buffer_positive, config.vector.buffer_negative, config.vector.tolerance) From a3ab21bc31ea1c75969bf2f5925ece7d61fc9ba3 Mon Sep 17 00:00:00 2001 From: mdupays Date: Tue, 7 May 2024 15:36:57 +0200 Subject: [PATCH 17/33] add test_clos_holes --- .../vectors/check_rectify_geometry.py | 38 +++-------- lidro/merge_mask_hydro/vectors/close_holes.py | 45 +++++++++++++ .../merge_mask_hydro/vectors/merge_vector.py | 64 ++++++++----------- lidro/merge_mask_hydro/vectors/remove_hole.py | 32 ---------- test/test_main_merge_mask.py | 12 ++-- test/vectors/test_check_rectify_geometry.py | 15 +---- ...est_remove_hole.py => test_close_holes.py} | 6 +- 7 files changed, 94 insertions(+), 118 deletions(-) create mode 100644 lidro/merge_mask_hydro/vectors/close_holes.py delete mode 100644 lidro/merge_mask_hydro/vectors/remove_hole.py rename test/vectors/{test_remove_hole.py => test_close_holes.py} (67%) diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py index e42f020e..11452120 100644 --- a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -2,43 +2,25 @@ """ Rectify and check geometry """ from shapely.geometry import CAP_STYLE -from shapely.validation import make_valid -def simplify_geometry(s, buffer_positive: float, buffer_negative: float): +def apply_buffers_to_geometry(hydro_mask, buffer_positive: float, buffer_negative: float): """Buffer geometry Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks Args: - s (GeoJSON): Hydro Mask geometry - buffer_positive (float): positive buffer from Mask Hydro - buffer_negative (float): negative buffer from Mask Hydro - + s (gpd.GeoDataFrame): geopandas dataframe with input geometry + buffer_positive (float): positive buffer to apply to the input geometry + buffer_negative (float): negative buffer to apply to the input geometry Returns: - GeoJSON: Hydro Mask geometry simplify + GeoJSON: updated geometry """ # Buffer - _geom = s.buffer(buffer_positive, cap_style=CAP_STYLE.square) + _geom = hydro_mask.buffer(buffer_positive, cap_style=CAP_STYLE.square) geom = _geom.buffer(buffer_negative, cap_style=CAP_STYLE.square) return geom -def fix_invalid_geometry(geometry): - """Set invalid geoemtries - - Args: - geometry (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry valid - """ - - if not geometry.is_valid: - return make_valid(geometry) - else: - return geometry - - def check_geometry(initial_gdf): """Check topology @@ -51,8 +33,8 @@ def check_geometry(initial_gdf): # Obtain simple geometries gdf_simple = initial_gdf.explode(ignore_index=True) # Delete duplicate geoemtries if any - gdf_without_duplicates = gdf_simple.drop_duplicates(ignore_index=True) - # Identify invalid geometries and keep only valid ones - gdf_valid = gdf_without_duplicates.copy() - gdf_valid.geometry = gdf_valid.geometry.apply(lambda geom: fix_invalid_geometry(geom)) + gdf_simple = gdf_simple.drop_duplicates(ignore_index=True) + # Identify invalid geometries, then returns a GeoSeries with valid geometrie + gdf_valid = gdf_simple.make_valid() + return gdf_valid diff --git a/lidro/merge_mask_hydro/vectors/close_holes.py b/lidro/merge_mask_hydro/vectors/close_holes.py new file mode 100644 index 00000000..6f43cf4c --- /dev/null +++ b/lidro/merge_mask_hydro/vectors/close_holes.py @@ -0,0 +1,45 @@ +# -*- coding: utf-8 -*- +""" Remove small holes """ +from shapely.geometry import MultiPolygon, Polygon +from shapely.ops import unary_union + + +def close_holes(polygon, min_hole_area): + """Remove small holes (surface < 100 m²) + + Args: + - polygon (GeoJSON): Hydro Mask geometry + - minhole_area (int): close holes in Mask Hydro : keep only holes with area bigger + than min_hole_area (> 100 m² by default) + + Returns: + GeoJSON: Hydro Mask geometry without holes (< 100 m²) + """ + # Exterior Hydro Mask + exterior_contours = [exterior for exterior in polygon.exterior] + # Create polygons from each exterior contour without holes + polygons = [Polygon(ring) for ring in exterior_contours] + # Merge all polygons into a single encompassing shape + unified_shape = unary_union(polygons) + + # Interior Hydro Mask + interiors_to_keep = [ + interior for poly in polygon for interior in poly.interiors if Polygon(interior).area >= min_hole_area + ] + interior_holes_filter = [[pt for pt in inner.coords] for inner in interiors_to_keep] + + # Determine the appropriate Polygon to use based on the type of unified_shape + if isinstance(unified_shape, Polygon): + final_polygon = Polygon( + shell=unified_shape.exterior, holes=interior_holes_filter if interior_holes_filter else [] + ) + elif isinstance(unified_shape, MultiPolygon): + # Explicitly find the largest polygon by area + largest_polygon = max(unified_shape.geoms, key=lambda p: p.area) + final_polygon = Polygon( + shell=largest_polygon.exterior, holes=interior_holes_filter if interior_holes_filter else [] + ) + else: + raise ValueError("The unified shape is neither a Polygon nor a MultiPolygon. Please handle accordingly.") + + return final_polygon diff --git a/lidro/merge_mask_hydro/vectors/merge_vector.py b/lidro/merge_mask_hydro/vectors/merge_vector.py index 6bf4409e..5f6779ed 100644 --- a/lidro/merge_mask_hydro/vectors/merge_vector.py +++ b/lidro/merge_mask_hydro/vectors/merge_vector.py @@ -4,76 +4,66 @@ import os import geopandas as gpd -from shapely.geometry import MultiPolygon from shapely.ops import unary_union from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( + apply_buffers_to_geometry, check_geometry, - simplify_geometry, ) -from lidro.merge_mask_hydro.vectors.remove_hole import remove_hole +from lidro.merge_mask_hydro.vectors.close_holes import close_holes def merge_geom( input_folder: str, output_folder: str, crs: str, - water_area: int, + min_water_area: int, buffer_positive: float, buffer_negative: float, tolerance: float, ): - """Merge severals masks of hydro surfaces from the points classification of the input LAS/LAZ file, - filter mask (keep only water's area > 150 m²) and save it as a GeoJSON file. + """Merge several masks of hydro surfaces from the points classification + of the input LAS/LAZ files from input_folder, + filter mask (keep only water's area > min_water_area) and save it as a GeoJSON file. Args: - input_folder (str): folder who contains severals Masks Hydro + input_folder (str): folder which contains several geojson geometries output_folder (str): output folder crs (str): a pyproj CRS object used to create the output GeoJSON file - water_area (int): filter Mask Hydro : keep only water's area (> 150 m² by default) - buffer_positive (int): positive buffer from Mask Hydro - buffer_negative (int): negative buffer from Mask Hydro + min_water_area (int): filter Mask Hydro : keep only water objects with area bigger + than min_water_area (> 150 m² by default) + buffer_positive (int): positive buffer to apply to the mask + buffer_negative (int): negative buffer to apply to the mask tolerance (float): All parts of a simplified geometry will be no more than tolerance distance from the original. """ - # List for stocking all GeoDataFrame - polys = [] - # Browse all files in folder - for fichier in os.listdir(input_folder): - if fichier.endswith(".GeoJSON"): - # Load each GeoJSON file as GeoDataFrame - geojson = gpd.read_file(os.path.join(input_folder, fichier)) - merge_geom = geojson["geometry"].unary_union - # Add the GeoDataFrame to the list - polys.append(merge_geom) + gdf = [ + gpd.read_file(os.path.join(input_folder, file), crs=crs).unary_union + for file in os.listdir(input_folder) + if file.endswith(".GeoJSON") + ] # Union geometry - mergedPolys = unary_union(polys) - - geometry = gpd.GeoSeries(mergedPolys, crs=crs).explode(index_parts=False) + gdf = unary_union(gdf) + gdf = gpd.GeoSeries(gdf, crs=crs).explode(index_parts=False) # keep only water's area (> 150 m² by default) - filter_geometry = [geom for geom in geometry if geom.area > water_area] - gdf_filter = gpd.GeoDataFrame(geometry=filter_geometry, crs=crs) + gdf = gdf[gdf.geometry.area > min_water_area] - # Check and rectify the invalid geometry - gdf = check_geometry(gdf_filter) + # Correct geometric errors: simplify certain shapes to make calculations easier + gdf = apply_buffers_to_geometry(gdf, buffer_positive, buffer_negative).explode(index_parts=False) + gdf = gdf.simplify(tolerance=tolerance, preserve_topology=True) - # Correcting geometric errors: simplifying certain shapes to make calculations easier - buffer_geom = simplify_geometry(gdf["geometry"], buffer_positive, buffer_negative).explode(index_parts=False) - simplify_geom = buffer_geom.simplify(tolerance=tolerance, preserve_topology=True) + # Check and rectify the invalid geometry + gdf = check_geometry(gdf) # Correction of holes (< 100m²) in Hydrological Masks - list_parts = [] - for poly in simplify_geom: - list_parts.append(poly) + gdf = close_holes(gdf, min_hole_area=100) + gdf = gpd.GeoSeries(gdf, crs=crs).explode(index_parts=False) - not_hole_geom = remove_hole(MultiPolygon(list_parts)) - geometry = gpd.GeoSeries(not_hole_geom.geoms, crs=crs).explode(index_parts=False) # keep only water's area (> 150 m² by default) : - filter_geometry_second = [geom for geom in geometry if geom.area > water_area] - gdf = gpd.GeoDataFrame(geometry=filter_geometry_second, crs=crs) + gdf = gdf[gdf.geometry.area > min_water_area] # save the result gdf.to_file(os.path.join(output_folder, "MaskHydro_merge.geojson"), driver="GeoJSON", crs=crs) diff --git a/lidro/merge_mask_hydro/vectors/remove_hole.py b/lidro/merge_mask_hydro/vectors/remove_hole.py deleted file mode 100644 index 420acbbd..00000000 --- a/lidro/merge_mask_hydro/vectors/remove_hole.py +++ /dev/null @@ -1,32 +0,0 @@ -# -*- coding: utf-8 -*- -""" Remove small holes """ -from shapely.geometry import MultiPolygon, Polygon - - -def remove_hole(multipoly): - """Remove small holes (surface < 100 m²) - - Args: - - multipoly (GeoJSON): Hydro Mask geometry - - Returns: - GeoJSON: Hydro Mask geometry without holes (< 100 m²) - """ - list_parts = [] - eps = 100 - - for polygon in multipoly.geoms: - list_interiors = [] - - for interior in polygon.interiors: - p = Polygon(interior) - - if p.area > eps: - list_interiors.append(interior) - - temp_pol = Polygon(polygon.exterior.coords, holes=list_interiors) - list_parts.append(temp_pol) - - new_multipolygon = MultiPolygon(list_parts) - - return new_multipolygon diff --git a/test/test_main_merge_mask.py b/test/test_main_merge_mask.py index 3f627249..a4d8beb3 100644 --- a/test/test_main_merge_mask.py +++ b/test/test_main_merge_mask.py @@ -27,7 +27,7 @@ def test_main_run_okay(): def test_main_lidro_fail_no_input_dir(): output_dir = OUTPUT_DIR / "main_no_input_dir" pixel_size = 1 - water_area = 150 + min_water_area = 150 buffer_positive = 0.5 buffer_negative = -1.5 tolerance = 1 @@ -39,7 +39,7 @@ def test_main_lidro_fail_no_input_dir(): overrides=[ f"io.output_dir={output_dir}", f"io.pixel_size={pixel_size}", - f"vector.water_area={water_area}", + f"vector.min_water_area={min_water_area}", f"vector.buffer_positive={buffer_positive}", f"vector.buffer_negative={buffer_negative}", f"vector.tolerance={tolerance}", @@ -53,7 +53,7 @@ def test_main_lidro_fail_no_input_dir(): def test_main_lidro_fail_wrong_input_dir(): output_dir = OUTPUT_DIR / "main_wrong_input_dir" pixel_size = 1 - water_area = 150 + min_water_area = 150 buffer_positive = 0.5 buffer_negative = -1.5 tolerance = 1 @@ -66,7 +66,7 @@ def test_main_lidro_fail_wrong_input_dir(): "io.input_dir=some_random_input_dir", f"io.output_dir={output_dir}", f"io.pixel_size={pixel_size}", - f"vector.water_area={water_area}", + f"vector.min_water_area={min_water_area}", f"vector.buffer_positive={buffer_positive}", f"vector.buffer_negative={buffer_negative}", f"vector.tolerance={tolerance}", @@ -80,7 +80,7 @@ def test_main_lidro_fail_wrong_input_dir(): def test_main_lidro_fail_no_output_dir(): input_dir = INPUT_DIR pixel_size = 1 - water_area = 150 + min_water_area = 150 buffer_positive = 0.5 buffer_negative = -1.5 tolerance = 1 @@ -92,7 +92,7 @@ def test_main_lidro_fail_no_output_dir(): overrides=[ f"io.input_dir={input_dir}", f"io.pixel_size={pixel_size}", - f"vector.water_area={water_area}", + f"vector.min_water_area={min_water_area}", f"vector.buffer_positive={buffer_positive}", f"vector.buffer_negative={buffer_negative}", f"vector.tolerance={tolerance}", diff --git a/test/vectors/test_check_rectify_geometry.py b/test/vectors/test_check_rectify_geometry.py index e010c657..3b6a77eb 100644 --- a/test/vectors/test_check_rectify_geometry.py +++ b/test/vectors/test_check_rectify_geometry.py @@ -2,19 +2,18 @@ from shapely.geometry import MultiPolygon from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( + apply_buffers_to_geometry, check_geometry, - fix_invalid_geometry, - simplify_geometry, ) input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" -def test_simplify_geometry_default(): +def test_apply_buffers_to_geometry_default(): # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) geometry = geojson["geometry"].unary_union - buffer = simplify_geometry(geometry, 1, -1) + buffer = apply_buffers_to_geometry(geometry, 1, -1) assert isinstance(buffer, MultiPolygon) @@ -23,11 +22,3 @@ def test_check_geometry_default(): geojson = gpd.read_file(input) check_geom = check_geometry(geojson) assert isinstance(check_geom.dtypes, object) - - -def test_fix_invalid_geometry_default(): - # Load each GeoJSON file as GeoDataFrame - geojson = gpd.read_file(input) - valid_geom = geojson.geometry.apply(lambda geom: fix_invalid_geometry(geom)) - assert isinstance(valid_geom.dtypes, object) - assert valid_geom[2] == geojson.geometry[2] diff --git a/test/vectors/test_remove_hole.py b/test/vectors/test_close_holes.py similarity index 67% rename from test/vectors/test_remove_hole.py rename to test/vectors/test_close_holes.py index 2c066927..b939320b 100644 --- a/test/vectors/test_remove_hole.py +++ b/test/vectors/test_close_holes.py @@ -1,15 +1,15 @@ import geopandas as gpd from shapely.geometry import MultiPolygon -from lidro.merge_mask_hydro.vectors.remove_hole import remove_hole +from lidro.merge_mask_hydro.vectors.close_holes import close_holes input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" -def test_remove_hole_default(): +def test_close_holes_default(): # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) geometry = geojson["geometry"].unary_union - mask_without_hole = remove_hole(geometry) + mask_without_hole = close_holes(geometry, 100) assert isinstance(mask_without_hole, MultiPolygon) From a2670f5a3e5842723186b4e1d58204ad3c8732cc Mon Sep 17 00:00:00 2001 From: mdupays Date: Tue, 7 May 2024 16:23:22 +0200 Subject: [PATCH 18/33] refacto test/vectors --- test/vectors/test_check_rectify_geometry.py | 14 ++++++ test/vectors/test_close_holes.py | 6 +-- test/vectors/test_convert_to_vector.py | 55 +++++++-------------- test/vectors/test_merge_vector.py | 54 +++++++------------- 4 files changed, 53 insertions(+), 76 deletions(-) diff --git a/test/vectors/test_check_rectify_geometry.py b/test/vectors/test_check_rectify_geometry.py index 3b6a77eb..907460d8 100644 --- a/test/vectors/test_check_rectify_geometry.py +++ b/test/vectors/test_check_rectify_geometry.py @@ -21,4 +21,18 @@ def test_check_geometry_default(): # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) check_geom = check_geometry(geojson) + assert isinstance(check_geom.dtypes, object) + + assert geojson.geometry.dtype == "geometry" + + +def test_check_geometry_integrity(): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(input) + + # duplicates in the data + assert not geojson.duplicated().any(), "There are duplicates in the data" + + # Check geometry + assert geojson["geometry"].is_valid.all(), "Geometry no-valid" diff --git a/test/vectors/test_close_holes.py b/test/vectors/test_close_holes.py index b939320b..8d4ffcd5 100644 --- a/test/vectors/test_close_holes.py +++ b/test/vectors/test_close_holes.py @@ -1,5 +1,5 @@ import geopandas as gpd -from shapely.geometry import MultiPolygon +from shapely.geometry import Polygon from lidro.merge_mask_hydro.vectors.close_holes import close_holes @@ -9,7 +9,7 @@ def test_close_holes_default(): # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) - geometry = geojson["geometry"].unary_union + geometry = geojson["geometry"] mask_without_hole = close_holes(geometry, 100) - assert isinstance(mask_without_hole, MultiPolygon) + assert isinstance(mask_without_hole, Polygon) diff --git a/test/vectors/test_convert_to_vector.py b/test/vectors/test_convert_to_vector.py index 8678c381..a196b834 100644 --- a/test/vectors/test_convert_to_vector.py +++ b/test/vectors/test_convert_to_vector.py @@ -1,15 +1,17 @@ -import json import os import shutil from pathlib import Path +import geopandas as gpd +from pyproj import CRS +from shapely.geometry import Polygon + from lidro.create_mask_hydro.vectors.convert_to_vector import create_hydro_vector_mask TMP_PATH = Path("./tmp/create_mask_hydro/vectors/convert_to_vector") las_file = "./data/pointcloud/Semis_2021_0830_6291_LA93_IGN69.laz" output = "./tmp/create_mask_hydro/vectors/convert_to_vector/MaskHydro_Semis_2021_0830_6291_LA93_IGN69.GeoJSON" -output_main = "./tmp/create_mask_hydro/main/main_lidro_default/MaskHydro_Semis_2021_0830_6291_LA93_IGN69.GeoJSON" def setup_module(module): @@ -19,42 +21,21 @@ def setup_module(module): def test_create_hydro_vector_mask_default(): - crs = 'PROJCS["RGF93 v1 / Lambert-93",GEOGCS["RGF93 v1",DATUM["Reseau_Geodesique_Francais_1993_v1",\ - SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],\ - AUTHORITY["EPSG","6171"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],\ - UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4171"]],\ - PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",46.5],\ - PARAMETER["central_meridian",3],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",44],\ - PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],\ - UNIT["metre",1,AUTHORITY["EPSG","9001"]],\ - AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2154"]]' - - create_hydro_vector_mask(las_file, output, 1, 1000, [0, 1, 2, 3, 4, 5, 6, 17, 66], crs, 3) - + # Parameters + pixel_size = 1 + tile_size = 1000 + classes = [0, 1, 2, 3, 4, 5, 6, 17, 66] + crs = CRS.from_epsg(2154) + dilatation_size = 3 + + create_hydro_vector_mask(las_file, output, pixel_size, tile_size, classes, crs, dilatation_size) assert Path(output).exists() + gdf = gpd.read_file(output) -def test_check_structure_default(): - # Output - with open(output, "r") as f: - geojson_data = json.load(f) - - with open(output_main, "r") as f: - geojson_data_main = json.load(f) - - # CHECK STRUCTURE - assert "type" in geojson_data - assert geojson_data["type"] == "FeatureCollection" - assert "features" in geojson_data - assert isinstance(geojson_data["features"], list) - - # CHECK POLYGON - for feature in geojson_data["features"]: - geometry = feature["geometry"] - coordinates = geometry["coordinates"] - - for feature in geojson_data_main["features"]: - geometry_main = feature["geometry"] - coordinates_main = geometry_main["coordinates"] + assert not gdf.empty # GeoDataFrame shouldn't empty + assert gdf.crs.to_string() == crs # CRS is identical + assert all(isinstance(geom, Polygon) for geom in gdf.geometry) # All geometry should Polygons - assert coordinates[0] == coordinates_main[0] + expected_number_of_geometries = 2820 + assert len(gdf) == expected_number_of_geometries # the number of geometries must be identical diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py index eee73225..4234a8ca 100644 --- a/test/vectors/test_merge_vector.py +++ b/test/vectors/test_merge_vector.py @@ -1,8 +1,11 @@ -import json import os import shutil from pathlib import Path +import geopandas as gpd +from pyproj import CRS +from shapely.geometry import Polygon + from lidro.merge_mask_hydro.vectors.merge_vector import merge_geom TMP_PATH = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro") @@ -19,42 +22,21 @@ def setup_module(module): def test_create_hydro_vector_mask_default(): - crs = 'PROJCS["RGF93 v1 / Lambert-93",GEOGCS["RGF93 v1",DATUM["Reseau_Geodesique_Francais_1993_v1",\ - SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],\ - AUTHORITY["EPSG","6171"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],\ - UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4171"]],\ - PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",46.5],\ - PARAMETER["central_meridian",3],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",44],\ - PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],\ - UNIT["metre",1,AUTHORITY["EPSG","9001"]],\ - AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","2154"]]' - - merge_geom(input_folder, TMP_PATH, crs, 150, 0.5, -1.5, 1) - + # Parameters + min_water_area = 150 + buffer_positive = 0.5 + buffer_negative = -1.5 + tolerance = 1 + crs = CRS.from_epsg(2154) + + merge_geom(input_folder, TMP_PATH, crs, min_water_area, buffer_positive, buffer_negative, tolerance) assert Path(output).exists() + gdf = gpd.read_file(output) -def test_check_structure_default(): - # Output - with open(output, "r") as f: - geojson_data = json.load(f) - - with open(output_main, "r") as f: - geojson_data_main = json.load(f) - - # CHECK STRUCTURE - assert "type" in geojson_data - assert geojson_data["type"] == "FeatureCollection" - assert "features" in geojson_data - assert isinstance(geojson_data["features"], list) - - # CHECK POLYGON - for feature in geojson_data["features"]: - geometry = feature["geometry"] - coordinates = geometry["coordinates"] - - for feature in geojson_data_main["features"]: - geometry_main = feature["geometry"] - coordinates_main = geometry_main["coordinates"] + assert not gdf.empty # GeoDataFrame shouldn't empty + assert gdf.crs.to_string() == crs # CRS is identical + assert all(isinstance(geom, Polygon) for geom in gdf.geometry) # All geometry should Polygons - assert coordinates[0] == coordinates_main[0] + expected_number_of_geometries = 1 + assert len(gdf) == expected_number_of_geometries # One geometry From 31e9bbd6b917835c3d83138e4cedd9f92941901c Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 13 May 2024 10:51:30 +0200 Subject: [PATCH 19/33] update configs for classes --- configs/configs_lidro.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 9963df7b..79a19ebb 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -32,7 +32,7 @@ vector: filter: # Classes to be considered as "non-water" - keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 65, 66, 67] + keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes hydra: output_subdir: null From c753c1aaa4a6b536ef9244207c7d96cbdecdb983 Mon Sep 17 00:00:00 2001 From: mdupaysign <162304206+mdupaysign@users.noreply.github.com> Date: Mon, 13 May 2024 10:52:37 +0200 Subject: [PATCH 20/33] Update lidro/merge_mask_hydro/vectors/check_rectify_geometry.py Co-authored-by: leavauchier <120112647+leavauchier@users.noreply.github.com> --- lidro/merge_mask_hydro/vectors/check_rectify_geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py index 11452120..4fe431e7 100644 --- a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -9,7 +9,7 @@ def apply_buffers_to_geometry(hydro_mask, buffer_positive: float, buffer_negativ Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks Args: - s (gpd.GeoDataFrame): geopandas dataframe with input geometry + hydro_mask (gpd.GeoDataFrame): geopandas dataframe with input geometry buffer_positive (float): positive buffer to apply to the input geometry buffer_negative (float): negative buffer to apply to the input geometry Returns: From 339866a786107303899d5e30b682cd4ee274fd9b Mon Sep 17 00:00:00 2001 From: mdupaysign <162304206+mdupaysign@users.noreply.github.com> Date: Mon, 13 May 2024 10:53:35 +0200 Subject: [PATCH 21/33] Update lidro/merge_mask_hydro/vectors/check_rectify_geometry.py Co-authored-by: leavauchier <120112647+leavauchier@users.noreply.github.com> --- lidro/merge_mask_hydro/vectors/check_rectify_geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py index 4fe431e7..d6c325d3 100644 --- a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -34,7 +34,7 @@ def check_geometry(initial_gdf): gdf_simple = initial_gdf.explode(ignore_index=True) # Delete duplicate geoemtries if any gdf_simple = gdf_simple.drop_duplicates(ignore_index=True) - # Identify invalid geometries, then returns a GeoSeries with valid geometrie + # Identify invalid geometries, then return a GeoSeries with valid geometries gdf_valid = gdf_simple.make_valid() return gdf_valid From c7eb87c609a17a27ca320577b21c3c20ac4bfec5 Mon Sep 17 00:00:00 2001 From: mdupaysign <162304206+mdupaysign@users.noreply.github.com> Date: Mon, 13 May 2024 11:39:37 +0200 Subject: [PATCH 22/33] Update lidro/merge_mask_hydro/vectors/close_holes.py Co-authored-by: leavauchier <120112647+leavauchier@users.noreply.github.com> --- lidro/merge_mask_hydro/vectors/close_holes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/merge_mask_hydro/vectors/close_holes.py b/lidro/merge_mask_hydro/vectors/close_holes.py index 6f43cf4c..3eb3f57b 100644 --- a/lidro/merge_mask_hydro/vectors/close_holes.py +++ b/lidro/merge_mask_hydro/vectors/close_holes.py @@ -18,7 +18,7 @@ def close_holes(polygon, min_hole_area): # Exterior Hydro Mask exterior_contours = [exterior for exterior in polygon.exterior] # Create polygons from each exterior contour without holes - polygons = [Polygon(ring) for ring in exterior_contours] + exterior_polygons = [Polygon(ring) for ring in exterior_contours] # Merge all polygons into a single encompassing shape unified_shape = unary_union(polygons) From 06b86df87d52fad0cda69f6c8b2586746feb5da3 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 13 May 2024 12:09:57 +0200 Subject: [PATCH 23/33] delete parameters of function merge_mask --- lidro/main_merge_mask.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/lidro/main_merge_mask.py b/lidro/main_merge_mask.py index be6dfcae..350625de 100644 --- a/lidro/main_merge_mask.py +++ b/lidro/main_merge_mask.py @@ -36,13 +36,6 @@ def main(config: DictConfig): os.makedirs(output_dir, exist_ok=True) - # Parameters for merging Mask Hydro - min_water_area = config.vector.min_water_area # keep only water's area (> 150 m² by default) - buffer_positive = config.vector.buffer_positive # positive buffer from Mask Hydro - buffer_negative = config.vector.buffer_negative # negative buffer from Mask Hydro - tolerance = config.vector.tolerance # Tolerance from Douglas-Peucker - crs = CRS.from_user_input(config.io.srid) - if os.path.isdir(input_dir): os.makedirs(output_dir, exist_ok=True) # Create folder "merge" # Merge all Mash Hydro From d956c289b1610da8c5b17325570c4ef3d716b950 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 13 May 2024 12:23:22 +0200 Subject: [PATCH 24/33] refacto check_rectify_geometry --- lidro/merge_mask_hydro/vectors/check_rectify_geometry.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py index d6c325d3..28fe815c 100644 --- a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -32,8 +32,9 @@ def check_geometry(initial_gdf): """ # Obtain simple geometries gdf_simple = initial_gdf.explode(ignore_index=True) - # Delete duplicate geoemtries if any + gdf_simple = gdf_simple.drop_duplicates(ignore_index=True) + # Identify invalid geometries, then return a GeoSeries with valid geometries gdf_valid = gdf_simple.make_valid() From d710577de823f8455409ee808475d0856a032045 Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 13 May 2024 12:23:59 +0200 Subject: [PATCH 25/33] refacto: delete output_main --- test/vectors/test_merge_vector.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py index 4234a8ca..7facd361 100644 --- a/test/vectors/test_merge_vector.py +++ b/test/vectors/test_merge_vector.py @@ -11,8 +11,7 @@ TMP_PATH = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro") input_folder = "./data/mask_hydro/" -output = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro/MaskHydro_merge.geojson") -output_main = "./data/merge_mask_hydro/MaskHydro_merge.geojson" +output = Path("./data/merge_mask_hydro/MaskHydro_merge.geojson") def setup_module(module): From de187614ef5f8f6b11fd37f120b47fd751312dfd Mon Sep 17 00:00:00 2001 From: mdupays Date: Mon, 13 May 2024 13:40:46 +0200 Subject: [PATCH 26/33] modify test close_holes and merge_vector --- lidro/merge_mask_hydro/vectors/close_holes.py | 2 +- test/vectors/test_merge_vector.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/lidro/merge_mask_hydro/vectors/close_holes.py b/lidro/merge_mask_hydro/vectors/close_holes.py index 3eb3f57b..cec8047c 100644 --- a/lidro/merge_mask_hydro/vectors/close_holes.py +++ b/lidro/merge_mask_hydro/vectors/close_holes.py @@ -20,7 +20,7 @@ def close_holes(polygon, min_hole_area): # Create polygons from each exterior contour without holes exterior_polygons = [Polygon(ring) for ring in exterior_contours] # Merge all polygons into a single encompassing shape - unified_shape = unary_union(polygons) + unified_shape = unary_union(exterior_polygons) # Interior Hydro Mask interiors_to_keep = [ diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py index 7facd361..5f160cf5 100644 --- a/test/vectors/test_merge_vector.py +++ b/test/vectors/test_merge_vector.py @@ -37,5 +37,5 @@ def test_create_hydro_vector_mask_default(): assert gdf.crs.to_string() == crs # CRS is identical assert all(isinstance(geom, Polygon) for geom in gdf.geometry) # All geometry should Polygons - expected_number_of_geometries = 1 + expected_number_of_geometries = 3 assert len(gdf) == expected_number_of_geometries # One geometry From 89a18152f0428e2fec4d47bde8862b42a836b0f7 Mon Sep 17 00:00:00 2001 From: mdupaysign <162304206+mdupaysign@users.noreply.github.com> Date: Tue, 14 May 2024 09:38:52 +0200 Subject: [PATCH 27/33] Update lidro/merge_mask_hydro/vectors/merge_vector.py Co-authored-by: leavauchier <120112647+leavauchier@users.noreply.github.com> --- lidro/merge_mask_hydro/vectors/merge_vector.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lidro/merge_mask_hydro/vectors/merge_vector.py b/lidro/merge_mask_hydro/vectors/merge_vector.py index 5f6779ed..80bf1d61 100644 --- a/lidro/merge_mask_hydro/vectors/merge_vector.py +++ b/lidro/merge_mask_hydro/vectors/merge_vector.py @@ -62,7 +62,7 @@ def merge_geom( gdf = close_holes(gdf, min_hole_area=100) gdf = gpd.GeoSeries(gdf, crs=crs).explode(index_parts=False) - # keep only water's area (> 150 m² by default) : + # filter out water area < min_water_area (150 m² by default) again to make sure that previous geometry updates did not generate new small water areas gdf = gdf[gdf.geometry.area > min_water_area] # save the result From df974059322c76b3ddff2f6d21479bb027f8c433 Mon Sep 17 00:00:00 2001 From: mdupays Date: Tue, 14 May 2024 14:55:52 +0200 Subject: [PATCH 28/33] refacto with Lea --- configs/configs_lidro.yaml | 2 +- .../vectors/check_rectify_geometry.py | 15 ++++++----- lidro/merge_mask_hydro/vectors/close_holes.py | 26 ++++--------------- .../merge_mask_hydro/vectors/merge_vector.py | 11 ++++---- test/vectors/test_check_rectify_geometry.py | 12 +++------ test/vectors/test_close_holes.py | 2 +- test/vectors/test_merge_vector.py | 2 +- 7 files changed, 26 insertions(+), 44 deletions(-) diff --git a/configs/configs_lidro.yaml b/configs/configs_lidro.yaml index 79a19ebb..4eee3a7b 100644 --- a/configs/configs_lidro.yaml +++ b/configs/configs_lidro.yaml @@ -26,7 +26,7 @@ vector: min_water_area: 150 # Parameters for buffer buffer_positive: 1 - buffer_negative: -0.5 + buffer_negative: -1.5 # negative buffer should be bigger than positive buffer to prevent protruding over the banks # Tolerance from Douglas-Peucker tolerance: 1 diff --git a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py index 28fe815c..12852b3a 100644 --- a/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py +++ b/lidro/merge_mask_hydro/vectors/check_rectify_geometry.py @@ -2,11 +2,12 @@ """ Rectify and check geometry """ from shapely.geometry import CAP_STYLE - +import geopandas as gpd def apply_buffers_to_geometry(hydro_mask, buffer_positive: float, buffer_negative: float): """Buffer geometry - Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks + Objective: create a HYDRO mask at the edge of the bank, not protruding over the banks. + In this case, negative buffer should be bigger than positive buffer Args: hydro_mask (gpd.GeoDataFrame): geopandas dataframe with input geometry @@ -21,15 +22,15 @@ def apply_buffers_to_geometry(hydro_mask, buffer_positive: float, buffer_negativ return geom -def check_geometry(initial_gdf): - """Check topology +def fix_topology(initial_gdf: gpd.GeoDataFrame)-> gpd.GeoDataFrame: + """Fix topology Args: - initial_gdf (GeoJSON): Hydro Mask geometry + initial_gdf (gpd.GeoDataFrame): Hydro Mask geometry Returns: - GeoJSON: Hydro Mask geometry valid - """ + gpd.GeoDataFrame: Hydro Mask geometry valid + """ # Obtain simple geometries gdf_simple = initial_gdf.explode(ignore_index=True) diff --git a/lidro/merge_mask_hydro/vectors/close_holes.py b/lidro/merge_mask_hydro/vectors/close_holes.py index cec8047c..26626c76 100644 --- a/lidro/merge_mask_hydro/vectors/close_holes.py +++ b/lidro/merge_mask_hydro/vectors/close_holes.py @@ -4,7 +4,7 @@ from shapely.ops import unary_union -def close_holes(polygon, min_hole_area): +def close_holes(polygon: Polygon, min_hole_area): """Remove small holes (surface < 100 m²) Args: @@ -15,31 +15,15 @@ def close_holes(polygon, min_hole_area): Returns: GeoJSON: Hydro Mask geometry without holes (< 100 m²) """ - # Exterior Hydro Mask - exterior_contours = [exterior for exterior in polygon.exterior] - # Create polygons from each exterior contour without holes - exterior_polygons = [Polygon(ring) for ring in exterior_contours] - # Merge all polygons into a single encompassing shape - unified_shape = unary_union(exterior_polygons) - # Interior Hydro Mask interiors_to_keep = [ - interior for poly in polygon for interior in poly.interiors if Polygon(interior).area >= min_hole_area + interior for interior in polygon.interiors if Polygon(interior).area >= min_hole_area ] interior_holes_filter = [[pt for pt in inner.coords] for inner in interiors_to_keep] # Determine the appropriate Polygon to use based on the type of unified_shape - if isinstance(unified_shape, Polygon): - final_polygon = Polygon( - shell=unified_shape.exterior, holes=interior_holes_filter if interior_holes_filter else [] - ) - elif isinstance(unified_shape, MultiPolygon): - # Explicitly find the largest polygon by area - largest_polygon = max(unified_shape.geoms, key=lambda p: p.area) - final_polygon = Polygon( - shell=largest_polygon.exterior, holes=interior_holes_filter if interior_holes_filter else [] - ) - else: - raise ValueError("The unified shape is neither a Polygon nor a MultiPolygon. Please handle accordingly.") + final_polygon = Polygon( + shell=polygon.exterior, holes=interior_holes_filter) + return final_polygon diff --git a/lidro/merge_mask_hydro/vectors/merge_vector.py b/lidro/merge_mask_hydro/vectors/merge_vector.py index 80bf1d61..88019203 100644 --- a/lidro/merge_mask_hydro/vectors/merge_vector.py +++ b/lidro/merge_mask_hydro/vectors/merge_vector.py @@ -8,7 +8,7 @@ from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( apply_buffers_to_geometry, - check_geometry, + fix_topology, ) from lidro.merge_mask_hydro.vectors.close_holes import close_holes @@ -56,13 +56,14 @@ def merge_geom( gdf = gdf.simplify(tolerance=tolerance, preserve_topology=True) # Check and rectify the invalid geometry - gdf = check_geometry(gdf) + gdf = fix_topology(gdf) # Correction of holes (< 100m²) in Hydrological Masks - gdf = close_holes(gdf, min_hole_area=100) - gdf = gpd.GeoSeries(gdf, crs=crs).explode(index_parts=False) + gs = gdf.geometry.apply(lambda p: close_holes(p, min_hole_area=100)) + gdf = gpd.GeoDataFrame(geometry=gs, crs=crs) - # filter out water area < min_water_area (150 m² by default) again to make sure that previous geometry updates did not generate new small water areas + # filter out water area < min_water_area (150 m² by default) again to make sure + # that previous geometry updates did not generate new small water areas gdf = gdf[gdf.geometry.area > min_water_area] # save the result diff --git a/test/vectors/test_check_rectify_geometry.py b/test/vectors/test_check_rectify_geometry.py index 907460d8..ff62075c 100644 --- a/test/vectors/test_check_rectify_geometry.py +++ b/test/vectors/test_check_rectify_geometry.py @@ -3,7 +3,7 @@ from lidro.merge_mask_hydro.vectors.check_rectify_geometry import ( apply_buffers_to_geometry, - check_geometry, + fix_topology, ) input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" @@ -17,22 +17,18 @@ def test_apply_buffers_to_geometry_default(): assert isinstance(buffer, MultiPolygon) -def test_check_geometry_default(): +def test_fix_topology_default(): # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) - check_geom = check_geometry(geojson) + check_geom = fix_topology(geojson) assert isinstance(check_geom.dtypes, object) assert geojson.geometry.dtype == "geometry" - -def test_check_geometry_integrity(): - # Load each GeoJSON file as GeoDataFrame - geojson = gpd.read_file(input) - # duplicates in the data assert not geojson.duplicated().any(), "There are duplicates in the data" # Check geometry assert geojson["geometry"].is_valid.all(), "Geometry no-valid" + diff --git a/test/vectors/test_close_holes.py b/test/vectors/test_close_holes.py index 8d4ffcd5..b77ebb4d 100644 --- a/test/vectors/test_close_holes.py +++ b/test/vectors/test_close_holes.py @@ -9,7 +9,7 @@ def test_close_holes_default(): # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) - geometry = geojson["geometry"] + geometry = geojson["geometry"][0] mask_without_hole = close_holes(geometry, 100) assert isinstance(mask_without_hole, Polygon) diff --git a/test/vectors/test_merge_vector.py b/test/vectors/test_merge_vector.py index 5f160cf5..b0599f1e 100644 --- a/test/vectors/test_merge_vector.py +++ b/test/vectors/test_merge_vector.py @@ -11,7 +11,7 @@ TMP_PATH = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro") input_folder = "./data/mask_hydro/" -output = Path("./data/merge_mask_hydro/MaskHydro_merge.geojson") +output = Path("./tmp/merge_mask_hydro/vectors/merge_mask_hydro/MaskHydro_merge.geojson") def setup_module(module): From bb1e35596f89e1aebc70f4f937a8aadcd33bf67a Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 15 May 2024 10:29:50 +0200 Subject: [PATCH 29/33] add test with geometry no valid --- test/vectors/test_check_rectify_geometry.py | 25 ++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/test/vectors/test_check_rectify_geometry.py b/test/vectors/test_check_rectify_geometry.py index ff62075c..c2dfa244 100644 --- a/test/vectors/test_check_rectify_geometry.py +++ b/test/vectors/test_check_rectify_geometry.py @@ -7,6 +7,7 @@ ) input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" +input_error = "./data/merge_mask_hydro/MaskHydro_merge_NoValid.geojson" def test_apply_buffers_to_geometry_default(): @@ -24,11 +25,25 @@ def test_fix_topology_default(): assert isinstance(check_geom.dtypes, object) - assert geojson.geometry.dtype == "geometry" + assert check_geom.geometry.dtype == "geometry" - # duplicates in the data - assert not geojson.duplicated().any(), "There are duplicates in the data" + # Check not duplicates in the data + assert not check_geom.duplicated().any(), "There are duplicates in the data" - # Check geometry - assert geojson["geometry"].is_valid.all(), "Geometry no-valid" + # # Check geometry + assert check_geom.geometry.is_valid.all(), "Geometry no-valid" +def test_fix_topology_error(): + # Load each GeoJSON file as GeoDataFrame + geojson = gpd.read_file(input_error) + check_geom = fix_topology(geojson) + + assert isinstance(check_geom.dtypes, object) + + assert check_geom.geometry.dtype == "geometry" + + # Check not duplicates in the data + assert not check_geom.duplicated().any(), "There are duplicates in the data" + + # # Check geometry + assert check_geom.geometry.is_valid.all(), "Geometry no-valid" From cdfdfcdee4d3959581080428648b3999fab185fc Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 15 May 2024 10:40:19 +0200 Subject: [PATCH 30/33] rectify test:check geometry with severals holes --- test/vectors/test_close_holes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/vectors/test_close_holes.py b/test/vectors/test_close_holes.py index b77ebb4d..1376cc04 100644 --- a/test/vectors/test_close_holes.py +++ b/test/vectors/test_close_holes.py @@ -9,7 +9,7 @@ def test_close_holes_default(): # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) - geometry = geojson["geometry"][0] + geometry = geojson["geometry"][1] # with contain severals holes mask_without_hole = close_holes(geometry, 100) assert isinstance(mask_without_hole, Polygon) From cd78e07ca05a4ab0bfac99e448164fbb06aa15c3 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 15 May 2024 10:52:26 +0200 Subject: [PATCH 31/33] lidro-data submodule updated --- data | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data b/data index 36d9aab8..f7ebd371 160000 --- a/data +++ b/data @@ -1 +1 @@ -Subproject commit 36d9aab8ae1e8eeb43449c89c2cbb136a0b1caa2 +Subproject commit f7ebd37180937ba4ceb9f231c0962234654c00e8 From fcf004c11df563be0797fd0103573fe11755c28c Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 15 May 2024 11:10:35 +0200 Subject: [PATCH 32/33] modify input on each test --- test/vectors/test_check_rectify_geometry.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/vectors/test_check_rectify_geometry.py b/test/vectors/test_check_rectify_geometry.py index c2dfa244..01719ffa 100644 --- a/test/vectors/test_check_rectify_geometry.py +++ b/test/vectors/test_check_rectify_geometry.py @@ -6,11 +6,8 @@ fix_topology, ) -input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" -input_error = "./data/merge_mask_hydro/MaskHydro_merge_NoValid.geojson" - - def test_apply_buffers_to_geometry_default(): + input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) geometry = geojson["geometry"].unary_union @@ -19,6 +16,7 @@ def test_apply_buffers_to_geometry_default(): def test_fix_topology_default(): + input = "./data/merge_mask_hydro/MaskHydro_merge.geojson" # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input) check_geom = fix_topology(geojson) @@ -34,6 +32,7 @@ def test_fix_topology_default(): assert check_geom.geometry.is_valid.all(), "Geometry no-valid" def test_fix_topology_error(): + input_error = "./data/merge_mask_hydro/MaskHydro_merge_NoValid.geojson" # Load each GeoJSON file as GeoDataFrame geojson = gpd.read_file(input_error) check_geom = fix_topology(geojson) From cea09f4e707be74b88a43f345c3a14f40ddb3371 Mon Sep 17 00:00:00 2001 From: mdupays Date: Wed, 15 May 2024 11:16:28 +0200 Subject: [PATCH 33/33] rectify docstring --- lidro/merge_mask_hydro/vectors/close_holes.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lidro/merge_mask_hydro/vectors/close_holes.py b/lidro/merge_mask_hydro/vectors/close_holes.py index 26626c76..7b88a910 100644 --- a/lidro/merge_mask_hydro/vectors/close_holes.py +++ b/lidro/merge_mask_hydro/vectors/close_holes.py @@ -1,19 +1,19 @@ # -*- coding: utf-8 -*- """ Remove small holes """ -from shapely.geometry import MultiPolygon, Polygon +from shapely.geometry import Polygon from shapely.ops import unary_union -def close_holes(polygon: Polygon, min_hole_area): +def close_holes(polygon: Polygon, min_hole_area)-> Polygon: """Remove small holes (surface < 100 m²) Args: - - polygon (GeoJSON): Hydro Mask geometry + - polygon (Polygon): Hydro Mask geometry - minhole_area (int): close holes in Mask Hydro : keep only holes with area bigger than min_hole_area (> 100 m² by default) Returns: - GeoJSON: Hydro Mask geometry without holes (< 100 m²) + Polygon: Hydro Mask geometry without holes (< 100 m²) """ # Interior Hydro Mask interiors_to_keep = [