diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_Api.ts b/D03_PlanAnalysis/D03_PlanAnalysis_Api.ts index bd9be9c4b..baff2fb76 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_Api.ts +++ b/D03_PlanAnalysis/D03_PlanAnalysis_Api.ts @@ -179,5 +179,27 @@ export interface PlanBand { collection: GeoJsonCollection; } +/** 선 겹(584-1) — 받은 범위 안 선 + 가닥 · 길이 m. */ +export interface PlanLines extends GeoJsonCollection { + summary: { count: number; length_m: number }; +} + +/** 단층 · 지질구조선(브이월드 WFS · 서버 `LINE_KINDS` 열쇠). */ +export const fetchPlanLines = (projectId: string, kind: string): Promise => + request(`${base(projectId)}/layers/lines?kind=${encodeURIComponent(kind)}`); + +/** 0차곡(584-1 · 지형 추정) — 위 집수 넓이 [min_m2, max_m2) 선 + 가닥 · 길이 m. */ +export interface PlanZeroOrder { + min_m2: number; + max_m2: number; + cell_m: number; + count: number; + length_m: number; + collection: GeoJsonCollection; +} + +export const fetchPlanZeroOrder = (projectId: string): Promise => + request(`${base(projectId)}/layers/zero-order`); + export const fetchPlanBands = (projectId: string): Promise<{ stream: PlanBand; road: PlanBand }> => request(`${base(projectId)}/layers/bands`); diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py b/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py index 6dbd2492b..25a991b2e 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py +++ b/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py @@ -13,6 +13,9 @@ - 산림입지도 = 브이월드 WFS `lt_c_fsdifrsts`(원자료 산림청 산림입지토양도 · PLAN 581-2) · 받은 범위로 잘라 토양형 코드 · 이름 · 넓이만 · 저장 = `processed/산림입지도_bounds.geojson` + 곁 열쇠 +- 선 겹 = 브이월드 WFS 단층 `lt_l_gimsfault` · 지질구조선 `lt_l_gimslinea`(PLAN 584-1) · + 받은 범위로 잘라 가닥 · 길이만 · 저장 = `processed/<이름>_bounds.geojson` + 곁 열쇠 +- 0차곡 = `D03_PlanAnalysis_ZeroOrder` - 계류 100 m · 도로 300 m 띠 = 도엽 하천중심선(세류 뺌) · 도로중심선(미분류 · 소로 · 면리간도로 뺌) 합친 띠 · 받은 범위로 자름 · 파생 저장본(열쇠 = 두 파일 지문) @@ -473,6 +476,67 @@ def ensure_soils(processed: Path, box: Box, epsg: int) -> dict[str, Any]: return {"path": path, "fetched": True} +# ── 선 겹(단층 · 지질구조선) ───────────────────────────── + +# 열쇠 → (이름, 브이월드 WFS 겹 이름) — 새 선 겹은 여기 한 줄(PLAN 584-1) +LINE_KINDS: dict[str, tuple[str, str]] = { + "faults": ("단층", "lt_l_gimsfault"), + "lineaments": ("지질구조선", "lt_l_gimslinea"), +} +LINE_VERSION = 1 + + +def clip_lines( + found: list[dict[str, Any]], box: Box, epsg: int +) -> tuple[list[dict[str, Any]], dict[str, Any]]: + """선 지물 → 받은 범위로 자른 선(길이 m · 원 속성 `legend`) + {가닥, 길이 m}.""" + from shapely.geometry import box as rect + from shapely.geometry import mapping, shape + from shapely.ops import transform + + to_m, to_ll = _transforms(epsg) + frame = rect(*box) + features: list[dict[str, Any]] = [] + total = 0.0 + for feature in dedup_parcels(found): + try: + part = transform(to_m, shape(feature["geometry"])).intersection(frame) + except (KeyError, TypeError, ValueError, AttributeError): + continue + if part.is_empty or part.length <= 0: + continue + total += part.length + props = feature.get("properties") or {} + features.append( + { + "type": "Feature", + "geometry": mapping(transform(to_ll, part)), + "properties": {"legend": props.get("legend"), "length_m": round(part.length, 1)}, + } + ) + return features, {"count": len(features), "length_m": round(total, 1)} + + +def ensure_lines(processed: Path, box: Box, epsg: int, kind: str) -> dict[str, Any]: + """선 겹 — 열쇠 상자가 덮으면 그대로 · 아니면 범위 통째로 다시 받음 · 실패 = 오류.""" + label, typename = LINE_KINDS[kind] + path, key_path = processed / f"{label}_bounds.geojson", processed / f"{label}_bounds.key.json" + key = _read_json(key_path) or {} + if path.is_file() and key.get("version") == LINE_VERSION and _covers(key.get("box"), box): + return {"path": path, "fetched": False} + got, failed = fetch_tile(box, epsg, 0, typename) + if failed: + raise RuntimeError(f"{label}을 받지 못했습니다 — 브이월드 WFS 실패") + features, summary = clip_lines(got, box, epsg) + processed.mkdir(parents=True, exist_ok=True) + atomic_write_json(path, {"type": "FeatureCollection", "features": features, "summary": summary}) + atomic_write_json( + key_path, {"version": LINE_VERSION, "box": [float(v) for v in box], "epsg": int(epsg)} + ) + logger.info("D %s 받기: 선 %d · %.0f m", label, summary["count"], summary["length_m"]) + return {"path": path, "fetched": True} + + # ── 계류 · 도로 띠 ─────────────────────────────────────── BAND_KIND = "plan_band" diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_Router.py b/D03_PlanAnalysis/D03_PlanAnalysis_Router.py index b97221240..00e016170 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_Router.py +++ b/D03_PlanAnalysis/D03_PlanAnalysis_Router.py @@ -4,7 +4,8 @@ - `POST …/plan/terrain` — 사업 전체 지형 한 벌 받기(뒤에서 · 누를 때만) · 끝나면 `surface_models` 행 - `GET …/plan/progress` — 일의 진행률 - `GET …/plan/layers/sheet?key=` · `…/cadastral` · `…/raster?kind=slope|hillshade` · - `…/zones` · `…/soils` · `…/bands` — 지도 겹(PLAN 569-1 · 581-1 · 581-2) + `…/zones` · `…/soils` · `…/bands` · `…/lines?kind=` · `…/zero-order` — 지도 겹 + (PLAN 569-1 · 581-1 · 581-2 · 584-1) 범위 = 전체 시점 · 종점(D02 정본 `load_base`) + 그은 점(`extra`) 상자 + 여유. 노선 정본(손 그은 선)은 @@ -30,6 +31,7 @@ from config.config_db import get_db_pool from D02_PlanBase import D02_PlanBase_Router as base from D03_PlanAnalysis import D03_PlanAnalysis_Layers as layers from D03_PlanAnalysis import D03_PlanAnalysis_Terrain as terrain +from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zero_order logger = logging.getLogger(__name__) router = APIRouter(prefix="/api/projects", tags=["D03 Plan Analysis"]) @@ -244,3 +246,36 @@ async def plan_layer_bands(project_id: UUID) -> dict[str, Any]: except FileNotFoundError as error: raise HTTPException(status_code=404, detail=str(error)) from error return {"status": "success", **value} + + +@router.get("/{project_id}/plan/layers/lines", response_model=None) +async def plan_layer_lines(project_id: UUID, kind: str, request: Request) -> Response: + """선 겹 GeoJSON(단층 · 지질구조선 · `summary` = 가닥 · 길이 m) — 범위 밖이면 받음.""" + if kind not in layers.LINE_KINDS: + raise HTTPException(status_code=400, detail=f"없는 선 겹: {kind}") + root, epsg, box = await _layer_root(project_id) + processed = root / terrain.STAGE / "processed" + try: + got = await asyncio.to_thread(layers.ensure_lines, processed, box, epsg, kind) + except RuntimeError as error: + raise HTTPException(status_code=502, detail=str(error)) from error + return cached_file_response(request, got["path"], "application/geo+json") + + +@router.get("/{project_id}/plan/layers/zero-order", response_model=None) +async def plan_layer_zero_order( + project_id: UUID, + min_m2: float = zero_order.ZERO_ORDER_MIN_M2, + max_m2: float = zero_order.ZERO_ORDER_MAX_M2, +) -> dict[str, Any]: + """0차곡 선(지형 추정 · 위 집수 `min_m2` ~ `max_m2`) — 격자 지문 · 임계 같으면 저장본.""" + if not 0 < min_m2 < max_m2: + raise HTTPException(status_code=400, detail="0 < min_m2 < max_m2 여야 합니다.") + root, epsg, _ = await _layer_root(project_id) + try: + value = await asyncio.to_thread( + zero_order.zero_order_layer, root, terrain.models_dir_of(root), epsg, min_m2, max_m2 + ) + except FileNotFoundError as error: + raise HTTPException(status_code=404, detail=str(error)) from error + return {"status": "success", **value} diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts b/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts index e5bde1e2c..b6709226e 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts +++ b/D03_PlanAnalysis/D03_PlanAnalysis_UI_LayerDefs.ts @@ -1,15 +1,19 @@ /* ============================================================================= * D03_PlanAnalysis_UI_LayerDefs.ts - * D03 바탕 겹 정의(PLAN 569-2 · 581-1 · 581-2) — 겹 단추 묶음 · 처음 값 · 보호구역 · 토양형 색 · 면 겹 채우기. + * D03 바탕 겹 정의(PLAN 569-2 · 581-1 · 581-2 · 584-1) — 겹 단추 묶음 · 처음 값 · 보호구역 · 토양형 색 · + * 선 겹(단층 · 지질구조선 · 0차곡) 꼴 · 면 겹 채우기 · 선 긋기 · 그림 겹 놓기. * 부르는 곳 = `D03_PlanAnalysis_UI_Layers.ts`(받기 · 그리기 차례 · 범례). * ========================================================================== */ import type { IconName } from "@ui/ui_template_icons"; import { + drawPreparedLayer, normalizedToScreen, type PreparedLayer, type ViewState, } from "../B04_PreProcess/B04_PreProcess_UI_MapRender"; +import type { PlanRaster } from "./D03_PlanAnalysis_Api"; +import type { ScreenMap } from "./D03_PlanAnalysis_UI_Map"; export type PlanLayerKey = | "satellite" @@ -23,6 +27,9 @@ export type PlanLayerKey = | "cadastral" | "zones" | "soils" + | "faults" + | "lineaments" + | "zeroOrder" | "bands" | "building" | "place"; @@ -67,6 +74,19 @@ export const GROUPS: readonly (readonly LayerDef[])[] = [ title: "산림입지도(토양형 · 산림청 산림입지토양도 · 참고용)", on: false, }, + { key: "faults", icon: "bolt", title: "단층(브이월드 지질 · 참고용)", on: false }, + { + key: "lineaments", + icon: "lineDashed", + title: "지질구조선(브이월드 지질 · 참고용)", + on: false, + }, + { + key: "zeroOrder", + icon: "arrowsJoin", + title: "0차곡(지형 추정 골짜기 머리 · 위 집수 0.1~1 ha · 곡두 통과 피함)", + on: false, + }, { key: "bands", icon: "borderOuter", @@ -89,6 +109,23 @@ export const ZONE_COLORS: Record = { wildlife: "--draw-map-zone-wildlife", greenbelt: "--draw-map-zone-greenbelt", }; +/** 선 겹 — 이름 · 그림 꼴(서버 `LINE_KINDS` 열쇠 = `kind` · 없으면 0차곡). 새 선 겹 = 여기 한 줄. */ +export type LineLayerKey = "faults" | "lineaments" | "zeroOrder"; +export const LINE_LAYERS: Record< + LineLayerKey, + { label: string; kind?: string; color: string; width: number; dash: number[] } +> = { + faults: { label: "단층", kind: "faults", color: "--draw-map-fault", width: 2.4, dash: [] }, + lineaments: { + label: "지질구조선", + kind: "lineaments", + color: "--draw-map-lineament", + width: 1.8, + dash: [8, 4], + }, + zeroOrder: { label: "0차곡", color: "--draw-map-zero-order", width: 1.8, dash: [4, 3] }, +}; +export const LINE_KEYS = Object.keys(LINE_LAYERS) as LineLayerKey[]; /** 토양형 색 — 받은 범위 안 코드 차례로 돌려 씀. */ export const SOIL_COLORS = [1, 2, 3, 4, 5, 6].map((n) => `--draw-map-soil-${n}`); /** 도로 중 기존 임도로 짐작하는 구분. */ @@ -142,3 +179,37 @@ export function fillLayer( } context.setLineDash([]); } + +export const loadImage = (src: string): Promise => + new Promise((done, fail) => { + const image = new Image(); + image.onload = () => done(image); + image.onerror = () => fail(new Error("그림을 읽지 못했습니다")); + image.src = src; + }); + +/** 선 겹 긋기(점 솎음 `dot` · `everyM` = 등고 간격 거름). */ +export function stroke( + context: CanvasRenderingContext2D, + map: ScreenMap, + layer: PreparedLayer, + color: string, + width: number, + everyM?: number, +): void { + context.strokeStyle = color; + context.lineWidth = width; + drawPreparedLayer(context, layer, map.view, "dot", everyM); +} + +/** 그림 겹(작업 좌표계 m 상자 · 위 = 북). */ +export function drawRaster( + context: CanvasRenderingContext2D, + map: ScreenMap, + raster: { raster: PlanRaster; image: HTMLImageElement }, +): void { + const [x0, y0, x1, y1] = raster.raster.box_m; + const [left, top] = map.toScreen([x0, y1]); + const [right, bottom] = map.toScreen([x1, y0]); + context.drawImage(raster.image, left, top, right - left, bottom - top); +} diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_UI_Layers.ts b/D03_PlanAnalysis/D03_PlanAnalysis_UI_Layers.ts index cee75dca3..cb65ef2ea 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_UI_Layers.ts +++ b/D03_PlanAnalysis/D03_PlanAnalysis_UI_Layers.ts @@ -1,15 +1,16 @@ /* ============================================================================= * D03_PlanAnalysis_UI_Layers.ts - * D03 평면 바탕 겹 켜고 끔(PLAN 569-2 · 581-1 · 581-2 · 요구 395-2 · 440-8 · 442) — 지도 왼쪽 위 아이콘 한 줄 + * D03 평면 바탕 겹 켜고 끔(PLAN 569-2 · 581-1 · 581-2 · 584-1 · 요구 395-2 · 440-8 · 442 · 443) — 지도 왼쪽 위 아이콘 한 줄 * (공용 지도 도구줄 `createMapTools` · 넘치면 가로 밀기 571-1 · 오른쪽 평면 도구줄과 따로). * * 겹 = 바탕(위성 · 음영기복 · 도엽 등고선(계곡선 25 m 굵게 · 수치) · 급경사 색) · * 물 · 길(세류선 · 하천 · 물 면(실폭하천 · 저수지) · 기존 도로 · 기존 임도 짐작(소로 · 미분류)) · - * 구역(지적도 · 보호 · 규제 구역 여덟 종 · 산림입지도(토양형) · 계류 100 m · 도로 300 m 띠 · 건물 · 지명). + * 구역(지적도 · 보호 · 규제 구역 여덟 종 · 산림입지도(토양형) · 선 겹(단층 · 지질구조선 · 0차곡) · + * 계류 100 m · 도로 300 m 띠 · 건물 · 지명). * 처음 켬 = 위성 · 등고선 · 세류 · 고른 값 = 화면 취향(`aislo:pref:d03.layers`) 기억. * 받기 = 켤 때만(한 번 받으면 다시 안 받음) · 지형을 안 받았으면 단추를 막음 · 받기 실패 = 그 겹만 꺼짐 + 알림. * 그리기 차례(아래 → 위) = 위성(`scene.image`) → 음영기복 → 급경사 → 산림입지도 → 보호구역 → 띠 → 지적 → 물 면 → - * 건물 → 등고선 → 도로 → 임도 짐작 → 세류 → 지명 → 그은 선. + * 건물 → 등고선 → 도로 → 임도 짐작 → 세류 → 선 겹 → 지명 → 그은 선. * 겹 모두 = 바탕 캔버스(`base` · 582) · 등고선 = 줌 단계 색인 · 글 폭 기억. * 범례 = 켠 겹의 받은 범위 안 수 · 넓이(왼쪽 아래). * ========================================================================== */ @@ -29,7 +30,6 @@ import { import { createNormalizer, drawPreparedLabels, - drawPreparedLayer, lonLatToScreen, pickLevelStep, type GeoJsonCollection, @@ -42,9 +42,11 @@ import { import { fetchPlanBands, fetchPlanCadastral, + fetchPlanLines, fetchPlanRaster, fetchPlanSheetLayer, fetchPlanSoils, + fetchPlanZeroOrder, fetchPlanZones, type PlanBand, type PlanRaster, @@ -55,20 +57,32 @@ import { FOREST_ROAD_KINDS, GROUPS, LAYERS, + LINE_KEYS, + LINE_LAYERS, SOIL_COLORS, WATER_KEYS, ZONE_COLORS, + drawRaster, fillLayer, hectare, + loadImage, + stroke, + type LineLayerKey, type PlanLayerKey, } from "./D03_PlanAnalysis_UI_LayerDefs"; -import type { PlanMap, PlanMapLayer, ScreenMap } from "./D03_PlanAnalysis_UI_Map"; +import type { PlanMap, PlanMapLayer } from "./D03_PlanAnalysis_UI_Map"; import type { Vertex } from "./logic/D03_PlanAnalysis_Plan"; const PREF_KEY = "aislo:pref:d03.layers"; /** 등고선을 받아 올 그은 선 둘레(m). */ const CONTOUR_PAD_M = 600; -interface Prepared { +/** 선 겹 — 그릴 선 · 받은 범위 안 가닥 · 길이 m. */ +interface PreparedLines { + layer: PreparedLayer; + count: number; + length_m: number; +} +type Prepared = Record & { contour: PreparedLayer; /** 세류 · 하천(소하천 · 지방하천 · 연결선). */ stream: [PreparedLayer, PreparedLayer]; @@ -84,7 +98,7 @@ interface Prepared { bands: { stream: PreparedLayer; road: PreparedLayer; info: { stream: PlanBand; road: PlanBand } }; building: PreparedLayer; place: Array<{ lon: number; lat: number; name: string }>; -} +}; type DataKey = keyof Prepared; /** 겹 → 받는 자료(임도 짐작 = 도로 자료). */ @@ -121,14 +135,6 @@ function readPref(): Record { return on; } -const loadImage = (src: string): Promise => - new Promise((done, fail) => { - const image = new Image(); - image.onload = () => done(image); - image.onerror = () => fail(new Error("그림을 읽지 못했습니다")); - image.src = src; - }); - export function createPlanLayers(opts: { projectId: string; map: PlanMap; @@ -220,6 +226,20 @@ export function createPlanLayers(opts: { }, ]); } + const lines = LINE_KEYS.flatMap((key) => { + const got = on[key] && data[key]; + if (!got) return []; + const { label, color, dash } = LINE_LAYERS[key]; + return [ + { + label: `${label} ${got.count}가닥 ${(got.length_m / 1000).toFixed(2)} km`, + color: `var(${color})`, + kind: dash.length ? ("dash" as const) : ("line" as const), + off: !got.count, + }, + ]; + }); + if (lines.length) rows.push(lines); const tail: LegendItem[] = []; if (on.water && data.water) tail.push({ label: `물 면 ${data.water.features.length}곳`, color: "var(--draw-map-water)" }); @@ -377,6 +397,28 @@ export function createPlanLayers(opts: { })), summary: soils.summary, }; + } else if (key === "faults" || key === "lineaments" || key === "zeroOrder") { + const kind = LINE_LAYERS[key].kind; + // 0차곡 첫 셈 = 몇 초 — 늦으면 알림 + const slow = window.setTimeout( + () => showToast(`${LINE_LAYERS[key].label} 셈 중 — 처음 한 번만`, "info"), + 600, + ); + try { + const got = kind + ? await fetchPlanLines(projectId, kind).then((lines) => ({ + collection: lines, + ...lines.summary, + })) + : await fetchPlanZeroOrder(projectId); + data[key] = { + layer: prepare(got.collection), + count: got.count, + length_m: got.length_m, + }; + } finally { + window.clearTimeout(slow); + } } else if (key === "bands") { const info = await fetchPlanBands(projectId); data.bands = { @@ -396,30 +438,6 @@ export function createPlanLayers(opts: { } } - const stroke = ( - context: CanvasRenderingContext2D, - map: ScreenMap, - layer: PreparedLayer, - color: string, - width: number, - everyM?: number, - ): void => { - context.strokeStyle = color; - context.lineWidth = width; - drawPreparedLayer(context, layer, map.view, "dot", everyM); - }; - - const drawRaster = ( - context: CanvasRenderingContext2D, - map: ScreenMap, - raster: { raster: PlanRaster; image: HTMLImageElement }, - ): void => { - const [x0, y0, x1, y1] = raster.raster.box_m; - const [left, top] = map.toScreen([x0, y1]); - const [right, bottom] = map.toScreen([x1, y0]); - context.drawImage(raster.image, left, top, right - left, bottom - top); - }; - /** 등고 글 폭 기억(글꼴 하나 · 582). */ const labelWidths = new Map(); const layers: PlanMapLayer[] = [ @@ -560,6 +578,18 @@ export function createPlanLayers(opts: { stroke(context, map, data.stream[1], color, 2.4); }, }, + ...LINE_KEYS.map((key): PlanMapLayer => ({ + base: true, + draw: (context, map) => { + const got = on[key] && data[key]; + if (!got) return; + const { color, width, dash } = LINE_LAYERS[key]; + context.setLineDash(dash); + stroke(context, map, got.layer, "rgba(0, 0, 0, 0.55)", width + 1.6); + stroke(context, map, got.layer, themeColor(color, "#f43f5e"), width); + context.setLineDash([]); + }, + })), { base: true, draw: (context, map) => { @@ -629,6 +659,17 @@ export function createPlanLayers(opts: { band_road_m2: data.bands.info.road.area_m2, }), ...(data.hillshade && { hillshade_px: data.hillshade.image.naturalWidth }), + ...Object.fromEntries( + LINE_KEYS.flatMap((key) => { + const got = data[key]; + return got + ? [ + [key, got.count], + [`${key}_m`, got.length_m], + ] + : []; + }), + ), }); return { diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py new file mode 100644 index 000000000..10b87ca10 --- /dev/null +++ b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py @@ -0,0 +1,183 @@ +"""D 노선분석 지도 0차곡 겹 — 지형 추정 골짜기 머리(PLAN 584-1 · 요구 443). + +0차곡 = 뚜렷한 유로가 없는 곡두의 집수지형(`노선_배치.md:31` · 교본 `3. 임도망의 배치 절차.md:156`) +→ 사업 DEM(평활 안 한 `models_plan/dtm_sheet_<방법>.npz`)을 2 m 로 솎아 D8 흐름누적 · +위 집수 넓이가 `min_m2` 이상 · `max_m2` 미만인 칸 = 0차곡(그보다 크면 1차 계류) → +세선화 · 선 뽑기(B05 `_mask_to_polylines`) · 경위도 선 + 가닥 · 길이 합. +임계값 = 지식DB 수치 없음 → 정함(브레인) 0.1 ha · 1 ha · 바꿈 = 요청 `min_m2` · `max_m2`. +파생 저장본(열쇠 = 격자 지문 · 임계 · 칸). + +계산 자리 = ③ 서버 단독. +""" + +from __future__ import annotations + +import math +from pathlib import Path +from typing import Any + +import numpy as np + +from common_util.common_util_derived_store import ( + derived_key, + file_stamp, + read_derived, + write_derived, +) +from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD + +ZERO_ORDER_KIND = "plan_zero_order" +ZERO_ORDER_VERSION = 1 +ZERO_ORDER_MIN_M2 = 1_000.0 # 집수 시작(0.1 ha) +ZERO_ORDER_MAX_M2 = 10_000.0 # 이 넓이부터 1차 계류(1 ha) +ZERO_ORDER_CELL_M = 2.0 # 흐름누적 칸(1 m 격자는 솎음) +ZERO_ORDER_MIN_LEN_M = 10.0 # 세선화 갈림목 토막 거름 +SIMPLIFY_M = 1.0 + +_D8 = ((-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)) + + +def _shift(grid: np.ndarray, dr: int, dc: int, fill: Any) -> np.ndarray: + """칸마다 (dr, dc) 이웃 값 — 격자 밖 = `fill`.""" + rows, cols = grid.shape + pad = np.pad(grid, 1, constant_values=fill) + return pad[1 + dr : 1 + dr + rows, 1 + dc : 1 + dc + cols] + + +def d8_accumulation(z: np.ndarray, valid: np.ndarray) -> np.ndarray: + """D8 흐름누적(칸 수 · 자기 칸 포함 · 빈 칸 0). + + 웅덩이 메움(형태 재구성 · 가장자리 · 빈 칸 곁 = 출구) → 가장 가파른 아래 이웃 하나로 흐름 → + 평지(등고선 칸 · 메운 웅덩이)는 이미 흐르는 같은 높이 이웃으로 이어 감 → 위에서 아래로 더함. + """ + from skimage.morphology import reconstruction + + rows, cols = z.shape + height = np.where(valid, np.asarray(z, float), np.nan) + floor = np.nanmin(height) if valid.any() else 0.0 + edge = valid & ~np.all([_shift(valid, dr, dc, False) for dr, dc in _D8], axis=0) + ground = np.where(valid, height, floor) + seed = np.where(edge | ~valid, ground, ground.max()) + filled = np.where(valid, reconstruction(seed, ground, method="erosion"), np.nan) + index = np.arange(rows * cols).reshape(rows, cols) + best = np.zeros(z.shape) + receiver = np.full(z.shape, -1, np.int64) + with np.errstate(invalid="ignore"): + for dr, dc in _D8: + drop = (filled - _shift(filled, dr, dc, np.nan)) / math.hypot(dr, dc) + better = drop > best + best = np.where(better, drop, best) + receiver = np.where(better, index + dr * cols + dc, receiver) + done = (receiver >= 0) | edge | ~valid + while True: + grew = False + for dr, dc in _D8: + take = ~done & _shift(done & valid, dr, dc, False) + take &= _shift(filled, dr, dc, np.nan) == filled + if take.any(): + receiver[take] = index[take] + dr * cols + dc + done |= take + grew = True + if not grew: + break + down = receiver.ravel().tolist() + acc = valid.astype(float).ravel().tolist() + upstream = np.bincount(receiver[receiver >= 0], minlength=rows * cols).tolist() + stack = [i for i, n in enumerate(upstream) if n == 0] + while stack: + i = stack.pop() + j = down[i] + if j >= 0: + acc[j] += acc[i] + upstream[j] -= 1 + if upstream[j] == 0: + stack.append(j) + return np.asarray(acc).reshape(rows, cols) + + +def zero_order_mask( + acc: np.ndarray, valid: np.ndarray, cell_m: float, min_m2: float, max_m2: float +) -> np.ndarray: + """위 집수 넓이(칸 수 × 칸 넓이)가 [min_m2, max_m2) 인 칸.""" + area = acc * cell_m * cell_m + return valid & (area >= min_m2) & (area < max_m2) + + +def zero_order_lines( + z: np.ndarray, + valid: np.ndarray, + x: np.ndarray, + y: np.ndarray, + cell_m: float, + min_m2: float = ZERO_ORDER_MIN_M2, + max_m2: float = ZERO_ORDER_MAX_M2, +) -> list[list[tuple[float, float]]]: + """격자(행 = 남 → 북) → 0차곡 선(작업 좌표계 m 점 목록 · `ZERO_ORDER_MIN_LEN_M` 미만 거름).""" + from B05_Profile.B05_Profile_Engine_Skeleton import _mask_to_polylines + + mask = zero_order_mask(d8_accumulation(z, valid), valid, cell_m, min_m2, max_m2) + lines = [[(p[0], p[1]) for p in line["polyline"]] for line in _mask_to_polylines(mask, x, y, z)] + return [ + line + for line in lines + if sum(math.dist(a, b) for a, b in zip(line, line[1:])) >= ZERO_ORDER_MIN_LEN_M + ] + + +def zero_order_layer( + project_root: Path, + models_dir: Path, + epsg: int, + min_m2: float = ZERO_ORDER_MIN_M2, + max_m2: float = ZERO_ORDER_MAX_M2, +) -> dict[str, Any]: + """0차곡 겹 — 지문(격자 파일 · 임계 · 칸) 같으면 저장본 · 경위도 선 + 가닥 · 길이 합 m.""" + from shapely.geometry import LineString, mapping + from shapely.ops import transform + + from D03_PlanAnalysis.D03_PlanAnalysis_Layers import _transforms + + npz = models_dir / f"dtm_sheet_{SHEET_SURFACE_DEFAULT_METHOD}.npz" + stamp = file_stamp(npz) + if stamp is None: + raise FileNotFoundError("사업 지형면이 없습니다 — 지형 받기 뒤 보입니다.") + key = derived_key( + ZERO_ORDER_KIND, + ZERO_ORDER_VERSION, + npz.name, + stamp, + int(epsg), + float(min_m2), + float(max_m2), + ZERO_ORDER_CELL_M, + ZERO_ORDER_MIN_LEN_M, + SIMPLIFY_M, + ) + saved = read_derived(project_root, ZERO_ORDER_KIND, key) + if saved is not None: + return saved["value"] + with np.load(npz) as dem: + x, y = np.asarray(dem["x"], float), np.asarray(dem["y"], float) + cell = float(dem["resolution"][0]) if "resolution" in dem.files else float(x[1] - x[0]) + step = max(1, round(ZERO_ORDER_CELL_M / cell)) + z = np.asarray(dem["z"], float)[::step, ::step] + valid = np.asarray(dem["valid_mask"], bool)[::step, ::step] + lines = zero_order_lines(z, valid, x[::step], y[::step], cell * step, min_m2, max_m2) + _, to_ll = _transforms(epsg) + shapes = [LineString(line).simplify(SIMPLIFY_M) for line in lines] + value = { + "min_m2": float(min_m2), + "max_m2": float(max_m2), + "cell_m": cell * step, + "count": len(shapes), + "length_m": round(sum(s.length for s in shapes), 1), + "collection": { + "type": "FeatureCollection", + "features": [ + {"type": "Feature", "geometry": mapping(transform(to_ll, s)), "properties": {}} + for s in shapes + ], + }, + } + write_derived(project_root, ZERO_ORDER_KIND, key, value) + return value diff --git a/resources/tester/test_553_4_de_parts.py b/resources/tester/test_553_4_de_parts.py index 0256d4cb3..5cb67c81a 100644 --- a/resources/tester/test_553_4_de_parts.py +++ b/resources/tester/test_553_4_de_parts.py @@ -67,6 +67,9 @@ def test_d03_parts(): ("cadastral", "false"), ("zones", "false"), ("soils", "false"), + ("faults", "false"), + ("lineaments", "false"), + ("zeroOrder", "false"), ("bands", "false"), ("building", "false"), ("place", "false"), diff --git a/resources/tester/test_584_1_plan_geo_lines.py b/resources/tester/test_584_1_plan_geo_lines.py new file mode 100644 index 000000000..7a742d9e6 --- /dev/null +++ b/resources/tester/test_584_1_plan_geo_lines.py @@ -0,0 +1,176 @@ +"""D 노선분석 GIS 겹 더함 2(PLAN 584-1) — 단층 · 지질구조선 잘라 가닥 · 길이 · 실패 = 오류 · +덮으면 다시 안 받음 · 0차곡 흐름누적(웅덩이 · 평지 넘김) · 임계 · 선 뽑기 · 저장본 · 화면 겹 정의. + +./venv/Scripts/python.exe -m pytest resources/tester/test_584_1_plan_geo_lines.py -q +""" + +import json +import math +import re +from pathlib import Path + +import numpy as np +import pytest + +from common_util.common_util_crs import tm_to_lonlat +from D03_PlanAnalysis import D03_PlanAnalysis_Layers as layers +from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zo + +EPSG = 5186 +BOX = (230000.0, 450000.0, 232000.0, 451000.0) # 2 × 1 km +ROOT = Path(__file__).resolve().parents[2] + + +def _line(gid: str, points: list[tuple[float, float]]) -> dict: + coords = [list(tm_to_lonlat(x, y, EPSG)) for x, y in points] + return {"id": gid, "geometry": {"type": "LineString", "coordinates": coords}, "properties": {}} + + +# ── 선 겹 ── + + +def test_clip_lines_length() -> None: + found = [ + _line("a", [(229000.0, 450500.0), (231000.0, 450500.0)]), # 반만 안 = 1,000 m + _line("a", [(229000.0, 450500.0), (231000.0, 450500.0)]), # 겹침 + _line("b", [(231500.0, 450100.0), (231500.0, 450400.0)]), # 300 m + _line("c", [(240000.0, 460000.0), (240100.0, 460000.0)]), # 밖 + ] + features, summary = layers.clip_lines(found, BOX, EPSG) + assert summary["count"] == 2 == len(features) + assert summary["length_m"] == pytest.approx(1300.0, rel=0.005) + assert features[0]["properties"]["length_m"] == pytest.approx(1000.0, rel=0.005) + + +class FakeWfs: + def __init__(self, fail: bool = False) -> None: + self.fail, self.calls = fail, [] + + def __call__(self, typename, filename, bounds, output_dir) -> None: + self.calls.append(typename) + if self.fail: + return + feats = [_line("a", [(230100.0, 450100.0), (230100.0, 450600.0)])] + body = {"type": "FeatureCollection", "features": feats} + (Path(output_dir) / filename).write_text(json.dumps(body), "utf-8") + + +@pytest.mark.parametrize("kind", list(layers.LINE_KINDS)) +def test_lines_fetch_once(tmp_path, monkeypatch, kind) -> None: + fake = FakeWfs() + monkeypatch.setattr(layers, "request_vworld_wfs", fake) + label, typename = layers.LINE_KINDS[kind] + assert layers.ensure_lines(tmp_path, BOX, EPSG, kind)["fetched"] and fake.calls == [typename] + saved = json.loads((tmp_path / f"{label}_bounds.geojson").read_text("utf-8")) + assert saved["summary"]["count"] == 1 + assert saved["summary"]["length_m"] == pytest.approx(500.0, rel=0.005) + inner = (230050.0, 450050.0, 231000.0, 450900.0) + assert not layers.ensure_lines(tmp_path, inner, EPSG, kind)["fetched"] # 덮으면 다시 안 받음 + assert layers.ensure_lines(tmp_path, (229000.0, *BOX[1:]), EPSG, kind)["fetched"] + + +def test_lines_fail_raises(tmp_path, monkeypatch) -> None: + monkeypatch.setattr(layers, "request_vworld_wfs", FakeWfs(fail=True)) + with pytest.raises(RuntimeError): + layers.ensure_lines(tmp_path, BOX, EPSG, "faults") + assert not list(tmp_path.iterdir()) + + +def test_line_kinds_vworld_names() -> None: + assert layers.LINE_KINDS["faults"][1] == "lt_l_gimsfault" + assert layers.LINE_KINDS["lineaments"][1] == "lt_l_gimslinea" + + +# ── 0차곡 ── + + +def _valley(n: int = 60) -> np.ndarray: + """가운데(열 n/2)로 모이고 북(행 끝)에서 남으로 내려가는 V 골짜기 — 행 = 남 → 북.""" + r, c = np.mgrid[0:n, 0:n].astype(float) + return 0.5 * r + 0.8 * np.abs(c - n // 2) + + +def test_accumulation_valley_floor() -> None: + z = _valley() + acc = zo.d8_accumulation(z, np.ones(z.shape, bool)) + assert acc[0, 30] == acc.max() # 골 바닥 남쪽 끝이 가장 큼 + assert acc[0, 30] > 0.4 * z.size # 반쯤 넘게 모임 + + +def test_accumulation_crosses_pit_and_flat() -> None: + z = _valley() + z[30, 30] -= 5.0 # 웅덩이 + z[20:24, 20:41] = z[22, 30] # 평지 띠(등고선 칸) + acc = zo.d8_accumulation(z, np.ones(z.shape, bool)) + assert acc[0, 30] > 0.4 * z.size # 웅덩이 · 평지에서 안 끊김 + + +def test_zero_order_mask_threshold() -> None: + acc = np.array([[10.0, 249.0, 250.0, 2499.0, 2500.0]]) + mask = zo.zero_order_mask(acc, np.ones(acc.shape, bool), 2.0, 1000.0, 10000.0) # 칸 4 m² + assert mask.tolist() == [[False, False, True, True, False]] + + +def test_zero_order_lines_head_only() -> None: + n = 60 + z = _valley(n) + x = y = np.arange(n) * 2.0 + lines = zo.zero_order_lines(z, np.ones(z.shape, bool), x, y, 2.0, min_m2=400.0, max_m2=2000.0) + assert lines + # 골 바닥(x = 60 m) 위 · 집수 400 ~ 2,000 m² 사이만(골 위쪽 머리) + floor = sorted(p[1] for line in lines for p in line if abs(p[0] - 60.0) < 1e-6) + assert floor and floor[0] > 0.0 and floor[-1] < y[-1] # 남 끝(1차 계류) · 북 끝 빠짐 + for line in lines: + assert sum(math.dist(a, b) for a, b in zip(line, line[1:])) >= zo.ZERO_ORDER_MIN_LEN_M + + +def test_zero_order_layer_saved_and_threshold_key(tmp_path, monkeypatch) -> None: + n = 60 + models = tmp_path / "models" + models.mkdir() + np.savez( + models / f"dtm_sheet_{zo.SHEET_SURFACE_DEFAULT_METHOD}.npz", + x=230000.0 + np.arange(n) * 2.0, + y=450000.0 + np.arange(n) * 2.0, + z=_valley(n), + valid_mask=np.ones((n, n), bool), + resolution=np.array([2.0, 2.0]), + ) + calls = [] + real = zo.zero_order_lines + monkeypatch.setattr(zo, "zero_order_lines", lambda *a, **k: calls.append(1) or real(*a, **k)) + first = zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 2000.0) + assert first["count"] >= 1 and first["length_m"] > 0 and first["cell_m"] == 2.0 + assert first["collection"]["features"][0]["geometry"]["type"] == "LineString" + again = zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 2000.0) + assert (again["count"], again["length_m"]) == (first["count"], first["length_m"]) + assert len(calls) == 1 # 두 번째 = 저장본 + zo.zero_order_layer(tmp_path, models, EPSG, 400.0, 3000.0) # 임계 바뀜 = 다시 셈 + assert len(calls) == 2 + + +def test_zero_order_no_dem(tmp_path) -> None: + with pytest.raises(FileNotFoundError): + zo.zero_order_layer(tmp_path, tmp_path, EPSG) + + +# ── 화면 ── + + +def test_screen_defs_icons_colors() -> None: + defs = (ROOT / "D03_PlanAnalysis" / "D03_PlanAnalysis_UI_LayerDefs.ts").read_text("utf-8") + icons = (ROOT / "ui_template" / "ui_template_icons.ts").read_text("utf-8") + theme = (ROOT / "ui_template" / "ui_template_theme.css").read_text("utf-8") + for key, color in ( + ("faults", "fault"), + ("lineaments", "lineament"), + ("zeroOrder", "zero-order"), + ): + entry = re.search(rf'key: "{key}",\s*icon: "(\w+)",[^}}]*on: (\w+)', defs) + assert entry and entry.group(2) == "false", key # 처음 값 끔 + assert re.search(rf"^\s+{entry.group(1)}:", icons, re.M), key + assert f"--draw-map-{color}:" in theme and f'"--draw-map-{color}"' in defs + for kind in layers.LINE_KINDS: + assert f'kind: "{kind}"' in defs # 서버 열쇠와 같음 + page = (ROOT / "D03_PlanAnalysis" / "D03_PlanAnalysis_UI_Layers.ts").read_text("utf-8") + assert page.count("base: true,") == page.count("draw: (context, map) =>") # 겹 모두 바탕(582) diff --git a/ui_template/ui_template_icons.ts b/ui_template/ui_template_icons.ts index bd09b386c..0b66c4480 100644 --- a/ui_template/ui_template_icons.ts +++ b/ui_template/ui_template_icons.ts @@ -144,6 +144,10 @@ const PATHS = { "M4 6a2 2 0 0 1 2 -2h12a2 2 0 0 1 2 2v12a2 2 0 0 1 -2 2h-12a2 2 0 0 1 -2 -2z|M12 8l0 .01|M8 12l0 .01|M12 12l0 .01|M16 12l0 .01|M12 16l0 .01", // 581-2 — 지도 겹(산림입지도) leaf: "M5 21c.5 -4.5 2.5 -8 7 -10|M9 18c6.218 0 10.5 -3.288 11 -12v-2h-4.014c-9 0 -11.986 4 -12 9c0 1 0 3 2 5h3z", + // 584-1 — 지도 겹(단층 · 지질구조선 · 0차곡) + bolt: "M13 3l0 7l6 0l-8 11l0 -7l-6 0l8 -11", + lineDashed: "M5 12h2|M17 12h2|M11 12h2", + arrowsJoin: "M3 7h5l3.5 5h9.5|M3 17h5l3.495 -5|M18 15l3 -3l-3 -3", } as const; export type IconName = keyof typeof PATHS; diff --git a/ui_template/ui_template_theme.css b/ui_template/ui_template_theme.css index 64ae0a506..21b83afa4 100644 --- a/ui_template/ui_template_theme.css +++ b/ui_template/ui_template_theme.css @@ -575,6 +575,9 @@ html.is-bm-route .m01lab-src__body * { --draw-map-soil-4: #bef264; --draw-map-soil-5: #86efac; --draw-map-soil-6: #e879f9; + --draw-map-fault: #f43f5e; /* 단층 */ + --draw-map-lineament: #c4b5fd; /* 지질구조선(점선) */ + --draw-map-zero-order: #7dd3fc; /* 0차곡(지형 추정 골짜기 머리 · 점선) */ } :root[data-theme="dark"] {