From 282809d0d820bf2109d3363c7b4bac3da37b740c Mon Sep 17 00:00:00 2001 From: "s0wlz (Matthias Puchstein)" Date: Tue, 10 Feb 2026 01:09:38 +0100 Subject: [PATCH] Fix building texture mapping and tree water-mask fallback --- geodata_pipeline/buildings.py | 90 ++++++++++++++++++++++----- geodata_pipeline/trees.py | 112 +++++++++++++++++++++++++++++++--- 2 files changed, 179 insertions(+), 23 deletions(-) diff --git a/geodata_pipeline/buildings.py b/geodata_pipeline/buildings.py index 4ae0ce9..0bab661 100644 --- a/geodata_pipeline/buildings.py +++ b/geodata_pipeline/buildings.py @@ -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 diff --git a/geodata_pipeline/trees.py b/geodata_pipeline/trees.py index cb8c19f..0494b4a 100644 --- a/geodata_pipeline/trees.py +++ b/geodata_pipeline/trees.py @@ -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)