Traitement LiDAR Drone avec Python : Du Vol au Jumeau Numérique

from dataclasses import dataclass

@dataclass
class FlightPlan:
    """Parametres de vol pour acquisition LiDAR drone."""
    altitude_m: float = 80.0       # Altitude AGL
    speed_ms: float = 5.0          # Vitesse de vol
    overlap_pct: float = 70.0      # Recouvrement lateral
    scan_rate_hz: int = 240000     # Frequence du LiDAR
    fov_deg: float = 70.0         # Champ de vision
    returns: int = 5              # Nombre de retours

    @property
    def swath_width(self):
        import math
        return 2 * self.altitude_m * math.tan(math.radians(self.fov_deg / 2))

    @property
    def point_density(self):
        """Densite theorique en pts/m2."""
        return self.scan_rate_hz / (self.speed_ms * self.swath_width)

plan = FlightPlan(altitude_m=60, speed_ms=4)
print(f"Fauchee : {plan.swath_width:.1f} m")
print(f"Densite : {plan.point_density:.0f} pts/m2")

import laspy
import numpy as np
from pathlib import Path

def ingest_drone_laz(laz_dir, output_path):
    """Fusionne et valide les dalles LAZ d'un vol drone."""
    laz_files = sorted(Path(laz_dir).glob("*.laz"))
    print(f"Dalles detectees : {len(laz_files)}")

    all_points = []
    total_count = 0

    for laz_file in laz_files:
        las = laspy.read(str(laz_file))

        # Controle qualite : verifier le CRS et les retours
        if las.header.point_count == 0:
            print(f"ATTENTION : {laz_file.name} est vide, ignore.")
            continue

        # Extraire les coordonnees avec precision float64
        coords = np.vstack((las.x, las.y, las.z)).T

        # Filtre de bruit : supprimer les points aberrants (z < 0 ou z > 500)
        valid = (coords[:, 2] > 0) & (coords[:, 2] < 500)
        coords = coords[valid]

        all_points.append(coords)
        total_count += len(coords)
        print(f"  {laz_file.name}: {len(coords):,} points valides")

    merged = np.vstack(all_points)
    print(f"Total fusionne : {total_count:,} points")
    return merged

points = ingest_drone_laz("./vol_2026_02/", "merged.laz")

import pdal
import json

def classify_ground(input_laz, output_laz):
    """Pipeline PDAL : classification sol + vegetation + batiments."""
    pipeline_json = {
        "pipeline": [
            input_laz,
            {
                "type": "filters.assign",
                "assignment": "Classification[:]=0"
            },
            {
                # Filtre morphologique pour extraire le sol
                "type": "filters.smrf",
                "slope": 0.15,
                "window": 18,
                "threshold": 0.5,
                "scalar": 1.2
            },
            {
                # Classification hauteur : vegetation basse / haute
                "type": "filters.hag_nn"
            },
            {
                "type": "filters.assign",
                "assignment": "Classification[0:3]=3",
                "where": "HeightAboveGround > 0.3 AND HeightAboveGround < 2.0"
            },
            {
                "type": "filters.assign",
                "assignment": "Classification[0:5]=5",
                "where": "HeightAboveGround >= 2.0 AND HeightAboveGround < 30.0"
            },
            {
                "type": "writers.las",
                "filename": output_laz,
                "compression": "laszip"
            }
        ]
    }

    pipeline = pdal.Pipeline(json.dumps(pipeline_json))
    count = pipeline.execute()
    print(f"Points traites : {count:,}")
    return pipeline.arrays[0]

classified = classify_ground("merged.laz", "classified.laz")

import numpy as np
import rasterio
from rasterio.transform import from_bounds
from scipy.interpolate import griddata

def generate_raster(points, values, resolution, output_path, crs="EPSG:2154"):
    """Genere un raster GeoTIFF a partir de points 3D."""
    x_min, y_min = points[:, :2].min(axis=0)
    x_max, y_max = points[:, :2].max(axis=0)

    # Grille reguliere
    cols = int((x_max - x_min) / resolution)
    rows = int((y_max - y_min) / resolution)
    xi = np.linspace(x_min, x_max, cols)
    yi = np.linspace(y_min, y_max, rows)
    grid_x, grid_y = np.meshgrid(xi, yi)

    # Interpolation IDW sur la grille
    grid_z = griddata(
        points[:, :2], values,
        (grid_x, grid_y),
        method="linear",
        fill_value=np.nan
    )

    # Ecriture GeoTIFF
    transform = from_bounds(x_min, y_min, x_max, y_max, cols, rows)
    with rasterio.open(
        output_path, "w",
        driver="GTiff",
        height=rows, width=cols,
        count=1, dtype="float32",
        crs=crs, transform=transform,
        compress="lzw"
    ) as dst:
        dst.write(grid_z.astype(np.float32), 1)

    print(f"Raster genere : {output_path} ({cols}x{rows} px)")

# DTM : uniquement les points sol (classe 2)
ground_mask = classified["Classification"] == 2
ground_pts = np.vstack((classified["X"], classified["Y"], classified["Z"])).T
generate_raster(ground_pts[ground_mask], ground_pts[ground_mask, 2], 0.5, "dtm_50cm.tif")

# DSM : tous les points (premier retour)
generate_raster(ground_pts, ground_pts[:, 2], 0.5, "dsm_50cm.tif")

import rasterio
import numpy as np
from scipy.ndimage import generic_filter

def compute_slope_aspect(dtm_path, output_slope, output_aspect):
    """Calcul de la pente et de l'aspect a partir du DTM."""
    with rasterio.open(dtm_path) as src:
        dem = src.read(1)
        res = src.res[0]  # Resolution en metres
        profile = src.profile.copy()

    # Gradient (derivees partielles)
    dy, dx = np.gradient(dem, res)

    # Pente en degres
    slope = np.degrees(np.arctan(np.sqrt(dx**2 + dy**2)))

    # Aspect (orientation) en degres
    aspect = np.degrees(np.arctan2(-dy, dx))
    aspect = (aspect + 360) % 360

    # Ecriture des rasters derivees
    for data, path in [(slope, output_slope), (aspect, output_aspect)]:
        with rasterio.open(path, "w", **profile) as dst:
            dst.write(data.astype(np.float32), 1)

    print(f"Pente max : {slope.max():.1f} degres")
    return slope, aspect

def compute_roughness(dtm_path, output_path, window=5):
    """Rugosite = ecart-type local de l'elevation."""
    with rasterio.open(dtm_path) as src:
        dem = src.read(1)
        profile = src.profile.copy()

    roughness = generic_filter(dem, np.std, size=window)

    with rasterio.open(output_path, "w", **profile) as dst:
        dst.write(roughness.astype(np.float32), 1)
    return roughness

slope, aspect = compute_slope_aspect("dtm_50cm.tif", "slope.tif", "aspect.tif")
roughness = compute_roughness("dtm_50cm.tif", "roughness.tif")

import ifcopenshell
import ifcopenshell.api as api
import numpy as np

def create_terrain_ifc(dtm_points, output_path):
    """Cree un fichier IFC avec le terrain en IfcSite."""
    ifc = ifcopenshell.file(schema="IFC4")

    # Entites de base
    project = api.run("root.create_entity", ifc, ifc_class="IfcProject", name="Drone Survey")
    site = api.run("root.create_entity", ifc, ifc_class="IfcSite", name="Terrain")

    # Sous-echantillonner pour le maillage IFC
    step = max(1, len(dtm_points) // 50000)
    sampled = dtm_points[::step]

    # Creer la representation geometrique du terrain
    # (Triangulation de Delaunay simplifiee)
    from scipy.spatial import Delaunay
    tri = Delaunay(sampled[:, :2])

    print(f"Terrain IFC : {len(sampled)} vertices, {len(tri.simplices)} triangles")

    ifc.write(output_path)
    print(f"Fichier IFC exporte : {output_path}")

create_terrain_ifc(ground_pts[ground_mask], "terrain_drone.ifc")

import numpy as np
import rasterio

def validate_dtm(dtm_path, gcp_file):
    """Compare le DTM avec les points de controle GNSS."""
    gcps = np.loadtxt(gcp_file, delimiter=",", skiprows=1)
    # Colonnes : X, Y, Z_gnss

    with rasterio.open(dtm_path) as src:
        errors = []
        for x, y, z_gnss in gcps:
            row, col = src.index(x, y)
            z_dtm = src.read(1)[row, col]
            if not np.isnan(z_dtm):
                errors.append(z_dtm - z_gnss)

    errors = np.array(errors)
    rmse = np.sqrt(np.mean(errors**2))
    mae = np.mean(np.abs(errors))
    bias = np.mean(errors)

    print(f"GCPs valides    : {len(errors)}")
    print(f"RMSE            : {rmse:.3f} m")
    print(f"MAE             : {mae:.3f} m")
    print(f"Biais           : {bias:.3f} m")
    return {"rmse": rmse, "mae": mae, "bias": bias}

metrics = validate_dtm("dtm_50cm.tif", "gcps_rtk.csv")