add pdal pipeline examples

macro/ex_filtering_points.py

@@ -0,0 +1,55 @@
+import argparse
+import pdal
+from macro import *
+
+"""
+This tool shows how to use function of macro in a pdal pipeline
+"""
+
+def parse_args():
+    parser = argparse.ArgumentParser("Tools for apply pdal pipelines")
+    parser.add_argument("--input", "-i", type=str, required=True, help="Input las file")
+    parser.add_argument("--output", "-o", type=str, required=True, help="Output las file")
+    return parser.parse_args()
+
+
+if __name__ == "__main__":
+    args = parse_args()
+
+    pipeline = pdal.Reader.las(args.input)
+
+    # step 1 a 9
+    pipeline = add_radius_search(pipeline, 1, False, "Classification==2", build_condition("Classification", [4,5]), "Classification=102")
+    pipeline = add_radius_search(pipeline, 1, False, "Classification==102", "Classification==2", "Classification=2")
+    pipeline = add_radius_search(pipeline, 1, False, "Classification==3", "Classification==5", "Classification=103")
+    pipeline = add_grid_decimation(pipeline, 0.75, "max", build_condition("Classification", [4,5,102,103]), "Classification=100")
+    pipeline |= pdal.Filter.assign(value="Classification=2", where="Classification==102")
+    pipeline |= pdal.Filter.assign(value="Classification=3", where="Classification==103")
+    pipeline = add_grid_decimation(pipeline, 0.5, "max", "Classification==2", "Classification=102")
+    pipeline = add_grid_decimation(pipeline, 0.5, "max", build_condition("Classification", [2,3,4,5,6,9,17,64,100]), "Classification=200")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==102", build_condition("Classification", [4,5,6,9,17,64,100]), "Classification=100")
+    pipeline |= pdal.Filter.assign(value="Classification=2", where="Classification==102")
+
+    # step 10
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==2", "Classification==17", "Classification=102")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==102", build_condition("Classification", [2,3,4,5]), "Classification=2")
+
+    # step 11
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==3", "Classification==17", "Classification=103")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==103", build_condition("Classification", [2,3,4,5]), "Classification=3")
+
+    # step 12
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==4", "Classification==17", "Classification=104")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==104", build_condition("Classification", [2,3,4,5]), "Classification=4")
+
+    # step 13
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==5", "Classification==17", "Classification=105")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==105", build_condition("Classification", [2,3,4,5]), "Classification=5")
+
+    # step 14
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==9", "Classification==17", "Classification=109")
+    pipeline = add_radius_search(pipeline, 1.5, False, "Classification==109", "Classification==9", "Classification=9")
+
+    pipeline |= pdal.Writer.las(extra_dims="all",minor_version=4,dataformat_id=6,filename=args.output)
+    pipeline.execute()
+

macro/macro.py

@@ -0,0 +1,81 @@
+import argparse
+import pdal
+
+"""
+Some useful filters combinations for complete pdal pipeline
+"""
+
+
+def add_radius_search(pipeline, radius, search_3d, condition_src, condition_ref, condition_out ):
+    """
+    search points from "condition_src" closed from "condition_ref", and reassign them to "condition_out"
+    This combination is equivalent to the CloseBy macro of TerraScan
+    radius : the search distance
+    search_3d : the distance reseach is in 3d if True
+    condition_src, condition_ref, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.ferry(dimensions=f"=>REF_DOMAIN, =>SRS_DOMAIN, =>radius_search")
+    pipeline |= pdal.Filter.assign(value=["SRS_DOMAIN = 0", f"SRS_DOMAIN = 1 WHERE {condition_src}",
+                                          "REF_DOMAIN = 0", f"REF_DOMAIN = 1 WHERE {condition_ref}",
+                                          "radius_search = 0"])
+    pipeline |= pdal.Filter.radius_search(radius=radius, src_domain="SRS_DOMAIN",reference_domain="REF_DOMAIN",
+                                          output_name_attribute="radius_search", search_3d=search_3d)
+    pipeline |= pdal.Filter.assign(value=condition_out,where="radius_search==1")
+    return pipeline
+
+
+
+def add_grid_decimation(pipeline, grid_resolution, output_type, condition, condition_out):
+    """
+    Select a points in a grid  from "condition"; points not selected are reassign to "condition_out"
+    This combination is equivalent to the Thin Points macro of TerraScan
+    grid_resolution : resolution of the grid
+    output_type : "max" or "min" (the highest or lower points of the grid)
+    condition, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.ferry(dimensions=f"=>grid,")
+    pipeline |= pdal.Filter.assign(value="grid = 0")
+    pipeline |= pdal.Filter.grid_decimation(resolution=grid_resolution, output_name_attribut="grid",
+                                            output_type=output_type, where=condition)
+    pipeline |= pdal.Filter.assign(value=condition_out,where=f"grid==0 && ({condition})")
+    return pipeline
+
+
+
+def classify_hgt_ground(pipeline, hmin, hmax, condition, condition_out):
+    """
+       reassign points from "condition" between "hmin" and "hmax" of the ground to "condition_out"
+       This combination is equivalent to the ClassifyHgtGrd macro of TerraScan
+       condition, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.hag_delaunay(allow_extrapolation=True)
+    condition_h = f"HeightAboveGround>{hmin} && HeightAboveGround<={hmax}"
+    condition_h += " && " + condition
+    pipeline |= pdal.Filter.assign(value=condition_out, where=condition_h)
+    return pipeline
+
+
+
+def keep_non_planar_pts(pipeline, condition, condition_out):
+    """
+       reassign points from "condition" who are planar to "condition_out"
+       This combination is equivalent to the ClassifyModelKey macro of TerraScan
+       condition, condition_out : a pdal condition as "Classification==2"
+    """
+    pipeline |= pdal.Filter.approximatecoplanar(knn=8,thresh1=25,thresh2=6,where=condition)
+    pipeline |= pdal.Filter.assign(value=condition_out,where=f"Coplanar==0 && ({condition})")
+    return pipeline
+
+
+
+
+def build_condition(key, values):
+    """
+       build 'key==values[0] || key==values[1] ...'
+    """
+    condition = ""
+    for v in values:
+        condition += key+"=="+str(v)
+        if v!=values[-1]:condition += " || "
+    return condition
+