Propagate metadata through pdal and laspy pipelines

.gitignore

@@ -7,6 +7,9 @@ outputs/
 tests/files/large/V0.5_792000_6272000.las
 tests/files/large/842_6521_invalid_band.las
 
+# test output
+tmp/
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

.pre-commit-config.yaml

@@ -7,6 +7,7 @@ repos:
       - id: end-of-file-fixer
       - id: check-yaml
       - id: check-added-large-files
+        args: ["--maxkb=1000"]
       - id: debug-statements
       - id: detect-private-key
 

CHANGELOG.md

@@ -1,5 +1,9 @@
 # main
 
+### 1.10.2
+- Add support for metadata propagation through compound pdal pipelines:
+  - fix epsg propagation
+
 ### 1.10.1
 - Fix edge case when BD uni does not have data for the requested bbox
 

environment.yml

@@ -43,3 +43,4 @@ dependencies:
       - sphinxnotes-mock==1.0.*
       - sphinx-argparse==0.4.*
       - sphinxcontrib-mermaid==0.9.*
+      - ign-pdal-tools>=1.5.2

lidar_prod/application.py

@@ -54,7 +54,7 @@ def get_list_las_path_from_src(src_path: str):
 def identify_vegetation_unclassified(config, src_las_path: str, dest_las_path: str):
     log.info(f"Identifying on {src_las_path}")
     data_format = config["data_format"]
-    las_data = get_las_data_from_las(src_las_path)
+    las_data = get_las_data_from_las(src_las_path, config.data_format.epsg)
 
     # add the necessary dimension to store the results
     cleaner: Cleaner = hydra.utils.instantiate(data_format.cleaning.input_vegetation_unclassified)
@@ -91,7 +91,7 @@ def just_clean(config, src_las_path: str, dest_las_path: str):
     avoid delays when doing the same operations over and over again )"""
     log.info(f"Cleaning {src_las_path}")
     data_format = config["data_format"]
-    las_data = get_las_data_from_las(src_las_path)
+    las_data = get_las_data_from_las(src_las_path, 2154)
 
     # remove unwanted dimensions
     cleaner = hydra.utils.instantiate(data_format.cleaning.input)
@@ -131,15 +131,15 @@ def apply_building_module(config: DictConfig, src_las_path: str, dest_las_path:
             thresholds=bv_cfg.thresholds,
             use_final_classification_codes=bv_cfg.use_final_classification_codes,
         )
-        bv.run(tmp_las_path)
+        las_metadata = bv.run(tmp_las_path)
 
         # Complete buildings with non-candidates that were nevertheless confirmed
         bc: BuildingCompletor = hydra.utils.instantiate(config.building_completion)
-        bc.run(bv.pipeline)
+        las_metadata = bc.run(bv.pipeline, las_metadata)
 
         # Define groups of confirmed building points among non-candidates
         bi: BuildingIdentifier = hydra.utils.instantiate(config.building_identification)
-        bi.run(bc.pipeline, tmp_las_path)
+        bi.run(bc.pipeline, tmp_las_path, las_metadata=las_metadata)
 
         # Remove unnecessary intermediary dimensions
         cl: Cleaner = hydra.utils.instantiate(config.data_format.cleaning.output_building)
@@ -166,7 +166,7 @@ def get_shapefile(config: DictConfig, src_las_path: str, dest_las_path: str):
             get_pipeline(
                 src_las_path,
                 config.data_format.epsg,
-            ),
+            )[0],
             buffer=config.building_validation.application.bd_uni_request.buffer,
         ),  # bbox
         config.data_format.epsg,

lidar_prod/tasks/building_completion.py

@@ -39,27 +39,34 @@ def __init__(
         self.data_format = data_format
         self.pipeline: pdal.pipeline.Pipeline = None
 
-    def run(self, input_values: Union[str, pdal.pipeline.Pipeline]):
+    def run(
+        self, input_values: Union[str, pdal.pipeline.Pipeline], las_metadata: dict = None
+    ) -> dict:
         """Application.
 
         Transform cloud at `src_las_path` following building completion logic
 
         Args:
             input_values (str|pdal.pipeline.Pipeline): path to either input LAS file or a pipeline
             target_las_path (str): path for saving updated LAS file.
+            las_metadata (dict): current pipeline metadata, used to propagate input metadata to the
+        application output las (epsg, las version, etc)
 
         Returns:
-            str: returns `target_las_path` for potential terminal piping.
-
+            str: returns `las_metadata`: metadata of the initial las, which contain
+            information to pass to the writer in order for the application to have an output
+            with the same header (las version, srs, ...) as the input
         """
         log.info(
             "Completion of building with relatively distant points that have high enough "
             + "probability"
         )
-        pipeline = get_pipeline(input_values, self.data_format.epsg)
+        pipeline, las_metadata = get_pipeline(input_values, self.data_format.epsg, las_metadata)
         self.prepare_for_building_completion(pipeline)
         self.update_classification()
 
+        return las_metadata
+
     def prepare_for_building_completion(self, pipeline: pdal.pipeline.Pipeline) -> None:
         """Prepare for building completion.
 
@@ -70,9 +77,7 @@ def prepare_for_building_completion(self, pipeline: pdal.pipeline.Pipeline) -> N
         the same building and they will be confirmed as well.
 
         Args:
-            src_las_path (pdal.pipeline.Pipeline): input LAS pipeline
-            target_las_path (str): output, prepared LAS.
-
+            pipeline (pdal.pipeline.Pipeline): input LAS pipeline
         """
 
         # Reset Cluster dim out of safety

lidar_prod/tasks/building_identification.py

@@ -34,22 +34,27 @@ def run(
         self,
         input_values: Union[str, pdal.pipeline.Pipeline],
         target_las_path: str = None,
-    ) -> str:
+        las_metadata: dict = None,
+    ) -> dict:
         """Identify potential buildings in a new channel, excluding former candidates as well as
         already confirmed building (confirmed by either Validation or Completion).
 
         Args:
             input_values (str | pdal.pipeline.Pipeline): path or pipeline to input LAS file with
             a building probability channel
             target_las_path (str): output LAS
+            las_metadata (dict): current pipeline metadata, used to propagate input metadata to the
+        application output las (epsg, las version, etc)
+
+        Returns: updated las_metadata
 
         """
         # aliases
         _cid = self.data_format.las_dimensions.cluster_id
         _completion_flag = self.data_format.las_dimensions.completion_non_candidate_flag
 
         log.info("Clustering of points with high building proba.")
-        pipeline = get_pipeline(input_values, self.data_format.epsg)
+        pipeline, las_metadata = get_pipeline(input_values, self.data_format.epsg, las_metadata)
 
         # Considered for identification:
         non_candidates = f"({self.data_format.las_dimensions.candidate_buildings_flag} == 0)"
@@ -78,10 +83,10 @@ def run(
             dimensions=f"{_cid}=>{self.data_format.las_dimensions.ai_building_identified}"
         )
         if target_las_path:
-            pipeline |= get_pdal_writer(target_las_path)
+            pipeline |= get_pdal_writer(target_las_path, las_metadata)
             os.makedirs(osp.dirname(target_las_path), exist_ok=True)
         pipeline.execute()
 
         self.pipeline = pipeline
 
-        return target_las_path
+        return las_metadata

lidar_prod/tasks/building_validation.py

@@ -81,7 +81,8 @@ def run(
         self,
         input_values: Union[str, pdal.pipeline.Pipeline],
         target_las_path: str = None,
-    ) -> str:
+        las_metadata: dict = None,
+    ) -> dict:
         """Runs application.
 
         Transforms cloud at `input_values` following building validation logic,
@@ -91,30 +92,35 @@ def run(
             input_values (str| pdal.pipeline.Pipeline): path or pipeline to input LAS file with
             a building probability channel
             target_las_path (str): path for saving updated LAS file.
+            las_metadata (dict): current pipeline metadata, used to propagate input metadata to the
+        application output las (epsg, las version, etc)
 
         Returns:
-            str: returns `target_las_path`
+            str: returns `las_metadata`: metadata of the input las, which contain
+            information to pass to the writer in order for the application to have an output
+            with the same header (las version, srs, ...) as the input
 
         """
-        self.pipeline = get_pipeline(input_values, self.data_format.epsg)
+        self.pipeline, las_metadata = get_pipeline(input_values, self.data_format.epsg)
         with TemporaryDirectory() as td:
             log.info("Preparation : Clustering of candidates buildings & Import vectors")
             if isinstance(input_values, str):
                 log.info(f"Applying Building Validation to file \n{input_values}")
                 temp_f = osp.join(td, osp.basename(input_values))
             else:
                 temp_f = ""
-            self.prepare(input_values, temp_f)
+            las_metadata = self.prepare(input_values, temp_f, las_metadata)
             log.info("Using AI and Databases to update cloud Classification")
-            self.update()
-        return target_las_path
+            las_metadata = self.update(target_las_path=target_las_path, las_metadata=las_metadata)
+        return las_metadata
 
     def prepare(
         self,
         input_values: Union[str, pdal.pipeline.Pipeline],
         prepared_las_path: str,
         save_result: bool = False,
-    ) -> None:
+        las_metadata: dict = None,
+    ) -> dict:
         f"""
         Prepare las for later decision process. .
         1. Cluster candidates points, in a new
@@ -135,6 +141,11 @@ def prepare(
             a building probability channel
             target_las_path (str): path for saving prepared LAS file.
             save_result (bool): True to save a las instead of propagating a pipeline
+            las_metadata (dict): current pipeline metadata, used to propagate input metadata to the
+        application output las (epsg, las version, etc)
+
+        Returns:
+            updated las metadata
 
         """
 
@@ -143,7 +154,7 @@ def prepare(
         dim_cluster_id_candidates = self.data_format.las_dimensions.ClusterID_candidate_building
         dim_overlay = self.data_format.las_dimensions.uni_db_overlay
 
-        self.pipeline = get_pipeline(input_values, self.data_format.epsg)
+        self.pipeline, las_metadata = get_pipeline(input_values, self.data_format.epsg)
         # Identify candidates buildings points with a boolean flag
         self.pipeline |= pdal.Filter.ferry(dimensions=f"=>{dim_candidate_flag}")
         _is_candidate_building = (
@@ -204,17 +215,21 @@ def prepare(
             )
 
         if save_result:
-            self.pipeline |= get_pdal_writer(prepared_las_path)
+            self.pipeline |= get_pdal_writer(prepared_las_path, las_metadata)
             os.makedirs(osp.dirname(prepared_las_path), exist_ok=True)
         self.pipeline.execute()
 
         if temp_dirpath:
             shutil.rmtree(temp_dirpath)
 
-    def update(self, src_las_path: str = None, target_las_path: str = None) -> None:
+        return las_metadata
+
+    def update(
+        self, src_las_path: str = None, target_las_path: str = None, las_metadata: dict = None
+    ) -> None:
         """Updates point cloud classification channel."""
         if src_las_path:
-            self.pipeline = get_pipeline(src_las_path, self.data_format.epsg)
+            self.pipeline, las_metadata = get_pipeline(src_las_path, self.data_format.epsg)
 
         points = self.pipeline.arrays[0]
 
@@ -250,10 +265,12 @@ def update(self, src_las_path: str = None, target_las_path: str = None) -> None:
         self.pipeline = pdal.Pipeline(arrays=[points])
 
         if target_las_path:
-            self.pipeline = get_pdal_writer(target_las_path).pipeline(points)
+            self.pipeline = get_pdal_writer(target_las_path, las_metadata).pipeline(points)
             os.makedirs(osp.dirname(target_las_path), exist_ok=True)
             self.pipeline.execute()
 
+        return las_metadata
+
     def _extract_cluster_info_by_idx(
         self, las: np.ndarray, pts_idx: np.ndarray
     ) -> BuildingValidationClusterInfo:

lidar_prod/tasks/building_validation_optimization.py

@@ -236,7 +236,7 @@ def _extract_clusters_from_las(
             candidate buildings
 
         """
-        las = pdal_read_las_array(prepared_las_path, self.bv.data_format.epsg)
+        las, _ = pdal_read_las_array(prepared_las_path, self.bv.data_format.epsg)
         # las: laspy.LasData = laspy.read(prepared_las_path)
         dim_cluster_id = las[self.bv.data_format.las_dimensions.ClusterID_candidate_building]
         dim_classification = las[self.bv.data_format.las_dimensions.classification]

lidar_prod/tasks/cleaning.py

@@ -60,15 +60,16 @@ def run(self, src_las_path: str, target_las_path: str, epsg: int | str):
         it from the las metadata)
 
         """
-        points = pdal_read_las_array(src_las_path, epsg)
+        points, metadata = pdal_read_las_array(src_las_path, epsg)
+
         # Check input dims to see what we can keep.
         input_dims = points.dtype.fields.keys()
         self.extra_dims_as_dict = {
             k: v for k, v in self.extra_dims_as_dict.items() if k in input_dims
         }
 
         pipeline = pdal.Pipeline(arrays=[points]) | get_pdal_writer(
-            target_las_path, extra_dims=self.get_extra_dims_as_str()
+            target_las_path, reader_metadata=metadata, extra_dims=self.get_extra_dims_as_str()
         )
         os.makedirs(osp.dirname(target_las_path), exist_ok=True)
         pipeline.execute()