Files
Aislo/B05_wf2_Route/B05_wf2_Route_Router_Drainage.py
T
eomsangdonandClaude Fable 5 638af3ecc6 feat(B05): 세류망 흐름 새김 + 32방위 화살표 + 확장 루프 분리
1. 세류선 위 셀이 파랑으로 나오던 문제
   원인: 등고선 TIN 보간면은 실제 물골(thalweg)을 재현하지 못해, 세류선 위
   셀인데도 D8이 옆 사면으로 흘려보내 도로에 닿지 못했다.
   조치 (2가지 함께):
   - burn_stream_flow: 확정된 상류 세류망을 따라 격자 흐름 방향을 강제로 새긴다.
     세류망은 이미 도로를 건너 하류로 빠지는 물길로 확정된 자료다.
   - 사슬 추적의 종결 조건에 세류망 셀을 추가. 물이 세류에 합류한 시점에
     도로 도달이 결정된다(사슬 끝이 도로 셀에 정확히 닿지 않아도 된다).
   - StreamSplit.upstream 을 물 흐름 방향(상류->하류)으로 정렬. 방향은
     _spread_network 확산 시 진입 끝점을 기록해 한 번에 정한다.

2. 화살표 8방위 -> 32방위
   descent_azimuth: 지표면 기울기에서 연속 최급강하 방위를 뽑아 32단계로 양자화.
   D8은 연결(도로 도달 판정)에만 쓰고 표시는 실제 지형 방위를 따른다.
   세류망 새김 셀은 확정된 물길 방향을 그대로 쓴다.
   응답 바이트: 하위 6비트=32방위(32=제자리, 33=표고없음), 0x80=도로 도달.

3. 확장 루프 분리
   Basin._solve_grid 안에 있던 루프를 Flow.expand_until_closed 로 분리.
   단계 검증 미리보기는 이 경로를 타지 않는다.

4. 도로 미도달 셀 화살표를 백색으로 변경(파랑 채움 + 백색 화살표).

합성 검증: 세류 새김 389셀 전부 도로 도달, 32방위 전부 등장,
지류가 본류 중간에 합류하는 갈래를 미연결로 오탐하던 경고 제거.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
2026-07-31 18:50:21 +09:00

451 lines
19 KiB
Python

"""배수유역도 API 라우터.
구조물 측점(관 매설) 후보 제안과 배수유역 산정을 제공한다. 지형 근거는 **도엽 등고선과
세류선 GeoJSON**뿐이며(표고점은 유효 데이터가 적어 2026-07-31 사용자 지시로 제외),
좌표는 사업지 CRS(m)에서 계산하고 응답만 WGS84로 바꿔 내보낸다.
"""
import asyncio
import base64
import json
import logging
from pathlib import Path
from typing import Any
from uuid import UUID
import numpy as np
from fastapi import APIRouter
from fastapi.responses import JSONResponse
from pyproj import Transformer
from shapely.geometry import LineString, Polygon, box
from B03_FileInput.B03_FileInput_Repository import get_project_storage_relative_path
from B05_wf2_Route.B05_wf2_Route_Engine_Drainage import (
StructureCandidate,
build_route_vertices,
propose_structure_stations,
)
from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Basin import (
build_drainage_watershed,
preview_stages,
)
from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Export import write_grid_arrays, write_stage
from B05_wf2_Route.B05_wf2_Route_Engine_Watershed_Grid import (
AZIMUTH_INVALID,
AZIMUTH_SINK,
AZIMUTH_STEPS,
mask_row_spans,
)
from B05_wf2_Route.B05_wf2_Route_Repository import (
get_latest_route,
get_route_points,
get_surface_crs_epsg,
)
from common_util.common_util_storage import resolve_stored_project_path
from config.config_db import get_db_pool
from config.config_system import DRAINAGE_CACHE_DIRNAME, DRAINAGE_CACHE_FILENAME
logger = logging.getLogger(__name__)
router = APIRouter(prefix="/api/projects", tags=["B05 Route Drainage"])
# 도엽 레이어 파일명 (B04 전처리 산출물과 동일 위치)
_CONTOUR_FILE = "도엽_등고선.geojson"
_STREAM_FILE = "도엽_하천중심선.geojson"
def _sheet_dir(stored_path: str) -> Path:
return Path(resolve_stored_project_path(stored_path)) / "B04_wf1_Surface" / "processed"
def _cache_path(stored_path: str) -> Path:
"""격자 해석 캐시(.npz) 경로. 관을 옮겨도 격자를 다시 풀지 않게 여기에 남긴다."""
root = Path(resolve_stored_project_path(stored_path)) / "B05_wf2_Route"
return root / DRAINAGE_CACHE_DIRNAME / DRAINAGE_CACHE_FILENAME
def _load_features(directory: Path, filename: str) -> list[dict[str, Any]]:
"""도엽 GeoJSON을 읽어 피처 목록만 돌려준다. 없으면 빈 목록."""
path = directory / filename
if not path.exists():
return []
try:
with path.open("r", encoding="utf-8") as file:
data = json.load(file)
except (OSError, json.JSONDecodeError):
logger.warning("도엽 GeoJSON을 읽지 못했습니다: %s", path)
return []
features = data.get("features")
return features if isinstance(features, list) else []
def _reproject_features(
features: list[dict[str, Any]],
transformer: Transformer | None,
) -> list[dict[str, Any]]:
"""WGS84 도엽 좌표를 사업지 CRS(m)로 바꾼다. 거리·면적을 미터로 계산하기 위함."""
if transformer is None:
return features
converted: list[dict[str, Any]] = []
for feature in features:
geometry = feature.get("geometry")
if not geometry:
continue
coordinates = _map_coordinates(geometry.get("coordinates"), transformer)
if coordinates is None:
continue
converted.append(
{
"type": "Feature",
"properties": feature.get("properties") or {},
"geometry": {"type": geometry.get("type"), "coordinates": coordinates},
}
)
return converted
def _map_coordinates(coordinates: Any, transformer: Transformer) -> Any:
"""중첩 좌표 배열을 재귀적으로 변환한다."""
if not isinstance(coordinates, list) or not coordinates:
return None
first = coordinates[0]
if isinstance(first, (int, float)):
x, y = transformer.transform(float(coordinates[0]), float(coordinates[1]))
return [x, y]
mapped = [_map_coordinates(item, transformer) for item in coordinates]
return [item for item in mapped if item is not None]
def _candidate_payload(
candidate: StructureCandidate,
to_lonlat: Any,
) -> dict[str, Any]:
lon, lat = to_lonlat(candidate.x, candidate.y)
return {
"chainage_m": round(candidate.chainage_m, 2),
"x": candidate.x,
"y": candidate.y,
"lon": lon,
"lat": lat,
"reason": candidate.reason,
"stream_name": candidate.stream_name,
}
async def _prepare(project_id: UUID) -> dict[str, Any] | JSONResponse:
"""노선 정점·도엽 피처·좌표 변환기를 한 번에 준비한다."""
pool = get_db_pool()
async with pool.acquire() as connection:
stored_path = await get_project_storage_relative_path(connection, project_id)
route = await get_latest_route(connection, project_id)
if not route:
return JSONResponse(
status_code=404,
content={"status": "error", "message": "확정된 경로가 없습니다."},
)
points = await get_route_points(connection, int(route["id"]))
surface_model_id = route.get("surface_model_id")
epsg = await get_surface_crs_epsg(
connection, project_id, int(surface_model_id) if surface_model_id else 0
)
vertices = build_route_vertices(points)
if len(vertices) < 2:
return JSONResponse(
status_code=400,
content={"status": "error", "message": "노선 좌표가 부족합니다."},
)
source_crs = f"EPSG:{epsg}" if epsg else "EPSG:5186"
to_lonlat_transformer = Transformer.from_crs(source_crs, "EPSG:4326", always_xy=True)
to_metric_transformer = Transformer.from_crs("EPSG:4326", source_crs, always_xy=True)
directory = _sheet_dir(stored_path)
streams = _reproject_features(_load_features(directory, _STREAM_FILE), to_metric_transformer)
contour_features = _reproject_features(
_load_features(directory, _CONTOUR_FILE), to_metric_transformer
)
return {
"route_id": int(route["id"]),
"vertices": vertices,
"route_line": LineString([(vertex.x, vertex.y) for vertex in vertices]),
"streams": streams,
"contours": contour_features,
"stored_path": stored_path,
"cache_path": _cache_path(stored_path),
"to_lonlat": lambda x, y: to_lonlat_transformer.transform(x, y),
}
@router.get("/{project_id}/drainage/candidates", response_model=None)
async def get_structure_candidates(project_id: UUID) -> dict[str, Any] | JSONResponse:
"""관 매설 구조물 측점 후보를 제안한다(세류 교차 + 300m 보충, 성토부 제외)."""
prepared = await _prepare(project_id)
if isinstance(prepared, JSONResponse):
return prepared
candidates = propose_structure_stations(prepared["vertices"], prepared["streams"])
to_lonlat = prepared["to_lonlat"]
return {
"status": "success",
"project_id": str(project_id),
"route_id": prepared["route_id"],
"candidates": [_candidate_payload(candidate, to_lonlat) for candidate in candidates],
}
@router.get("/{project_id}/drainage/primary-region", response_model=None)
async def get_primary_region(project_id: UUID) -> dict[str, Any] | JSONResponse:
"""1차 배수유역 근거를 돌려준다 — 단계 검증용, TIN·흐름 계산은 하지 않는다.
도로 교차점 상류로 이어진 세류망, 제외된 하류망, 그 상류망을 반경 버퍼한 1차 영역,
그 bbox로 잡은 격자 정보를 함께 준다. 같은 내용을 영구저장소에 GeoJSON으로도 남겨
QGIS 등으로 직접 열어 대조할 수 있게 한다.
"""
prepared = await _prepare(project_id)
if isinstance(prepared, JSONResponse):
return prepared
preview = await asyncio.to_thread(
preview_stages,
prepared["vertices"],
prepared["contours"],
prepared["streams"],
)
if preview is None:
return JSONResponse(
status_code=400,
content={"status": "error", "message": "1차 배수유역을 정할 등고선이 없습니다."},
)
region = preview.region
to_lonlat = prepared["to_lonlat"]
spec = region.spec
payload = {
"status": "success",
"project_id": str(project_id),
"route_id": prepared["route_id"],
"radius_m": region.radius_m,
# 채택된 상류 세류망 = 1차 영역의 기준선.
"upstream_lines": [_line_lonlat(line, to_lonlat) for line in region.split.upstream],
# 도로 아래로 이어진 하류망 — 판정이 맞는지 눈으로 대조하기 위해 함께 준다.
"downstream_lines": [_line_lonlat(line, to_lonlat) for line in region.split.downstream],
"no_contact_count": region.split.no_contact,
# 노선이 1차 영역 밖으로 나간 길이(m). 크면 반경을 올려야 한다는 신호.
"road_outside_m": round(region.road_outside_m, 1),
# 1차 영역(버퍼 합집합) 외곽 링 목록.
"region_rings": _polygon_rings(region.area, to_lonlat),
"grid": {
"cell_m": spec.cell_m,
"rows": spec.n_rows,
"cols": spec.n_cols,
# bbox 전체 셀 수와, 1차 영역에 걸쳐 실제로 생성된 셀 수.
"bbox_cells": spec.size,
"cells": region.active_cells,
"width_m": round(spec.n_cols * spec.cell_m, 1),
"height_m": round(spec.n_rows * spec.cell_m, 1),
# 격자 bbox 링. 프론트는 이 사각형을 rows×cols로 나눠 행·열 좌표를 얻는다.
"bbox_lonlat": _grid_bbox_lonlat(spec, to_lonlat),
# 실제 생성된 셀을 행별 연속 구간 [행, 시작열, 끝열]으로 압축해 보낸다.
# 셀을 낱개로 보내면 수십만 건이라 응답이 감당되지 않는다.
"row_spans": [list(span) for span in mask_row_spans(region.cell_mask)]
if region.cell_mask is not None
else [],
},
# 셀별 흐름 방향과 도로 도달 여부. row_spans 순서(행 → 구간 → 열 오름차순)로 1바이트씩.
"flow": _flow_payload(preview, region),
}
# 단계 산출물을 영구저장소에 남긴다 — 기능을 붙일 때마다 여기에 단계가 하나씩 늘어난다.
payload["saved_to"] = write_stage(
prepared["stored_path"],
"primary_region",
{
"primary_region": _as_polygons(region.area),
"upstream": region.split.upstream,
"downstream": region.split.downstream,
"route": [prepared["route_line"]],
"grid_bbox": [_grid_bbox_polygon(spec)],
},
{
"radius_m": region.radius_m,
"road_outside_m": payload["road_outside_m"],
"no_contact_count": region.split.no_contact,
"grid": {key: value for key, value in payload["grid"].items() if key != "bbox_lonlat"},
},
to_lonlat,
)
_write_stage_arrays(prepared["stored_path"], preview, region, spec)
return payload
def _write_stage_arrays(stored_path: str, preview: Any, region: Any, spec: Any) -> None:
"""격자 규모 배열(셀 마스크·흐름 방향·도달 여부)을 단계별 `.npz`로 남긴다."""
if region.cell_mask is not None:
write_grid_arrays(
stored_path,
"primary_region",
spec,
{"mask": region.cell_mask},
{"cells": region.active_cells, "bbox_cells": spec.size},
)
flow = preview.flow
if flow is None:
return
arrays = {
"direction": flow.direction.reshape(spec.n_rows, spec.n_cols),
"reaches_road": flow.reaches_road.reshape(spec.n_rows, spec.n_cols),
"analyzed": flow.analyzed.reshape(spec.n_rows, spec.n_cols),
}
if flow.burned is not None:
arrays["burned"] = flow.burned.reshape(spec.n_rows, spec.n_cols)
write_grid_arrays(
stored_path,
"flow_direction",
spec,
arrays,
{
"azimuth_steps": AZIMUTH_STEPS,
"sink_code": AZIMUTH_SINK,
"invalid_code": AZIMUTH_INVALID,
"analyzed": int(flow.analyzed.sum()),
"reaches_road": int((flow.reaches_road & flow.analyzed).sum()),
"no_road": int((~flow.reaches_road & flow.analyzed).sum()),
"burned": 0 if flow.burned is None else int(flow.burned.sum()),
"outer_seeds": flow.outer_seeds,
"interior_seeds": flow.interior_seeds,
},
)
def _flow_payload(preview: Any, region: Any) -> dict[str, Any] | None:
"""셀별 흐름 방향·도로 도달 여부를 바이트 배열로 압축한다.
셀이 수십만 개라 JSON 객체로는 못 보낸다. 셀 하나당 1바이트로 줄이고 base64로 싣는다:
하위 6비트(0x3F) = 32방위 코드(0~31, 0=화면 오른쪽·시계방향), 32=제자리, 33=표고 없음
최상위 비트(0x80) = 도로 도달(적색). 꺼져 있으면 미도달(백색 화살표).
바이트 순서는 `grid.row_spans`를 행 → 구간 → 열 오름차순으로 훑은 순서와 같다.
"""
flow = preview.flow
if flow is None or region.cell_mask is None:
return None
order = np.flatnonzero(region.cell_mask.reshape(-1))
analyzed = flow.analyzed[order]
reaches = flow.reaches_road[order]
packed = np.clip(flow.direction[order], 0, AZIMUTH_INVALID).astype(np.uint8)
packed |= np.where(reaches, 0x80, 0).astype(np.uint8)
burned = flow.burned
return {
"encoding": "base64-uint8",
"azimuth_steps": AZIMUTH_STEPS,
"sink_code": AZIMUTH_SINK,
"invalid_code": AZIMUTH_INVALID,
"cells": int(order.size),
"reaches_road": int((reaches & analyzed).sum()),
"no_road": int((~reaches & analyzed).sum()),
# 격자에는 있으나 등고선 TIN 밖이라 표고가 없어 판정 못한 셀.
"unanalyzed": int((~analyzed).sum()),
# 확정된 상류 세류망을 따라 흐름을 강제로 새긴 셀 수.
"burned": 0 if burned is None else int(burned[order].sum()),
"outer_seeds": flow.outer_seeds,
"interior_seeds": flow.interior_seeds,
"data": base64.b64encode(packed.tobytes()).decode("ascii"),
}
def _as_polygons(geometry: Any) -> list[Any]:
if geometry is None or geometry.is_empty:
return []
return list(geometry.geoms) if geometry.geom_type == "MultiPolygon" else [geometry]
def _grid_bbox_polygon(spec: Any) -> Polygon:
x_max = spec.x_min + spec.n_cols * spec.cell_m
y_min = spec.y_max - spec.n_rows * spec.cell_m
return box(spec.x_min, y_min, x_max, spec.y_max)
def _line_lonlat(line: Any, to_lonlat: Any) -> list[list[float]]:
return [list(to_lonlat(x, y)) for x, y in line.coords]
def _polygon_rings(geometry: Any, to_lonlat: Any) -> list[list[list[float]]]:
"""폴리곤/멀티폴리곤의 외곽 링만 뽑아 lonlat으로 바꾼다."""
if geometry is None or geometry.is_empty:
return []
parts = geometry.geoms if geometry.geom_type == "MultiPolygon" else [geometry]
return [[list(to_lonlat(x, y)) for x, y in part.exterior.coords] for part in parts]
def _grid_bbox_lonlat(spec: Any, to_lonlat: Any) -> list[list[float]]:
x_min = spec.x_min
x_max = spec.x_min + spec.n_cols * spec.cell_m
y_max = spec.y_max
y_min = spec.y_max - spec.n_rows * spec.cell_m
corners = ((x_min, y_min), (x_min, y_max), (x_max, y_max), (x_max, y_min), (x_min, y_min))
return [list(to_lonlat(x, y)) for x, y in corners]
@router.post("/{project_id}/drainage/basins", response_model=None)
async def post_drainage_basins(
project_id: UUID,
payload: dict[str, Any] | None = None,
) -> dict[str, Any] | JSONResponse:
"""격자 흐름 해석으로 배수유역과 관 배치를 산정한다.
payload에 `chainages`(누가거리 목록)를 주면 그 위치로 관을 확정하고, 없으면 세류
교차 + 최소 보충으로 자동 배치한다. 격자 해석은 `.npz` 캐시를 재사용하므로 관만
옮기는 재요청은 세부유역 분할만 다시 돈다.
"""
prepared = await _prepare(project_id)
if isinstance(prepared, JSONResponse):
return prepared
raw_chainages = (payload or {}).get("chainages")
confirmed = _parse_chainages(raw_chainages) if isinstance(raw_chainages, list) else []
# 격자 해석은 수백만 셀 numpy 연산이라 이벤트 루프를 막지 않도록 스레드로 뺀다.
result = await asyncio.to_thread(
build_drainage_watershed,
prepared["vertices"],
prepared["contours"],
prepared["streams"],
confirmed,
prepared["cache_path"],
)
to_lonlat = prepared["to_lonlat"]
return {
"status": "success",
"project_id": str(project_id),
"route_id": prepared["route_id"],
# 계획선 위 배관(관 매설) 지점 — 유역이 없는 관도 마커로 표시해야 하므로 별도 목록.
"pipes": [_candidate_payload(candidate, to_lonlat) for candidate in result.pipes],
# 2차 전체 배수유역 외곽선 = 분수령. 세부유역은 전부 이 안쪽에 들어간다.
"main_polygon_lonlat": [list(to_lonlat(x, y)) for x, y in result.main_boundary_xy],
# 도로 위 흐름 강도 — [누가거리 m, 그 지점으로 모이는 상류 면적 ㎡].
"strength_profile": [
[round(chainage, 1), round(area, 1)] for chainage, area in result.strength_profile
],
"grid_cell_m": result.grid_cell_m,
"basins": [
{
"index": basin.index,
"chainage_m": round(basin.chainage_m, 2),
# 관(배관) 매설 지점 좌표 — 계획선 위 마커 렌더용.
"outlet_lonlat": list(to_lonlat(basin.outlet_x, basin.outlet_y)),
"polygon_lonlat": [list(to_lonlat(x, y)) for x, y in basin.boundary_xy],
"area_m2": round(basin.area_m2, 1),
"relief_m": round(basin.relief_m, 2),
"flow_length_m": round(basin.flow_length_m, 1),
# 관경 수식 미확정 — 산정 함수가 None을 돌려주면 프론트가 "미정"으로 표기한다.
"pipe_diameter_mm": basin.pipe_diameter_mm,
}
for basin in result.basins
],
}
def _parse_chainages(values: list[Any]) -> list[float]:
"""사용자가 확정·편집한 누가거리 목록을 숫자로 정리한다."""
parsed: list[float] = []
for value in values:
try:
parsed.append(float(value))
except (TypeError, ValueError):
continue
return parsed