Update script mnx (#19)

* new version algo * un coup de black * update code * fix mnx script --------- Co-authored-by: Lea Vauchier <[email protected]>
IGNF · Jul 18, 2024 · 08747ad · 08747ad
1 parent 0dc6d3a
commit 08747ad
Showing 1 changed file with 62 additions and 26 deletions.
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
@@ -89,14 +89,20 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
     pipeline = pdal.Pipeline() | pdal.Reader.las(input_las)
 
     # 0 - ajout de dimensions temporaires et de sortie
-    temporary_dimensions = ["PT_VEG_DSM", "PT_ON_BRIDGE", "PT_ON_BUILDING", "PT_ON_VEGET"]
+    temporary_dimensions = [
+        "PT_VEG_DSM",
+        "PT_ON_BRIDGE",
+        "PT_ON_BUILDING",
+        "PT_ON_VEGET",
+        "PT_ON_SOL",
+    ]
     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) proche de la végétation
+    #     vegetation (3) proches de la végétation
     #     pour le calcul du DSM
 
     pipeline |= pdal.Filter.assign(
@@ -116,11 +122,11 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         pipeline,
         1,
         False,
-        condition_src=macro.build_condition("Classification", [6, 17]),
+        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)
+        max2d_below=-1,  # prendre tous les points qui sont en dessous des points pont (condition_ref)
     )
     pipeline = macro.add_radius_assign(
         pipeline,
@@ -134,8 +140,8 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         pipeline,
         1,
         False,
-        condition_src="PT_ON_VEGET==1 && Classification==6",
-        condition_ref="Classification==6 && PT_ON_VEGET==0",
+        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)
@@ -150,6 +156,16 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         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==9",
+        condition_ref="Classification==9 && 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)
+    )
 
     # selection des points de veget basse proche de la veget haute
     pipeline = macro.add_radius_assign(
@@ -169,13 +185,37 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
 
     # 2 - sélection des points pour DTM et DSM
 
+    ###########################################################################################################################################
+    # Voir pour l'eau sous le sol dans les dévers corse
+    pipeline = macro.add_radius_assign(
+        pipeline,
+        1.25,
+        False,
+        condition_src="Classification==9",
+        condition_ref="Classification==2",
+        condition_out="PT_ON_SOL=1",
+        max2d_above=0,
+        max2d_below=-1,
+    )
+    pipeline = macro.add_radius_assign(
+        pipeline,
+        1,
+        False,
+        condition_src="PT_ON_SOL==1",
+        condition_ref="PT_ON_SOL==0 && Classification==9",
+        condition_out="PT_ON_SOL=0",
+        max2d_above=0.5,
+        max2d_below=0.5,
+    )
+    ###########################################################################################################################################
+
     # selection de points DTM (max) 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.5,
         output_dimension=dtm_dimension,
         output_type="max",
-        where="Classification==2",
+        where="(Classification==2 || PT_ON_SOL==0 && Classification==9)",
     )
 
     # selection de points DSM (max) sur une grille régulière
@@ -184,18 +224,18 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         resolution=0.5,
         output_dimension=dsm_dimension,
         output_type="max",
-        where="(PT_ON_VEGET==0 && ("
+        where="(PT_ON_VEGET==0 && "
         + macro.build_condition("Classification", [6, 9, 17, 64])
-        + f") || {dsm_dimension}==1)",
+        + f" || {dsm_dimension}==1)",
     )
 
-    # assigne des points sol sélectionnés : les points proches de la végétation, des ponts, de l'eau, 64
+    # assigne des points sol sélectionnés : les points proches de la végétation, des ponts, de l'eau, 67
     pipeline = macro.add_radius_assign(
         pipeline,
         1.5,
         False,
         condition_src=f"{dtm_dimension}==1",
-        condition_ref=macro.build_condition("Classification", [4, 5, 6, 9, 17, 64]),
+        condition_ref=macro.build_condition("Classification", [4, 5, 6, 17, 67]),
         condition_out=f"{dsm_dimension}=0",
     )
     # Test proximité batiment
@@ -217,46 +257,42 @@ def define_marking_pipeline(input_las, output_las, dsm_dimension, dtm_dimension)
         condition_out=f"{dsm_dimension}=1",
     )
     # 3 - gestion des ponts
-    #     bouche trou : on filtre les points (2,3,4,5,9) au milieu du pont en les mettant à PT_ON_BRIDGE=1
+    #     bouche trou : on filtre les points au milieu du pont en les mettant à PT_ON_BRIDGE=1
     pipeline = macro.add_radius_assign(
         pipeline,
         1.5,
         False,
-        condition_src=macro.build_condition("Classification", [2, 3, 4, 5, 6, 9]),
+        condition_src=macro.build_condition("Classification", [2, 3, 4, 5, 6, 9, 67]),
         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=900,  # prendre tous les points qui sont en dessous 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.25,
         False,
-        # condition_ref=macro.build_condition("Classification", [2, 3, 4, 5]),
         condition_src="PT_ON_BRIDGE==1",
-        condition_ref="PT_ON_BRIDGE==0 && ( "
-        + macro.build_condition("Classification", [2, 3, 4, 5, 6, 9])
-        + " )",
+        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 && NOT(Classification==17))"])
-    pipeline |= pdal.Filter.assign(
-        value=[f"{dsm_dimension}=0 WHERE PT_ON_BRIDGE==1"]
-        # value=["dsm_marker=0 WHERE (PT_ON_BRIDGE==1 AND ( " + macro.build_condition("Classification", [2,3,4,5,6,9]) + " ))"]
-    )
+    pipeline |= pdal.Filter.assign(value=[f"{dsm_dimension}=0 WHERE PT_ON_BRIDGE==1"])
 
     # 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 )"
+            f"{dsm_dimension}=1 WHERE ({dtm_dimension}==1 && PT_VEG_DSM==0 && PT_ON_BRIDGE==0 && PT_ON_BUILDING==0 && PT_ON_VEGET==0)"
         ]
     )
-    # ERREUR EN 4!###############################################################################################!
-    # 5 - export du nuage et des DSM
 
+    # 5 - Ajout de la classe 66 pts virtuels dans DTM et DSM
+    pipeline |= pdal.Filter.assign(value=[f"{dtm_dimension}=1 WHERE (Classification==66)"])
+
+    # 6 - export du nuage et des DSM
     pipeline |= pdal.Writer.las(extra_dims="all", forward="all", filename=output_las)
 
     return pipeline, temporary_dimensions