Remove PMF, fix NaN in gradient visualizations, fix pos_open/neg_open shared param
- Remove PMF from ground classification options (PDAL recommends SMRF over PMF) - Auto-detection now uses CSF for urban/complex terrain instead of PMF - Add z_std > 30m heuristic to auto-select CSF for complex terrain - Fix pos_open/neg_open lambda missing 'shared' parameter (NameError in workers) - Fix NaN mask not restored in hillshade, slope, aspect, curvature (gradient-based products computed on filled DEM lost NaN transparency) - Add nan_mask parameter to _save_tif for centralized NaN restoration - DTM TIF kept by default (no longer deleted after WebP conversion) Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
This commit is contained in:
Jacquin Antoine
@ -118,15 +118,15 @@ Exemples:
|
||||
help="Reclassifier le sol même si le fichier .las existe déjà"
|
||||
)
|
||||
parser.add_argument(
|
||||
"--no-keep-tif",
|
||||
"--keep-tif",
|
||||
action="store_true",
|
||||
help="Supprimer les fichiers TIFF intermédiaires après conversion WebP (par défaut: conservés)"
|
||||
help="Conserver les fichiers TIFF (DTM + visualisations) pour pouvoir régénérer les WebP sans recalculer"
|
||||
)
|
||||
parser.add_argument(
|
||||
"--ground-classification",
|
||||
choices=["auto", "smrf", "pmf", "csf"],
|
||||
choices=["auto", "smrf", "csf"],
|
||||
default="auto",
|
||||
help="Méthode de classification du sol : auto (détection), smrf, pmf, csf (défaut: auto)"
|
||||
help="Méthode de classification du sol : auto (détection), smrf, csf (défaut: auto)"
|
||||
)
|
||||
parser.add_argument(
|
||||
"--file",
|
||||
@ -176,7 +176,7 @@ Exemples:
|
||||
force=args.force,
|
||||
ground_method=args.ground_classification,
|
||||
force_classify=args.force_classification,
|
||||
keep_tif=not args.no_keep_tif,
|
||||
keep_tif=args.keep_tif,
|
||||
)
|
||||
|
||||
# If --file is specified, process only matching files
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
"""DTM generation from classified LiDAR point clouds.
|
||||
|
||||
Handles ground classification via PDAL/SMRF and DTM rasterisation
|
||||
Handles ground classification via PDAL (SMRF or CSF) and DTM rasterisation
|
||||
using scipy binned_statistic_2d. Zones without LiDAR data remain as NaN.
|
||||
"""
|
||||
|
||||
@ -27,13 +27,13 @@ def _create_ground_pipeline(input_laz, output_las, method):
|
||||
1. Reset Classification to 0
|
||||
2. ELM (Extended Local Minimum) — mark low outliers as noise (Classification=7)
|
||||
3. Statistical outlier removal
|
||||
4. Ground classification (SMRF/PMF/CSF)
|
||||
4. Ground classification (SMRF or CSF)
|
||||
5. Extract ground points (Classification=2)
|
||||
|
||||
Args:
|
||||
input_laz: Path to input LAZ/LAS file.
|
||||
output_las: Path to output classified LAS file.
|
||||
method: Ground classification method ('smrf', 'pmf', or 'csf').
|
||||
method: Ground classification method ('smrf' or 'csf').
|
||||
|
||||
Returns:
|
||||
JSON string of the PDAL pipeline.
|
||||
@ -84,15 +84,6 @@ def _create_ground_pipeline(input_laz, output_las, method):
|
||||
"threshold": 0.5,
|
||||
"scalar": 1.25
|
||||
}
|
||||
elif method == 'pmf':
|
||||
ground_step = {
|
||||
"type": "filters.pmf",
|
||||
"max_window": 33,
|
||||
"slope": 0.15,
|
||||
"max_distance": 2.5,
|
||||
"initial_distance": 0.15,
|
||||
"cell_size": 1.0
|
||||
}
|
||||
elif method == 'csf':
|
||||
ground_step = {
|
||||
"type": "filters.csf",
|
||||
@ -128,11 +119,6 @@ def create_smrf_pipeline(input_laz, output_las):
|
||||
return _create_ground_pipeline(input_laz, output_las, 'smrf')
|
||||
|
||||
|
||||
def create_pmf_pipeline(input_laz, output_las):
|
||||
"""Create a PDAL pipeline JSON for PMF ground classification."""
|
||||
return _create_ground_pipeline(input_laz, output_las, 'pmf')
|
||||
|
||||
|
||||
def create_csf_pipeline(input_laz, output_las):
|
||||
"""Create a PDAL pipeline JSON for CSF ground classification."""
|
||||
return _create_ground_pipeline(input_laz, output_las, 'csf')
|
||||
@ -141,9 +127,9 @@ def create_csf_pipeline(input_laz, output_las):
|
||||
def detect_ground_method(laz_file):
|
||||
"""Detect the best ground classification method based on point cloud statistics.
|
||||
|
||||
Auto-selects between SMRF (natural terrain) and PMF (urban) only.
|
||||
CSF is available only via --ground-classification csf (slow but robust
|
||||
on complex terrain).
|
||||
Auto-selects between SMRF and CSF:
|
||||
- SMRF: fast, robust for most natural terrain (PDAL recommended default)
|
||||
- CSF: cloth simulation, better for complex/urban terrain
|
||||
|
||||
Falls back to SMRF if the file cannot be read or attributes are missing.
|
||||
|
||||
@ -151,7 +137,7 @@ def detect_ground_method(laz_file):
|
||||
laz_file: Path to input LAZ/LAS file.
|
||||
|
||||
Returns:
|
||||
String: 'smrf', 'pmf', or 'csf'
|
||||
String: 'smrf' or 'csf'
|
||||
"""
|
||||
import laspy
|
||||
|
||||
@ -182,13 +168,16 @@ def detect_ground_method(laz_file):
|
||||
f"ratio_retours_uniques={single_return_ratio:.2f}, "
|
||||
f"écart_Z={z_std:.1f}m, amplitude_Z={z_range:.1f}m")
|
||||
|
||||
# Decision logic (auto selects between SMRF and PMF only):
|
||||
# - High single-return ratio (>0.6) → urban (buildings, roads) → PMF
|
||||
# Decision logic:
|
||||
# - High single-return ratio (>0.6) → urban (buildings, roads) → CSF (cloth simulation)
|
||||
# - High elevation variance (>30m) → complex/mountainous terrain → CSF
|
||||
# - Default → SMRF (fast, robust for most natural terrain)
|
||||
# Note: CSF is available only via --ground-classification csf (slow but robust on complex terrain)
|
||||
if single_return_ratio > 0.6:
|
||||
method = 'pmf'
|
||||
method = 'csf'
|
||||
reason = f"ratio retours uniques={single_return_ratio:.2f} > 0.6 → milieu urbain"
|
||||
elif z_std > 30:
|
||||
method = 'csf'
|
||||
reason = f"écart_Z={z_std:.1f}m > 30m → terrain complexe"
|
||||
else:
|
||||
method = 'smrf'
|
||||
reason = f"terrain naturel standard"
|
||||
@ -203,7 +192,7 @@ def classify_ground(laz_file, temp_dir, method='auto', force=False):
|
||||
Args:
|
||||
laz_file: Path to input LAZ/LAS file.
|
||||
temp_dir: Directory for temporary files (pipeline.json, ground.las).
|
||||
method: Ground classification method ('auto', 'smrf', 'pmf', or 'csf').
|
||||
method: Ground classification method ('auto', 'smrf', or 'csf').
|
||||
force: If True, reclassify even if output file already exists.
|
||||
|
||||
Returns:
|
||||
|
||||
@ -78,8 +78,8 @@ VIZ_STEPS = [
|
||||
('curvature', generate_curvature),
|
||||
('svf', generate_svf),
|
||||
('lrm', generate_lrm),
|
||||
('pos_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=True)),
|
||||
('neg_open', lambda d, b, v, r: generate_openness(d, b, v, r, positive=False)),
|
||||
('pos_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=True, shared=shared)),
|
||||
('neg_open', lambda d, b, v, r, shared=None: generate_openness(d, b, v, r, positive=False, shared=shared)),
|
||||
('mslrm', generate_mslrm),
|
||||
('tpi', generate_tpi),
|
||||
('sailore', generate_sailore),
|
||||
@ -323,9 +323,6 @@ class LidarArchaeoPipeline:
|
||||
t_pdf = time.time() - t4
|
||||
logger.info(f" ✓ Rapport PDF terminé ({t_pdf:.1f}s)")
|
||||
|
||||
# Step 5: Keep DTM TIF for reuse (regenerating WebPs skips classification)
|
||||
# Use --no-keep-tif to delete DTM files after processing
|
||||
|
||||
t_total = time.time() - t_start
|
||||
logger.info(f"✓ {basename} terminé en {t_total:.1f}s")
|
||||
logger.debug(f" Détails: classification={t_classif:.1f}s, DTM={t_dtm:.1f}s, PDF={t_pdf:.1f}s")
|
||||
|
||||
@ -100,8 +100,18 @@ def _filter_nanaware_from_filled(shared, filter_func, *args, **kwargs):
|
||||
return result
|
||||
|
||||
|
||||
def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None):
|
||||
"""Helper to save a 2D or 3D array as GeoTIFF."""
|
||||
def _save_tif(output_path, data, transform, crs, dtype='float32', count=1, nodata=None, nan_mask=None):
|
||||
"""Helper to save a 2D or 3D array as GeoTIFF.
|
||||
|
||||
Args:
|
||||
nan_mask: Optional boolean mask (True=NaN) to apply before saving.
|
||||
Restores NaN zones in gradient-derived products that were
|
||||
computed on the filled DEM.
|
||||
"""
|
||||
if nan_mask is not None:
|
||||
data = np.array(data, dtype=dtype, copy=True)
|
||||
data[nan_mask] = np.nan
|
||||
|
||||
# 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')
|
||||
@ -238,7 +248,9 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
combined_hillshade = xp.mean(xp.array(hillshades), axis=0)
|
||||
slope_shaded = cos_slope
|
||||
combined = 0.7 * combined_hillshade + 0.3 * slope_shaded
|
||||
_save_tif(output, to_cpu(combined), transform, crs)
|
||||
|
||||
nan_mask = shared.nan_mask if shared else np.isnan(to_cpu(dem_np))
|
||||
_save_tif(output, to_cpu(combined), transform, crs, nan_mask=nan_mask)
|
||||
logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s){gpu_tag}")
|
||||
return output
|
||||
except Exception as e:
|
||||
@ -258,6 +270,7 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
transform = shared.transform
|
||||
crs = shared.crs
|
||||
slope = shared.slope_deg
|
||||
nan_mask = shared.nan_mask
|
||||
if HAS_GPU:
|
||||
slope = to_gpu(slope)
|
||||
else:
|
||||
@ -265,7 +278,8 @@ def generate_slope(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
dem = to_gpu(dem_np)
|
||||
dy, dx = xp.gradient(dem)
|
||||
slope = xp.arctan(xp.sqrt(dx**2 + dy**2)) * 180 / xp.pi
|
||||
_save_tif(output, to_cpu(slope) if HAS_GPU else slope, transform, crs)
|
||||
nan_mask = np.isnan(dem_np)
|
||||
_save_tif(output, to_cpu(slope) if HAS_GPU else slope, transform, crs, nan_mask=nan_mask)
|
||||
logger.info(f" ✓ Pente terminée ({time.time()-t0:.1f}s){gpu_tag}")
|
||||
return output
|
||||
except Exception as e:
|
||||
@ -285,6 +299,7 @@ def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
transform = shared.transform
|
||||
crs = shared.crs
|
||||
aspect = shared.aspect
|
||||
nan_mask = shared.nan_mask
|
||||
if HAS_GPU:
|
||||
aspect = to_gpu(aspect)
|
||||
else:
|
||||
@ -293,7 +308,8 @@ def generate_aspect(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
dy, dx = xp.gradient(dem)
|
||||
aspect = xp.arctan2(dy, dx) * 180 / xp.pi
|
||||
aspect = xp.mod(aspect, 360)
|
||||
_save_tif(output, to_cpu(aspect) if HAS_GPU else aspect, transform, crs)
|
||||
nan_mask = np.isnan(dem_np)
|
||||
_save_tif(output, to_cpu(aspect) if HAS_GPU else aspect, transform, crs, nan_mask=nan_mask)
|
||||
logger.info(f" ✓ Aspect terminé ({time.time()-t0:.1f}s){gpu_tag}")
|
||||
return output
|
||||
except Exception as e:
|
||||
@ -314,6 +330,7 @@ def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
crs = shared.crs
|
||||
dx = shared.dx
|
||||
dy = shared.dy
|
||||
nan_mask = shared.nan_mask
|
||||
if HAS_GPU:
|
||||
dx = to_gpu(dx)
|
||||
dy = to_gpu(dy)
|
||||
@ -321,10 +338,11 @@ def generate_curvature(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
dem_np, transform, crs = _read_dem(dem_file)
|
||||
dem = to_gpu(dem_np)
|
||||
dy, dx = xp.gradient(dem)
|
||||
nan_mask = np.isnan(dem_np)
|
||||
d2z_dx2 = xp.gradient(dx, axis=1)
|
||||
d2z_dy2 = xp.gradient(dy, axis=0)
|
||||
curvature = (d2z_dx2 + d2z_dy2) / 2
|
||||
_save_tif(output, to_cpu(curvature), transform, crs)
|
||||
_save_tif(output, to_cpu(curvature), transform, crs, nan_mask=nan_mask)
|
||||
logger.info(f" ✓ Courbure terminée ({time.time()-t0:.1f}s){gpu_tag}")
|
||||
return output
|
||||
except Exception as e:
|
||||
|
||||
Reference in New Issue
Block a user