def build_cost_surface(
*,
wv3_path: str | Path,
chm_path: str | Path,
slic_path: str | Path | None = None,
weights: tuple[float, float, float, float] = (0.4, 0.3, 0.3, 0),
nodata_value: float = -9999.0,
):
"""Build cost surface."""
w = CostSurfaceWeights(*weights)
with rasterio.open(wv3_path) as src:
C, B, G, Y, R, RE, N1, N2 = src.read(masked=True).astype(np.float32)
profile = src.profile.copy()
profile["bounds"] = rasterio.transform.array_bounds(
profile["height"], profile["width"], profile["transform"]
)
chm, _ = _read_band(chm_path)
grad = _chm_gradient(chm)
gap = _normalise(1 - _ndvi(R, N1))
tex = _texture_entropy(C)
if slic_path:
if str(slic_path).lower().endswith(".gpkg"):
slic_lab = _rasterise_slic_gpkg(slic_path, profile)
else: # assume raster label image
slic_lab, _ = _read_band(slic_path)
edge = _slic_edge(slic_lab)
w_slic = w.slic
else:
edge = 0.0
s = w.grad + w.gap + w.tex
w.grad, w.gap, w.tex = (w.grad / s, w.gap / s, w.tex / s)
w_slic = 0.0
warnings.warn("No SLIC provided – cost built from 3 terms only.")
cost = (
w.grad * grad +
w.gap * gap +
w.tex * tex +
w_slic * edge
)
cost = np.clip(cost, 0, 1).astype(np.float32)
cost[np.isnan(cost)] = nodata_value
profile.update(count=1, dtype="float32",
compress="deflate",
nodata=nodata_value)
return cost, profile