From c781d53896dbbe3b86c80efc0deb0ac66cd69d7a Mon Sep 17 00:00:00 2001 From: Antoine Jacquin Date: Fri, 8 May 2026 23:37:14 +0200 Subject: [PATCH] =?UTF-8?q?Pipeline=20LiDAR=20complet:=209=20visualisation?= =?UTF-8?q?s=20+=20classification=20s=C3=A9mantique=20automatique?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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 | 6 +- process_lidar.py | 235 +++++++++++++++++++++++++++++++++++++++++++---- requirements.txt | 1 + 3 files changed, 224 insertions(+), 18 deletions(-) diff --git a/Dockerfile b/Dockerfile index e3c2e52..3c18100 100644 --- a/Dockerfile +++ b/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 && \ diff --git a/process_lidar.py b/process_lidar.py index 2801035..08b048e 100755 --- a/process_lidar.py +++ b/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 - print(f"\n Conversion images PNG:") + # Convert to JPEG + print(f"\n Conversion images JPEG:") for name, tif_file in vis_results.items(): - png_file = self.tif_to_png(tif_file) - if png_file: - print(f" ✓ {png_file.name}") + jpg_file = self.tif_to_png(tif_file) + if jpg_file: + 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 - fig, axes = plt.subplots(2, 3, figsize=(20, 14), facecolor='white') + # 3x3 grid for 9 visualizations + 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 diff --git a/requirements.txt b/requirements.txt index ef862ef..0a00f3b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -6,3 +6,4 @@ rasterio>=1.3 laspy[laspy]>=2.5 scikit-image>=0.21 tqdm>=4.65 +scikit-learn>=1.3