379 lines
16 KiB
Python
379 lines
16 KiB
Python
"""Survey-coordinate map layers embedded in ground_model.json."""
|
|
from __future__ import annotations
|
|
|
|
import base64
|
|
import io
|
|
import math
|
|
import urllib.error
|
|
import urllib.request
|
|
|
|
from PIL import Image, UnidentifiedImageError
|
|
from pyproj import Transformer
|
|
|
|
from .pipeline import GroundModel
|
|
|
|
_WORLD = 20037508.342789244
|
|
_TILE_SIZE = 256
|
|
_MAX_IMAGE_SIDE = 2048
|
|
_MAX_TILES = 64
|
|
_USER_AGENT = "PierzainerGround/0.1 (ground_model map export)"
|
|
_SURVEY_TO_WEB = Transformer.from_crs("EPSG:5186", "EPSG:3857", always_xy=True)
|
|
_SURVEY_TO_WGS84 = Transformer.from_crs("EPSG:5186", "EPSG:4326", always_xy=True)
|
|
_WEB_TO_WGS84 = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)
|
|
_WEB_TO_SURVEY = Transformer.from_crs("EPSG:3857", "EPSG:5186", always_xy=True)
|
|
|
|
_PROVIDERS = {
|
|
"ground": {
|
|
"name": "2D 지도",
|
|
"format": "PNG",
|
|
"mime_type": "image/png",
|
|
"attribution": "© OpenStreetMap contributors",
|
|
"url": "https://tile.openstreetmap.org/{z}/{x}/{y}.png",
|
|
},
|
|
"satellite": {
|
|
"name": "위성사진",
|
|
"format": "JPEG",
|
|
"mime_type": "image/jpeg",
|
|
"attribution": "Source: Esri, Vantor, Earthstar Geographics, and the GIS User Community",
|
|
"url": "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
|
|
},
|
|
}
|
|
|
|
|
|
def _tile_extent_for_survey_bounds(
|
|
northing0: float,
|
|
easting0: float,
|
|
northing1: float,
|
|
easting1: float,
|
|
) -> tuple[int, int, int, int, int]:
|
|
projected = [
|
|
_SURVEY_TO_WEB.transform(easting, northing)
|
|
for northing, easting in (
|
|
(northing0, easting0),
|
|
(northing0, easting1),
|
|
(northing1, easting0),
|
|
(northing1, easting1),
|
|
)
|
|
]
|
|
if not all(math.isfinite(value) for point in projected for value in point):
|
|
raise ValueError("교각 좌표를 지도 좌표로 변환하지 못했습니다.")
|
|
|
|
west = min(point[0] for point in projected)
|
|
east = max(point[0] for point in projected)
|
|
south = min(point[1] for point in projected)
|
|
north = max(point[1] for point in projected)
|
|
|
|
for zoom in range(18, 9, -1):
|
|
scale = _TILE_SIZE * 2**zoom / (2 * _WORLD)
|
|
left = math.floor((west + _WORLD) * scale)
|
|
right = math.ceil((east + _WORLD) * scale)
|
|
top = math.floor((_WORLD - north) * scale)
|
|
bottom = math.ceil((_WORLD - south) * scale)
|
|
width, height = right - left, bottom - top
|
|
tile_columns = (right - 1) // _TILE_SIZE - left // _TILE_SIZE + 1
|
|
tile_rows = (bottom - 1) // _TILE_SIZE - top // _TILE_SIZE + 1
|
|
if max(width, height) <= _MAX_IMAGE_SIDE and tile_columns * tile_rows <= _MAX_TILES:
|
|
return zoom, left, top, right, bottom
|
|
raise ValueError("지도 범위가 너무 큽니다. 교각 좌표 범위를 줄여 다시 시도하세요.")
|
|
|
|
|
|
def _tile_extent(grid: dict) -> tuple[int, int, int, int, int]:
|
|
northing0, easting0 = grid["origin"]
|
|
northing1 = northing0 + (grid["nx"] - 1) * grid["step"]
|
|
easting1 = easting0 + (grid["ny"] - 1) * grid["step"]
|
|
return _tile_extent_for_survey_bounds(northing0, easting0, northing1, easting1)
|
|
|
|
|
|
def _read_tile(provider: str, zoom: int, x: int, y: int) -> Image.Image:
|
|
config = _PROVIDERS[provider]
|
|
url = config["url"].format(z=zoom, x=x, y=y)
|
|
request = urllib.request.Request(url, headers={"User-Agent": _USER_AGENT})
|
|
try:
|
|
with urllib.request.urlopen(request, timeout=30) as response:
|
|
content = response.read(10_000_001)
|
|
except (urllib.error.URLError, TimeoutError, OSError) as exc:
|
|
raise RuntimeError(f"{config['name']} 타일을 가져오지 못했습니다: {url}") from exc
|
|
if len(content) > 10_000_000:
|
|
raise RuntimeError(f"{config['name']} 타일 응답이 너무 큽니다: {url}")
|
|
|
|
try:
|
|
with Image.open(io.BytesIO(content)) as image:
|
|
image.load()
|
|
if image.size != (_TILE_SIZE, _TILE_SIZE):
|
|
raise RuntimeError(f"{config['name']} 타일 크기가 올바르지 않습니다: {image.size}")
|
|
return image.convert("RGB")
|
|
except (UnidentifiedImageError, OSError) as exc:
|
|
raise RuntimeError(f"{config['name']} 타일 이미지가 올바르지 않습니다: {url}") from exc
|
|
|
|
|
|
def _mosaic(provider: str, extent: tuple[int, int, int, int, int]) -> tuple[bytes, tuple[float, float, float, float]]:
|
|
zoom, left, top, right, bottom = extent
|
|
first_x, last_x = left // _TILE_SIZE, (right - 1) // _TILE_SIZE
|
|
first_y, last_y = top // _TILE_SIZE, (bottom - 1) // _TILE_SIZE
|
|
if first_x < 0 or first_y < 0 or last_x >= 2**zoom or last_y >= 2**zoom:
|
|
raise ValueError("지도 좌표가 Web Mercator 범위를 벗어났습니다.")
|
|
|
|
mosaic = Image.new("RGB", ((last_x - first_x + 1) * _TILE_SIZE, (last_y - first_y + 1) * _TILE_SIZE))
|
|
for tile_y in range(first_y, last_y + 1):
|
|
for tile_x in range(first_x, last_x + 1):
|
|
tile = _read_tile(provider, zoom, tile_x, tile_y)
|
|
mosaic.paste(tile, ((tile_x - first_x) * _TILE_SIZE, (tile_y - first_y) * _TILE_SIZE))
|
|
|
|
image = mosaic.crop((left - first_x * _TILE_SIZE, top - first_y * _TILE_SIZE,
|
|
right - first_x * _TILE_SIZE, bottom - first_y * _TILE_SIZE))
|
|
output = io.BytesIO()
|
|
image.save(output, format=_PROVIDERS[provider]["format"], optimize=True, quality=88)
|
|
|
|
scale = _TILE_SIZE * 2**zoom / (2 * _WORLD)
|
|
bounds = (
|
|
left / scale - _WORLD,
|
|
_WORLD - bottom / scale,
|
|
right / scale - _WORLD,
|
|
_WORLD - top / scale,
|
|
)
|
|
return output.getvalue(), bounds
|
|
|
|
|
|
def _image_layer(provider: str, extent: tuple[int, int, int, int, int]) -> tuple[dict, tuple[float, float, float, float]]:
|
|
config = _PROVIDERS[provider]
|
|
image, bounds = _mosaic(provider, extent)
|
|
west, south, east, north = bounds
|
|
southwest = _WEB_TO_WGS84.transform(west, south)
|
|
northeast = _WEB_TO_WGS84.transform(east, north)
|
|
layer_bounds = [[southwest[1], southwest[0]], [northeast[1], northeast[0]]]
|
|
west, south, east, north = bounds
|
|
survey_corners = [
|
|
_WEB_TO_SURVEY.transform(west, south),
|
|
_WEB_TO_SURVEY.transform(east, south),
|
|
_WEB_TO_SURVEY.transform(west, north),
|
|
_WEB_TO_SURVEY.transform(east, north),
|
|
]
|
|
layer = {
|
|
"id": provider,
|
|
"name": config["name"],
|
|
"type": "image",
|
|
"bounds": layer_bounds,
|
|
"survey_corners": [[northing, easting] for easting, northing in survey_corners],
|
|
"mime_type": config["mime_type"],
|
|
"image": f"data:{config['mime_type']};base64,{base64.b64encode(image).decode('ascii')}",
|
|
"attribution": config["attribution"],
|
|
}
|
|
return layer, bounds
|
|
|
|
|
|
def _features(model: GroundModel) -> dict:
|
|
features = []
|
|
|
|
def coordinate(northing: float, easting: float) -> list[float]:
|
|
longitude, latitude = _SURVEY_TO_WGS84.transform(easting, northing)
|
|
if not math.isfinite(longitude) or not math.isfinite(latitude):
|
|
raise ValueError("교각 또는 시추공 좌표를 위경도로 변환하지 못했습니다.")
|
|
return [longitude, latitude]
|
|
|
|
for result in model.results:
|
|
pier = result.pier
|
|
features.append({
|
|
"type": "Feature",
|
|
"geometry": {"type": "Point", "coordinates": coordinate(pier.X, pier.Y)},
|
|
"properties": {"id": pier.id, "kind": "pier", "ground_el": result.ground_center},
|
|
})
|
|
corners = pier.corners()
|
|
if corners:
|
|
ring = [coordinate(northing, easting) for northing, easting in corners]
|
|
ring.append(ring[0])
|
|
features.append({
|
|
"type": "Feature",
|
|
"geometry": {"type": "Polygon", "coordinates": [ring]},
|
|
"properties": {"id": pier.id, "kind": "footing"},
|
|
})
|
|
|
|
for borehole in model.boreholes:
|
|
if borehole.X is None or borehole.Y is None:
|
|
continue
|
|
features.append({
|
|
"type": "Feature",
|
|
"geometry": {"type": "Point", "coordinates": coordinate(borehole.X, borehole.Y)},
|
|
"properties": {"id": borehole.id, "kind": "borehole"},
|
|
})
|
|
return {"type": "FeatureCollection", "features": features}
|
|
|
|
|
|
def _features_from_json(data: dict) -> dict:
|
|
features = []
|
|
|
|
def coordinate(northing: float, easting: float) -> list[float]:
|
|
longitude, latitude = _SURVEY_TO_WGS84.transform(easting, northing)
|
|
if not math.isfinite(longitude) or not math.isfinite(latitude):
|
|
raise ValueError("교각 또는 시추공 좌표를 위경도로 변환하지 못했습니다.")
|
|
return [longitude, latitude]
|
|
|
|
for pier in data.get("piers", []):
|
|
x, y = pier.get("X"), pier.get("Y")
|
|
if x is None or y is None:
|
|
raise ValueError(f"{pier.get('id', '교각')}: 지도에 필요한 X/Y 좌표가 없습니다.")
|
|
features.append({
|
|
"type": "Feature",
|
|
"geometry": {"type": "Point", "coordinates": coordinate(float(x), float(y))},
|
|
"properties": {
|
|
"id": pier["id"],
|
|
"kind": "pier",
|
|
"ground_el": (pier.get("ground") or {}).get("center"),
|
|
},
|
|
})
|
|
|
|
source_features = (data.get("visualization") or {}).get("features") or {}
|
|
existing_footings = [
|
|
feature for feature in source_features.get("features", [])
|
|
if feature.get("properties", {}).get("kind") == "footing"
|
|
]
|
|
features.extend(existing_footings)
|
|
|
|
for borehole in data.get("boreholes_dxf", []):
|
|
x, y = borehole.get("X"), borehole.get("Y")
|
|
if x is None or y is None:
|
|
continue
|
|
features.append({
|
|
"type": "Feature",
|
|
"geometry": {"type": "Point", "coordinates": coordinate(float(x), float(y))},
|
|
"properties": {"id": borehole["id"], "kind": "borehole"},
|
|
})
|
|
return {"type": "FeatureCollection", "features": features}
|
|
|
|
|
|
def build_visualization_from_json(data: dict, margin_m: float = 50.0) -> dict:
|
|
"""Create map layers for an existing ground_model.json that may not contain a terrain grid."""
|
|
if not math.isfinite(margin_m) or margin_m < 0:
|
|
raise ValueError("지도 여백은 0 이상이어야 합니다.")
|
|
|
|
points = [
|
|
(float(item["X"]), float(item["Y"]))
|
|
for item in [*data.get("piers", []), *data.get("boreholes_dxf", [])]
|
|
if item.get("X") is not None and item.get("Y") is not None
|
|
]
|
|
if not points:
|
|
raise ValueError("ground_model.json에 지도 범위를 계산할 X/Y 좌표가 없습니다.")
|
|
|
|
northings = [point[0] for point in points]
|
|
eastings = [point[1] for point in points]
|
|
if not all(math.isfinite(value) for value in [*northings, *eastings]):
|
|
raise ValueError("ground_model.json의 지도 좌표에 유효하지 않은 값이 있습니다.")
|
|
extent = _tile_extent_for_survey_bounds(
|
|
min(northings) - margin_m,
|
|
min(eastings) - margin_m,
|
|
max(northings) + margin_m,
|
|
max(eastings) + margin_m,
|
|
)
|
|
features = _features_from_json(data)
|
|
existing_layers = (data.get("visualization") or {}).get("layers") or []
|
|
plan_layer = next((layer for layer in existing_layers if layer.get("id") == "plan"), None)
|
|
if plan_layer is None:
|
|
plan_layer = next(
|
|
(
|
|
layer for layer in existing_layers
|
|
if layer.get("id") == "ground" and layer.get("name") != "2D 지도"
|
|
),
|
|
None,
|
|
)
|
|
if plan_layer is not None:
|
|
plan_layer = {
|
|
**plan_layer,
|
|
"id": "plan",
|
|
"name": "종평면도 지반·등고선",
|
|
}
|
|
missing_boreholes = sum(
|
|
borehole.get("X") is None or borehole.get("Y") is None
|
|
for borehole in data.get("boreholes_dxf", [])
|
|
)
|
|
return _visualization(extent, features, missing_boreholes, plan_layer)
|
|
|
|
|
|
def add_maps_to_ground_json(data: dict, margin_m: float = 50.0) -> dict:
|
|
"""Return a map-enabled copy of an existing ground_model.json without changing calculations."""
|
|
updated = dict(data)
|
|
updated.setdefault("schema", 2)
|
|
updated.setdefault("crs", "EPSG:5186")
|
|
updated.setdefault("crs_assumed", True)
|
|
updated["visualization"] = build_visualization_from_json(data, margin_m)
|
|
return updated
|
|
|
|
|
|
def _visualization(
|
|
extent: tuple[int, int, int, int, int],
|
|
features: dict,
|
|
missing_boreholes: int = 0,
|
|
plan_layer: dict | None = None,
|
|
) -> dict:
|
|
street, map_bounds = _image_layer("ground", extent)
|
|
satellite, satellite_bounds = _image_layer("satellite", extent)
|
|
if map_bounds != satellite_bounds:
|
|
raise RuntimeError("2D 지도와 위성사진의 좌표 범위가 일치하지 않습니다.")
|
|
|
|
west, south, east, north = map_bounds
|
|
southwest = _WEB_TO_WGS84.transform(west, south)
|
|
northeast = _WEB_TO_WGS84.transform(east, north)
|
|
bounds = [[southwest[1], southwest[0]], [northeast[1], northeast[0]]]
|
|
notes = [
|
|
"2D 지도와 위성사진을 이미지 데이터로 포함하여 원본 이미지 파일 없이 불러올 수 있습니다.",
|
|
"지도 위치는 측량 X(북)·Y(동) 좌표에 맞췄으며 지반고 Z/EL은 사용하지 않았습니다.",
|
|
"입력 좌표계는 EPSG:5186으로 가정했으므로 원도면 좌표계 표기를 확인하세요.",
|
|
]
|
|
if missing_boreholes:
|
|
notes.append(f"좌표가 없는 시추공 {missing_boreholes}개는 지도 표식에서 제외했습니다.")
|
|
return {
|
|
"schema_version": 1,
|
|
"type": "georeferenced_raster_layers",
|
|
"source_crs": "EPSG:5186",
|
|
"source_crs_status": "assumed",
|
|
"display_crs": "EPSG:3857",
|
|
"bounds_crs": "EPSG:4326",
|
|
"bounds_order": "southwest_northeast_lat_lon",
|
|
"georeferencing": {
|
|
"bounds": bounds,
|
|
"crs": "EPSG:5186 (가정)",
|
|
"controlPoints": 4,
|
|
"pixelResidualMax": 0.0,
|
|
"pixelResidualRms": 0.0,
|
|
"source": "ground_model.json 교각·시추공 측량 좌표",
|
|
"sourceDate": None,
|
|
"method": "측량 X(북)·Y(동) 좌표로 범위를 계산. 지반고 Z/EL은 지도 정합에 사용하지 않음.",
|
|
},
|
|
"layers": [satellite, street, *([plan_layer] if plan_layer else [])],
|
|
"features": features,
|
|
"view": {
|
|
"base": "satellite",
|
|
"ground_visible": plan_layer is not None,
|
|
"ground_opacity": 0.55,
|
|
"piers": True,
|
|
"bores": False,
|
|
"footing": True,
|
|
"alignment": False,
|
|
},
|
|
"notes": notes,
|
|
}
|
|
|
|
|
|
def build_visualization(model: GroundModel, terrain_grid: dict) -> dict:
|
|
"""Build portable street and satellite images using survey X/Y, not ground elevations."""
|
|
missing_boreholes = sum(borehole.X is None or borehole.Y is None for borehole in model.boreholes)
|
|
plan_layer = None
|
|
existing_layers = (getattr(model, "visualization", None) or {}).get("layers", [])
|
|
for layer in existing_layers:
|
|
if layer.get("id") == "plan":
|
|
plan_layer = layer
|
|
break
|
|
if layer.get("id") == "ground" and layer.get("name") != "2D 지도":
|
|
plan_layer = {**layer, "id": "plan", "name": "종평면도 지반·등고선"}
|
|
return _visualization(_tile_extent(terrain_grid), _features(model), missing_boreholes, plan_layer)
|
|
|
|
|
|
def validate_embedded_image(layer: dict) -> None:
|
|
"""Raise if an embedded image URI is malformed; used by tests and downstream checks."""
|
|
header, separator, encoded = layer["image"].partition(",")
|
|
if not separator or header != f"data:{layer['mime_type']};base64":
|
|
raise ValueError(f"{layer.get('name', '지도')} 이미지 데이터 URI가 올바르지 않습니다.")
|
|
try:
|
|
with Image.open(io.BytesIO(base64.b64decode(encoded, validate=True))) as image:
|
|
image.verify()
|
|
except (ValueError, UnidentifiedImageError, OSError) as exc:
|
|
raise ValueError(f"{layer.get('name', '지도')} 이미지 데이터가 손상되었습니다.") from exc
|