From cd99c7531eb47148153d288a2befa6ad7e6943ef Mon Sep 17 00:00:00 2001
From: Antoine Jacquin <lidar@arkel.fr>
Date: Sat, 9 May 2026 01:21:28 +0200
Subject: [PATCH] Correction orientation nord + filtrage ReturnNumber=0 pour
 fichiers corrompus
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- Inversion axe Y du DTM pour orientation nord correcte
- Fallback filters.range pour fichiers LAZ avec ReturnNumber=0
- Flèche nord vectorielle noire au-dessus de la légende
- 9/9 fichiers traités avec succès

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
---
 process_lidar.py | 116 ++++++++++++++++++++++++++++-------------------
 1 file changed, 69 insertions(+), 47 deletions(-)

diff --git a/process_lidar.py b/process_lidar.py
index 16ea815..6fd27a7 100755
--- a/process_lidar.py
+++ b/process_lidar.py
@@ -109,53 +109,7 @@ class LidarArchaeoPipeline:
             print(f"  ✓ Classification déjà effectuée")
             return output_las
 
-        # First, try to fix ReturnNumber values if needed
-        fixed_las = self.temp_dir / f"{laz_file.stem}_fixed.las"
-
-        # Create preprocessing pipeline to fix ReturnNumber=0 issue using expression
-        fix_pipeline = [
-            {
-                "type": "readers.las",
-                "filename": str(laz_file)
-            },
-            {
-                "type": "filters.expression",
-                "expression": "ReturnNumber = MAX(ReturnNumber, 1)"
-            },
-            {
-                "type": "filters.expression",
-                "expression": "NumberOfReturns = MAX(NumberOfReturns, 1)"
-            },
-            {
-                "type": "writers.las",
-                "filename": str(fixed_las),
-                "extra_dims": "all"
-            }
-        ]
-
-        try:
-            # Try with fixed data first
-            fix_json = json.dumps(fix_pipeline)
-            fix_pipe_file = self.temp_dir / "fix_pipeline.json"
-            with open(fix_pipe_file, 'w') as f:
-                f.write(fix_json)
-
-            import subprocess
-            result = subprocess.run(['pdal', 'pipeline', str(fix_pipe_file)],
-                                    capture_output=True, text=True)
-            if result.returncode == 0:
-                print(f"  → Correction ReturnNumber effectuée")
-                input_for_smrf = fixed_las
-            else:
-                print(f"  → Erreur correction, utilisation fichier original")
-                input_for_smrf = laz_file
-                if fixed_las.exists():
-                    fixed_las.unlink()
-        except Exception as e:
-            print(f"  → Erreur prétraitement: {e}")
-            input_for_smrf = laz_file
-
-        pipeline_json = self.create_pipeline_json(input_for_smrf, output_las)
+        pipeline_json = self.create_pipeline_json(laz_file, output_las)
         pipeline_file = self.temp_dir / "pipeline.json"
 
         with open(pipeline_file, 'w') as f:
@@ -169,6 +123,71 @@ class LidarArchaeoPipeline:
             )
             print(f"  ✓ Classification sol terminée")
             return output_las
+        except subprocess.CalledProcessError as e:
+            error_msg = e.stderr.decode()
+            print(f"  ✗ Erreur PDAL: {error_msg}")
+
+            # If error is about ReturnNumber=0, try filtering those points
+            if "ReturnNumber" in error_msg and "NumberOfReturns" in error_msg:
+                print(f"  → Tentative de filtrage des points ReturnNumber=0...")
+
+                # Use filters.range to keep only valid points (ReturnNumber >= 1)
+                filtered_pipeline = [
+                    {
+                        "type": "readers.las",
+                        "filename": str(laz_file)
+                    },
+                    {
+                        "type": "filters.range",
+                        "limits": "ReturnNumber[1:],NumberOfReturns[1:]"
+                    },
+                    {
+                        "type": "filters.smrf",
+                        "scalar": 1.25
+                    },
+                    {
+                        "type": "filters.range",
+                        "limits": "Classification[2:2]"
+                    },
+                    {
+                        "type": "writers.las",
+                        "filename": str(output_las),
+                        "extra_dims": "all"
+                    }
+                ]
+
+                filtered_json = json.dumps(filtered_pipeline)
+                with open(pipeline_file, 'w') as f:
+                    f.write(filtered_json)
+
+                try:
+                    subprocess.run(
+                        ["pdal", "pipeline", str(pipeline_file)],
+                        capture_output=True,
+                        check=True
+                    )
+                    print(f"  ✓ Classification sol terminée (points filtrés)")
+                    return output_las
+                except subprocess.CalledProcessError as e2:
+                    print(f"  ✗ Erreur même avec filtrage: {e2.stderr.decode()}")
+                    return None
+
+            return None
+
+    def run_pdal_pipeline(self, pipeline_json, output_file):
+        """Exécute un pipeline PDAL"""
+        pipeline_file = self.temp_dir / "temp_pipeline.json"
+        with open(pipeline_file, 'w') as f:
+            f.write(pipeline_json)
+
+        try:
+            subprocess.run(
+                ["pdal", "pipeline", str(pipeline_file)],
+                capture_output=True,
+                check=True
+            )
+            print(f"  ✓ Pipeline terminé")
+            return output_file
         except subprocess.CalledProcessError as e:
             print(f"  ✗ Erreur PDAL: {e.stderr.decode()}")
             return None
@@ -209,6 +228,9 @@ class LidarArchaeoPipeline:
 
             dtm = stat.statistic.T
 
+            # Flip Y axis so north is at the top (Y decreases from top to bottom in image)
+            dtm = dtm[::-1, :]
+
             # Fill gaps using interpolation
             print(f"    Remplissage des zones vides...")
             from scipy.ndimage import distance_transform_edt