Grigo/RadioPropagationApi
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

added conopy

Browse Source
This commit is contained in:
grigo
2026-06-24 09:01:11 +03:00
parent d9c262cee2
commit 7ec7a51648
18 changed files with 652 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
api/app/core/viewshed.py
+98 -31
View File
@@ -8,17 +8,20 @@ from typing import Any
import numpy as np
import rasterio
from pyproj import Geod, Transformer
from rasterio.crs import CRS
from rasterio.features import shapes
from rasterio.merge import merge
from rasterio.warp import (
Resampling,
calculate_default_transform,
reproject,
transform,
transform_bounds,
)
from pyproj import Transformer
from shapely.geometry import mapping, shape
from shapely.ops import transform as shapely_transform, unary_union
from shapely.ops import transform as shapely_transform
from shapely.ops import unary_union
from app.core.dem import _dem_files, _sample_dataset
from app.core.geo import GeoPoint
@@ -31,11 +34,31 @@ class RasterResult:
metadata: dict[str, Any]
_GEOD = Geod(ellps="WGS84")
def _utm_epsg(lon: float, lat: float) -> int:
zone = int((lon + 180) // 6) + 1
return 32600 + zone if lat >= 0 else 32700 + zone
def _intersects(
a: tuple[float, float, float, float],
b: tuple[float, float, float, float],
) -> bool:
return not (a[2] < b[0] or a[0] > b[2] or a[3] < b[1] or a[1] > b[3])
def _radius_bbox(observer: GeoPoint, radius_m: float) -> tuple[float, float, float, float]:
coords = [
_GEOD.fwd(observer.lon, observer.lat, azimuth, radius_m)[:2]
for azimuth in range(0, 360, 45)
]
lons = [observer.lon, *(coord[0] for coord in coords)]
lats = [observer.lat, *(coord[1] for coord in coords)]
return min(lons), min(lats), max(lons), max(lats)
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:
@@ -44,14 +67,48 @@ def _find_dem_tile(lat: float, lon: float, dem_path: Path) -> Path | None:
return None
def _find_dem_tiles(
dem_path: Path,
bbox_wgs84: tuple[float, float, float, float],
) -> list[Path]:
matches: list[Path] = []
for path in _dem_files(dem_path):
with rasterio.open(path) as dataset:
bounds = dataset.bounds
tile_bbox = (
bounds.left,
bounds.bottom,
bounds.right,
bounds.top,
)
if dataset.crs is not None and dataset.crs != CRS.from_epsg(4326):
tile_bbox = transform_bounds(dataset.crs, CRS.from_epsg(4326), *tile_bbox)
if _intersects(tile_bbox, bbox_wgs84):
matches.append(path)
return matches
def _prepare_metric_dem(
source_path: Path,
source_paths: list[Path],
observer: GeoPoint,
radius_m: float,
output_path: Path,
) -> tuple[Path, float, float]:
) -> tuple[Path, float, float, int]:
if not source_paths:
raise FileNotFoundError("No DEM tiles found for viewshed radius")
dst_crs = CRS.from_epsg(_utm_epsg(observer.lon, observer.lat))
with rasterio.open(source_path) as src:
bbox_wgs84 = _radius_bbox(observer, radius_m)
datasets = [rasterio.open(path) for path in source_paths]
try:
src_crs = datasets[0].crs or CRS.from_epsg(4326)
merge_bounds = bbox_wgs84
if src_crs != CRS.from_epsg(4326):
merge_bounds = transform_bounds(CRS.from_epsg(4326), src_crs, *bbox_wgs84)
mosaic, src_transform = merge(datasets, bounds=merge_bounds)
source_data = mosaic[0].astype(np.float32)
height, width = source_data.shape
src_bounds = rasterio.transform.array_bounds(height, width, src_transform)
observer_x, observer_y = transform(
CRS.from_epsg(4326),
dst_crs,
@@ -59,37 +116,38 @@ def _prepare_metric_dem(
[observer.lat],
)
transform_affine, width, height = calculate_default_transform(
src.crs,
src_crs,
dst_crs,
src.width,
src.height,
*src.bounds,
width,
height,
*src_bounds,
)
data = np.empty((height, width), dtype=np.float32)
reproject(
source=rasterio.band(src, 1),
source=source_data,
destination=data,
src_transform=src.transform,
src_crs=src.crs,
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,
}
)
profile = {
"driver": "GTiff",
"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]
finally:
for dataset in datasets:
dataset.close()
return output_path, observer_x[0], observer_y[0], len(source_paths)
def _run_gdal_viewshed(
@@ -187,19 +245,24 @@ def compute_viewshed(
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:
calc_radius_m = request.calc_radius_m or request.radius_m
bbox_wgs84 = _radius_bbox(observer, calc_radius_m)
source_tiles = _find_dem_tiles(dem_root, bbox_wgs84)
if source_tile is not None and source_tile not in source_tiles:
source_tiles.append(source_tile)
if not source_tiles:
raise FileNotFoundError(
f"No DEM tile covers observer at lat={observer.lat}, lon={observer.lon}"
f"No DEM tile covers viewshed radius 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,
metric_path, observer_x, observer_y, dem_tiles_count = _prepare_metric_dem(
source_tiles,
observer,
request.radius_m,
calc_radius_m,
metric_dem,
)
_run_gdal_viewshed(
@@ -208,7 +271,7 @@ def compute_viewshed(
observer_y,
request.observer.height_agl,
request.target_height_agl,
request.radius_m,
calc_radius_m,
viewshed_raster,
gdal_viewshed_bin,
)
@@ -226,6 +289,8 @@ def compute_viewshed(
"format": "geojson",
"surface": request.surface,
"radius_m": request.radius_m,
"calc_radius_m": calc_radius_m,
"dem_tiles_count": dem_tiles_count,
"feature_count": len(result.get("features", [])),
},
)
@@ -238,5 +303,7 @@ def compute_viewshed(
"format": request.format,
"surface": request.surface,
"radius_m": request.radius_m,
"calc_radius_m": calc_radius_m,
"dem_tiles_count": dem_tiles_count,
},
)