diff --git a/B04_PreProcess/B04_PreProcess_Router_Watershed.py b/B04_PreProcess/B04_PreProcess_Router_Watershed.py index c6faa675..1ac02733 100644 --- a/B04_PreProcess/B04_PreProcess_Router_Watershed.py +++ b/B04_PreProcess/B04_PreProcess_Router_Watershed.py @@ -8,7 +8,6 @@ """ import asyncio -import base64 import json import logging import math @@ -16,33 +15,39 @@ 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 Point, Polygon, box +from shapely.geometry import Point from B03_FileInput.B03_FileInput_Repository import get_project_storage_relative_path from B04_PreProcess.B04_PreProcess_Engine_Watershed_Analyze import preview_stages from B04_PreProcess.B04_PreProcess_Engine_Watershed_Export import ( drainage_dir, - write_grid_arrays, write_stage, ) from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( - AZIMUTH_INVALID, - AZIMUTH_SINK, - AZIMUTH_STEPS, mask_row_spans, ) +from B04_PreProcess.B04_PreProcess_Router_Watershed_Output import ( + _as_polygons, + _boundary_geometry, + _flow_payload, + _grid_bbox_lonlat, + _grid_bbox_polygon, + _line_lonlat, + _polygon_rings, + _write_road_routing, + _write_stage_arrays, +) from B05_Profile.B05_Profile_Repository import get_surface_crs_epsg +from common_util.common_util_crs import resolve_project_crs from common_util.common_util_route_geometry import ( StructureCandidate, find_planned_route_file, load_design_route, read_planned_route, ) -from common_util.common_util_crs import resolve_project_crs from common_util.common_util_storage import resolve_stored_project_path from common_util.common_util_surface_confirmation import get_surface_confirmation_params from common_util.common_util_wamis_rainfall import ( @@ -511,204 +516,3 @@ async def get_primary_region( # 응답 자체를 캐시로 남긴다 — 다음 조회는 배열을 재조립하지 않고 이 파일을 그대로 준다. _save_response(prepared["stored_path"], payload) return payload - - -def _write_road_routing( - stored_path: str, preview: Any, spec: Any, route_line: Any, to_lonlat: Any -) -> None: - """B05가 세부유역을 나눌 때 쓸 최소 산출물을 남긴다. - - B05는 일반 사용자용이라 가벼워야 한다. 화살표(방향 코드)·밴드 표고 같은 확인용 배열은 - 빼고, **셀 → 도로 셀 귀속**과 도로 셀 제원만 담는다. 여기에 표고를 함께 넣는 이유는 - 유역 낙차를 내려면 셀 표고가 필요해서다(2026-07-31 사용자 지시). - """ - routing = preview.routing - road = preview.road - if routing is None or road is None or road.count == 0: - return - write_grid_arrays( - stored_path, - "road_routing", - spec, - { - "road_slot": routing.road_slot, - "path_length": routing.path_length, - "strength": routing.strength, - "road_cell_index": road.cell_index, - "road_chainage": road.chainage, - "elevation": preview.terrain.elevation.reshape(-1), - }, - { - "road_cells": road.count, - "reached_cells": int((routing.road_slot >= 0).sum()), - "basin_area_m2": round(preview.basin_area_m2, 1), - "pipe_count": len(preview.pipes), - }, - ) - # B05가 그대로 그릴 기하 — 계획도로선 · 기본 배관 · 2차 전체 배수유역, 이 셋뿐이다. - write_stage( - stored_path, - "road_routing", - { - "route": [route_line], - "basin_boundary": _boundary_geometry(preview.basin_boundary_xy), - "pipe": [ - ( - Point(pipe.x, pipe.y), - {"chainage_m": round(pipe.chainage_m, 2), "reason": pipe.reason}, - ) - for pipe in preview.pipes - ], - # 평균 흐름 화살표 — B05도 같은 그림을 그려야 하므로 여기 함께 남긴다. - "flow_arrow": [ - ( - Point(x, y), - { - # B05 화면은 사업지 CRS(m)로 그리므로 미터 좌표도 함께 남긴다. - "x": round(x, 2), - "y": round(y, 2), - "azimuth_deg": round(math.degrees(angle), 1), - "reaches_road": reaches, - "cells": cells, - }, - ) - for x, y, angle, reaches, cells in preview.flow_arrows - ], - }, - { - "basin_area_m2": round(preview.basin_area_m2, 1), - "pipe_count": len(preview.pipes), - "route_length_m": round(route_line.length, 1), - "arrow_count": len(preview.flow_arrows), - "arrow_spacing_m": DRAINAGE_ARROW_SPACING_M, - }, - to_lonlat, - ) - - -def _write_stage_arrays(stored_path: str, preview: Any, domain: Any, spec: Any) -> None: - """격자 규모 배열(셀 마스크·흐름 방향·도달 여부)을 단계별 `.npz`로 남긴다.""" - if domain is not None: - write_grid_arrays( - stored_path, - "primary_region", - spec, - {"mask": domain}, - { - "cells": int(domain.sum()), - "bbox_cells": spec.size, - "expand_rounds": preview.expand_rounds, - "expand_closed": preview.expand_closed, - }, - ) - 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), - # 사후 진단용 — 수신 셀이 있어야 사슬을 다시 따라가 볼 수 있다. - "receiver": preview.terrain.receiver.reshape(spec.n_rows, spec.n_cols), - } - if flow.burned is not None: - arrays["burned"] = flow.burned.reshape(spec.n_rows, spec.n_cols) - if preview.descent is not None: - arrays["band_elevation"] = preview.descent.band_elevation - # ⑥ 흐름 강도 곡선 — 기하가 아니라 수치 곡선이라 GeoJSON이 아닌 여기에 함께 담는다. - if preview.strength_profile: - curve = np.asarray(preview.strength_profile, dtype=np.float64) - arrays["strength_chainage_m"] = curve[:, 0] - arrays["strength_area_m2"] = curve[:, 1] - 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, - "strength_points": len(preview.strength_profile), - "strength_total_m2": round(sum(area for _, area in preview.strength_profile), 1), - }, - ) - - -def _flow_payload(preview: Any, domain: 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 domain is None: - return None - order = np.flatnonzero(domain.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 _boundary_geometry(ring: list[tuple[float, float]]) -> list[Any]: - """2차 유역 외곽 링을 저장용 폴리곤으로 만든다(정점 3개 미만이면 비운다).""" - return [Polygon(ring)] if len(ring) >= 4 else [] - - -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] diff --git a/B04_PreProcess/B04_PreProcess_Router_Watershed_Output.py b/B04_PreProcess/B04_PreProcess_Router_Watershed_Output.py new file mode 100644 index 00000000..87affee5 --- /dev/null +++ b/B04_PreProcess/B04_PreProcess_Router_Watershed_Output.py @@ -0,0 +1,225 @@ +"""배수유역 분석 **산출물 기록·응답 기하** 헬퍼. + +`B04_PreProcess_Router_Watershed` 가 700줄을 넘겨 떼어낸 조각이다(2026-09-04). +B05 가 읽을 `.npz`·스테이지 기하를 남기는 쓰기 함수와, 응답용 좌표 변환 헬퍼만 모았다. +함수 본문·값은 옮기기 전 그대로이고, 라우터가 이 모듈에서 가져다 쓴다. +""" + +import base64 +import math +from typing import Any + +import numpy as np +from shapely.geometry import Point, Polygon, box + +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Export import ( + write_grid_arrays, + write_stage, +) +from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( + AZIMUTH_INVALID, + AZIMUTH_SINK, + AZIMUTH_STEPS, +) +from config.config_system import DRAINAGE_ARROW_SPACING_M + + +def _write_road_routing( + stored_path: str, preview: Any, spec: Any, route_line: Any, to_lonlat: Any +) -> None: + """B05가 세부유역을 나눌 때 쓸 최소 산출물을 남긴다. + + B05는 일반 사용자용이라 가벼워야 한다. 화살표(방향 코드)·밴드 표고 같은 확인용 배열은 + 빼고, **셀 → 도로 셀 귀속**과 도로 셀 제원만 담는다. 여기에 표고를 함께 넣는 이유는 + 유역 낙차를 내려면 셀 표고가 필요해서다(2026-07-31 사용자 지시). + """ + routing = preview.routing + road = preview.road + if routing is None or road is None or road.count == 0: + return + write_grid_arrays( + stored_path, + "road_routing", + spec, + { + "road_slot": routing.road_slot, + "path_length": routing.path_length, + "strength": routing.strength, + "road_cell_index": road.cell_index, + "road_chainage": road.chainage, + "elevation": preview.terrain.elevation.reshape(-1), + }, + { + "road_cells": road.count, + "reached_cells": int((routing.road_slot >= 0).sum()), + "basin_area_m2": round(preview.basin_area_m2, 1), + "pipe_count": len(preview.pipes), + }, + ) + # B05가 그대로 그릴 기하 — 계획도로선 · 기본 배관 · 2차 전체 배수유역, 이 셋뿐이다. + write_stage( + stored_path, + "road_routing", + { + "route": [route_line], + "basin_boundary": _boundary_geometry(preview.basin_boundary_xy), + "pipe": [ + ( + Point(pipe.x, pipe.y), + {"chainage_m": round(pipe.chainage_m, 2), "reason": pipe.reason}, + ) + for pipe in preview.pipes + ], + # 평균 흐름 화살표 — B05도 같은 그림을 그려야 하므로 여기 함께 남긴다. + "flow_arrow": [ + ( + Point(x, y), + { + # B05 화면은 사업지 CRS(m)로 그리므로 미터 좌표도 함께 남긴다. + "x": round(x, 2), + "y": round(y, 2), + "azimuth_deg": round(math.degrees(angle), 1), + "reaches_road": reaches, + "cells": cells, + }, + ) + for x, y, angle, reaches, cells in preview.flow_arrows + ], + }, + { + "basin_area_m2": round(preview.basin_area_m2, 1), + "pipe_count": len(preview.pipes), + "route_length_m": round(route_line.length, 1), + "arrow_count": len(preview.flow_arrows), + "arrow_spacing_m": DRAINAGE_ARROW_SPACING_M, + }, + to_lonlat, + ) + + +def _write_stage_arrays(stored_path: str, preview: Any, domain: Any, spec: Any) -> None: + """격자 규모 배열(셀 마스크·흐름 방향·도달 여부)을 단계별 `.npz`로 남긴다.""" + if domain is not None: + write_grid_arrays( + stored_path, + "primary_region", + spec, + {"mask": domain}, + { + "cells": int(domain.sum()), + "bbox_cells": spec.size, + "expand_rounds": preview.expand_rounds, + "expand_closed": preview.expand_closed, + }, + ) + 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), + # 사후 진단용 — 수신 셀이 있어야 사슬을 다시 따라가 볼 수 있다. + "receiver": preview.terrain.receiver.reshape(spec.n_rows, spec.n_cols), + } + if flow.burned is not None: + arrays["burned"] = flow.burned.reshape(spec.n_rows, spec.n_cols) + if preview.descent is not None: + arrays["band_elevation"] = preview.descent.band_elevation + # ⑥ 흐름 강도 곡선 — 기하가 아니라 수치 곡선이라 GeoJSON이 아닌 여기에 함께 담는다. + if preview.strength_profile: + curve = np.asarray(preview.strength_profile, dtype=np.float64) + arrays["strength_chainage_m"] = curve[:, 0] + arrays["strength_area_m2"] = curve[:, 1] + 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, + "strength_points": len(preview.strength_profile), + "strength_total_m2": round(sum(area for _, area in preview.strength_profile), 1), + }, + ) + + +def _flow_payload(preview: Any, domain: 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 domain is None: + return None + order = np.flatnonzero(domain.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 _boundary_geometry(ring: list[tuple[float, float]]) -> list[Any]: + """2차 유역 외곽 링을 저장용 폴리곤으로 만든다(정점 3개 미만이면 비운다).""" + return [Polygon(ring)] if len(ring) >= 4 else [] + + +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]