added covarage
This commit is contained in:
+222
-8
@@ -1,6 +1,28 @@
|
||||
from __future__ import annotations
|
||||
|
||||
import subprocess
|
||||
import tempfile
|
||||
from dataclasses import dataclass
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import rasterio
|
||||
from rasterio.crs import CRS
|
||||
from rasterio.features import shapes
|
||||
from rasterio.warp import (
|
||||
Resampling,
|
||||
calculate_default_transform,
|
||||
reproject,
|
||||
transform,
|
||||
)
|
||||
from shapely.geometry import mapping, shape
|
||||
from shapely.ops import unary_union
|
||||
|
||||
from app.core.dem import _dem_files, _sample_dataset
|
||||
from app.core.geo import GeoPoint
|
||||
from app.models.viewshed import ViewshedRequest
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class RasterResult:
|
||||
@@ -8,13 +30,205 @@ class RasterResult:
|
||||
metadata: dict[str, Any]
|
||||
|
||||
|
||||
def compute_viewshed(
|
||||
observer: object,
|
||||
def _utm_epsg(lon: float, lat: float) -> int:
|
||||
zone = int((lon + 180) // 6) + 1
|
||||
return 32600 + zone if lat >= 0 else 32700 + zone
|
||||
|
||||
|
||||
def _find_dem_tile(lat: float, lon: float, dem_path: Path) -> Path | None:
|
||||
for path in _dem_files(dem_path):
|
||||
with rasterio.open(path) as dataset:
|
||||
if _sample_dataset(dataset, lat, lon) is not None:
|
||||
return path
|
||||
return None
|
||||
|
||||
|
||||
def _prepare_metric_dem(
|
||||
source_path: Path,
|
||||
observer: GeoPoint,
|
||||
radius_m: float,
|
||||
target_h: float,
|
||||
surface: str,
|
||||
k: float,
|
||||
output_path: Path,
|
||||
) -> tuple[Path, float, float]:
|
||||
dst_crs = CRS.from_epsg(_utm_epsg(observer.lon, observer.lat))
|
||||
with rasterio.open(source_path) as src:
|
||||
observer_x, observer_y = transform(
|
||||
CRS.from_epsg(4326),
|
||||
dst_crs,
|
||||
[observer.lon],
|
||||
[observer.lat],
|
||||
)
|
||||
transform_affine, width, height = calculate_default_transform(
|
||||
src.crs,
|
||||
dst_crs,
|
||||
src.width,
|
||||
src.height,
|
||||
*src.bounds,
|
||||
)
|
||||
data = np.empty((height, width), dtype=np.float32)
|
||||
reproject(
|
||||
source=rasterio.band(src, 1),
|
||||
destination=data,
|
||||
src_transform=src.transform,
|
||||
src_crs=src.crs,
|
||||
dst_transform=transform_affine,
|
||||
dst_crs=dst_crs,
|
||||
resampling=Resampling.bilinear,
|
||||
)
|
||||
profile = src.profile.copy()
|
||||
profile.update(
|
||||
{
|
||||
"crs": dst_crs,
|
||||
"transform": transform_affine,
|
||||
"width": width,
|
||||
"height": height,
|
||||
"dtype": "float32",
|
||||
"count": 1,
|
||||
"nodata": -9999,
|
||||
}
|
||||
)
|
||||
with rasterio.open(output_path, "w", **profile) as dst:
|
||||
dst.write(data, 1)
|
||||
return output_path, observer_x[0], observer_y[0]
|
||||
|
||||
|
||||
def _run_gdal_viewshed(
|
||||
dem_path: Path,
|
||||
observer_x: float,
|
||||
observer_y: float,
|
||||
observer_height: float,
|
||||
target_height: float,
|
||||
output_path: Path,
|
||||
gdal_bin: str,
|
||||
) -> None:
|
||||
command = [
|
||||
gdal_bin,
|
||||
"-ox",
|
||||
str(observer_x),
|
||||
"-oy",
|
||||
str(observer_y),
|
||||
"-oz",
|
||||
str(observer_height),
|
||||
"-tz",
|
||||
str(target_height),
|
||||
"-vv",
|
||||
"255",
|
||||
"-iv",
|
||||
"0",
|
||||
"-ov",
|
||||
"0",
|
||||
str(dem_path),
|
||||
str(output_path),
|
||||
]
|
||||
subprocess.run(command, check=True, capture_output=True, text=True)
|
||||
|
||||
|
||||
def _viewshed_geojson(viewshed_path: Path, observer: GeoPoint) -> dict[str, Any]:
|
||||
features: list[dict[str, Any]] = []
|
||||
with rasterio.open(viewshed_path) as dataset:
|
||||
data = dataset.read(1, masked=True)
|
||||
mask = (data == 255).astype(np.uint8)
|
||||
for geom, value in shapes(mask, mask=mask.astype(bool), transform=dataset.transform):
|
||||
if int(value) != 1:
|
||||
continue
|
||||
polygon = shape(geom)
|
||||
if polygon.is_empty:
|
||||
continue
|
||||
features.append(
|
||||
{
|
||||
"type": "Feature",
|
||||
"properties": {"visible": True},
|
||||
"geometry": mapping(polygon),
|
||||
}
|
||||
)
|
||||
if not features:
|
||||
return {
|
||||
"type": "FeatureCollection",
|
||||
"features": [],
|
||||
"properties": {"observer": {"lat": observer.lat, "lon": observer.lon}},
|
||||
}
|
||||
merged = unary_union([shape(feature["geometry"]) for feature in features])
|
||||
return {
|
||||
"type": "FeatureCollection",
|
||||
"features": [
|
||||
{
|
||||
"type": "Feature",
|
||||
"properties": {"visible": True},
|
||||
"geometry": mapping(merged),
|
||||
}
|
||||
],
|
||||
"properties": {"observer": {"lat": observer.lat, "lon": observer.lon}},
|
||||
}
|
||||
|
||||
|
||||
def compute_viewshed(
|
||||
request: ViewshedRequest,
|
||||
*,
|
||||
dem_path: str | Path,
|
||||
output_dir: str | Path,
|
||||
gdal_viewshed_bin: str = "gdal_viewshed",
|
||||
) -> RasterResult:
|
||||
raise NotImplementedError(
|
||||
"Viewshed requires WhiteboxTools or gdal_viewshed integration in a later stage"
|
||||
)
|
||||
dem_root = Path(dem_path)
|
||||
output_root = Path(output_dir)
|
||||
output_root.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
observer = GeoPoint(lat=request.observer.lat, lon=request.observer.lon)
|
||||
source_tile = _find_dem_tile(observer.lat, observer.lon, dem_root)
|
||||
if source_tile is None:
|
||||
raise FileNotFoundError(
|
||||
f"No DEM tile covers observer at lat={observer.lat}, lon={observer.lon}"
|
||||
)
|
||||
|
||||
with tempfile.TemporaryDirectory(prefix="viewshed-") as tmp_dir:
|
||||
tmp = Path(tmp_dir)
|
||||
metric_dem = tmp / "metric_dem.tif"
|
||||
viewshed_raster = tmp / "viewshed.tif"
|
||||
metric_path, observer_x, observer_y = _prepare_metric_dem(
|
||||
source_tile,
|
||||
observer,
|
||||
request.radius_m,
|
||||
metric_dem,
|
||||
)
|
||||
with rasterio.open(source_tile) as source_dataset:
|
||||
observer_elevation = _sample_dataset(
|
||||
source_dataset,
|
||||
observer.lat,
|
||||
observer.lon,
|
||||
)
|
||||
observer_abs = float(observer_elevation or 0.0) + request.observer.height_agl
|
||||
_run_gdal_viewshed(
|
||||
metric_path,
|
||||
observer_x,
|
||||
observer_y,
|
||||
observer_abs,
|
||||
request.target_height_agl,
|
||||
viewshed_raster,
|
||||
gdal_viewshed_bin,
|
||||
)
|
||||
|
||||
if request.format == "geojson":
|
||||
result = _viewshed_geojson(viewshed_raster, observer)
|
||||
output_file = output_root / f"viewshed_{observer.lat:.5f}_{observer.lon:.5f}.json"
|
||||
output_file.write_text(
|
||||
__import__("json").dumps(result),
|
||||
encoding="utf-8",
|
||||
)
|
||||
return RasterResult(
|
||||
uri=str(output_file),
|
||||
metadata={
|
||||
"format": "geojson",
|
||||
"surface": request.surface,
|
||||
"radius_m": request.radius_m,
|
||||
"feature_count": len(result.get("features", [])),
|
||||
},
|
||||
)
|
||||
|
||||
output_file = output_root / f"viewshed_{observer.lat:.5f}_{observer.lon:.5f}.tif"
|
||||
output_file.write_bytes(viewshed_raster.read_bytes())
|
||||
return RasterResult(
|
||||
uri=str(output_file),
|
||||
metadata={
|
||||
"format": request.format,
|
||||
"surface": request.surface,
|
||||
"radius_m": request.radius_m,
|
||||
},
|
||||
)
|
||||
|
||||
Reference in New Issue
Block a user