diff --git a/lidar_pipeline/dtm.py b/lidar_pipeline/dtm.py
index b5129df..18d3b4d 100644
--- a/lidar_pipeline/dtm.py
+++ b/lidar_pipeline/dtm.py
@@ -1,7 +1,7 @@
 """DTM generation from classified LiDAR point clouds.
 
 Handles ground classification via PDAL/SMRF and DTM rasterisation
-using scipy binned_statistic_2d with gap-filling interpolation.
+using scipy binned_statistic_2d. Zones without LiDAR data remain as NaN.
 """
 
 import json
@@ -12,10 +12,7 @@ from pathlib import Path
 import numpy as np
 import rasterio
 from rasterio.transform import from_bounds
-from scipy import ndimage as scipy_ndimage
-from scipy.ndimage import distance_transform_edt, gaussian_filter
 from scipy.stats import binned_statistic_2d
-from scipy import interpolate
 
 logger = logging.getLogger("lidar")
 
@@ -137,34 +134,10 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution):
         dtm = stat.statistic.T
         dtm = dtm[::-1, :]  # Flip Y so north is at top
 
-        # Fill gaps using interpolation
-        mask = np.isnan(dtm)
-
-        if np.any(mask):
-            logger.info("    Remplissage des zones vides...")
-
-            dist_to_nearest = distance_transform_edt(mask)
-            max_fill_distance = max(20, int(10 / resolution))
-
-            y_coords, x_coords = np.where(~mask)
-            z_values = dtm[~mask]
-
-            if len(y_coords) > 100:
-                interp = interpolate.NearestNDInterpolator(
-                    np.column_stack((y_coords, x_coords)),
-                    z_values
-                )
-
-                y_missing, x_missing = np.where(mask)
-                dtm[y_missing, x_missing] = interp(y_missing, x_missing)
-
-                dtm_smooth = gaussian_filter(dtm, sigma=2.0)
-
-                fill_mask = mask & (dist_to_nearest <= max_fill_distance)
-                dtm[fill_mask] = dtm_smooth[fill_mask]
-
-                far_mask = mask & (dist_to_nearest > max_fill_distance)
-                dtm[far_mask] = np.nan
+        # No interpolation: keep NaN for zones without LiDAR data
+        nan_count = np.count_nonzero(np.isnan(dtm))
+        if nan_count > 0:
+            logger.info(f"    {nan_count:,} pixels sans données (conservés en NaN)")
 
         # Save as GeoTIFF
         output_tif = dtm_dir / f"{basename}_dtm.tif"
@@ -175,6 +148,7 @@ def create_dtm_fast(las_file, basename, dtm_dir, resolution):
             driver='GTiff', height=height, width=width,
             count=1, dtype='float32',
             crs='EPSG:2154', transform=transform,
+            nodata=float('nan'),
             compress='lzw'
         ) as dst:
             dst.write(dtm.astype('float32'), 1)
diff --git a/lidar_pipeline/visualizations.py b/lidar_pipeline/visualizations.py
index 1b99ae5..406b601 100644
--- a/lidar_pipeline/visualizations.py
+++ b/lidar_pipeline/visualizations.py
@@ -25,14 +25,19 @@ else:
     xp = np
 
 
-def _save_tif(output_path, data, transform, crs, dtype='float32', count=1):
+def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None):
     """Helper to save a 2D or 3D array as GeoTIFF."""
+    # Auto-detect nodata for float types with NaN
+    if nodata is None and dtype.startswith('float') and np.any(np.isnan(data)):
+        nodata = float('nan')
+
     if data.ndim == 2:
         height, width = data.shape
         with rasterio.open(
             output_path, 'w', driver='GTiff',
             height=height, width=width, count=count,
-            dtype=dtype, crs=crs, transform=transform, compress='lzw'
+            dtype=dtype, crs=crs, transform=transform,
+            compress='lzw', nodata=nodata
         ) as dst:
             dst.write(data.astype(dtype), 1)
     elif data.ndim == 3:
@@ -40,7 +45,8 @@ def _save_tif(output_path, data, transform, crs, dtype='float32', count=1):
         with rasterio.open(
             output_path, 'w', driver='GTiff',
             height=height, width=width, count=bands,
-            dtype=dtype, crs=crs, transform=transform, compress='lzw'
+            dtype=dtype, crs=crs, transform=transform,
+            compress='lzw', nodata=nodata
         ) as dst:
             for i in range(bands):
                 dst.write(data[i].astype(dtype), i + 1)