diff --git a/B04_PreProcess/B04_PreProcess_Engine_Smooth.py b/B04_PreProcess/B04_PreProcess_Engine_Smooth.py index d06858f3f..03e5b60fa 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Smooth.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Smooth.py @@ -90,7 +90,7 @@ def smooth_dtm( y=y_coords, z=z_smooth, valid_mask=valid_grid, - resolution=np.array([preview_res], np.float32), + resolution=np.array([x_coords[1] - x_coords[0]], np.float32), # 넓힌 실제 칸(633) ) write_glb(preview_path, vertices_compact, faces_compact, context.bounds) return { diff --git a/B07_DesignDetail/B07_DesignDetail_Router_Support_Basin.py b/B07_DesignDetail/B07_DesignDetail_Router_Support_Basin.py index 82f38e41f..9062e89d2 100644 --- a/B07_DesignDetail/B07_DesignDetail_Router_Support_Basin.py +++ b/B07_DesignDetail/B07_DesignDetail_Router_Support_Basin.py @@ -25,6 +25,7 @@ from B07_DesignDetail.B07_DesignDetail_Engine_Cad_Plan import ( plan_drawing_label, ) from common_util.common_util_drainage_pipes import detail_basins_path +from common_util.common_util_grid_cell import grid_cell_m from config.config_system import ( BASIN_CONTOUR_INTERVAL_M, DRAWING_SCALE_BASIN, @@ -490,7 +491,7 @@ def _hillshade_for_box( grid_y = np.asarray(data["y"], dtype=np.float64) grid_z = np.asarray(data["z"], dtype=np.float64) valid = np.asarray(data["valid_mask"], dtype=bool) - resolution = float(np.asarray(data["resolution"]).reshape(-1)[0]) + resolution = grid_cell_m(data, grid_x) min_x, min_y, max_x, max_y = box columns = np.where((grid_x >= min_x) & (grid_x <= max_x))[0] diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py b/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py index a0d249f54..f1f027bb1 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py +++ b/D03_PlanAnalysis/D03_PlanAnalysis_Layers.py @@ -47,6 +47,7 @@ from common_util.common_util_derived_store import ( read_derived, write_derived, ) +from common_util.common_util_grid_cell import grid_cell_m from common_util.common_util_json import atomic_write_json from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD @@ -64,9 +65,9 @@ WFS_WORKERS = 4 SLOPE_KIND = "plan_slope" HILLSHADE_KIND = "plan_hillshade" -HILLSHADE_VERSION = 2 # 빈 칸 투명 +HILLSHADE_VERSION = 3 # 칸 = 좌표 간격(633) RASTER_KINDS = ("slope", "hillshade") -SLOPE_VERSION = 2 # 값 그림 + 1 % 칸 수(PLAN 603) +SLOPE_VERSION = 3 # 칸 = 좌표 간격(633) SLOPE_MAX_PCT = 250 # 칸 값 0~250 = 기울기 % 내림(250 = 250 % 이상) SLOPE_EMPTY = 255 # 빈 칸 MAX_SIDE_PX = 4096 # ponytail: 한 변 넘으면 솎음 · 넓은 사업은 조각 그림으로 @@ -273,7 +274,7 @@ def slope_layer( 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]) + cell = grid_cell_m(dem, x) value = build(dem["z"], np.asarray(dem["valid_mask"], bool), cell) half = cell / 2 box = (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half) diff --git a/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py index 885b91935..7a3b75c95 100644 --- a/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py +++ b/D03_PlanAnalysis/D03_PlanAnalysis_ZeroOrder.py @@ -41,6 +41,7 @@ from common_util.common_util_derived_store import ( read_derived, write_derived, ) +from common_util.common_util_grid_cell import grid_cell_m from common_util.common_util_json import atomic_write_json from config.config_system_terrain import SHEET_SURFACE_DEFAULT_METHOD @@ -420,7 +421,7 @@ def _project_box(npz: Path) -> tuple[float, float, float, float]: """사업 DEM 범위(칸 모서리) — 창고 격자에서 이 상자 안 칸을 씀.""" 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]) + cell = grid_cell_m(dem, x) half = cell / 2 return (x[0] - half, y[0] - half, x[-1] + half, y[-1] + half) @@ -500,6 +501,7 @@ def zero_order_layer( DEM_ZERO_ORDER_KIND = "plan_zero_order_dem" +DEM_ZERO_ORDER_VERSION = 4 # 칸 = 좌표 간격(633) · 도엽 길 판과 따로 def dem_zero_order_layer( @@ -519,7 +521,7 @@ def dem_zero_order_layer( raise FileNotFoundError("지형 격자가 없습니다.") key = derived_key( DEM_ZERO_ORDER_KIND, - ZERO_ORDER_VERSION, + DEM_ZERO_ORDER_VERSION, npz.name, stamp, int(epsg), @@ -534,7 +536,7 @@ def dem_zero_order_layer( return saved["value"] if saved is not None else None 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]) + cell = grid_cell_m(dem, x) 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] diff --git a/common_util/common_util_grid_cell.py b/common_util/common_util_grid_cell.py new file mode 100644 index 000000000..6172c61ef --- /dev/null +++ b/common_util/common_util_grid_cell.py @@ -0,0 +1,27 @@ +"""격자 칸 크기(m) — 적힌 칸(`resolution`)보다 좌표 간격을 믿음(633). + +B04 평활 격자는 꼭짓점이 많으면 칸을 넓히고 적힌 칸은 요청값 그대로라 둘이 어긋남. +둘이 맞으면(1 % 안) 적힌 값을 그대로 돌려줘 옛 결과와 같음. +""" + +from __future__ import annotations + +import logging +import math +from typing import Any + +import numpy as np + +logger = logging.getLogger(__name__) + + +def grid_cell_m(dem: Any, x: Any) -> float: + """npz 격자(`dem`)의 칸 — `x` = 열 좌표.""" + spacing = abs(float(x[1] - x[0])) + if "resolution" not in dem.files: + return spacing + written = float(np.asarray(dem["resolution"]).reshape(-1)[0]) + if math.isclose(written, spacing, rel_tol=0.01): + return written + logger.warning("격자 칸 %.3f m 로 적힘 · 좌표 간격 %.3f m — 좌표 간격을 씀", written, spacing) + return spacing diff --git a/resources/tester/test_633_grid_cell.py b/resources/tester/test_633_grid_cell.py new file mode 100644 index 000000000..74ca22587 --- /dev/null +++ b/resources/tester/test_633_grid_cell.py @@ -0,0 +1,110 @@ +"""격자 칸 = 좌표 간격(요구 538 · PLAN 633). + +- 적힌 칸(0.5) ≠ 좌표 간격(1.0) → 급경사 · 음영 · 0차곡이 좌표 간격으로 셈 +- 둘이 맞는 격자(D 수치지형도) → 적힌 값 그대로 · 결과 전과 같음 + +./venv/Scripts/python.exe -m pytest resources/tester/test_633_grid_cell.py -q +""" + +import logging + +import numpy as np + +from common_util.common_util_grid_cell import grid_cell_m +from D03_PlanAnalysis import D03_PlanAnalysis_Layers as layers +from D03_PlanAnalysis import D03_PlanAnalysis_ZeroOrder as zero_order + +EPSG = 5186 +X0, Y0 = 200000.0, 450000.0 + + +def _npz(path, cell: float, written: float | None, n: int = 60, tilt: float = 0.5): + x, y = X0 + np.arange(n) * cell, Y0 + np.arange(n) * cell + z = np.add.outer(np.zeros(n), x - X0) * tilt # 동쪽으로 tilt(= 50 %) + extra = {} if written is None else {"resolution": np.array([written], np.float32)} + path.parent.mkdir(parents=True, exist_ok=True) + np.savez(path, x=x, y=y, z=z, valid_mask=np.ones(z.shape, bool), **extra) + return path + + +def test_grid_cell_prefers_spacing(tmp_path, caplog) -> None: + with np.load(_npz(tmp_path / "a.npz", 2.0, 2.0)) as dem: + assert grid_cell_m(dem, dem["x"]) == 2.0 + with np.load(_npz(tmp_path / "b.npz", 1.5, None)) as dem: + assert grid_cell_m(dem, dem["x"]) == 1.5 + with caplog.at_level(logging.WARNING), np.load(_npz(tmp_path / "c.npz", 1.078, 0.5)) as dem: + assert abs(grid_cell_m(dem, dem["x"]) - 1.078) < 1e-9 + assert "좌표 간격" in caplog.text + + +def test_slope_uses_real_cell(tmp_path) -> None: + npz = _npz(tmp_path / "models" / "dtm_csf_smooth.npz", 1.0, 0.5) + value = layers.slope_layer(tmp_path, tmp_path / "models", EPSG, npz=npz) + counts = np.asarray(value["counts"]) + assert int(counts.argmax()) == 50 # 0.5 를 믿으면 100 + assert value["box_m"][0] == X0 - 0.5 + + +def test_sheet_slope_unchanged(tmp_path) -> None: + """적힌 칸 = 좌표 간격 → 옛 셈(적힌 칸)과 같은 그림.""" + npz = _npz(tmp_path / "models" / "dtm_sheet_x.npz", 2.0, 2.0) + value = layers.slope_layer(tmp_path, tmp_path / "models", EPSG, npz=npz) + with np.load(npz) as dem: + old = layers.slope_values_png(dem["z"], dem["valid_mask"], 2.0) + assert value["image"] == old["image"] and value["counts"] == old["counts"] + + +def test_hillshade_uses_real_cell(tmp_path) -> None: + npz = _npz(tmp_path / "models" / "dtm_csf_smooth.npz", 1.0, 0.5) + value = layers.slope_layer(tmp_path, tmp_path / "models", EPSG, kind="hillshade", npz=npz) + with np.load(npz) as dem: + real = layers.hillshade_layer_png(dem["z"], dem["valid_mask"], 1.0) + assert value["image"] == real["image"] + + +def test_zero_order_dem_uses_real_cell(tmp_path, monkeypatch) -> None: + npz = _npz(tmp_path / "models" / "dtm_csf_smooth.npz", 1.0, 0.5, n=40) + seen = {} + + def spy(z, valid, x, y, cell, *rest): + seen.update(shape=z.shape, dx=float(x[1] - x[0]), cell=cell) + return [] + + monkeypatch.setattr(zero_order, "zero_order_lines", spy) + value = zero_order.dem_zero_order_layer(tmp_path, npz, EPSG) + # 0.5 를 믿으면 4 칸 건너(10×10 · 실제 4 m 를 2 m 로 셈) + assert seen == {"shape": (20, 20), "dx": 2.0, "cell": 2.0} and value["cell_m"] == 2.0 + assert zero_order._project_box(npz)[0] == X0 - 0.5 + + +def test_b07_frame_hillshade_uses_real_cell(tmp_path) -> None: + """B07 도곽 음영 — 옛 평활 파일(0.5 적힘)도 좌표 간격으로.""" + from types import SimpleNamespace + + from B07_DesignDetail.B07_DesignDetail_Engine_Cad_Lidar import hillshade_png + from B07_DesignDetail.B07_DesignDetail_Router_Support_Basin import _hillshade_for_box + + npz = _npz(tmp_path / "B04_PreProcess" / "models" / "dtm_csf_smooth.npz", 1.0, 0.5) + context = SimpleNamespace(project_root=tmp_path, surface_params={}) + url, _ = _hillshade_for_box(context, (X0, Y0, X0 + 100, Y0 + 100)) + with np.load(npz) as dem: + assert url == hillshade_png(dem["z"], dem["valid_mask"], 1.0)[0] + + +def test_smooth_writes_real_cell(tmp_path) -> None: + """B04 평활 격자 — 칸을 넓히면(꼭짓점 상한) 넓힌 칸을 적음.""" + from types import SimpleNamespace + + from B04_PreProcess.B04_PreProcess_Engine_Smooth import smooth_dtm + + src = _npz(tmp_path / "dtm_csf.npz", 1.0, 1.0, n=20) + wide = X0 + np.arange(10) * 2.0 # 요청 0.5 → 상한 때문에 2.0 + context = SimpleNamespace( + config={"dtm_grid_resolution_meters": 1.0}, + preview_grid=lambda res: (wide, Y0 + np.arange(10) * 2.0, None), + contains_xy=lambda xs, ys: np.ones(np.shape(xs), bool), + bounds=np.array([[X0, X0 + 18.0], [Y0, Y0 + 18.0], [0.0, 10.0]]), + ) + smooth_dtm(context, tmp_path, "dtm_csf", src) + with np.load(tmp_path / "dtm_csf_smooth.npz") as dem: + assert float(dem["resolution"][0]) == 2.0