DTTrierFlood/GeoData
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fix building texture mapping and tree water-mask fallback

Browse Source
This commit is contained in:
Matthias Puchstein
2026-02-10 01:09:38 +01:00
parent c63f21cf81
commit 282809d0d8
2 changed files with 179 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

90
geodata_pipeline/buildings.py
View File

@@ -364,6 +364,38 @@ def _bounds_from_row(row: Dict[str, str]) -> Tuple[float, float, float, float]:
return float(row["xmin"]), float(row["ymin"]), float(row["xmax"]), float(row["ymax"])
def _texture_bounds_from_geometry(
source_xy: np.ndarray,
manifest_bounds: Tuple[float, float, float, float],
tile_size_m: float,
) -> Tuple[float, float, float, float]:
"""Expand texture sampling bounds to the snapped extent of the tile geometry."""
xmin, ymin, xmax, ymax = manifest_bounds
if source_xy.size == 0 or tile_size_m <= 0:
return manifest_bounds
world_x = source_xy[:, 0] + xmin
world_y = source_xy[:, 1] + ymin
if world_x.size == 0 or world_y.size == 0:
return manifest_bounds
min_x = float(np.min(world_x))
min_y = float(np.min(world_y))
max_x = float(np.max(world_x))
max_y = float(np.max(world_y))
if not (math.isfinite(min_x) and math.isfinite(min_y) and math.isfinite(max_x) and math.isfinite(max_y)):
return manifest_bounds
tex_xmin = math.floor(min_x / tile_size_m) * tile_size_m
tex_ymin = math.floor(min_y / tile_size_m) * tile_size_m
tex_xmax = math.ceil(max_x / tile_size_m) * tile_size_m
tex_ymax = math.ceil(max_y / tile_size_m) * tile_size_m
if tex_xmax <= tex_xmin or tex_ymax <= tex_ymin:
return manifest_bounds
return tex_xmin, tex_ymin, tex_xmax, tex_ymax
def _run(cmd: list[str], desc: str) -> bool:
# If the first element is a command like "uv run cjio", we split it.
# Subsequent elements (arguments) are kept as-is to preserve paths with spaces.
@@ -602,6 +634,9 @@ def export_buildings(cfg: Config) -> int:
xmin, ymin, xmax, ymax = bounds
tile_size_m = float(getattr(cfg.heightmap, "tile_size_m", 1000.0) or 1000.0)
if tile_size_m <= 0:
tile_size_m = 1000.0
w = xmax - xmin
h = ymax - ymin
if w == 0 or h == 0:
@@ -609,6 +644,19 @@ def export_buildings(cfg: Config) -> int:
continue
# Convert to glTF-friendly axes: x=east, y=height (z), z=-north
source_xy = vertices[:, :2].copy()
world_x = source_xy[:, 0] + xmin
world_y = source_xy[:, 1] + ymin
texture_bounds = _texture_bounds_from_geometry(source_xy, bounds, tile_size_m)
tex_xmin, tex_ymin, tex_xmax, tex_ymax = texture_bounds
tex_w = tex_xmax - tex_xmin
tex_h = tex_ymax - tex_ymin
if tex_w <= 0 or tex_h <= 0:
print(f"[buildings] invalid texture bounds for {tile_id}, falling back to manifest bounds")
texture_bounds = bounds
tex_xmin, tex_ymin, tex_xmax, tex_ymax = texture_bounds
tex_w = tex_xmax - tex_xmin
tex_h = tex_ymax - tex_ymin
gltf_vertices = np.zeros_like(vertices)
gltf_vertices[:, 0] = vertices[:, 0]
gltf_vertices[:, 1] = vertices[:, 2]
@@ -617,8 +665,9 @@ def export_buildings(cfg: Config) -> int:
roof_normals = _compute_normals(gltf_vertices, roof_faces) if roof_faces.size else np.zeros_like(gltf_vertices)
wall_normals = _compute_normals(gltf_vertices, wall_faces) if wall_faces.size else np.zeros_like(gltf_vertices)
uv = np.zeros((len(vertices), 2), dtype=np.float32)
uv[:, 0] = source_xy[:, 0] / w
uv[:, 1] = 1.0 - (source_xy[:, 1] / h)
uv[:, 0] = (world_x - tex_xmin) / tex_w
uv[:, 1] = 1.0 - ((world_y - tex_ymin) / tex_h)
np.clip(uv, 0.0, 1.0, out=uv)
# Wall colors sampled from ortho if available (fallback constant)
wall_color = np.zeros((len(vertices), 3), dtype=np.float32) + 0.75
@@ -627,20 +676,33 @@ def export_buildings(cfg: Config) -> int:
ortho_bytes = None
if ortho_vrt_ds:
try:
# Calculate resolution based on tile size relative to 1km base
tile_w_km = (xmax - xmin) / 1000.0
target_res = int(round(cfg.ortho.out_res * tile_w_km))
ortho_bytes = _extract_texture_from_vrt(bounds, ortho_vrt_ds, width_px=target_res)
coverage_tiles = tex_w / tile_size_m
target_res = int(round(cfg.ortho.out_res * coverage_tiles))
target_res = max(cfg.ortho.out_res, target_res)
max_res = max(cfg.ortho.out_res, cfg.ortho.out_res * 8)
if target_res > max_res:
print(
f"[buildings] clamping texture resolution for {tile_id}: "
f"{target_res} -> {max_res}"
)
target_res = max_res
print(
f"[buildings] texture map for {tile_id}: "
f"manifest=({xmin:.2f},{ymin:.2f},{xmax:.2f},{ymax:.2f}) "
f"sample=({tex_xmin:.2f},{tex_ymin:.2f},{tex_xmax:.2f},{tex_ymax:.2f}) "
f"res={target_res}"
)
ortho_bytes = _extract_texture_from_vrt(texture_bounds, ortho_vrt_ds, width_px=target_res)
if ortho_bytes:
# Also sample wall colors from the VRT (already implemented using windowed read)
gt_o = ortho_vrt_ds.GetGeoTransform()
# Calculate srcwin for the tile
xoff_o = int((xmin - gt_o[0]) / gt_o[1])
yoff_o = int((ymax - gt_o[3]) / gt_o[5])
xsize_o = int((xmax - xmin) / gt_o[1]) + 1
ysize_o = int((ymax - ymin) / abs(gt_o[5])) + 1
# Calculate srcwin for the sampled texture bounds
xoff_o = int((tex_xmin - gt_o[0]) / gt_o[1])
yoff_o = int((tex_ymax - gt_o[3]) / gt_o[5])
xsize_o = int((tex_xmax - tex_xmin) / gt_o[1]) + 1
ysize_o = int((tex_ymax - tex_ymin) / abs(gt_o[5])) + 1
# Clamp to raster dimensions
xoff_o_clamped = max(0, min(xoff_o, ortho_vrt_ds.RasterXSize - 1))
@@ -653,9 +715,9 @@ def export_buildings(cfg: Config) -> int:
for i in range(min(3, ortho_vrt_ds.RasterCount))
]
for idx, (vx, vy, _) in enumerate(vertices):
wx = vx + xmin
wy = vy + ymin
for idx in range(len(vertices)):
wx = world_x[idx]
wy = world_y[idx]
col = int((wx - gt_o[0]) / gt_o[1]) - xoff_o_clamped
row = int((wy - gt_o[3]) / gt_o[5]) - yoff_o_clamped

112
geodata_pipeline/trees.py
View File

@@ -143,20 +143,108 @@ def _bridge_mask_from_chm(
return bridge.astype(np.uint8) if np.any(bridge) else None
def _water_mask(tile_id: str, bounds: Tuple[float, float, float, float], like_ds: gdal.Dataset, cfg: Config) -> np.ndarray | None:
def _water_mask_candidates(tile_id: str) -> List[str]:
base_id = tile_id.replace("_1_rp", "_rp")
raw = [
f"{tile_id}.png",
f"{tile_id}_mask.png",
f"{tile_id}_mask_viz.png",
f"{tile_id}_viz.png",
]
if base_id != tile_id:
raw.extend(
[
f"{base_id}.png",
f"{base_id}_mask.png",
f"{base_id}_mask_viz.png",
f"{base_id}_viz.png",
]
)
# Preserve order while removing accidental duplicates.
return list(dict.fromkeys(raw))
def _water_mask_from_dir(tile_id: str, like_ds: gdal.Dataset, search_dir: str) -> np.ndarray | None:
if not os.path.isdir(search_dir):
return None
mask_path = None
for candidate in _water_mask_candidates(tile_id):
candidate_path = os.path.join(search_dir, candidate)
if os.path.exists(candidate_path):
mask_path = candidate_path
break
if not mask_path:
return None
src_ds = gdal.Open(mask_path)
if src_ds is None:
return None
like_gt = like_ds.GetGeoTransform()
width = like_ds.RasterXSize
height = like_ds.RasterYSize
xmin = like_gt[0]
ymax = like_gt[3]
xmax = xmin + like_gt[1] * width
ymin = ymax + like_gt[5] * height
like_proj = like_ds.GetProjection() or ""
src_proj = src_ds.GetProjection() or ""
warp_kwargs = {
"format": "MEM",
"outputBounds": (xmin, ymin, xmax, ymax),
"width": width,
"height": height,
"resampleAlg": "near",
"dstNodata": 0,
}
if like_proj:
warp_kwargs["dstSRS"] = like_proj
if not src_proj:
warp_kwargs["srcSRS"] = like_proj
warped = gdal.Warp("", src_ds, options=gdal.WarpOptions(**warp_kwargs))
if warped is None or warped.RasterCount == 0:
return None
if warped.RasterCount >= 3:
blue = warped.GetRasterBand(3).ReadAsArray()
mask = (blue > 0).astype(np.uint8)
else:
band = warped.GetRasterBand(1).ReadAsArray()
mask = (band > 0).astype(np.uint8)
return mask if np.any(mask) else None
def _water_mask(
tile_id: str,
bounds: Tuple[float, float, float, float],
like_ds: gdal.Dataset,
cfg: Config,
) -> Tuple[np.ndarray | None, str]:
# Preferred source order:
# 1) curated raw masks, 2) generated river masks, 3) LiDAR classification fallback.
for search_dir, label in (("raw/water_masks", "raw"), ("work/river_masks", "river")):
mask = _water_mask_from_dir(tile_id, like_ds, search_dir)
if mask is not None:
dilate_px = max(1, int(round(1.5 / max(cfg.trees.grid_res_m, 0.1))))
mask = _dilate_mask(mask, dilate_px)
return mask, label
lidar_cfg = getattr(cfg.river_erosion, "lidar", None)
source_dir = getattr(lidar_cfg, "source_dir", "raw/bdom20rgbi")
water_class = getattr(lidar_cfg, "classification_water", 9)
parts = tile_id.split("_")
if len(parts) < 6:
return None
return None, "none"
x_idx = parts[2]
y_idx = parts[3]
laz_name = f"bdom20rgbi_32_{x_idx}_{y_idx}_2_rp.laz"
laz_path = os.path.join(source_dir, laz_name)
if not os.path.exists(laz_path):
return None
return None, "none"
pipeline_json = [
{"type": "readers.las", "filename": laz_path},
@@ -167,14 +255,14 @@ def _water_mask(tile_id: str, bounds: Tuple[float, float, float, float], like_ds
pipeline = pdal.Pipeline(json.dumps(pipeline_json))
count = pipeline.execute()
except RuntimeError:
return None
return None, "none"
if count == 0:
return None
return None, "none"
arrays = pipeline.arrays
if not arrays:
return None
return None, "none"
xs = arrays[0]["X"]
ys = arrays[0]["Y"]
@@ -190,14 +278,14 @@ def _water_mask(tile_id: str, bounds: Tuple[float, float, float, float], like_ds
row = ((ymax - ys) / yres).astype(int)
valid = (col >= 0) & (col < width) & (row >= 0) & (row < height)
if not np.any(valid):
return None
return None, "none"
mask = np.zeros((height, width), dtype=np.uint8)
mask[row[valid], col[valid]] = 1
dilate_px = max(1, int(round(1.5 / max(cfg.trees.grid_res_m, 0.1))))
mask = _dilate_mask(mask, dilate_px)
return mask
return mask, "lidar"
def _dilate_mask(mask: np.ndarray, radius_px: int) -> np.ndarray:
@@ -785,12 +873,18 @@ def export_trees(cfg: Config, *, force_vrt: bool = False) -> int:
bmask = _dilate_mask(bmask, px_buffer)
chm = np.where(bmask == 1, np.nan, chm)
wmask = _water_mask(tile_id, bounds, dtm_ds, cfg)
wmask, wmask_source = _water_mask(tile_id, bounds, dtm_ds, cfg)
brmask = _bridge_mask_from_chm(chm, wmask, cfg, cfg.trees.grid_res_m)
if wmask is not None:
chm = np.where(wmask == 1, np.nan, chm)
if brmask is not None:
chm = np.where(brmask == 1, np.nan, chm)
water_pixels = int(np.sum(wmask > 0)) if wmask is not None else 0
bridge_pixels = int(np.sum(brmask > 0)) if brmask is not None else 0
print(
f"[trees] {tile_id}: water_mask_source={wmask_source}, "
f"water_pixels={water_pixels}, bridge_pixels={bridge_pixels}"
)
spacing_px = max(2, int(math.ceil((cfg.trees.grid_res_m * 2.5) / cfg.trees.grid_res_m)))
maxima = _local_maxima(chm, cfg.trees.min_height_m, spacing_px)