Préserver les zones sans données: NaN-aware filtering dans les visualisations

- DTM: plus d'interpolation, les zones sans LiDAR restent NaN
- Ajout _fill_nans() et _filter_nanaware(): remplissent les NaN par
  nearest-neighbor avant filtrage, puis restaurent le masque NaN
- Toutes les visualisations avec filtres (LRM, MSLRM, TPI, SAILORE,
  roughness, anomalies, wavelet) utilisent _filter_nanaware pour
  éviter l'érosion des bords de données
- _save_tif() écrit nodata=float('nan') quand le tableau contient des NaN
- Les zones sans données restent vides dans les visualisations
- Les calculs ne sont pas faussés par des valeurs interpolées

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

lidar_pipeline/visualizations.py

@@ -58,6 +58,59 @@ def _read_dem(dem_file):
         return src.read(1), src.transform, src.crs
 
 
+def _fill_nans(arr):
+    """Fill NaN values using nearest-neighbor interpolation.
+
+    Returns (filled_array, nan_mask) so the caller can restore NaN after filtering.
+    """
+    from scipy.interpolate import NearestNDInterpolator
+    nan_mask = np.isnan(arr)
+    if not np.any(nan_mask):
+        return arr, nan_mask
+    valid = ~nan_mask
+    y_coords, x_coords = np.where(valid)
+    if len(y_coords) == 0:
+        return np.zeros_like(arr), nan_mask
+    z_values = arr[valid]
+    interp = NearestNDInterpolator(
+        np.column_stack((y_coords, x_coords)), z_values
+    )
+    y_missing, x_missing = np.where(nan_mask)
+    filled = arr.copy()
+    filled[y_missing, x_missing] = interp(y_missing, x_missing)
+    return filled, nan_mask
+
+
+def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
+    """Apply a filter to an array while preserving NaN zones.
+
+    1. Fill NaN with nearest-neighbor interpolation
+    2. Apply the filter
+    3. Restore original NaN mask on the result
+
+    Args:
+        arr: Input array (numpy or cupy).
+        filter_func: Function that takes (array, *args, **kwargs) and returns filtered array.
+        use_gpu: If True, apply filter on GPU (send filled array to GPU first).
+    Returns:
+        Filtered array with original NaN positions preserved.
+    """
+    is_gpu_arr = HAS_GPU and isinstance(arr, cp.ndarray)
+    arr_np = to_cpu(arr) if is_gpu_arr else arr
+
+    filled, nan_mask = _fill_nans(arr_np)
+
+    if use_gpu and HAS_GPU:
+        filled_gpu = to_gpu(filled)
+        result_gpu = filter_func(filled_gpu, *args, **kwargs)
+        result = to_cpu(result_gpu)
+    else:
+        result = filter_func(filled, *args, **kwargs)
+
+    result[nan_mask] = np.nan
+    return result
+
+
 # ============================================================
 # Core terrain visualizations
 # ============================================================
@@ -179,11 +232,13 @@ def generate_lrm(dem_file, basename, vis_dir, resolution):
 
     try:
         dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np)
+        nan_mask = np.isnan(dem_np)
-        local_mean = xp_gaussian_filter(dem, sigma=15/resolution)
+        local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15/resolution)
-        lrm = dem - local_mean
+        lrm = dem_np - local_mean
-        lrm_np = to_cpu(lrm).astype(np.float32)
+        lrm[nan_mask] = np.nan
-        _save_tif(output, lrm_np, transform, crs)
+        _save_tif(output, lrm.astype(np.float32), transform, crs)
+        logger.info(f"  ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
+        return output
         logger.info(f"  ✓ LRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
         return output
     except Exception as e:
@@ -320,21 +375,25 @@ def generate_mslrm(dem_file, basename, vis_dir, resolution):
 
     try:
         dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np)
+        nan_mask = np.isnan(dem_np)
 
         sigmas = [5, 10, 25, 50, 100]
         lrm_stack = []
 
         for sigma in sigmas:
             sigma_px = sigma / resolution
-            local_mean = xp_gaussian_filter(dem, sigma=sigma_px)
+            local_mean = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_px)
-            lrm = dem - local_mean
+            lrm = dem_np - local_mean
-            lrm_norm = lrm / max(float(xp.nanstd(lrm)), 0.01)
+            lrm[nan_mask] = np.nan
-            lrm_stack.append(lrm_norm)
+            valid_lrm = lrm[~nan_mask]
+            lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
+            lrm_norm = lrm / lrm_std
+            lrm_stack.append(lrm_norm.astype(np.float32))
 
-        mslrm = xp.sqrt(xp.mean(xp.array(lrm_stack) ** 2, axis=0))
+        lrm_array = np.array(lrm_stack)
-        mslrm_np = to_cpu(mslrm).astype(np.float32)
+        mslrm = np.sqrt(np.nanmean(lrm_array ** 2, axis=0))
-        _save_tif(output, mslrm_np, transform, crs)
+        mslrm[nan_mask] = np.nan
+        _save_tif(output, mslrm.astype(np.float32), transform, crs)
         logger.info(f"  ✓ MSRM terminé ({time.time()-t0:.1f}s){gpu_tag}")
         return output
     except Exception as e:
@@ -355,24 +414,28 @@ def generate_tpi(dem_file, basename, vis_dir, resolution):
 
     try:
         dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np)
+        nan_mask = np.isnan(dem_np)
 
         fine_size = int(5 / resolution)
         if fine_size % 2 == 0:
             fine_size += 1
-        tpi_fine = dem - xp_uniform_filter(dem, size=fine_size)
+        fine_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=fine_size)
+        tpi_fine = dem_np - fine_mean
+        tpi_fine[nan_mask] = np.nan
 
         broad_size = int(100 / resolution)
         if broad_size % 2 == 0:
             broad_size += 1
-        tpi_broad = dem - xp_uniform_filter(dem, size=broad_size)
+        broad_mean = _filter_nanaware(dem_np, xp_uniform_filter, size=broad_size)
+        tpi_broad = dem_np - broad_mean
+        tpi_broad[nan_mask] = np.nan
 
-        fine_std = max(float(xp.nanstd(tpi_fine)), 0.01)
+        fine_std = max(np.nanstd(tpi_fine[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
-        broad_std = max(float(xp.nanstd(tpi_broad)), 0.01)
+        broad_std = max(np.nanstd(tpi_broad[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
         tpi_combined = 0.6 * (tpi_fine / fine_std) + 0.4 * (tpi_broad / broad_std)
-        tpi_np = to_cpu(tpi_combined).astype(np.float32)
+        tpi_combined[nan_mask] = np.nan
 
-        _save_tif(output, tpi_np, transform, crs)
+        _save_tif(output, tpi_combined.astype(np.float32), transform, crs)
         logger.info(f"  ✓ TPI terminé ({time.time()-t0:.1f}s){gpu_tag}")
         return output
     except Exception as e:
@@ -448,31 +511,34 @@ def generate_sailore(dem_file, basename, vis_dir, resolution):
 
     try:
         dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np)
+        nan_mask = np.isnan(dem_np)
 
-        gy, gx = xp.gradient(dem, resolution)
+        gy, gx = np.gradient(dem_np, resolution)
-        slope = xp.arctan(xp.sqrt(gx**2 + gy**2))
+        slope = np.arctan(np.sqrt(gx**2 + gy**2))
-        slope_deg = xp.degrees(slope)
+        slope_deg = np.degrees(slope)
+        slope_deg[nan_mask] = np.nan
 
         sigma_min = 2.0 / resolution
         sigma_max = 25.0 / resolution
-        slope_norm = xp.clip(slope_deg / 30.0, 0, 1)
+        slope_norm = np.clip(slope_deg / 30.0, 0, 1)
-        adaptive_sigma = sigma_max - slope_norm * (sigma_max - sigma_min)
 
-        lrm_fine = dem - xp_gaussian_filter(dem, sigma=sigma_min)
+        lrm_fine = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_min)
-        lrm_medium = dem - xp_gaussian_filter(dem, sigma=(sigma_min + sigma_max) / 2)
+        lrm_fine[nan_mask] = np.nan
-        lrm_coarse = dem - xp_gaussian_filter(dem, sigma=sigma_max)
+        lrm_medium = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=(sigma_min + sigma_max) / 2)
+        lrm_medium[nan_mask] = np.nan
+        lrm_coarse = dem_np - _filter_nanaware(dem_np, xp_gaussian_filter, sigma=sigma_max)
+        lrm_coarse[nan_mask] = np.nan
 
         w_fine = slope_norm
-        w_medium = 1 - 2 * xp.abs(slope_norm - 0.5)
+        w_medium = 1 - 2 * np.abs(slope_norm - 0.5)
         w_coarse = 1 - slope_norm
         w_total = w_fine + w_medium + w_coarse
         w_total[w_total == 0] = 1
 
         sailore = (w_fine * lrm_fine + w_medium * lrm_medium + w_coarse * lrm_coarse) / w_total
-        sailore_np = to_cpu(sailore).astype(np.float32)
+        sailore[nan_mask] = np.nan
 
-        _save_tif(output, sailore_np, transform, crs)
+        _save_tif(output, sailore.astype(np.float32), transform, crs)
         logger.info(f"  ✓ SAILORE terminé ({time.time()-t0:.1f}s){gpu_tag}")
         return output
     except Exception as e:
@@ -493,15 +559,16 @@ def generate_roughness(dem_file, basename, vis_dir, resolution):
 
     try:
         dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np.astype(np.float64))
+        nan_mask = np.isnan(dem_np)
         window_size = int(5 / resolution)
         if window_size % 2 == 0:
             window_size += 1
 
-        # Vectorized std: sqrt(E[X²] - (E[X])²) via uniform_filter (GPU-accelerated)
+        # Vectorized std: sqrt(E[X²] - (E[X])²) via uniform_filter (NaN-aware)
-        local_mean = xp_uniform_filter(dem, size=window_size)
+        local_mean = _filter_nanaware(dem_np.astype(np.float64), xp_uniform_filter, size=window_size)
-        local_mean_sq = xp_uniform_filter(dem * dem, size=window_size)
+        local_mean_sq = _filter_nanaware(dem_np.astype(np.float64)**2, xp_uniform_filter, size=window_size)
-        roughness = xp.sqrt(local_mean_sq - local_mean * local_mean)
+        roughness = np.sqrt(np.maximum(local_mean_sq - local_mean * local_mean, 0))
+        roughness[nan_mask] = np.nan
 
         roughness = to_cpu(roughness)
         _save_tif(output, roughness, transform, crs)
@@ -525,24 +592,28 @@ def generate_anomalies(dem_file, basename, vis_dir, resolution):
 
     try:
         dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np)
+        nan_mask = np.isnan(dem_np)
 
-        lrm = dem - xp_gaussian_filter(dem, sigma=15 / resolution)
+        lrm_mean_val = _filter_nanaware(dem_np, xp_gaussian_filter, sigma=15 / resolution)
-        lrm_mean = xp.nanmean(lrm)
+        lrm = dem_np - lrm_mean_val
-        lrm_std = max(float(xp.nanstd(lrm)), 0.01)
+        lrm[nan_mask] = np.nan
+        valid_lrm = lrm[~nan_mask]
+        lrm_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0
+        lrm_std = max(np.nanstd(valid_lrm), 0.01) if len(valid_lrm) > 0 else 0.01
         z_score = (lrm - lrm_mean) / lrm_std
 
         window = int(10 / resolution)
         if window % 2 == 0:
             window += 1
 
-        local_mean = xp_uniform_filter(z_score, size=window)
+        local_mean = _filter_nanaware(z_score, xp_uniform_filter, size=window)
-        z_mean = xp.nanmean(z_score)
+        z_mean = np.nanmean(valid_lrm) if len(valid_lrm) > 0 else 0
-        z_std = max(float(xp.nanstd(z_score)), 0.01)
+        z_std = max(np.nanstd(z_score[~nan_mask]), 0.01) if np.any(~nan_mask) else 0.01
         morans_i = z_score * (local_mean - z_mean) / z_std
-        anomaly_score = xp.abs(z_score) * xp.sign(morans_i)
+        anomaly_score = np.abs(z_score) * np.sign(morans_i)
+        anomaly_score[nan_mask] = np.nan
 
-        _save_tif(output, to_cpu(anomaly_score), transform, crs)
+        _save_tif(output, anomaly_score.astype(np.float32), transform, crs)
         logger.info(f"  ✓ Anomalies terminé ({time.time()-t0:.1f}s){gpu_tag}")
         return output
     except Exception as e:
@@ -566,26 +637,31 @@ def generate_wavelet(dem_file, basename, vis_dir, resolution):
 
     try:
         dem_np, transform, crs = _read_dem(dem_file)
-        dem = to_gpu(dem_np)
+        nan_mask = np.isnan(dem_np)
 
         scales = [2, 5, 10, 20, 50]
         wavelet_stack = []
 
         for scale_m in scales:
             sigma_px = scale_m / resolution
+            filled, _ = _fill_nans(dem_np.astype(np.float64))
             if HAS_GPU:
                 from cupyx.scipy.ndimage import gaussian_laplace as gpu_gaussian_laplace
-                response = -gpu_gaussian_laplace(dem, sigma=sigma_px)
+                response = -gpu_gaussian_laplace(to_gpu(filled), sigma=sigma_px)
+                response = to_cpu(response)
             else:
                 from scipy.ndimage import gaussian_laplace
-                response = to_gpu(-gaussian_laplace(to_cpu(dem).astype(np.float64), sigma=sigma_px))
+                response = -gaussian_laplace(filled, sigma=sigma_px)
-            response /= max(float(xp.nanstd(response)), 0.01)
+            response[nan_mask] = np.nan
+            valid = response[~nan_mask]
+            std_val = max(np.nanstd(valid), 0.01) if len(valid) > 0 else 0.01
+            response = response / std_val
             wavelet_stack.append(response)
 
-        combined = xp.sqrt(xp.mean(xp.array(wavelet_stack) ** 2, axis=0))
+        combined = np.sqrt(np.nanmean(np.array(wavelet_stack) ** 2, axis=0))
-        combined_np = to_cpu(combined).astype(np.float32)
+        combined[nan_mask] = np.nan
 
-        _save_tif(output, combined_np, transform, crs)
+        _save_tif(output, combined.astype(np.float32), transform, crs)
         logger.info(f"  ✓ Ondelette terminée ({time.time()-t0:.1f}s){gpu_tag}")
         return output
     except Exception as e: