Remove RRIM and Multi-Hillshade RGB, fix DTM resolution reuse bug, add --init to docker run
- Remove generate_rrim, generate_multi_hillshade, _compute_openness_both - Remove corresponding VIZ_STEPS entries, COLORMAPS, RGB_LEGENDS, and tests - Fix DTM resolution mismatch: existing DTM at different resolution is now regenerated instead of silently reused - Propagate actual DTM resolution to visualizations and rendering - Add --init to docker run commands for proper signal handling on Ctrl+C - Add .playwright-mcp/ to .gitignore Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
This commit is contained in:
Jacquin Antoine
@ -201,11 +201,11 @@ def _filter_nanaware(arr, filter_func, *args, use_gpu=True, **kwargs):
|
||||
# ============================================================
|
||||
|
||||
def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
"""Generate multi-directional hillshade with slope shading — GPU if available.
|
||||
"""Generate multi-directional hillshade with contrast enhancement — GPU if available.
|
||||
|
||||
Combines 4-direction hillshade (NW, NE, SW, SE) with slope shading
|
||||
for improved micro-relief visibility on flat terrain.
|
||||
Result = 0.7 * hillshade + 0.3 * cos(slope).
|
||||
Combines 4-direction hillshade (NW, NE, SW, SE) with slope shading.
|
||||
Applies percentile normalization and gamma correction to restore
|
||||
contrast lost by averaging multiple azimuths.
|
||||
"""
|
||||
gpu_tag = " [GPU]" if HAS_GPU else ""
|
||||
logger.info(f" → Hillshade multidirectionnel{gpu_tag}...")
|
||||
@ -249,8 +249,20 @@ def generate_hillshade(dem_file, basename, vis_dir, resolution, shared=None):
|
||||
slope_shaded = cos_slope
|
||||
combined = 0.7 * combined_hillshade + 0.3 * slope_shaded
|
||||
|
||||
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)
|
||||
# Contrast enhancement: percentile stretch + gamma
|
||||
# Averaging 4 azimuths flattens contrast — this restores it
|
||||
combined_np = to_cpu(combined)
|
||||
nan_mask = shared.nan_mask if shared else np.isnan(to_cpu(dem_np) if HAS_GPU else dem_np)
|
||||
valid = combined_np[~nan_mask]
|
||||
if len(valid) > 0:
|
||||
p2, p98 = np.percentile(valid, 2), np.percentile(valid, 98)
|
||||
if p98 - p2 > 0.01:
|
||||
combined_np = np.clip((combined_np - p2) / (p98 - p2), 0, 1)
|
||||
# Gamma correction to enhance shadows
|
||||
gamma = 0.8
|
||||
combined_np = np.power(combined_np, gamma)
|
||||
|
||||
_save_tif(output, combined_np.astype(np.float32), transform, crs, nan_mask=nan_mask)
|
||||
logger.info(f" ✓ Hillshade terminé ({time.time()-t0:.1f}s){gpu_tag}")
|
||||
return output
|
||||
except Exception as e:
|
||||
@ -528,194 +540,6 @@ def generate_openness(dem_file, basename, vis_dir, resolution, positive=True, sh
|
||||
return None
|
||||
|
||||
|
||||
def _compute_openness_both(dem, resolution, nan_mask, n_dirs=8, radius=50):
|
||||
"""Compute positive and negative openness in one ray-tracing pass.
|
||||
|
||||
Traces rays in n_dirs directions up to radius pixels, measuring:
|
||||
- positive openness: max angle above horizontal to visible terrain
|
||||
- negative openness: max angle below horizontal to visible terrain
|
||||
|
||||
Returns (pos_open, neg_open) as float32 arrays in degrees.
|
||||
NaN mask is applied after computation.
|
||||
"""
|
||||
rows, cols = dem.shape
|
||||
res = resolution
|
||||
max_dist = int(radius / res)
|
||||
|
||||
angles = np.linspace(0, 2 * np.pi, n_dirs, endpoint=False)
|
||||
dx_dir = np.cos(angles)
|
||||
dy_dir = np.sin(angles)
|
||||
|
||||
padded = np.pad(dem, max_dist, mode='constant', constant_values=np.nanmax(dem[~nan_mask]) + 10000 if np.any(~nan_mask) else 0)
|
||||
|
||||
pos_sum = np.zeros_like(dem)
|
||||
neg_sum = np.zeros_like(dem)
|
||||
|
||||
for d_idx in range(n_dirs):
|
||||
ddx, ddy = dx_dir[d_idx], dy_dir[d_idx]
|
||||
max_pos_angle = np.zeros_like(dem)
|
||||
max_neg_angle = np.zeros_like(dem)
|
||||
|
||||
for step in range(1, max_dist + 1):
|
||||
px = int(round(ddx * step))
|
||||
py = int(round(ddy * step))
|
||||
dist_m = np.sqrt((ddx * step * res) ** 2 + (ddy * step * res) ** 2)
|
||||
if dist_m < res * 0.5:
|
||||
continue
|
||||
|
||||
elev_diff = padded[max_dist + py:max_dist + py + rows,
|
||||
max_dist + px:max_dist + px + cols] - dem
|
||||
|
||||
pos_angle = np.arctan2(np.maximum(elev_diff, 0), dist_m)
|
||||
neg_angle = np.arctan2(np.maximum(-elev_diff, 0), dist_m)
|
||||
|
||||
valid = ~np.isnan(elev_diff)
|
||||
max_pos_angle[valid] = np.maximum(max_pos_angle[valid], pos_angle[valid])
|
||||
max_neg_angle[valid] = np.maximum(max_neg_angle[valid], neg_angle[valid])
|
||||
|
||||
pos_sum += max_pos_angle
|
||||
neg_sum += max_neg_angle
|
||||
|
||||
pos_open = np.degrees(pos_sum / n_dirs).astype(np.float32)
|
||||
neg_open = np.degrees(neg_sum / n_dirs).astype(np.float32)
|
||||
pos_open[nan_mask] = np.nan
|
||||
neg_open[nan_mask] = np.nan
|
||||
return pos_open, neg_open
|
||||
|
||||
|
||||
def generate_rrim(dem_file, basename, vis_dir, resolution, shared=None,
|
||||
n_dirs=8, radius=50, pmin=2, pmax=98, contrast=1.5):
|
||||
"""Red Relief Image Map — RGB composite for archaeological prospection.
|
||||
|
||||
Combines slope, positive openness, and negative openness into a single
|
||||
false-color image where:
|
||||
Red = positive openness (ridges, elevated features)
|
||||
Green = inverted slope (flat = bright, steep = dark)
|
||||
Blue = negative openness (depressions, ditches)
|
||||
|
||||
Each channel is normalized via percentiles and enhanced with a gamma curve.
|
||||
"""
|
||||
gpu_tag = " [GPU]" if HAS_GPU else ""
|
||||
logger.info(f" → RRIM (Red Relief Image){gpu_tag}...")
|
||||
t0 = time.time()
|
||||
output = vis_dir / f"{basename}_rrim.tif"
|
||||
|
||||
try:
|
||||
if shared:
|
||||
transform = shared.transform
|
||||
crs = shared.crs
|
||||
dem_np = shared.dem_np
|
||||
nan_mask = shared.nan_mask
|
||||
slope_rad = shared.slope_rad
|
||||
dem_for_ray = to_gpu(shared.filled) if HAS_GPU else shared.filled
|
||||
else:
|
||||
dem_np, transform, crs = _read_dem(dem_file)
|
||||
nan_mask = np.isnan(dem_np)
|
||||
filled, _ = _fill_nans(dem_np)
|
||||
dem_for_ray = to_gpu(filled) if HAS_GPU else filled
|
||||
dy, dx = np.gradient(filled, resolution)
|
||||
slope_rad = np.arctan(np.sqrt(dx**2 + dy**2))
|
||||
|
||||
# Compute both openness values (ray-tracing on filled DEM)
|
||||
pos_open, neg_open = _compute_openness_both(
|
||||
to_cpu(dem_for_ray) if HAS_GPU else dem_for_ray,
|
||||
resolution, nan_mask, n_dirs=n_dirs, radius=radius
|
||||
)
|
||||
|
||||
# Normalize each component to [0, 1] using percentiles
|
||||
slope_deg = np.degrees(slope_rad)
|
||||
slope_deg[nan_mask] = np.nan
|
||||
|
||||
def _normalize(arr, lo, hi):
|
||||
valid = arr[~np.isnan(arr)]
|
||||
if len(valid) == 0:
|
||||
return np.zeros_like(arr, dtype=np.float32)
|
||||
vlo = np.percentile(valid, lo)
|
||||
vhi = np.percentile(valid, hi)
|
||||
if vhi - vlo < 1e-6:
|
||||
return np.full_like(arr, 0.5, dtype=np.float32)
|
||||
norm = np.clip((arr - vlo) / (vhi - vlo), 0, 1)
|
||||
# Apply gamma for contrast
|
||||
norm = np.power(norm, 1.0 / contrast)
|
||||
return norm.astype(np.float32)
|
||||
|
||||
r = _normalize(pos_open, pmin, pmax) # Red: positive openness (ridges)
|
||||
g = _normalize(90 - slope_deg, pmin, pmax) # Green: inverted slope (flat=bright)
|
||||
g[nan_mask] = np.nan
|
||||
b = _normalize(neg_open, pmin, pmax) # Blue: negative openness (ditches)
|
||||
|
||||
# Assemble RGB (uint8)
|
||||
rgb = np.stack([r, g, b], axis=0) # (3, H, W)
|
||||
rgb = np.nan_to_num(rgb, nan=0.0)
|
||||
rgb_uint8 = (np.clip(rgb, 0, 1) * 255).astype(np.uint8)
|
||||
|
||||
_save_tif(output, rgb_uint8, transform, crs, dtype='uint8', count=3)
|
||||
logger.info(f" ✓ RRIM terminé ({time.time()-t0:.1f}s){gpu_tag}")
|
||||
return output
|
||||
except Exception as e:
|
||||
logger.error(f" ✗ Erreur RRIM: {e}", exc_info=True)
|
||||
return None
|
||||
|
||||
|
||||
def generate_multi_hillshade(dem_file, basename, vis_dir, resolution, shared=None,
|
||||
azimuths=(315, 135, 45), altitude=30, blend_slope=0.3):
|
||||
"""Multi-directional hillshade RGB composite — 3 azimuths mapped to R/G/B.
|
||||
|
||||
Each azimuth produces a hillshade mapped to a color channel:
|
||||
Red = azimuth 315° (NW illumination)
|
||||
Green = azimuth 135° (SE illumination)
|
||||
Blue = azimuth 45° (NE illumination)
|
||||
|
||||
Shadow direction reveals structure orientation through color.
|
||||
"""
|
||||
gpu_tag = " [GPU]" if HAS_GPU else ""
|
||||
logger.info(f" → Hillshade Composite RGB{gpu_tag}...")
|
||||
t0 = time.time()
|
||||
output = vis_dir / f"{basename}_multi_hillshade.tif"
|
||||
|
||||
try:
|
||||
if shared:
|
||||
transform = shared.transform
|
||||
crs = shared.crs
|
||||
nan_mask = shared.nan_mask
|
||||
slope_rad = to_gpu(shared.slope_rad) if HAS_GPU else shared.slope_rad
|
||||
aspect = to_gpu(shared.aspect) if HAS_GPU else shared.aspect
|
||||
else:
|
||||
dem_np, transform, crs = _read_dem(dem_file)
|
||||
nan_mask = np.isnan(dem_np)
|
||||
filled, _ = _fill_nans(dem_np)
|
||||
dem = to_gpu(filled) if HAS_GPU else filled
|
||||
dy, dx = xp.gradient(dem, resolution)
|
||||
slope_rad = xp.arctan(xp.sqrt(dx**2 + dy**2))
|
||||
aspect = xp.arctan2(dy, dx)
|
||||
|
||||
alt_rad = xp.radians(xp.array(altitude))
|
||||
sin_alt = xp.sin(alt_rad)
|
||||
cos_alt = xp.cos(alt_rad)
|
||||
cos_slope = xp.cos(slope_rad)
|
||||
|
||||
channels = []
|
||||
for az in azimuths:
|
||||
az_rad = xp.radians(xp.array(az))
|
||||
hs = sin_alt * xp.sin(slope_rad) + cos_alt * cos_slope * xp.cos(az_rad - aspect)
|
||||
blended = (1 - blend_slope) * xp.clip(hs, 0, 1) + blend_slope * cos_slope
|
||||
channels.append(to_cpu(blended).astype(np.float32))
|
||||
|
||||
gpu_cleanup()
|
||||
|
||||
# Assemble RGB (uint8)
|
||||
rgb = np.stack(channels, axis=0) # (3, H, W)
|
||||
rgb[:, nan_mask] = 0.0
|
||||
rgb_uint8 = (np.clip(rgb, 0, 1) * 255).astype(np.uint8)
|
||||
|
||||
_save_tif(output, rgb_uint8, transform, crs, dtype='uint8', count=3)
|
||||
logger.info(f" ✓ Hillshade Composite RGB terminé ({time.time()-t0:.1f}s){gpu_tag}")
|
||||
return output
|
||||
except Exception as e:
|
||||
logger.error(f" ✗ Erreur multi_hillshade: {e}", exc_info=True)
|
||||
return None
|
||||
|
||||
|
||||
def generate_local_dominance(dem_file, basename, vis_dir, resolution, shared=None,
|
||||
radius=15, pmin=2, pmax=98):
|
||||
"""Local Dominance — proportion of neighborhood below center point.
|
||||
|
||||
Reference in New Issue
Block a user