bjkim1
21ebb16fb7
- SPipe_server_2026-10-30.zip 내용을 풀어서 등록 (압축 파일은 올리지 않음) - README에 데모 v2 서버 업데이트 순서 추가, 저장소에 없는 문서 목차 정리 - 패치내역을 docs/로 이동, _worklog.md 추가 - 임시 파일(backend/temp.*)과 API 키가 든 파일은 제외
234 lines
10 KiB
Python
234 lines
10 KiB
Python
"""GIS layer operations: export in any CRS/format, geometry ops, Voronoi, grids.
|
|
|
|
Layers travel as GeoJSON FeatureCollections in WGS84 (what the map holds).
|
|
Metric operations run in a planar CRS given by the caller (the model's CRS),
|
|
and results come back in WGS84.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import io
|
|
import os
|
|
import tempfile
|
|
import zipfile
|
|
|
|
import numpy as np
|
|
import shapely
|
|
from pyogrio import raw as ogr_raw
|
|
from pyproj import CRS, Transformer
|
|
from shapely.geometry import mapping, shape
|
|
|
|
from .batch import BatchError
|
|
|
|
EXPORT_FORMATS = {
|
|
"shp": ("ESRI Shapefile", ".shp"),
|
|
"gpkg": ("GPKG", ".gpkg"),
|
|
"dxf": ("DXF", ".dxf"),
|
|
"geojson": ("GeoJSON", ".geojson"),
|
|
"kml": ("KML", ".kml"),
|
|
}
|
|
WGS84 = CRS.from_epsg(4326)
|
|
MAX_FEATURES = 200_000
|
|
MAX_GRID_CELLS = 200_000
|
|
|
|
|
|
def _crs(user: str | None, what: str = "좌표계") -> CRS:
|
|
if not user:
|
|
raise BatchError(f"{what}가 필요합니다.")
|
|
try:
|
|
return CRS.from_user_input(user)
|
|
except Exception as e:
|
|
raise BatchError(f"{what}를 읽을 수 없습니다: {str(user)[:60]}") from e
|
|
|
|
|
|
def _planar(user: str | None) -> CRS:
|
|
c = _crs(user, "평면 좌표계")
|
|
if c.is_geographic:
|
|
raise BatchError("거리 계산에는 평면 좌표계(예: EPSG:5186)가 필요합니다.")
|
|
return c
|
|
|
|
|
|
def _project(geoms: np.ndarray, src: CRS, dst: CRS) -> np.ndarray:
|
|
if src.equals(dst):
|
|
return geoms
|
|
tr = Transformer.from_crs(src, dst, always_xy=True)
|
|
return shapely.transform(geoms, lambda xy: np.column_stack(tr.transform(xy[:, 0], xy[:, 1])))
|
|
|
|
|
|
def _read_fc(fc: dict) -> tuple[np.ndarray, list[dict]]:
|
|
feats = (fc or {}).get("features") or []
|
|
if len(feats) > MAX_FEATURES:
|
|
raise BatchError(f"도형이 {MAX_FEATURES:,}개를 넘습니다.")
|
|
geoms, props = [], []
|
|
for f in feats:
|
|
g = f.get("geometry")
|
|
try:
|
|
geoms.append(shape(g) if g else None)
|
|
except Exception:
|
|
geoms.append(None)
|
|
props.append(dict(f.get("properties") or {}))
|
|
return np.array(geoms, dtype=object), props
|
|
|
|
|
|
def _to_fc(geoms, props: list[dict]) -> dict:
|
|
out = []
|
|
for g, p in zip(geoms, props):
|
|
if g is None or g.is_empty:
|
|
continue
|
|
out.append({"type": "Feature", "geometry": mapping(g), "properties": p})
|
|
return {"type": "FeatureCollection", "features": out}
|
|
|
|
|
|
# ---------------------------------------------------------------- export
|
|
|
|
|
|
def _field_arrays(props: list[dict], fmt: str):
|
|
names: list[str] = []
|
|
for p in props:
|
|
for k in p:
|
|
if k not in names:
|
|
names.append(k)
|
|
if fmt == "shp":
|
|
# dBASE: 10-char unique names
|
|
seen, short = set(), []
|
|
for n in names:
|
|
base = str(n)[:10] or "F"
|
|
s, k = base, 1
|
|
while s in seen:
|
|
k += 1
|
|
s = f"{base[:10 - len(str(k))]}{k}"
|
|
seen.add(s)
|
|
short.append(s)
|
|
else:
|
|
short = [str(n) for n in names]
|
|
arrays = []
|
|
for n in names:
|
|
vals = [p.get(n) for p in props]
|
|
nums = [v for v in vals if v is not None and v != ""]
|
|
if nums and all(isinstance(v, (int, float)) and not isinstance(v, bool) for v in nums):
|
|
arrays.append(np.array([float(v) if v is not None and v != "" else np.nan for v in vals], dtype="float64"))
|
|
else:
|
|
arrays.append(np.array(["" if v is None else str(v) for v in vals], dtype=object))
|
|
return short, arrays
|
|
|
|
|
|
def write_vector(fc: dict, target_crs: str | None, fmt: str, name: str) -> tuple[bytes, str, str]:
|
|
"""FeatureCollection (WGS84) -> file bytes in `fmt` and `target_crs` (shapefile zipped)."""
|
|
if fmt not in EXPORT_FORMATS:
|
|
raise BatchError(f"지원하지 않는 형식입니다: {fmt}")
|
|
driver, ext = EXPORT_FORMATS[fmt]
|
|
geoms, props = _read_fc(fc)
|
|
keep = [i for i, g in enumerate(geoms) if g is not None and not g.is_empty]
|
|
if not keep:
|
|
raise BatchError("내보낼 도형이 없습니다.")
|
|
geoms, props = geoms[keep], [props[i] for i in keep]
|
|
dst = _crs(target_crs) if target_crs else WGS84
|
|
if fmt == "kml":
|
|
dst = WGS84
|
|
geoms = _project(geoms, WGS84, dst)
|
|
safe = "".join(ch for ch in (name or "layer") if ch.isalnum() or ch in "-_ ").strip() or "layer"
|
|
kinds = {g.geom_type.replace("Multi", "") for g in geoms}
|
|
if fmt == "shp" and len(kinds) > 1:
|
|
raise BatchError("Shapefile은 한 가지 도형(점/선/면)만 담을 수 있습니다. GeoPackage로 내보내세요.")
|
|
geometry_type = "Unknown" if len(kinds) > 1 else ("Multi" + kinds.pop() if any(g.geom_type.startswith("Multi") for g in geoms) else geoms[0].geom_type)
|
|
# Promote singles to multi when the layer mixes them (pyogrio needs one type).
|
|
if geometry_type.startswith("Multi"):
|
|
geoms = np.array([g if g.geom_type.startswith("Multi") else getattr(shapely, geometry_type)([g]) for g in geoms], dtype=object)
|
|
fields, arrays = _field_arrays(props, fmt) if fmt != "dxf" else ([], [])
|
|
with tempfile.TemporaryDirectory() as tmp:
|
|
path = os.path.join(tmp, safe + ext)
|
|
kwargs = {"encoding": "UTF-8"} if fmt == "shp" else {}
|
|
ogr_raw.write(path, shapely.to_wkb(geoms), arrays, fields, driver=driver, geometry_type=geometry_type,
|
|
crs=dst.to_wkt() if fmt != "dxf" else None, **kwargs)
|
|
if fmt == "shp":
|
|
buf = io.BytesIO()
|
|
with zipfile.ZipFile(buf, "w", zipfile.ZIP_DEFLATED) as z:
|
|
for fn in sorted(os.listdir(tmp)):
|
|
z.write(os.path.join(tmp, fn), fn)
|
|
return buf.getvalue(), safe + ".zip", "application/zip"
|
|
with open(path, "rb") as f:
|
|
data = f.read()
|
|
media = {"gpkg": "application/geopackage+sqlite3", "dxf": "application/dxf", "geojson": "application/geo+json", "kml": "application/vnd.google-earth.kml+xml"}[fmt]
|
|
return data, safe + ext, media
|
|
|
|
|
|
# ---------------------------------------------------------------- operations
|
|
|
|
|
|
def reproject_fc(fc: dict, src: str, dst: str) -> dict:
|
|
"""Coordinates that were read as `src` but are really `dst`: fix the layer (WGS84 in/out)."""
|
|
a, b = _crs(src, "원래 좌표계"), _crs(dst, "새 좌표계")
|
|
geoms, props = _read_fc(fc)
|
|
raw = _project(geoms, WGS84, a) # back to the numbers in the file
|
|
return _to_fc(_project(raw, b, WGS84), props)
|
|
|
|
|
|
def vector_op(op: str, params: dict) -> dict:
|
|
planar = _planar(params.get("crs"))
|
|
if op == "voronoi":
|
|
pts = np.asarray(params.get("points") or [], dtype="float64").reshape(-1, 2) # in planar CRS
|
|
ids = [str(i) for i in params.get("ids") or []]
|
|
if len(pts) < 3 or len(ids) != len(pts):
|
|
raise BatchError("점이 3개 이상 필요합니다.")
|
|
margin = float(params.get("margin") or 0)
|
|
hull = shapely.MultiPoint(pts).convex_hull
|
|
clip = hull.buffer(margin if margin > 0 else max(1.0, 0.05 * np.sqrt(max(hull.area, 1.0))))
|
|
cells = shapely.get_parts(shapely.voronoi_polygons(shapely.MultiPoint(pts), extend_to=clip.envelope))
|
|
tree = shapely.STRtree(shapely.points(pts))
|
|
out_g, out_p = [], []
|
|
for c in cells:
|
|
hit = tree.query(c, predicate="contains")
|
|
if len(hit) == 0:
|
|
continue
|
|
g = c.intersection(clip)
|
|
if g.is_empty:
|
|
continue
|
|
k = int(hit[0])
|
|
out_g.append(g)
|
|
out_p.append({"node": ids[k], "area_m2": round(float(g.area), 2)})
|
|
return _to_fc(_project(np.array(out_g, dtype=object), planar, WGS84), out_p)
|
|
|
|
if op == "grid":
|
|
b = params.get("bounds") # planar [minx, miny, maxx, maxy]
|
|
cell = float(params.get("cell") or 0)
|
|
if not b or cell <= 0:
|
|
raise BatchError("범위와 격자 크기가 필요합니다.")
|
|
x0, y0, x1, y1 = map(float, b)
|
|
nx, ny = int(np.ceil((x1 - x0) / cell)), int(np.ceil((y1 - y0) / cell))
|
|
if nx * ny > MAX_GRID_CELLS:
|
|
raise BatchError(f"격자가 {nx * ny:,}칸입니다. 격자 크기를 키우세요 (최대 {MAX_GRID_CELLS:,}).")
|
|
gs, ps = [], []
|
|
for j in range(ny):
|
|
for i in range(nx):
|
|
gs.append(shapely.box(x0 + i * cell, y0 + j * cell, x0 + (i + 1) * cell, y0 + (j + 1) * cell))
|
|
ps.append({"col": i, "row": j, "cell": f"{i}_{j}"})
|
|
return _to_fc(_project(np.array(gs, dtype=object), planar, WGS84), ps)
|
|
|
|
geoms, props = _read_fc(params.get("layer") or {})
|
|
if not len(geoms):
|
|
raise BatchError("대상 레이어가 비어 있습니다.")
|
|
g = _project(geoms, WGS84, planar)
|
|
if op == "buffer":
|
|
d = float(params.get("distance") or 0)
|
|
if d == 0:
|
|
raise BatchError("버퍼 거리를 입력하세요.")
|
|
res = shapely.buffer(g, d)
|
|
return _to_fc(_project(res, planar, WGS84), props)
|
|
if op == "dissolve":
|
|
field = params.get("field")
|
|
groups: dict = {}
|
|
for geom, p in zip(g, props):
|
|
if geom is None:
|
|
continue
|
|
groups.setdefault(p.get(field) if field else None, []).append(geom)
|
|
out_g = [shapely.union_all(v) for v in groups.values()]
|
|
out_p = [{field: k, "count": len(v)} if field else {"count": len(v)} for k, v in groups.items()]
|
|
return _to_fc(_project(np.array(out_g, dtype=object), planar, WGS84), out_p)
|
|
if op in ("intersection", "difference", "clip"):
|
|
other_g, _ = _read_fc(params.get("other") or {})
|
|
if not len(other_g):
|
|
raise BatchError("두 번째 레이어가 비어 있습니다.")
|
|
mask = shapely.union_all(_project(other_g, WGS84, planar))
|
|
res = shapely.difference(g, mask) if op == "difference" else shapely.intersection(g, mask)
|
|
return _to_fc(_project(res, planar, WGS84), props)
|
|
raise BatchError(f"알 수 없는 연산입니다: {op}")
|