244 lines
14 KiB
Python
244 lines
14 KiB
Python
"""Resolve terrain and coordinate axes from the drawing and pier input alone."""
|
||
from __future__ import annotations
|
||
|
||
import hashlib
|
||
import json
|
||
import re
|
||
from pathlib import Path
|
||
|
||
import numpy as np
|
||
from ezdxf import recover
|
||
from pyproj import CRS
|
||
from scipy.spatial import QhullError
|
||
|
||
from .config import GroundConfig
|
||
from .dtm import DTM
|
||
from .dxf_terrain import TerrainData, convert_dwg_to_dxf, read_terrain_dxf, read_source_document
|
||
from .piers import pier_from_dict, read_piers
|
||
from .project_coordinates import workbook_piers
|
||
from .survey_terrain import read_survey_block
|
||
|
||
|
||
def pier_controls(path):
|
||
piers=([pier_from_dict(p) for p in workbook_piers(path)]
|
||
if Path(path).suffix.lower()=='.xlsx' else read_piers(path))
|
||
if not piers:raise ValueError('교각 Excel에서 IP1/IP3 좌표를 찾지 못했습니다.')
|
||
return piers
|
||
|
||
|
||
def referenced_profile(source,root):
|
||
"""Reuse established source conventions only for an identical original file."""
|
||
digest=hashlib.sha256(source.read_bytes()).hexdigest()
|
||
for folder in sorted((root/'data/input').iterdir()):
|
||
config=folder/'project.json'
|
||
if not config.is_file():continue
|
||
profile=json.loads(config.read_text(encoding='utf-8'))
|
||
if not profile.get('terrain_source_available'):continue
|
||
for original in folder.iterdir():
|
||
if original.suffix.lower()!=source.suffix.lower() or original.stat().st_size!=source.stat().st_size:continue
|
||
if hashlib.sha256(original.read_bytes()).hexdigest()==digest:
|
||
return profile,original.name
|
||
return None,None
|
||
|
||
|
||
def referenced_coordinate_choice(source,workbook,root):
|
||
"""Retain a user's earlier CRS choice for exactly the same two inputs."""
|
||
path=root/'data/input/coordinate_references.json'
|
||
if not path.is_file():return None
|
||
drawing_hash=hashlib.sha256(source.read_bytes()).hexdigest()
|
||
workbook_hash=hashlib.sha256(Path(workbook).read_bytes()).hexdigest()
|
||
for reference in json.loads(path.read_text(encoding='utf-8')).get('references',[]):
|
||
if reference.get('drawing_sha256')==drawing_hash and reference.get('pier_sha256')==workbook_hash:
|
||
crs=CRS.from_user_input(reference['crs'])
|
||
if crs.is_projected and all(abs(a.unit_conversion_factor-1)<1e-8 for a in crs.axis_info[:2]):return reference
|
||
return None
|
||
|
||
|
||
def explicit_crs(doc,workbook):
|
||
"""Only explicit EPSG/WKT metadata qualifies; coordinates alone do not."""
|
||
found=set()
|
||
texts=[]
|
||
for obj in doc.objects:
|
||
if obj.dxftype()=='GEODATA':
|
||
definition=getattr(obj,'coordinate_system_definition','')
|
||
if definition:
|
||
texts.append(definition)
|
||
try:
|
||
epsg=CRS.from_wkt(definition).to_epsg()
|
||
if epsg:found.add(epsg)
|
||
except Exception:pass
|
||
for block in doc.blocks:
|
||
for entity in block:
|
||
if entity.dxftype() in ('TEXT','MTEXT'):
|
||
text=entity.plain_text() if entity.dxftype()=='MTEXT' else entity.dxf.text
|
||
if 'epsg' in text.lower():texts.append(text)
|
||
if Path(workbook).suffix.lower()=='.xlsx':
|
||
import openpyxl
|
||
book=openpyxl.load_workbook(workbook,data_only=True,read_only=True)
|
||
try:
|
||
for sheet in book:
|
||
for row in sheet.iter_rows(values_only=True):
|
||
texts.extend(v for v in row if isinstance(v,str) and 'epsg' in v.lower())
|
||
finally:book.close()
|
||
for text in texts:
|
||
found.update(int(v) for v in re.findall(r'EPSG\s*[:=：]?\s*(\d{4,6})',text,re.I))
|
||
valid=[]
|
||
for code in sorted(found):
|
||
try:
|
||
crs=CRS.from_epsg(code)
|
||
if crs.is_projected and all(abs(a.unit_conversion_factor-1)<1e-8 for a in crs.axis_info[:2]):valid.append(crs.to_string())
|
||
except Exception:pass
|
||
return (valid[0] if len(valid)==1 else ''),valid
|
||
|
||
|
||
def read_auto_inputs(job,work,log):
|
||
from .drawing_job import ROOT
|
||
source=Path(job.drawing)
|
||
piers=pier_controls(job.piers)
|
||
centers=np.array([[p.X,p.Y] for p in piers])
|
||
if not np.isfinite(centers).all():raise ValueError('교각 Excel 좌표에 유효하지 않은 수가 있습니다.')
|
||
profile,profile_source=referenced_profile(source,ROOT)
|
||
coordinate_choice=referenced_coordinate_choice(source,job.piers,ROOT)
|
||
if source.suffix.lower()=='.dwg':
|
||
log('DWG를 DXF로 변환하는 중... (원본이 같으면 변환 캐시 재사용)')
|
||
work.mkdir(parents=True,exist_ok=True)
|
||
source=convert_dwg_to_dxf(source,work)
|
||
log('도면의 측량 블록·실제 표고점과 Excel 교각 중심을 대조하는 중...')
|
||
document,_=read_source_document(source,work)
|
||
crs,crs_candidates=explicit_crs(document,job.piers)
|
||
job.crs=crs;job.terrain_crs=crs
|
||
job.terrain_block='';job.layers='';job.use_3d_lines=False
|
||
job.z_min=-500.;job.z_max=9000.
|
||
status='explicit_input_metadata' if crs else 'unresolved_no_input_metadata'
|
||
if not crs and coordinate_choice:
|
||
job.crs=coordinate_choice['crs'];job.terrain_crs=coordinate_choice.get('terrain_crs',job.crs)
|
||
status=coordinate_choice['status']
|
||
log(f'동일한 두 입력의 기존 사용자 지정 좌표계 적용: {job.crs} (공인 기준점 검증 전)')
|
||
if profile:
|
||
job.crs=profile['crs'];job.terrain_crs=profile['terrain_crs']
|
||
job.z_min=0.;job.z_max=200.
|
||
status=profile.get('crs_status','existing_source_reference')
|
||
log(f'동일 원본의 기존 좌표 처리 규약 재사용: {profile_source}')
|
||
default_applied=False
|
||
if not job.crs:
|
||
defaults_path=ROOT/'data/input/map_defaults.json'
|
||
if defaults_path.is_file():
|
||
defaults=json.loads(defaults_path.read_text(encoding='utf-8'))
|
||
default_crs=CRS.from_user_input(defaults['crs'])
|
||
if not default_crs.is_projected or any(abs(a.unit_conversion_factor-1)>1e-8 for a in default_crs.axis_info[:2]):
|
||
raise ValueError('프로젝트 기본 지도 좌표계는 m 단위 평면좌표계여야 합니다.')
|
||
job.crs=default_crs.to_string();job.terrain_crs=job.crs
|
||
status=defaults['status'];default_applied=True
|
||
log(f'원본에 좌표계 표기 없음 → 프로젝트 기본 {job.crs}로 지도 포함 생성 (좌표계 가정)')
|
||
job.auto_detection={'crs_status':status,'crs_candidates':crs_candidates,
|
||
'crs_assumed':bool(profile) or bool(coordinate_choice and not crs) or default_applied,
|
||
'project_default_applied':default_applied,'source_reference':profile_source,
|
||
'method':'terrain_TIN_coverage_at_Excel_IP1_IP3_centers',
|
||
'arbitrary_offset_applied':False}
|
||
cfg=GroundConfig(use_profile_points=False,z_min=job.z_min,z_max=job.z_max,
|
||
tin_max_edge=job.max_edge,use_3d_lines=False)
|
||
names=[]
|
||
for block in document.blocks:
|
||
if block.name.lower() in ('*paper_space','*model_space'):
|
||
if block.name.lower()!='*model_space':continue
|
||
elif block.name.startswith('*'):continue
|
||
measured=0
|
||
for e in block:
|
||
if e.dxftype()=='POINT' and 1e-6<abs(e.dxf.location.z)<9000:measured+=1
|
||
elif e.dxftype()=='LWPOLYLINE' and 1e-6<abs(e.dxf.elevation)<9000:measured+=1
|
||
elif e.dxftype()=='POLYLINE':
|
||
vertices=list(e.points_in_wcs())
|
||
if len(vertices)>1 and 1e-6<abs(vertices[0].z)<9000 and max(p.z for p in vertices)-min(p.z for p in vertices)<1e-6:
|
||
measured+=1
|
||
if measured>=3:break
|
||
if measured>=3:names.append(block.name)
|
||
candidates=[]
|
||
lower=centers.min(0)-job.margin-job.max_edge
|
||
upper=centers.max(0)+job.margin+job.max_edge
|
||
for name in names:
|
||
for axes in ('YX','XY'):
|
||
if profile and name!=profile['terrain_block']:continue
|
||
try:
|
||
points,counts,contours=read_survey_block(source,{
|
||
'terrain_block':name,'terrain_crs':job.terrain_crs,'crs':job.crs,
|
||
'cad_axes':axes,'contour_step_m':cfg.contour_step,'dedup_m':cfg.dedup_xy,
|
||
'z_min':job.z_min,'z_max':job.z_max,'exclude_zero':True},doc=document)
|
||
local=points[np.all((points[:,:2]>=lower)&(points[:,:2]<=upper),axis=1)]
|
||
if len(local)<3:continue
|
||
dtm=DTM(local,job.max_edge); values=dtm.z_many(centers)
|
||
covered=int(np.isfinite(values).sum())
|
||
if not covered:continue
|
||
candidates.append({'name':name,'axes':axes,'points':points,'counts':counts,
|
||
'contours':contours,'covered':covered,'values':values,'nearby':len(local)})
|
||
log(f'지형 후보 {name} / CAD {axes}: 교각 {covered}/{len(piers)}개에서 실측 TIN 확인')
|
||
except (ValueError,KeyError,QhullError):continue
|
||
if not candidates:
|
||
# Preserve the established text/contour layer reader as a fallback.
|
||
for axes in ('YX','XY'):
|
||
cfg.cad_axes=axes
|
||
terrain=read_terrain_dxf(source,cfg,work,document=document)
|
||
points=terrain.points
|
||
if len(points)<3:continue
|
||
points=points[(np.abs(points[:,2])>1e-6)&np.all((points[:,:2]>=lower)&(points[:,:2]<=upper),axis=1)]
|
||
if len(points)<3:continue
|
||
try:values=DTM(points,job.max_edge).z_many(centers)
|
||
except (ValueError,QhullError):continue
|
||
covered=int(np.isfinite(values).sum())
|
||
if covered:candidates.append({'name':'*Model_Space','axes':axes,'points':points,
|
||
'counts':{},'contours':len(terrain.contours),'covered':covered,'values':values,'nearby':len(points)})
|
||
if not candidates:
|
||
# A plan may have been rotated/translated for its printed sheet. Recover
|
||
# that frame only when original pier rectangles identify every Excel
|
||
# midpoint and axis, with independent hold-out residuals <= 3 cm.
|
||
# Unknown numeric CAD lines are never used as terrain in this branch.
|
||
from .plan_registration import resolve_plan_registration
|
||
registration=(resolve_plan_registration(document,workbook_piers(job.piers))
|
||
if len(piers)>=3 and Path(job.piers).suffix.lower()=='.xlsx' else {'status':'no_verified_frame'})
|
||
if registration['status']=='registered':
|
||
cfg.cad_axes='YX'
|
||
layers=tuple(name for name in (*cfg.layers.spot_block,*cfg.layers.spot_text,*cfg.layers.contour) if name in document.layers)
|
||
if layers:
|
||
cfg.allowed_terrain_layers=layers
|
||
original=read_terrain_dxf(source,cfg,work,document=document)
|
||
transform=registration['solutions'][0]
|
||
rotation=np.array(transform['rotation']);translation=np.array(transform['translation'])
|
||
xy=original.points[:,[1,0]]@rotation+translation
|
||
points=np.column_stack([xy[:,1],xy[:,0],original.points[:,2]])
|
||
local=points[np.all((points[:,:2]>=lower)&(points[:,:2]<=upper),axis=1)]
|
||
if len(local)>=3:
|
||
values=DTM(local,job.max_edge).z_many(centers)
|
||
if np.isfinite(values).all():
|
||
job.cad_axes='YX';job.terrain_block='*Model_Space';job.layers=','.join(layers)
|
||
registration.update(source_sha256=hashlib.sha256(Path(job.drawing).read_bytes()).hexdigest(),
|
||
workbook_sha256=hashlib.sha256(Path(job.piers).read_bytes()).hexdigest(),terrain_layers=list(layers))
|
||
job.auto_detection.update(plan_registration=registration,covered_piers=len(piers),total_piers=len(piers),
|
||
terrain_block=job.terrain_block,cad_axes=job.cad_axes)
|
||
original.points=points;original.contours=[];original.bridge_lines=[];original.pier_labels=[]
|
||
for bore in original.boreholes:
|
||
east,north=np.array([bore['Y'],bore['X']])@rotation+translation
|
||
bore.update(X=float(north),Y=float(east))
|
||
original.stats.update(source_block='*Model_Space',block_selection='verified_plan_pier_geometry_registration',
|
||
plan_registration=registration)
|
||
log(f"원본 평면 교각 {len(piers)}개로 배치 좌표 복구: 최대 잔차 {max(transform['residual_m']):.9f} m. 표고 변경 없음.")
|
||
return original,cfg
|
||
if not candidates:
|
||
raise ValueError('두 입력 파일에서 교각 위치와 겹치는 실제 지형을 찾지 못했습니다. 원본 도면에 측량 지형/XREF가 포함되어 있는지 확인이 필요합니다. 임의 좌표 보정이나 지형 생성은 하지 않았습니다.')
|
||
candidates.sort(key=lambda c:(c['covered'],c['nearby']),reverse=True)
|
||
best=candidates[0]
|
||
for other in candidates[1:]:
|
||
if other['covered']!=best['covered']:continue
|
||
shared=np.isfinite(other['values'])&np.isfinite(best['values'])
|
||
if np.any(shared) and np.max(np.abs(other['values'][shared]-best['values'][shared]))>.05:
|
||
raise ValueError('교각 위치에서 서로 다른 표고를 가진 지형 원본이 중복됩니다. 원도면의 최신 측량 지형을 확인해야 합니다: '+best['name']+', '+other['name'])
|
||
job.terrain_block=best['name'];job.cad_axes=best['axes'];cfg.cad_axes=best['axes']
|
||
job.auto_detection.update(terrain_block=best['name'],cad_axes=best['axes'],
|
||
covered_piers=best['covered'],total_piers=len(piers))
|
||
log(f"자동 선택: {best['name']} / CAD {best['axes']} / 교각 {best['covered']}/{len(piers)}개")
|
||
if not job.crs:log('원본에 EPSG/지도 좌표계 표기가 없습니다. 지반면은 Excel과 일치하는 측량 좌표로 생성합니다.')
|
||
terrain=TerrainData(best['points'],np.full(len(best['points']),'survey'),stats={
|
||
'file':source.name,'points_used':len(best['points']),'removed':0,
|
||
'source_block':best['name'],'source_counts':best['counts'],'contours':best['contours'],
|
||
'block_selection':'automatic_Excel_coordinate_TIN_coverage',
|
||
'source_metadata_repair':getattr(document,'_ground_source_repair',None)})
|
||
return terrain,cfg
|