Add HydroRIVERS river erosion post-process
This commit is contained in:
@@ -27,6 +27,44 @@ resample = "bilinear"
|
||||
tile_size_m = 1000
|
||||
use_tile_minmax = true
|
||||
|
||||
[river_erosion]
|
||||
enabled = true
|
||||
source_path = "raw/HydroRIVERS_v10_eu_shp/HydroRIVERS_v10_eu.shp"
|
||||
source_layer = "HydroRIVERS_v10_eu"
|
||||
output_dir = "export_unity/height_png16_river"
|
||||
manifest_vr = "export_unity/tile_index_river_vr.csv"
|
||||
manifest_server = "export_unity/tile_index_river_server.csv"
|
||||
|
||||
[river_erosion.vr]
|
||||
order_field = "ORD_STRA"
|
||||
invert_order = false
|
||||
invert_max_order = 10
|
||||
min_order = 3
|
||||
width_base_m = 8.0
|
||||
width_per_order_m = 6.0
|
||||
width_min_m = 4.0
|
||||
width_max_m = 90.0
|
||||
depth_base_m = 0.7
|
||||
depth_per_order_m = 0.5
|
||||
depth_min_m = 0.3
|
||||
depth_max_m = 6.0
|
||||
smooth_sigma_m = 3.0
|
||||
|
||||
[river_erosion.server]
|
||||
order_field = "ORD_FLOW"
|
||||
invert_order = true
|
||||
invert_max_order = 10
|
||||
min_order = 2
|
||||
width_base_m = 5.0
|
||||
width_per_order_m = 7.0
|
||||
width_min_m = 3.0
|
||||
width_max_m = 60.0
|
||||
depth_base_m = 0.4
|
||||
depth_per_order_m = 0.3
|
||||
depth_min_m = 0.2
|
||||
depth_max_m = 4.0
|
||||
smooth_sigma_m = 1.5
|
||||
|
||||
[tile_key]
|
||||
tile_size_x = 1000.0
|
||||
tile_size_y = 1000.0
|
||||
|
||||
@@ -25,6 +25,44 @@ out_res = 1025
|
||||
resample = "bilinear"
|
||||
tile_size_m = 1000
|
||||
|
||||
[river_erosion]
|
||||
enabled = true
|
||||
source_path = "raw/HydroRIVERS_v10_eu_shp/HydroRIVERS_v10_eu.shp"
|
||||
source_layer = "HydroRIVERS_v10_eu"
|
||||
output_dir = "export_unity/height_png16_river"
|
||||
manifest_vr = "export_unity/tile_index_river_vr.csv"
|
||||
manifest_server = "export_unity/tile_index_river_server.csv"
|
||||
|
||||
[river_erosion.vr]
|
||||
order_field = "ORD_STRA"
|
||||
invert_order = false
|
||||
invert_max_order = 10
|
||||
min_order = 3
|
||||
width_base_m = 8.0
|
||||
width_per_order_m = 6.0
|
||||
width_min_m = 4.0
|
||||
width_max_m = 90.0
|
||||
depth_base_m = 0.7
|
||||
depth_per_order_m = 0.5
|
||||
depth_min_m = 0.3
|
||||
depth_max_m = 6.0
|
||||
smooth_sigma_m = 3.0
|
||||
|
||||
[river_erosion.server]
|
||||
order_field = "ORD_FLOW"
|
||||
invert_order = true
|
||||
invert_max_order = 10
|
||||
min_order = 2
|
||||
width_base_m = 5.0
|
||||
width_per_order_m = 7.0
|
||||
width_min_m = 3.0
|
||||
width_max_m = 60.0
|
||||
depth_base_m = 0.4
|
||||
depth_per_order_m = 0.3
|
||||
depth_min_m = 0.2
|
||||
depth_max_m = 4.0
|
||||
smooth_sigma_m = 1.5
|
||||
|
||||
[tile_key]
|
||||
tile_size_x = 1000.0
|
||||
tile_size_y = 1000.0
|
||||
|
||||
@@ -49,6 +49,71 @@ class HeightmapConfig:
|
||||
use_tile_minmax: bool = True
|
||||
|
||||
|
||||
@dataclass
|
||||
class RiverErosionProfileConfig:
|
||||
order_field: str = "ORD_STRA"
|
||||
invert_order: bool = False
|
||||
invert_max_order: int = 10
|
||||
min_order: int = 3
|
||||
width_base_m: float = 8.0
|
||||
width_per_order_m: float = 6.0
|
||||
width_min_m: float = 4.0
|
||||
width_max_m: float = 90.0
|
||||
depth_base_m: float = 0.7
|
||||
depth_per_order_m: float = 0.5
|
||||
depth_min_m: float = 0.3
|
||||
depth_max_m: float = 6.0
|
||||
smooth_sigma_m: float = 3.0
|
||||
|
||||
|
||||
def _default_river_profile_vr() -> "RiverErosionProfileConfig":
|
||||
return RiverErosionProfileConfig(
|
||||
order_field="ORD_STRA",
|
||||
invert_order=False,
|
||||
invert_max_order=10,
|
||||
min_order=3,
|
||||
width_base_m=8.0,
|
||||
width_per_order_m=6.0,
|
||||
width_min_m=4.0,
|
||||
width_max_m=90.0,
|
||||
depth_base_m=0.7,
|
||||
depth_per_order_m=0.5,
|
||||
depth_min_m=0.3,
|
||||
depth_max_m=6.0,
|
||||
smooth_sigma_m=3.0,
|
||||
)
|
||||
|
||||
|
||||
def _default_river_profile_server() -> "RiverErosionProfileConfig":
|
||||
return RiverErosionProfileConfig(
|
||||
order_field="ORD_FLOW",
|
||||
invert_order=True,
|
||||
invert_max_order=10,
|
||||
min_order=2,
|
||||
width_base_m=5.0,
|
||||
width_per_order_m=7.0,
|
||||
width_min_m=3.0,
|
||||
width_max_m=60.0,
|
||||
depth_base_m=0.4,
|
||||
depth_per_order_m=0.3,
|
||||
depth_min_m=0.2,
|
||||
depth_max_m=4.0,
|
||||
smooth_sigma_m=1.5,
|
||||
)
|
||||
|
||||
|
||||
@dataclass
|
||||
class RiverErosionConfig:
|
||||
enabled: bool = True
|
||||
source_path: str = "raw/HydroRIVERS_v10_eu_shp/HydroRIVERS_v10_eu.shp"
|
||||
source_layer: str = "HydroRIVERS_v10_eu"
|
||||
output_dir: str = "export_unity/height_png16_river"
|
||||
manifest_vr: str = "export_unity/tile_index_river_vr.csv"
|
||||
manifest_server: str = "export_unity/tile_index_river_server.csv"
|
||||
vr: RiverErosionProfileConfig = field(default_factory=_default_river_profile_vr)
|
||||
server: RiverErosionProfileConfig = field(default_factory=_default_river_profile_server)
|
||||
|
||||
|
||||
@dataclass
|
||||
class OrthoConfig:
|
||||
out_res: int = 2048
|
||||
@@ -160,6 +225,7 @@ class Config:
|
||||
work: WorkConfig = field(default_factory=WorkConfig)
|
||||
export: ExportConfig = field(default_factory=ExportConfig)
|
||||
heightmap: HeightmapConfig = field(default_factory=HeightmapConfig)
|
||||
river_erosion: RiverErosionConfig = field(default_factory=RiverErosionConfig)
|
||||
ortho: OrthoConfig = field(default_factory=OrthoConfig)
|
||||
tile_key: TileKeyConfig = field(default_factory=TileKeyConfig)
|
||||
buildings: BuildingConfig = field(default_factory=BuildingConfig)
|
||||
@@ -189,6 +255,7 @@ class Config:
|
||||
work=WorkConfig(**_filter_kwargs(WorkConfig, data.get("work", {}))),
|
||||
export=ExportConfig(**_filter_kwargs(ExportConfig, data.get("export", {}))),
|
||||
heightmap=HeightmapConfig(**_filter_kwargs(HeightmapConfig, data.get("heightmap", {}))),
|
||||
river_erosion=_river_erosion_from_dict(data.get("river_erosion", {})),
|
||||
ortho=OrthoConfig(**_filter_kwargs(OrthoConfig, data.get("ortho", {}))),
|
||||
tile_key=TileKeyConfig(**_filter_kwargs(TileKeyConfig, data.get("tile_key", {}))),
|
||||
buildings=BuildingConfig(**_filter_kwargs(BuildingConfig, data.get("buildings", {}))),
|
||||
@@ -254,3 +321,14 @@ def _coerce_ranges(data: Dict[str, Any]) -> Dict[str, Any]:
|
||||
if isinstance(value, list):
|
||||
out[key] = tuple(value)
|
||||
return out
|
||||
|
||||
|
||||
def _river_erosion_from_dict(data: Dict[str, Any]) -> RiverErosionConfig:
|
||||
if not isinstance(data, dict):
|
||||
return RiverErosionConfig()
|
||||
base = _filter_kwargs(RiverErosionConfig, data)
|
||||
base.pop("vr", None)
|
||||
base.pop("server", None)
|
||||
vr_cfg = RiverErosionProfileConfig(**_filter_kwargs(RiverErosionProfileConfig, data.get("vr", {})))
|
||||
server_cfg = RiverErosionProfileConfig(**_filter_kwargs(RiverErosionProfileConfig, data.get("server", {})))
|
||||
return RiverErosionConfig(**base, vr=vr_cfg, server=server_cfg)
|
||||
|
||||
451
geodata_pipeline/river_erosion.py
Normal file
451
geodata_pipeline/river_erosion.py
Normal file
@@ -0,0 +1,451 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import csv
|
||||
import glob
|
||||
import math
|
||||
import os
|
||||
from typing import Dict, Iterable, List, Tuple
|
||||
|
||||
import numpy as np
|
||||
from osgeo import gdal, ogr, osr
|
||||
from scipy import ndimage
|
||||
|
||||
from .config import Config, RiverErosionProfileConfig
|
||||
from .gdal_utils import build_vrt, ensure_dir, ensure_parent, open_dataset
|
||||
|
||||
gdal.UseExceptions()
|
||||
|
||||
|
||||
def _coerce_int(value) -> int | None:
|
||||
if value is None:
|
||||
return None
|
||||
try:
|
||||
return int(value)
|
||||
except (TypeError, ValueError):
|
||||
return None
|
||||
|
||||
|
||||
def _normalize_order(value, profile: RiverErosionProfileConfig) -> int | None:
|
||||
order = _coerce_int(value)
|
||||
if order is None:
|
||||
return None
|
||||
if profile.invert_order:
|
||||
max_order = max(1, int(profile.invert_max_order))
|
||||
if order <= 0:
|
||||
order = max_order
|
||||
order = max(1, min(max_order, order))
|
||||
order = max_order + 1 - order
|
||||
return order
|
||||
|
||||
|
||||
def _order_width_depth(order: int, profile: RiverErosionProfileConfig) -> Tuple[float, float] | None:
|
||||
if order < profile.min_order:
|
||||
return None
|
||||
delta = order - profile.min_order
|
||||
width = profile.width_base_m + profile.width_per_order_m * delta
|
||||
depth = profile.depth_base_m + profile.depth_per_order_m * delta
|
||||
width = min(max(width, profile.width_min_m), profile.width_max_m)
|
||||
depth = min(max(depth, profile.depth_min_m), profile.depth_max_m)
|
||||
return width, depth
|
||||
|
||||
|
||||
def _tile_geotransform(bounds: Tuple[float, float, float, float], out_res: int) -> Tuple[float, float, float, float, float, float]:
|
||||
xmin, ymin, xmax, ymax = bounds
|
||||
xres = (xmax - xmin) / float(out_res)
|
||||
yres = (ymax - ymin) / float(out_res)
|
||||
return xmin, xres, 0.0, ymax, 0.0, -yres
|
||||
|
||||
|
||||
def _pixel_size_m(gt: Tuple[float, ...]) -> float:
|
||||
return (abs(gt[1]) + abs(gt[5])) / 2.0
|
||||
|
||||
|
||||
def _bounds_to_layer(
|
||||
bounds: Tuple[float, float, float, float],
|
||||
buffer_m: float,
|
||||
to_layer: osr.CoordinateTransformation,
|
||||
) -> Tuple[float, float, float, float]:
|
||||
xmin, ymin, xmax, ymax = bounds
|
||||
xmin -= buffer_m
|
||||
ymin -= buffer_m
|
||||
xmax += buffer_m
|
||||
ymax += buffer_m
|
||||
corners = [(xmin, ymin), (xmin, ymax), (xmax, ymin), (xmax, ymax)]
|
||||
xs: list[float] = []
|
||||
ys: list[float] = []
|
||||
for x, y in corners:
|
||||
lon, lat, _ = to_layer.TransformPoint(x, y)
|
||||
xs.append(lon)
|
||||
ys.append(lat)
|
||||
return min(xs), min(ys), max(xs), max(ys)
|
||||
|
||||
|
||||
def _compute_tile_minmax_from_vrt(
|
||||
vrt_ds: gdal.Dataset,
|
||||
bounds: Tuple[float, float, float, float],
|
||||
out_res: int,
|
||||
resample: str,
|
||||
) -> Tuple[float | None, float | None]:
|
||||
opts = gdal.WarpOptions(
|
||||
outputBounds=bounds,
|
||||
width=out_res,
|
||||
height=out_res,
|
||||
resampleAlg=resample,
|
||||
srcNodata=-9999,
|
||||
dstNodata=-9999,
|
||||
format="MEM",
|
||||
)
|
||||
tmp = gdal.Warp("", vrt_ds, options=opts)
|
||||
if tmp is None:
|
||||
return None, None
|
||||
band = tmp.GetRasterBand(1)
|
||||
nodata = band.GetNoDataValue()
|
||||
data = band.ReadAsArray()
|
||||
tmp = None
|
||||
if data is None or data.size == 0:
|
||||
return None, None
|
||||
valid = np.isfinite(data)
|
||||
if nodata is not None and not math.isnan(nodata):
|
||||
valid &= data != nodata
|
||||
if not np.any(valid):
|
||||
return None, None
|
||||
return float(data[valid].min()), float(data[valid].max())
|
||||
|
||||
|
||||
def _build_depth_raster(
|
||||
layer: ogr.Layer,
|
||||
bounds: Tuple[float, float, float, float],
|
||||
out_res: int,
|
||||
profile: RiverErosionProfileConfig,
|
||||
to_tile: osr.CoordinateTransformation,
|
||||
to_layer: osr.CoordinateTransformation,
|
||||
tile_srs: osr.SpatialReference,
|
||||
) -> Tuple[np.ndarray, Tuple[float, float, float, float, float, float]]:
|
||||
gt = _tile_geotransform(bounds, out_res)
|
||||
driver = gdal.GetDriverByName("MEM")
|
||||
depth_ds = driver.Create("", out_res, out_res, 1, gdal.GDT_Float32)
|
||||
depth_ds.SetGeoTransform(gt)
|
||||
depth_ds.SetProjection(tile_srs.ExportToWkt())
|
||||
band = depth_ds.GetRasterBand(1)
|
||||
band.Fill(0.0)
|
||||
|
||||
buffer_m = (profile.width_max_m * 0.5) + (profile.smooth_sigma_m * 3.0)
|
||||
minx, miny, maxx, maxy = _bounds_to_layer(bounds, buffer_m, to_layer)
|
||||
layer.SetSpatialFilterRect(minx, miny, maxx, maxy)
|
||||
layer.ResetReading()
|
||||
|
||||
for feat in layer:
|
||||
order = _normalize_order(feat.GetField(profile.order_field), profile)
|
||||
if order is None:
|
||||
continue
|
||||
width_depth = _order_width_depth(order, profile)
|
||||
if width_depth is None:
|
||||
continue
|
||||
width_m, depth_m = width_depth
|
||||
geom = feat.GetGeometryRef()
|
||||
if geom is None:
|
||||
continue
|
||||
geom = geom.Clone()
|
||||
geom.Transform(to_tile)
|
||||
if geom.IsEmpty():
|
||||
continue
|
||||
buffered = geom.Buffer(width_m * 0.5)
|
||||
if buffered is None or buffered.IsEmpty():
|
||||
continue
|
||||
gdal.RasterizeGeometries(
|
||||
depth_ds,
|
||||
[1],
|
||||
[buffered],
|
||||
burn_values=[depth_m],
|
||||
options=["MERGE_ALG=MAX"],
|
||||
)
|
||||
|
||||
depth = band.ReadAsArray().astype(np.float32)
|
||||
if profile.smooth_sigma_m > 0.0:
|
||||
sigma_px = profile.smooth_sigma_m / max(_pixel_size_m(gt), 1e-6)
|
||||
if sigma_px > 0.0:
|
||||
depth = ndimage.gaussian_filter(depth, sigma=sigma_px, mode="nearest")
|
||||
|
||||
return depth, gt
|
||||
|
||||
|
||||
def _write_manifest(
|
||||
manifest_path: str,
|
||||
rows: Iterable[Dict[str, str]],
|
||||
global_min: float,
|
||||
global_max: float,
|
||||
) -> None:
|
||||
ensure_parent(manifest_path)
|
||||
header = [
|
||||
"tile_id",
|
||||
"xmin",
|
||||
"ymin",
|
||||
"xmax",
|
||||
"ymax",
|
||||
"global_min",
|
||||
"global_max",
|
||||
"out_res",
|
||||
"tile_key",
|
||||
"tile_min",
|
||||
"tile_max",
|
||||
]
|
||||
with open(manifest_path, "w", encoding="utf-8", newline="") as f:
|
||||
writer = csv.DictWriter(f, fieldnames=header)
|
||||
writer.writeheader()
|
||||
for row in rows:
|
||||
out = dict(row)
|
||||
out["global_min"] = f"{global_min:.6f}"
|
||||
out["global_max"] = f"{global_max:.6f}"
|
||||
writer.writerow(out)
|
||||
|
||||
|
||||
def _process_profile(
|
||||
name: str,
|
||||
profile: RiverErosionProfileConfig,
|
||||
tiles: List[Dict[str, object]],
|
||||
layer: ogr.Layer,
|
||||
to_tile: osr.CoordinateTransformation,
|
||||
to_layer: osr.CoordinateTransformation,
|
||||
tile_srs: osr.SpatialReference,
|
||||
output_dir: str,
|
||||
manifest_path: str,
|
||||
) -> int:
|
||||
ensure_dir(output_dir)
|
||||
rows_out: List[Dict[str, str]] = []
|
||||
global_min = math.inf
|
||||
global_max = -math.inf
|
||||
skipped = 0
|
||||
|
||||
for tile in tiles:
|
||||
tile_id = tile["tile_id"]
|
||||
png_path = tile["png_path"]
|
||||
if not os.path.exists(png_path):
|
||||
print(f"[river_erosion:{name}] Missing PNG for {tile_id}: {png_path}")
|
||||
skipped += 1
|
||||
continue
|
||||
|
||||
bounds = tile["bounds"]
|
||||
out_res = tile["out_res"]
|
||||
scale_min = tile["scale_min"]
|
||||
scale_max = tile["scale_max"]
|
||||
|
||||
ds = open_dataset(png_path, f"[river_erosion:{name}] Could not open {png_path}.")
|
||||
band = ds.GetRasterBand(1)
|
||||
raw = band.ReadAsArray()
|
||||
ds = None
|
||||
if raw is None:
|
||||
print(f"[river_erosion:{name}] Empty raster for {tile_id}, skipping.")
|
||||
skipped += 1
|
||||
continue
|
||||
|
||||
if raw.shape[0] != out_res or raw.shape[1] != out_res:
|
||||
print(
|
||||
f"[river_erosion:{name}] Warning: {tile_id} PNG is {raw.shape[1]}x{raw.shape[0]} "
|
||||
f"but manifest says {out_res}; using PNG resolution."
|
||||
)
|
||||
out_res = raw.shape[0]
|
||||
|
||||
depth, gt = _build_depth_raster(
|
||||
layer,
|
||||
bounds,
|
||||
out_res,
|
||||
profile,
|
||||
to_tile,
|
||||
to_layer,
|
||||
tile_srs,
|
||||
)
|
||||
|
||||
data = raw.astype(np.float32)
|
||||
valid = data > 0
|
||||
height = np.full(data.shape, np.nan, dtype=np.float32)
|
||||
if scale_max <= scale_min:
|
||||
scale_max = scale_min + 1e-3
|
||||
height[valid] = scale_min + (data[valid] / 65535.0) * (scale_max - scale_min)
|
||||
height[valid] = height[valid] - depth[valid]
|
||||
|
||||
if not np.any(valid):
|
||||
print(f"[river_erosion:{name}] No valid samples for {tile_id}, copying source.")
|
||||
out_u16 = raw.astype(np.uint16, copy=False)
|
||||
tile_min = scale_min
|
||||
tile_max = scale_max
|
||||
else:
|
||||
tile_min = float(np.nanmin(height))
|
||||
tile_max = float(np.nanmax(height))
|
||||
if tile_max <= tile_min:
|
||||
tile_max = tile_min + 1e-3
|
||||
scaled = (height[valid] - tile_min) / (tile_max - tile_min)
|
||||
out_u16 = np.zeros(data.shape, dtype=np.uint16)
|
||||
out_u16[valid] = np.clip(scaled * 65535.0, 0.0, 65535.0).astype(np.uint16)
|
||||
|
||||
global_min = min(global_min, tile_min)
|
||||
global_max = max(global_max, tile_max)
|
||||
|
||||
out_path = os.path.join(output_dir, f"{tile_id}.png")
|
||||
driver = gdal.GetDriverByName("PNG")
|
||||
out_ds = driver.Create(out_path, out_u16.shape[1], out_u16.shape[0], 1, gdal.GDT_UInt16, options=["WORLDFILE=YES"])
|
||||
out_ds.SetGeoTransform(gt)
|
||||
out_ds.SetProjection(tile_srs.ExportToWkt())
|
||||
out_ds.GetRasterBand(1).WriteArray(out_u16)
|
||||
out_ds = None
|
||||
|
||||
rows_out.append(
|
||||
{
|
||||
"tile_id": tile_id,
|
||||
"xmin": f"{bounds[0]:.6f}",
|
||||
"ymin": f"{bounds[1]:.6f}",
|
||||
"xmax": f"{bounds[2]:.6f}",
|
||||
"ymax": f"{bounds[3]:.6f}",
|
||||
"out_res": str(out_res),
|
||||
"tile_key": tile["tile_key"],
|
||||
"tile_min": f"{tile_min:.6f}",
|
||||
"tile_max": f"{tile_max:.6f}",
|
||||
}
|
||||
)
|
||||
|
||||
print(f"[river_erosion:{name}] Wrote {out_path}")
|
||||
|
||||
if not rows_out:
|
||||
print(f"[river_erosion:{name}] No tiles processed; nothing written.")
|
||||
return 1
|
||||
|
||||
_write_manifest(manifest_path, rows_out, global_min, global_max)
|
||||
print(f"[river_erosion:{name}] Manifest: {manifest_path}")
|
||||
print(f"[river_erosion:{name}] Summary: wrote {len(rows_out)} tiles; skipped {skipped}.")
|
||||
return 1 if skipped else 0
|
||||
|
||||
|
||||
def erode_rivers(cfg: Config) -> int:
|
||||
re_cfg = cfg.river_erosion
|
||||
if not os.path.exists(cfg.export.manifest_path):
|
||||
print(f"[river_erosion] Missing manifest: {cfg.export.manifest_path}. Run heightmap export first.")
|
||||
return 1
|
||||
if not os.path.exists(re_cfg.source_path):
|
||||
print(f"[river_erosion] Missing HydroRIVERS source: {re_cfg.source_path}")
|
||||
return 1
|
||||
if not re_cfg.enabled:
|
||||
print("[river_erosion] Warning: river_erosion.enabled=false; proceeding because --erode-rivers was requested.")
|
||||
|
||||
ds = ogr.Open(re_cfg.source_path)
|
||||
if ds is None:
|
||||
print(f"[river_erosion] Could not open {re_cfg.source_path}")
|
||||
return 1
|
||||
layer = ds.GetLayerByName(re_cfg.source_layer) if re_cfg.source_layer else None
|
||||
if layer is None:
|
||||
layer = ds.GetLayer(0)
|
||||
if layer is None:
|
||||
print(f"[river_erosion] No layers found in {re_cfg.source_path}")
|
||||
return 1
|
||||
|
||||
layer_srs = layer.GetSpatialRef()
|
||||
if layer_srs is None:
|
||||
layer_srs = osr.SpatialReference()
|
||||
layer_srs.ImportFromEPSG(4326)
|
||||
tile_srs = osr.SpatialReference()
|
||||
tile_srs.ImportFromEPSG(25832)
|
||||
layer_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
|
||||
tile_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
|
||||
to_tile = osr.CoordinateTransformation(layer_srs, tile_srs)
|
||||
to_layer = osr.CoordinateTransformation(tile_srs, layer_srs)
|
||||
|
||||
with open(cfg.export.manifest_path, newline="", encoding="utf-8") as f:
|
||||
reader = csv.DictReader(f)
|
||||
if not reader.fieldnames:
|
||||
print(f"[river_erosion] Empty manifest: {cfg.export.manifest_path}")
|
||||
return 1
|
||||
required = {"tile_id", "xmin", "ymin", "xmax", "ymax", "global_min", "global_max", "out_res"}
|
||||
missing = required - set(reader.fieldnames)
|
||||
if missing:
|
||||
print(f"[river_erosion] Manifest missing columns: {', '.join(sorted(missing))}")
|
||||
return 1
|
||||
rows = list(reader)
|
||||
|
||||
if not rows:
|
||||
print(f"[river_erosion] No tiles in manifest: {cfg.export.manifest_path}")
|
||||
return 1
|
||||
|
||||
has_tile_minmax = "tile_min" in (reader.fieldnames or []) and "tile_max" in (reader.fieldnames or [])
|
||||
need_tile_minmax = not has_tile_minmax or any(
|
||||
not row.get("tile_min") or not row.get("tile_max") for row in rows
|
||||
)
|
||||
|
||||
vrt_ds = None
|
||||
if need_tile_minmax:
|
||||
if not os.path.isdir(cfg.raw.dgm1_dir):
|
||||
print(f"[river_erosion] Missing DGM1 directory: {cfg.raw.dgm1_dir}")
|
||||
return 1
|
||||
tif_paths = sorted(glob.glob(os.path.join(cfg.raw.dgm1_dir, "*.tif")))
|
||||
if not tif_paths:
|
||||
print(f"[river_erosion] No DGM1 TIFFs found in {cfg.raw.dgm1_dir}")
|
||||
return 1
|
||||
build_vrt(cfg.work.heightmap_vrt, tif_paths, force=False)
|
||||
vrt_ds = open_dataset(cfg.work.heightmap_vrt, f"[river_erosion] Could not open {cfg.work.heightmap_vrt}.")
|
||||
|
||||
tiles: List[Dict[str, object]] = []
|
||||
for row in rows:
|
||||
tile_id = row["tile_id"].strip()
|
||||
bounds = (
|
||||
float(row["xmin"]),
|
||||
float(row["ymin"]),
|
||||
float(row["xmax"]),
|
||||
float(row["ymax"]),
|
||||
)
|
||||
out_res = int(row["out_res"])
|
||||
tile_key = row.get("tile_key", "").strip()
|
||||
scale_min = None
|
||||
scale_max = None
|
||||
if row.get("tile_min") and row.get("tile_max"):
|
||||
scale_min = float(row["tile_min"])
|
||||
scale_max = float(row["tile_max"])
|
||||
if scale_min is None or scale_max is None:
|
||||
if vrt_ds is None:
|
||||
print(f"[river_erosion] Missing tile_min/tile_max for {tile_id} and no VRT available.")
|
||||
return 1
|
||||
scale_min, scale_max = _compute_tile_minmax_from_vrt(
|
||||
vrt_ds,
|
||||
bounds,
|
||||
out_res,
|
||||
cfg.heightmap.resample,
|
||||
)
|
||||
if scale_min is None or scale_max is None:
|
||||
print(f"[river_erosion] Could not compute tile min/max for {tile_id}, using global range.")
|
||||
scale_min = float(row["global_min"])
|
||||
scale_max = float(row["global_max"])
|
||||
|
||||
tiles.append(
|
||||
{
|
||||
"tile_id": tile_id,
|
||||
"bounds": bounds,
|
||||
"out_res": out_res,
|
||||
"tile_key": tile_key,
|
||||
"scale_min": float(scale_min),
|
||||
"scale_max": float(scale_max),
|
||||
"png_path": os.path.join(cfg.export.heightmap_dir, f"{tile_id}.png"),
|
||||
}
|
||||
)
|
||||
|
||||
output_vr = os.path.join(re_cfg.output_dir, "vr")
|
||||
output_server = os.path.join(re_cfg.output_dir, "server")
|
||||
|
||||
vr_exit = _process_profile(
|
||||
"vr",
|
||||
re_cfg.vr,
|
||||
tiles,
|
||||
layer,
|
||||
to_tile,
|
||||
to_layer,
|
||||
tile_srs,
|
||||
output_vr,
|
||||
re_cfg.manifest_vr,
|
||||
)
|
||||
server_exit = _process_profile(
|
||||
"server",
|
||||
re_cfg.server,
|
||||
tiles,
|
||||
layer,
|
||||
to_tile,
|
||||
to_layer,
|
||||
tile_srs,
|
||||
output_server,
|
||||
re_cfg.manifest_server,
|
||||
)
|
||||
return max(vr_exit, server_exit)
|
||||
@@ -21,6 +21,7 @@ from geodata_pipeline.heightmaps import export_heightmaps
|
||||
from geodata_pipeline.heightmaps_enhanced import export_heightmaps_enhanced
|
||||
from geodata_pipeline.lpolpg_split import split_lpolpg
|
||||
from geodata_pipeline.orthophotos import export_orthophotos
|
||||
from geodata_pipeline.river_erosion import erode_rivers
|
||||
from geodata_pipeline.setup_helpers import ensure_directories, materialize_archives
|
||||
from geodata_pipeline.street_furniture import export_street_furniture
|
||||
from geodata_pipeline.trees import export_trees
|
||||
@@ -123,6 +124,11 @@ def parse_args(argv: Iterable[str] | None = None) -> argparse.Namespace:
|
||||
action="store_true",
|
||||
help="Delete lpolpg source files after a successful split.",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--erode-rivers",
|
||||
action="store_true",
|
||||
help="Erode heightmap tiles using HydroRIVERS after heightmap export (requires tile_index.csv).",
|
||||
)
|
||||
return parser.parse_args(argv)
|
||||
|
||||
|
||||
@@ -164,8 +170,11 @@ def main(argv: Iterable[str] | None = None) -> int:
|
||||
args = parse_args(argv)
|
||||
cfg = load_config(args)
|
||||
target_export = None
|
||||
if args.export is not None or not (args.download or args.split_lpolpg):
|
||||
target_export = args.export or "all"
|
||||
action_flags = args.download or args.split_lpolpg or args.erode_rivers
|
||||
if args.export is not None:
|
||||
target_export = args.export
|
||||
elif not action_flags:
|
||||
target_export = "all"
|
||||
|
||||
if args.setup:
|
||||
ensure_directories(cfg)
|
||||
@@ -214,34 +223,35 @@ def main(argv: Iterable[str] | None = None) -> int:
|
||||
if split_exit != 0:
|
||||
return split_exit
|
||||
|
||||
if target_export is None:
|
||||
return 0
|
||||
|
||||
exit_codes = []
|
||||
|
||||
# Standard exports
|
||||
if target_export in ("heightmap", "all"):
|
||||
exit_codes.append(export_heightmaps(cfg, force_vrt=args.force_vrt))
|
||||
if target_export in ("textures", "all"):
|
||||
exit_codes.append(export_orthophotos(cfg, force_vrt=args.force_vrt))
|
||||
if target_export in ("buildings", "all"):
|
||||
exit_codes.append(export_buildings(cfg))
|
||||
if target_export in ("trees", "all"):
|
||||
exit_codes.append(export_trees(cfg, force_vrt=args.force_vrt))
|
||||
if target_export is not None:
|
||||
# Standard exports
|
||||
if target_export in ("heightmap", "all"):
|
||||
exit_codes.append(export_heightmaps(cfg, force_vrt=args.force_vrt))
|
||||
if target_export in ("textures", "all"):
|
||||
exit_codes.append(export_orthophotos(cfg, force_vrt=args.force_vrt))
|
||||
if target_export in ("buildings", "all"):
|
||||
exit_codes.append(export_buildings(cfg))
|
||||
if target_export in ("trees", "all"):
|
||||
exit_codes.append(export_trees(cfg, force_vrt=args.force_vrt))
|
||||
|
||||
# Enhanced exports (use point cloud data)
|
||||
# Order matters: heightmap-enhanced creates tile_index.csv needed by others
|
||||
# street-furniture must run before trees-enhanced (for exclusion mask)
|
||||
if target_export in ("heightmap-enhanced", "all-enhanced"):
|
||||
exit_codes.append(export_heightmaps_enhanced(cfg))
|
||||
if target_export in ("all-enhanced",):
|
||||
exit_codes.append(export_orthophotos(cfg, force_vrt=args.force_vrt))
|
||||
if target_export in ("street-furniture", "all-enhanced"):
|
||||
exit_codes.append(export_street_furniture(cfg))
|
||||
if target_export in ("buildings-enhanced", "all-enhanced"):
|
||||
exit_codes.append(export_buildings_enhanced(cfg))
|
||||
if target_export in ("trees-enhanced", "all-enhanced"):
|
||||
exit_codes.append(export_trees_enhanced(cfg))
|
||||
# Enhanced exports (use point cloud data)
|
||||
# Order matters: heightmap-enhanced creates tile_index.csv needed by others
|
||||
# street-furniture must run before trees-enhanced (for exclusion mask)
|
||||
if target_export in ("heightmap-enhanced", "all-enhanced"):
|
||||
exit_codes.append(export_heightmaps_enhanced(cfg))
|
||||
if target_export in ("all-enhanced",):
|
||||
exit_codes.append(export_orthophotos(cfg, force_vrt=args.force_vrt))
|
||||
if target_export in ("street-furniture", "all-enhanced"):
|
||||
exit_codes.append(export_street_furniture(cfg))
|
||||
if target_export in ("buildings-enhanced", "all-enhanced"):
|
||||
exit_codes.append(export_buildings_enhanced(cfg))
|
||||
if target_export in ("trees-enhanced", "all-enhanced"):
|
||||
exit_codes.append(export_trees_enhanced(cfg))
|
||||
|
||||
if args.erode_rivers:
|
||||
exit_codes.append(erode_rivers(cfg))
|
||||
|
||||
return max(exit_codes) if exit_codes else 0
|
||||
|
||||
|
||||
Reference in New Issue
Block a user