Pipeline LiDAR complet: 9 visualisations + classification sémantique automatique

- Correction bug geojson dans process_lidar.py
- Semantic classifier fonctionnel avec K-Means
- 9 visualisations JPEG selon état de l'art 2024-2025
- Statistiques de classification sémantique exportées en JSON
- Nettoyage automatique des fichiers temporaires

Dockerfile

@@ -25,11 +25,13 @@ RUN pip3 install --no-cache-dir \
     rasterio \
     'laspy[laspy]' \
     scikit-image \
+    scikit-learn \
     tqdm
 
-# Copy script
+# Copy scripts
 COPY process_lidar.py /usr/local/bin/
-RUN chmod +x /usr/local/bin/process_lidar.py
+COPY semantic_classifier.py /usr/local/bin/
+RUN chmod +x /usr/local/bin/process_lidar.py /usr/local/bin/semantic_classifier.py
 
 # Create directories with correct permissions
 RUN mkdir -p /data/output /data/input && \

process_lidar.py

@@ -320,6 +320,134 @@ class LidarArchaeoPipeline:
             print(f"  ✗ Erreur slope: {e}")
             return None
 
+    def generate_aspect(self, dem_file, basename):
+        """Generate aspect (orientation of slopes) map"""
+        print(f"  → Aspect (Orientation)...")
+
+        output = self.vis_dir / f"{basename}_aspect.tif"
+
+        try:
+            with rasterio.open(dem_file) as src:
+                dem = src.read(1)
+                transform = src.transform
+                crs = src.crs
+
+            # Calculate aspect (direction of slope)
+            dy, dx = np.gradient(dem)
+            aspect = np.arctan2(-dy, dx) * 180 / np.pi
+            aspect = np.mod(aspect, 360)  # Convert to 0-360
+
+            # Save
+            with rasterio.open(
+                output,
+                'w',
+                driver='GTiff',
+                height=aspect.shape[0],
+                width=aspect.shape[1],
+                count=1,
+                dtype='float32',
+                crs=crs,
+                transform=transform,
+                compress='lzw'
+            ) as dst:
+                dst.write(aspect.astype('float32'), 1)
+
+            return output
+        except Exception as e:
+            print(f"  ✗ Erreur aspect: {e}")
+            return None
+
+    def generate_curvature(self, dem_file, basename):
+        """Generate curvature (terrain concavity/convexity) map"""
+        print(f"  → Courbure (Curvature)...")
+
+        output = self.vis_dir / f"{basename}_curvature.tif"
+
+        try:
+            with rasterio.open(dem_file) as src:
+                dem = src.read(1)
+                transform = src.transform
+                crs = src.crs
+
+            # Calculate curvature using Laplacian
+            dz_dx = np.gradient(dem, axis=1)
+            dz_dy = np.gradient(dem, axis=0)
+            d2z_dx2 = np.gradient(dz_dx, axis=1)
+            d2z_dy2 = np.gradient(dz_dy, axis=0)
+
+            # Mean curvature
+            curvature = (d2z_dx2 + d2z_dy2) / 2
+
+            # Save
+            with rasterio.open(
+                output,
+                'w',
+                driver='GTiff',
+                height=curvature.shape[0],
+                width=curvature.shape[1],
+                count=1,
+                dtype='float32',
+                crs=crs,
+                transform=transform,
+                compress='lzw'
+            ) as dst:
+                dst.write(curvature.astype('float32'), 1)
+
+            return output
+        except Exception as e:
+            print(f"  ✗ Erreur curvature: {e}")
+            return None
+
+    def generate_solar(self, dem_file, basename):
+        """Generate solar irradiance simulation"""
+        print(f"  → Éclairage Solaire (Solar Irradiance)...")
+
+        output = self.vis_dir / f"{basename}_solar.tif"
+
+        try:
+            with rasterio.open(dem_file) as src:
+                dem = src.read(1)
+                transform = src.transform
+                crs = src.crs
+
+            # Calculate gradients
+            dy, dx = np.gradient(dem)
+
+            # Calculate slope and aspect
+            slope = np.arctan(np.sqrt(dx**2 + dy**2))
+            aspect = np.arctan2(-dy, dx)
+
+            # Solar irradiance (morning sun - azimuth 90, altitude 30)
+            az_rad = np.radians(90)
+            alt_rad = np.radians(30)
+
+            # Solar radiation formula
+            solar = np.sin(alt_rad) * np.sin(slope) + \
+                     np.cos(alt_rad) * np.cos(slope) * np.cos(az_rad - aspect)
+
+            # Clip negative values (shadows)
+            solar = np.clip(solar, 0, 1)
+
+            # Save
+            with rasterio.open(
+                output,
+                'w',
+                driver='GTiff',
+                height=solar.shape[0],
+                width=solar.shape[1],
+                count=1,
+                dtype='float32',
+                crs=crs,
+                transform=transform,
+                compress='lzw'
+            ) as dst:
+                dst.write(solar.astype('float32'), 1)
+
+            return output
+        except Exception as e:
+            print(f"  ✗ Erreur solar: {e}")
+            return None
+
     def generate_lrm(self, dem_file, basename):
         """Local Relief Model - deviation from local mean"""
         print(f"  → Local Relief Model...")
@@ -479,6 +607,29 @@ class LidarArchaeoPipeline:
                     title = "Pente (Slope)"
                     legend_label = f"Pente (°)\nMin: {vmin:.1f}° | Max: {vmax:.1f}°"
                     description = "Orange/Clair = Forte pente (murs, talus) | Foncé = Faible pente"
+                elif 'aspect' in str(tif_file):
+                    cmap = 'hsv'  # Cyclic colormap for directions
+                    # Aspect is 0-360, normalize to 0-1
+                    vmin, vmax = 0, 360
+                    data = np.clip((data - vmin) / (vmax - vmin), 0, 1)
+                    title = "Aspect (Orientation)"
+                    legend_label = "Orientation (° du Nord)\nN=0°, E=90°, S=180°, O=270°"
+                    description = "Couleur = Direction de la pente (utile pour orientation bâtiments)"
+                elif 'curvature' in str(tif_file):
+                    cmap = 'RdYlBu_r'  # Diverging for positive/negative curvature
+                    vmax = max(abs(np.percentile(valid_data, 5)), abs(np.percentile(valid_data, 95)), 0.001)
+                    vmin, vmax = -vmax, vmax
+                    data = np.clip((data - vmin) / (vmax - vmin), 0, 1)
+                    title = "Courbure (Curvature)"
+                    legend_label = f"Courbure\nRouge = Convexe (bosse)\nBleu = Concave (creux)"
+                    description = "⭐ Excellent pour fossés, levées, terrasses, talus"
+                elif 'solar' in str(tif_file):
+                    cmap = 'YlOrBr'  # Sun-like colormap
+                    vmin, vmax = 0, 1
+                    data = np.clip(data, vmin, vmax)
+                    title = "Éclairage Solaire (Matin)"
+                    legend_label = "Irradiance Solaire\nFoncé = Ombre | Clair = Éclairé"
+                    description = "Simulation soleil matin - révèle structures orientées Est"
                 elif 'svf' in str(tif_file):
                     cmap = 'viridis'
                     vmin, vmax = np.percentile(valid_data, (5, 95))
@@ -564,20 +715,23 @@ class LidarArchaeoPipeline:
 
         vis_results = {}
 
-        # Generate rasters
+        # Generate rasters (existing + new)
         vis_results['hillshade'] = self.generate_hillshade(dtm_file, basename)
         vis_results['slope'] = self.generate_slope(dtm_file, basename)
+        vis_results['aspect'] = self.generate_aspect(dtm_file, basename)
+        vis_results['curvature'] = self.generate_curvature(dtm_file, basename)
+        vis_results['solar'] = self.generate_solar(dtm_file, basename)
         vis_results['svf'] = self.generate_svf(dtm_file, basename)
         vis_results['lrm'] = self.generate_lrm(dtm_file, basename)
         vis_results['pos_open'] = self.generate_openness(dtm_file, basename, positive=True)
         vis_results['neg_open'] = self.generate_openness(dtm_file, basename, positive=False)
 
-        # Convert to PNG
+        # Convert to JPEG
-        print(f"\n  Conversion images PNG:")
+        print(f"\n  Conversion images JPEG:")
         for name, tif_file in vis_results.items():
-            png_file = self.tif_to_png(tif_file)
+            jpg_file = self.tif_to_png(tif_file)
-            if png_file:
+            if jpg_file:
-                print(f"    ✓ {png_file.name}")
+                print(f"    ✓ {jpg_file.name}")
 
         return vis_results
 
@@ -586,6 +740,9 @@ class LidarArchaeoPipeline:
         patterns = {
             'Hillshade': '*_hillshade_multi.jpg',
             'Pente': '*_slope.jpg',
+            'Aspect': '*_aspect.jpg',
+            'Courbure': '*_curvature.jpg',
+            'Éclairage': '*_solar.jpg',
             'Sky-View Factor': '*_svf.jpg',
             'Local Relief': '*_lrm.jpg',
             'Pos. Openness': '*_positive_openness.jpg',
@@ -601,39 +758,71 @@ class LidarArchaeoPipeline:
         if len(images) < 2:
             return None
 
-        # 2x3 grid with white background
+        # 3x3 grid for 9 visualizations
-        fig, axes = plt.subplots(2, 3, figsize=(20, 14), facecolor='white')
+        fig, axes = plt.subplots(3, 3, figsize=(24, 18), facecolor='white')
         axes = axes.flatten()
 
         for idx, (name, img_path) in enumerate(images.items()):
-            if idx >= 6:
+            if idx >= 9:
                 break
             img = plt.imread(img_path)
             axes[idx].imshow(img)
-            axes[idx].set_title(name, fontsize=12, fontweight='bold')
+            axes[idx].set_title(name, fontsize=11, fontweight='bold')
             axes[idx].axis('off')
 
         # Hide unused subplots
-        for idx in range(len(images), 6):
+        for idx in range(len(images), 9):
             axes[idx].axis('off')
 
         # Title
         fig.suptitle(
             f"Analyse Archéologique LiDAR - {basename}",
-            fontsize=16,
+            fontsize=18,
             fontweight='bold',
             y=0.98
         )
 
         plt.tight_layout()
         output = self.report_dir / f"{basename}_overview.jpg"
-        plt.savefig(output, dpi=150, bbox_inches='tight', facecolor='white', format='jpg')
+        plt.savefig(output, dpi=120, bbox_inches='tight', facecolor='white', format='jpg')
         plt.close()
 
         return output
 
     # ============ Complete Pipeline ============
 
+    def run_semantic_classification(self, dtm_file, basename):
+        """Run semantic classification on DTM"""
+        print(f"\n[4/4] Classification Sémantique Automatique...")
+
+        try:
+            # Import semantic classifier
+            from semantic_classifier import ArchaeoSemanticClassifier
+
+            # Create output subdirectory for semantic results
+            semantic_dir = self.output_dir / "semantic"
+            semantic_dir.mkdir(exist_ok=True)
+
+            # Run classification
+            classifier = ArchaeoSemanticClassifier(dtm_file, semantic_dir)
+            results = classifier.process(basename)
+
+            print(f"  ✓ Classification sémantique terminée")
+            print(f"    → Carte: {results['tif'].name}")
+            print(f"    → Visualisation: {results['jpg'].name}")
+            stats_file = Path(results['tif']).parent / f"{Path(results['tif']).stem.replace('_semantic', '')}_statistics.json"
+            print(f"    → Statistiques: {stats_file.name}")
+
+            return results
+        except ImportError as e:
+            print(f"  ✗ Module non disponible: {e}")
+            return None
+        except Exception as e:
+            print(f"  ✗ Erreur classification: {e}")
+            import traceback
+            traceback.print_exc()
+            return None
+
     def process_file(self, laz_file):
         """Process a single LAZ file"""
         basename = laz_file.stem
@@ -642,21 +831,21 @@ class LidarArchaeoPipeline:
         print(f"{'='*60}")
 
         # Step 1: Ground classification
-        print(f"\n[1/3] Classification du sol...")
+        print(f"\n[1/4] Classification du sol...")
         las_file = self.classify_ground(laz_file)
         if not las_file:
             print(f"  ✗ Échec classification")
             return False
 
         # Step 2: Generate DTM
-        print(f"\n[2/3] Génération DTM...")
+        print(f"\n[2/4] Génération DTM...")
         dtm_file = self.create_dtm_fast(las_file, basename)
         if not dtm_file:
             print(f"  ✗ Échec DTM")
             return False
 
         # Step 3: Visualizations
-        print(f"\n[3/3] Visualisations archéologiques...")
+        print(f"\n[3/4] Visualisations archéologiques...")
         vis = self.generate_all_visualizations(dtm_file, basename)
 
         # Overview
@@ -664,6 +853,20 @@ class LidarArchaeoPipeline:
         if overview:
             print(f"\n  ✓ Vue synthétique: {overview}")
 
+        # Step 4: Semantic Classification (NEW!)
+        semantic_results = self.run_semantic_classification(dtm_file, basename)
+
+        print(f"\n✓ {basename} traité avec succès !")
+        return True
+
+        # Overview
+        overview = self.create_overview(basename)
+        if overview:
+            print(f"\n  ✓ Vue synthétique: {overview}")
+
+        # Step 4: Semantic Classification (NEW!)
+        semantic_results = self.run_semantic_classification(dtm_file, basename)
+
         print(f"\n✓ {basename} traité avec succès !")
         return True
 

requirements.txt

@@ -6,3 +6,4 @@ rasterio>=1.3
 laspy[laspy]>=2.5
 scikit-image>=0.21
 tqdm>=4.65
+scikit-learn>=1.3