diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 1ff0699..4dc9a34 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -1,3 +1,7 @@ +### 0.3.0 + +Update algorithm of [mark_points_to_use_for_digital_models_with_new_dimension](pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py), + ### 0.2.1 Fix (and test) arguments parsing in [mark_points_to_use_for_digital_models_with_new_dimension](pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py) diff --git a/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py b/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py index 0fe1e19..cdc1417 100755 --- a/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py +++ b/pdal_ign_macro/mark_points_to_use_for_digital_models_with_new_dimension.py @@ -92,21 +92,26 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) "PT_ON_BUILDING", "PT_ON_VEGET", "PT_ON_SOL", + "PT_ON_VIRT", ] added_dimensions = [dtm_dimension, dsm_dimension] + temporary_dimensions pipeline |= pdal.Filter.ferry(dimensions="=>" + ", =>".join(added_dimensions)) - # 1 - recherche des points max de végétation (4,5) sur une grille régulière, avec prise en - # compte des points sol (2) et basse - # vegetation (3) proches de la végétation - # pour le calcul du DSM + ################################################################################################################### + # 1 - Gestion de la végétation pour le calcul du DSM: + ################################################################################################################### + # Recherche des points max de végétation (4,5) sur une grille régulière = PT_VEG_DSM=1, + # Prise en compte des points sol (2) proche de la végétation = PT_VEG_DSM=1 + # Isolement des classes (hors sol) sous la végetation comme les bâtis, l'eau, les ponts et les divers-bâtis + # Prise en compte des points basse vegetation (3) proche de la végétation = PT_VEG_DSM=1 + # Passage des premiers points "veget" PT_VEG_DSM=1 dans le taguage pour le MNS dsm_dimension=1 + # 1.1 Point veget max pipeline |= pdal.Filter.assign( value=["PT_VEG_DSM = 1 WHERE " + macro.build_condition("Classification", [4, 5])] ) - - # bouche trou : assigne les points sol à l'intérieur de la veget (4,5) + # 1.2 bouche trou : assigne les points sol à l'intérieur de la veget (4,5) pipeline = macro.add_radius_assign( pipeline, 1, @@ -115,16 +120,6 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref=macro.build_condition("Classification", [4, 5]), condition_out="PT_VEG_DSM=1", ) - pipeline = macro.add_radius_assign( - pipeline, - 1, - False, - condition_src=macro.build_condition("Classification", [6, 9, 17, 67]), - condition_ref=macro.build_condition("Classification", [4, 5]), - condition_out="PT_ON_VEGET=1", - max2d_above=0, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=-1, # prendre tous les points qui sont en dessous des points pont (condition_ref) - ) pipeline = macro.add_radius_assign( pipeline, 1, @@ -133,38 +128,32 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref="Classification==2 && PT_VEG_DSM==0", condition_out="PT_VEG_DSM=0", ) + # 1.3 Isolement en PT_ON_VEGET=1 des éléments sous la végétation (hors sol) pipeline = macro.add_radius_assign( pipeline, 1, False, - condition_src="PT_ON_VEGET==1 && " + macro.build_condition("Classification", [6, 67]), - condition_ref="PT_ON_VEGET==0 &&" + macro.build_condition("Classification", [6, 67]), - condition_out="PT_ON_VEGET=0", - max2d_above=0.5, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) - ) - pipeline = macro.add_radius_assign( - pipeline, - 1, - False, - condition_src="PT_ON_VEGET==1 && Classification==17", - condition_ref="Classification==17 && PT_ON_VEGET==0", - condition_out="PT_ON_VEGET=0", - max2d_above=0.5, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) + condition_src=macro.build_condition("Classification", [6, 9, 17, 67]), + condition_ref=macro.build_condition("Classification", [4, 5]), + condition_out="PT_ON_VEGET=1", + max2d_above=0, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) + max2d_below=900, # prendre tous les points qui sont en dessous des points pont (condition_ref) ) pipeline = macro.add_radius_assign( pipeline, 1, False, - condition_src="PT_ON_VEGET==1 && Classification==9", - condition_ref="Classification==9 && PT_ON_VEGET==0", + condition_src="PT_ON_VEGET==1 && ( " + + macro.build_condition("Classification", [6, 9, 17, 67]) + + " )", + condition_ref="PT_ON_VEGET==0 && ( " + + macro.build_condition("Classification", [6, 9, 17, 67]) + + " )", condition_out="PT_ON_VEGET=0", max2d_above=0.5, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) ) - - # selection des points de veget basse proche de la veget haute + # 1.4 selection des points de veget basse proche de la veget haute pipeline = macro.add_radius_assign( pipeline, 1, @@ -173,17 +162,19 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref="Classification==5", condition_out="PT_VEG_DSM=1", ) - + # 1.5 Premiers points tagués pour le MNS # max des points de veget (PT_VEG_DSM==1) sur une grille régulière : # TODO: remplacer par GridDecimation une fois le correctif mergé dans PDAL pipeline |= pdal.Filter.grid_decimation_deprecated( resolution=0.75, output_dimension=dsm_dimension, output_type="max", where="PT_VEG_DSM==1" ) + ################################################################################################################### + # 2 - Gestion de l'eau + ################################################################################################################### + # L'eau sous la roche en dévers (typique Corse) + # Gestion de l'eau sur les masques hydro qui ne doivent pas se supperposer aux points virtuels 66 - # 2 - sélection des points pour DTM et DSM - - ########################################################################################################################################### - # Voir pour l'eau sous le sol dans les dévers corse + # 2.1 L'eau sous la roche pipeline = macro.add_radius_assign( pipeline, 1.25, @@ -192,7 +183,7 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref="Classification==2", condition_out="PT_ON_SOL=1", max2d_above=0, - max2d_below=-1, + max2d_below=900, ) pipeline = macro.add_radius_assign( pipeline, @@ -204,29 +195,67 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) max2d_above=0.5, max2d_below=0.5, ) - ########################################################################################################################################### + # 2.2 Gestion de l'eau sur les masques hydro + pipeline = macro.add_radius_assign( + pipeline, + 1, + False, + condition_src="Classification==9", + condition_ref="Classification==66", + condition_out="PT_ON_VIRT=1", + ) + pipeline = macro.add_radius_assign( + pipeline, + 1, + False, + condition_src="PT_ON_VIRT==1 && Classification==9", + condition_ref="PT_ON_VIRT==0 && Classification==9", + condition_out="PT_ON_VIRT=0", + ) + ################################################################################################################### + # 3 - sélection des premiers points pour MNT et MNS + ################################################################################################################### + # Initialisation pour le MNT + # Initialisation pour le MNS - # selection de points DTM (max) sur une grille régulière + # 3.1 Pour le MNT (le point sol max sur uen grille de 50cm) # TODO: remplacer par GridDecimation une fois le correctif mergé dans PDAL pipeline |= pdal.Filter.grid_decimation_deprecated( resolution=0.5, output_dimension=dtm_dimension, output_type="max", - where="(Classification==2 || PT_ON_SOL==0 && Classification==9)", + where="(Classification==2)", ) - - # selection de points DSM (max) sur une grille régulière + # 3.2 Pour les MNS (Pour le moment: Les bâtis, ponts, veget) # TODO: remplacer par GridDecimation une fois le correctif mergé dans PDAL pipeline |= pdal.Filter.grid_decimation_deprecated( resolution=0.5, output_dimension=dsm_dimension, output_type="max", - where="(PT_ON_VEGET==0 && " - + macro.build_condition("Classification", [6, 9, 17, 64]) - + f" || {dsm_dimension}==1)", + where="(PT_ON_VEGET==0 && (" + + macro.build_condition("Classification", [6, 17, 67]) + + f") || {dsm_dimension}==1)", + ) + # 3.3 Pour les points "eau" on prendra le point le plus bas de la grille de 50cm et qui ne sont ni sous la roche ni près de pts virtuels + pipeline |= pdal.Filter.grid_decimation_deprecated( + resolution=0.5, + output_dimension=dtm_dimension, + output_type="min", + where="(PT_ON_SOL==0 && PT_ON_VIRT==0 && Classification==9)", + ) + pipeline |= pdal.Filter.grid_decimation_deprecated( + resolution=0.5, + output_dimension=dsm_dimension, + output_type="min", + where="(PT_ON_VEGET==0 && PT_ON_SOL==0 && PT_ON_VIRT==0 && Classification==9)", ) + ################################################################################################################### + # 4 - Gestion des points sol sous la veget,bâtis et ponts pour le MNS + ################################################################################################################### + # On enlève les points sol sous la véget, le bati et les ponts du taguage pour les MNS + # Particularité de reprise des points sol au plus près des bâtiments - # assigne des points sol sélectionnés : les points proches de la végétation, des ponts, de l'eau, 67 + # 4.1 Isolement des points sols sous la véget, le bati et les ponts pipeline = macro.add_radius_assign( pipeline, 1.5, @@ -235,26 +264,27 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref=macro.build_condition("Classification", [4, 5, 6, 17, 67]), condition_out=f"{dsm_dimension}=0", ) - # Test proximité batiment + # 4.2 Particularité de reprise des points sol au plus près des bâtiments pipeline = macro.add_radius_assign( pipeline, 1.25, False, condition_src="Classification==2 && PT_VEG_DSM==0", - condition_ref="Classification==6", + condition_ref=macro.build_condition("Classification", [6, 67]), condition_out="PT_ON_BUILDING=1", ) - # BUFFER INVERSE Se mettre pipeline = macro.add_radius_assign( pipeline, 1, False, condition_src=f"Classification==2 && {dsm_dimension}==0 && PT_ON_BUILDING==1 && {dtm_dimension}==1", - condition_ref="Classification==2 && PT_ON_BUILDING==0 && PT_VEG_DSM==0", + condition_ref=f"Classification==2 && PT_ON_BUILDING==0 && PT_VEG_DSM==0", condition_out=f"{dsm_dimension}=1", ) - # 3 - gestion des ponts - # bouche trou : on filtre les points au milieu du pont en les mettant à PT_ON_BRIDGE=1 + ################################################################################################################### + # 5 - Gestion des classes sous les ponts pour être détaguées pour le MNS dsm_dimension=0 + ################################################################################################################### + pipeline = macro.add_radius_assign( pipeline, 1.5, @@ -263,34 +293,78 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension) condition_ref="Classification==17", condition_out="PT_ON_BRIDGE=1", max2d_above=0, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) - max2d_below=-1, # prendre tous les points qui sont en dessous des points pont (condition_ref) + max2d_below=900, # prendre tous les points qui sont en dessous des points pont (condition_ref) ) pipeline = macro.add_radius_assign( pipeline, 1.25, False, condition_src="PT_ON_BRIDGE==1", - condition_ref="PT_ON_BRIDGE==0 && " - + macro.build_condition("Classification", [2, 3, 4, 5, 6, 9, 67]), + condition_ref="PT_ON_BRIDGE==0 && ( " + + macro.build_condition("Classification", [2, 3, 4, 5, 6, 9, 67]) + + " )", condition_out="PT_ON_BRIDGE=0", max2d_above=0.5, # ne pas prendre les points qui sont au dessus des points pont (condition_ref) max2d_below=0.5, # prendre tous les points qui sont en dessous des points pont (condition_ref) ) pipeline |= pdal.Filter.assign(value=[f"{dsm_dimension}=0 WHERE PT_ON_BRIDGE==1"]) + ################################################################################################################### + # 6 - Ajout des point pour MNT (sol) qui servent au MNS également + ################################################################################################################### - # 4 - point pour DTM servent au DSM également - # HOMOGENEISER L UTILISATION DE PT_VEG_DSM POUR LES POINT SOL SOUS VEGET AVEC PT_ON_VEGET pipeline |= pdal.Filter.assign( value=[ f"{dsm_dimension}=1 WHERE ({dtm_dimension}==1 && PT_VEG_DSM==0 && PT_ON_BRIDGE==0 && PT_ON_BUILDING==0 && PT_ON_VEGET==0)" ] ) - # 5 - Ajout de la classe 66 pts virtuels dans DTM et DSM - pipeline |= pdal.Filter.assign(value=[f"{dtm_dimension}=1 WHERE (Classification==66)"]) + ################################################################################################################### + # 7 - Gestion des pts virtuels 66 pont et eau / MNT et MNS + ################################################################################################################### + # Taguage pour les MNT de tous les points virtuels + # Gestion des pts virtuels qui sont sous la végétation ou autres pour le MNS + # Taguage pour les MNS des points vituels "eau" seulement - # 6 - export du nuage et des DSM - pipeline |= pdal.Writer.las(extra_dims="all", forward="all", filename=output_las) + # 7.1 Taguade pour les MNT des points virtuels ponts et eau + pipeline |= pdal.Filter.assign(value=[f"{dtm_dimension}=1 WHERE Classification==66"]) + # 7.2 gestion des pts 66 "eau" sous le sursol + pipeline = macro.add_radius_assign( + pipeline, + 0.5, + False, + condition_src="Classification==66", + condition_ref=macro.build_condition("Classification", [4, 5, 6, 17, 67]), + condition_out="PT_ON_VEGET=1", + ) + pipeline = macro.add_radius_assign( + pipeline, + 0.5, + False, + condition_src="PT_ON_VEGET==1 && Classification==66", + condition_ref="PT_ON_VEGET==0 && Classification==66", + condition_out="PT_ON_VEGET=0", + ) + # 7.3 Taguage pour les MNS des points virtuels eau seulement + pipeline = macro.add_radius_assign( + pipeline, + 0.5, + False, + condition_src="Classification==66", + condition_ref="Classification==17", + condition_out="PT_ON_BRIDGE=1", + ) + pipeline |= pdal.Filter.assign( + value=[ + f"{dsm_dimension}=1 WHERE (Classification==66 && PT_ON_VEGET==0 && PT_ON_BRIDGE==0)" + ] + ) + + ################################################################################################################## + # 8 - export du nuage et des DSM + ################################################################################################################### + pipeline |= pdal.Writer.las( + extra_dims="all", forward="all", filename=output_las, minor_version="4" + ) return pipeline, temporary_dimensions diff --git a/pdal_ign_macro/version.py b/pdal_ign_macro/version.py index 19f10f4..74bd34b 100755 --- a/pdal_ign_macro/version.py +++ b/pdal_ign_macro/version.py @@ -1,4 +1,4 @@ -__version__ = "0.2.1" +__version__ = "0.3.0" if __name__ == "__main__": diff --git a/test/data/mnx/input/crop_1.laz b/test/data/mnx/input/crop_1.laz new file mode 100644 index 0000000..3478412 Binary files /dev/null and b/test/data/mnx/input/crop_1.laz differ diff --git a/test/data/mnx/input/crop_2.laz b/test/data/mnx/input/crop_2.laz new file mode 100644 index 0000000..c070343 Binary files /dev/null and b/test/data/mnx/input/crop_2.laz differ diff --git a/test/data/mnx/input/crop_3.laz b/test/data/mnx/input/crop_3.laz new file mode 100644 index 0000000..d7348b7 Binary files /dev/null and b/test/data/mnx/input/crop_3.laz differ diff --git a/test/data/mnx/reference/crop_1.laz b/test/data/mnx/reference/crop_1.laz new file mode 100644 index 0000000..38048ad Binary files /dev/null and b/test/data/mnx/reference/crop_1.laz differ diff --git a/test/data/mnx/reference/crop_2.laz b/test/data/mnx/reference/crop_2.laz new file mode 100644 index 0000000..44d5fef Binary files /dev/null and b/test/data/mnx/reference/crop_2.laz differ diff --git a/test/data/mnx/reference/crop_3.laz b/test/data/mnx/reference/crop_3.laz new file mode 100644 index 0000000..fa6a21a Binary files /dev/null and b/test/data/mnx/reference/crop_3.laz differ diff --git a/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py b/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py index 07ff379..947a122 100644 --- a/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py +++ b/test/pdal_ign_macro/test_mark_points_to_use_for_digital_models_with_new_dimension.py @@ -4,6 +4,11 @@ import numpy as np import pdal +import pytest + +from pdaltools.las_remove_dimensions import remove_dimensions_from_las + + from pdal_ign_macro.mark_points_to_use_for_digital_models_with_new_dimension import ( main, mark_points_to_use_for_digital_models_with_new_dimension, @@ -131,3 +136,62 @@ def test_parse_args(): parsed_args_keys = args.__dict__.keys() main_parameters = inspect.signature(main).parameters.keys() assert parsed_args_keys == main_parameters + + +@pytest.mark.parametrize( + "crop", + [ + "crop_1.laz", + "crop_2.laz", + "crop_3.laz", + ], +) +def test_algo_mark_points_for_dm_with_reference(crop): + ini_las = "test/data/mnx/input/" + crop + dsm_dimension = "dsm_marker" + dtm_dimension = "dtm_marker" + with tempfile.NamedTemporaryFile(suffix="_after.las") as las_output: + main( + ini_las, + las_output.name, + dsm_dimension, + dtm_dimension, + "", + "", + keep_temporary_dims=False, + skip_buffer=True, + ) + pipeline = pdal.Pipeline() + pipeline.execute() + + def sort_points(points): + sorted_index = np.lexsort((points["Y"], points["X"], points["Z"], points["GpsTime"])) + np.set_printoptions(precision=30) + sorted_points = points[sorted_index] + return sorted_points + + pipeline_calc = pdal.Pipeline() | pdal.Reader.las(filename=las_output.name) + pipeline_calc.execute() + arr_result = pipeline_calc.arrays[0] + + # read reference : + ref_las = "test/data/mnx/reference/" + crop + + pipeline_ref = pdal.Pipeline() | pdal.Reader.las(filename=ref_las) + pipeline_ref.execute() + arr_reference = pipeline_ref.arrays[0] + + assert len(arr_result) == len(arr_reference) + + arr_result = sort_points(arr_result) + arr_reference = sort_points(arr_reference) + + arr_result_dimensions = list(arr_result.dtype.fields.keys()) + arr_ref_dimensions = list(arr_reference.dtype.fields.keys()) + + assert arr_result_dimensions == arr_ref_dimensions + + for dim in arr_result_dimensions: + diff_mask = np.where(arr_result[dim] == arr_reference[dim], 0, 1) + nb_pts_incorrect = np.count_nonzero(diff_mask) + assert nb_pts_incorrect == 0