353 lines
12 KiB
Python
353 lines
12 KiB
Python
from __future__ import annotations
|
|
|
|
import numpy as np
|
|
from pyproj import Transformer
|
|
from shapely.geometry import LineString, Point, mapping
|
|
from shapely.ops import transform as shapely_transform
|
|
from shapely.ops import unary_union
|
|
from sqlalchemy.exc import SQLAlchemyError
|
|
from sqlalchemy.orm import Session
|
|
|
|
from app.config import get_settings
|
|
from app.core import dem
|
|
from app.core.diffraction import deygout
|
|
from app.core.fresnel import LosSample, los_analysis
|
|
from app.core.geo import GeoPoint, PathPoint, linestring_geojson, sample_path
|
|
from app.core.raster_sampling import raster_files
|
|
from app.core.surface import SurfaceProfile, build_surface_profile
|
|
from app.models.common import DataSources
|
|
from app.models.terrain import (
|
|
ElevationResponse,
|
|
FresnelProfileSample,
|
|
FresnelSliceRequest,
|
|
FresnelSliceResponse,
|
|
LosRequest,
|
|
LosResponse,
|
|
Obstruction,
|
|
ProfileSample,
|
|
TerrainProfileRequest,
|
|
TerrainProfileResponse,
|
|
)
|
|
from app.services import buildings as buildings_service
|
|
from app.services.vegetation import vegetation_loss_along
|
|
|
|
|
|
def elevation_at(lat: float, lon: float, surface: str) -> ElevationResponse:
|
|
settings = get_settings()
|
|
elevation_m = dem.elevation_at(lat, lon, surface=surface, dem_path=settings.dem_path)
|
|
return ElevationResponse(
|
|
lat=lat,
|
|
lon=lon,
|
|
elevation_m=elevation_m,
|
|
surface=surface,
|
|
data_sources=DataSources(dem=True),
|
|
)
|
|
|
|
|
|
def terrain_profile(
|
|
request: TerrainProfileRequest,
|
|
db: Session | None = None,
|
|
) -> TerrainProfileResponse:
|
|
points = sample_path(
|
|
GeoPoint(lat=request.start.lat, lon=request.start.lon),
|
|
GeoPoint(lat=request.end.lat, lon=request.end.lon),
|
|
request.samples,
|
|
)
|
|
profile = surface_profile_from_points(
|
|
points,
|
|
include_buildings=request.include_buildings,
|
|
include_canopy=request.include_canopy,
|
|
db=db,
|
|
)
|
|
return TerrainProfileResponse(
|
|
distance_m=profile.distance_m,
|
|
samples=[ProfileSample(**sample.__dict__) for sample in profile.samples],
|
|
path_geojson=linestring_geojson(points),
|
|
data_sources=_data_sources(
|
|
dem_used=bool(raster_files(get_settings().dem_path)),
|
|
include_buildings=request.include_buildings,
|
|
db=db,
|
|
),
|
|
)
|
|
|
|
|
|
def los(
|
|
request: LosRequest,
|
|
db: Session | None = None,
|
|
) -> LosResponse:
|
|
profile_request = TerrainProfileRequest(
|
|
start=request.tx,
|
|
end=request.rx,
|
|
samples=request.samples,
|
|
include_buildings=request.include_buildings,
|
|
include_canopy=request.include_canopy,
|
|
k_factor=request.k_factor,
|
|
)
|
|
points = sample_path(
|
|
GeoPoint(lat=request.tx.lat, lon=request.tx.lon),
|
|
GeoPoint(lat=request.rx.lat, lon=request.rx.lon),
|
|
request.samples,
|
|
)
|
|
surface_profile = surface_profile_from_points(
|
|
points,
|
|
include_buildings=request.include_buildings,
|
|
include_canopy=request.include_canopy,
|
|
db=db,
|
|
)
|
|
settings = get_settings()
|
|
result = los_analysis(
|
|
surface_profile,
|
|
request.tx.height_agl,
|
|
request.rx.height_agl,
|
|
request.frequency_mhz * 1_000_000,
|
|
request.fresnel_clearance,
|
|
request.k_factor,
|
|
)
|
|
diffraction_loss = deygout(
|
|
surface_profile,
|
|
request.tx.height_agl,
|
|
request.rx.height_agl,
|
|
request.frequency_mhz * 1_000_000,
|
|
)
|
|
vegetation_loss = vegetation_loss_along(
|
|
points,
|
|
freq_hz=request.frequency_mhz * 1_000_000,
|
|
include_vegetation=request.include_vegetation,
|
|
)
|
|
|
|
def convert(sample: LosSample) -> Obstruction:
|
|
return Obstruction(
|
|
distance_m=sample.distance_m,
|
|
clearance_m=sample.clearance_m,
|
|
fresnel_radius_m=sample.fresnel_radius_m,
|
|
required_clearance_m=sample.required_clearance_m,
|
|
type=sample.obstruction_type,
|
|
)
|
|
|
|
def convert_fresnel_sample(i: int, sample: LosSample) -> FresnelProfileSample:
|
|
return FresnelProfileSample(
|
|
i=i,
|
|
lat=sample.lat,
|
|
lon=sample.lon,
|
|
distance_m=sample.distance_m,
|
|
surface_m=sample.surface_m,
|
|
path_height_m=sample.path_height_m,
|
|
obstacle_height_m=sample.obstacle_height_m,
|
|
clearance_m=sample.clearance_m,
|
|
fresnel_radius_m=sample.fresnel_radius_m,
|
|
required_clearance_m=sample.required_clearance_m,
|
|
fresnel_upper_60pct_m=sample.path_height_m + sample.required_clearance_m,
|
|
fresnel_upper_100pct_m=sample.path_height_m + sample.fresnel_radius_m,
|
|
fresnel_lower_60pct_m=sample.path_height_m - sample.required_clearance_m,
|
|
fresnel_lower_100pct_m=sample.path_height_m - sample.fresnel_radius_m,
|
|
type=sample.obstruction_type,
|
|
)
|
|
|
|
geometric_obstructions = [
|
|
sample for sample in result.obstructions if sample.clearance_m < 0
|
|
]
|
|
building_obstructions = [
|
|
sample
|
|
for sample in geometric_obstructions
|
|
if sample.obstruction_type == "building"
|
|
]
|
|
|
|
return LosResponse(
|
|
los_clear=result.los_clear,
|
|
geometric_los=result.geometric_los,
|
|
first_fresnel_clearance_pct=result.first_fresnel_clearance_pct,
|
|
worst_obstruction=convert(result.worst_sample) if result.worst_sample else None,
|
|
obstructions=[convert(sample) for sample in result.obstructions],
|
|
geometric_obstructions=[convert(sample) for sample in geometric_obstructions],
|
|
fresnel_profile=[
|
|
convert_fresnel_sample(i, sample) for i, sample in enumerate(result.samples)
|
|
],
|
|
fresnel_violations_count=len(result.obstructions),
|
|
geometric_obstructions_count=len(geometric_obstructions),
|
|
building_obstructions_count=len(building_obstructions),
|
|
diffraction_loss_db=diffraction_loss,
|
|
vegetation_loss_db=vegetation_loss,
|
|
data_sources=_data_sources(
|
|
dem_used=bool(raster_files(settings.dem_path)),
|
|
include_buildings=request.include_buildings,
|
|
db=db,
|
|
landcover_used=request.include_vegetation
|
|
and bool(raster_files(settings.landcover_path)),
|
|
canopy_used=request.include_canopy
|
|
and bool(raster_files(settings.canopy_path)),
|
|
),
|
|
profile_ref=f"inline:{profile_request.samples}",
|
|
)
|
|
|
|
|
|
def fresnel_slice(
|
|
request: FresnelSliceRequest,
|
|
db: Session | None = None,
|
|
) -> FresnelSliceResponse:
|
|
points = sample_path(
|
|
GeoPoint(lat=request.tx.lat, lon=request.tx.lon),
|
|
GeoPoint(lat=request.rx.lat, lon=request.rx.lon),
|
|
request.samples,
|
|
)
|
|
surface_profile = surface_profile_from_points(
|
|
points,
|
|
include_buildings=request.include_buildings,
|
|
include_canopy=request.include_canopy,
|
|
db=db,
|
|
)
|
|
settings = get_settings()
|
|
result = los_analysis(
|
|
surface_profile,
|
|
request.tx.height_agl,
|
|
request.rx.height_agl,
|
|
request.frequency_mhz * 1_000_000,
|
|
request.fresnel_fraction,
|
|
request.k_factor,
|
|
)
|
|
fresnel_geometry = _fresnel_slice_geojson(result.samples, request.fresnel_fraction)
|
|
buildings_geojson = {"type": "FeatureCollection", "features": []}
|
|
if request.include_buildings and db is not None:
|
|
try:
|
|
buildings_geojson = buildings_service.query_building_intersections(
|
|
fresnel_geometry["geometry"],
|
|
db,
|
|
request.building_limit,
|
|
).model_dump()
|
|
except (NotImplementedError, SQLAlchemyError):
|
|
buildings_geojson = {"type": "FeatureCollection", "features": []}
|
|
|
|
return FresnelSliceResponse(
|
|
path_geojson=linestring_geojson(points),
|
|
fresnel_geojson=fresnel_geometry,
|
|
buildings_geojson=buildings_geojson,
|
|
profile=[
|
|
_convert_fresnel_sample(index, sample)
|
|
for index, sample in enumerate(result.samples)
|
|
],
|
|
data_sources=_data_sources(
|
|
dem_used=bool(raster_files(settings.dem_path)),
|
|
include_buildings=request.include_buildings,
|
|
db=db,
|
|
canopy_used=request.include_canopy
|
|
and bool(raster_files(settings.canopy_path)),
|
|
),
|
|
)
|
|
|
|
|
|
def _building_heights(
|
|
points: list[PathPoint],
|
|
include_buildings: bool,
|
|
db: Session | None,
|
|
) -> list[float] | None:
|
|
if not include_buildings or db is None:
|
|
return None
|
|
try:
|
|
return buildings_service.building_heights_along(points, db)
|
|
except SQLAlchemyError:
|
|
return None
|
|
|
|
|
|
def surface_profile_from_points(
|
|
points: list[PathPoint],
|
|
include_buildings: bool,
|
|
include_canopy: bool,
|
|
db: Session | None = None,
|
|
) -> SurfaceProfile:
|
|
try:
|
|
ground_elevations = dem.elevations_along(points, dem_path=get_settings().dem_path).tolist()
|
|
except dem.DemNotConfiguredError:
|
|
ground_elevations = np.zeros(len(points), dtype=float).tolist()
|
|
|
|
building_heights = _building_heights(points, include_buildings, db)
|
|
return build_surface_profile(
|
|
points,
|
|
ground_elevations=ground_elevations,
|
|
building_heights=building_heights,
|
|
include_buildings=include_buildings,
|
|
include_canopy=include_canopy,
|
|
)
|
|
|
|
|
|
def _utm_epsg(lon: float, lat: float) -> int:
|
|
zone = int((lon + 180) // 6) + 1
|
|
return 32600 + zone if lat >= 0 else 32700 + zone
|
|
|
|
|
|
def _convert_fresnel_sample(i: int, sample: LosSample) -> FresnelProfileSample:
|
|
return FresnelProfileSample(
|
|
i=i,
|
|
lat=sample.lat,
|
|
lon=sample.lon,
|
|
distance_m=sample.distance_m,
|
|
surface_m=sample.surface_m,
|
|
path_height_m=sample.path_height_m,
|
|
obstacle_height_m=sample.obstacle_height_m,
|
|
clearance_m=sample.clearance_m,
|
|
fresnel_radius_m=sample.fresnel_radius_m,
|
|
required_clearance_m=sample.required_clearance_m,
|
|
fresnel_upper_60pct_m=sample.path_height_m + sample.required_clearance_m,
|
|
fresnel_upper_100pct_m=sample.path_height_m + sample.fresnel_radius_m,
|
|
fresnel_lower_60pct_m=sample.path_height_m - sample.required_clearance_m,
|
|
fresnel_lower_100pct_m=sample.path_height_m - sample.fresnel_radius_m,
|
|
type=sample.obstruction_type,
|
|
)
|
|
|
|
|
|
def _fresnel_slice_geojson(
|
|
samples: list[LosSample],
|
|
fresnel_fraction: float,
|
|
) -> dict[str, object]:
|
|
if not samples:
|
|
return {
|
|
"type": "Feature",
|
|
"properties": {"fresnel_fraction": fresnel_fraction},
|
|
"geometry": {"type": "Polygon", "coordinates": []},
|
|
}
|
|
|
|
midpoint = samples[len(samples) // 2]
|
|
epsg = _utm_epsg(midpoint.lon, midpoint.lat)
|
|
to_metric = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True)
|
|
to_wgs84 = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True)
|
|
metric_points = [
|
|
Point(*to_metric.transform(sample.lon, sample.lat)) for sample in samples
|
|
]
|
|
buffers = []
|
|
for index, point in enumerate(metric_points):
|
|
radius = max(samples[index].fresnel_radius_m * fresnel_fraction, 0.1)
|
|
buffers.append(point.buffer(radius, resolution=12))
|
|
if index + 1 < len(metric_points):
|
|
next_radius = max(samples[index + 1].fresnel_radius_m * fresnel_fraction, 0.1)
|
|
segment_radius = max(radius, next_radius)
|
|
buffers.append(
|
|
LineString([point, metric_points[index + 1]]).buffer(
|
|
segment_radius,
|
|
resolution=12,
|
|
)
|
|
)
|
|
|
|
geometry = unary_union(buffers)
|
|
geometry_wgs84 = shapely_transform(to_wgs84.transform, geometry)
|
|
return {
|
|
"type": "Feature",
|
|
"properties": {
|
|
"fresnel_fraction": fresnel_fraction,
|
|
"samples": len(samples),
|
|
},
|
|
"geometry": mapping(geometry_wgs84),
|
|
}
|
|
|
|
|
|
def _data_sources(
|
|
dem_used: bool,
|
|
include_buildings: bool,
|
|
db: Session | None,
|
|
landcover_used: bool = False,
|
|
canopy_used: bool = False,
|
|
) -> DataSources:
|
|
return DataSources(
|
|
dem=dem_used,
|
|
buildings=include_buildings and db is not None,
|
|
landcover=landcover_used,
|
|
canopy=canopy_used,
|
|
)
|