fix(ground): 측량 레이어로 운곡교 지형 추출 대상 제한

scripts/review-ungok-terrain.py

@@ -0,0 +1,102 @@
+"""Read-only terrain diagnostics: grid slopes, isolated peaks, and source TIN."""
+import json
+import sys
+from pathlib import Path
+import numpy as np
+from collections import Counter
+
+ROOT=Path(__file__).resolve().parents[1]
+sys.path.insert(0,str(ROOT/'pierzainer_ground'))
+from pierzainer_ground.survey_terrain import read_survey_block
+from pierzainer_ground.dxf_terrain import read_source_document, read_terrain_dxf
+from pierzainer_ground.dtm import DTM
+from pierzainer_ground.config import GroundConfig
+
+out=ROOT/'.verify/ungok-terrain-review-20261110'
+out.mkdir(parents=True,exist_ok=True)
+folder=ROOT/'pierzainer_ground/data/output/app_b614afc0-6829-4aef-96e3-3ff1365ff5db'
+data=json.loads((folder/'ground_model_3d_with_maps.json').read_text())
+settings=json.loads((folder/'generation_settings.json').read_text())
+g=data['terrain_grid']; z=np.array([np.nan if v is None else v for v in g['z']]).reshape(g['ny'],g['nx'])
+edges=[]
+for j in range(g['ny']):
+    for i in range(g['nx']):
+        for dj,di in [(0,1),(1,0)]:
+            if j+dj<g['ny'] and i+di<g['nx'] and np.isfinite(z[j,i]) and np.isfinite(z[j+dj,i+di]):
+                edges.append({'row':j,'col':i,'to_row':j+dj,'to_col':i+di,'z':float(z[j,i]),'to_z':float(z[j+dj,i+di]),'dz':float(abs(z[j+dj,i+di]-z[j,i]))})
+edges.sort(key=lambda e:-e['dz'])
+summary={'file':str(folder),'terrain':data['terrain'],'provenance':data.get('surface_provenance'),'grid':{k:v for k,v in g.items() if k!='z'},'valid_nodes':int(np.isfinite(z).sum()),'z_range':[float(np.nanmin(z)),float(np.nanmax(z))],'adjacent_edge_count':len(edges),'edges_over_5m':sum(e['dz']>5 for e in edges),'edges_over_8m':sum(e['dz']>8 for e in edges),'worst_edges':edges[:12]}
+source=next((folder/'cache').glob('*.dxf'))
+doc,_=read_source_document(source)
+cfg=GroundConfig()
+points,counts,_=read_survey_block(source,{'terrain_block':settings['terrain_block'],'terrain_crs':settings['terrain_crs'],'crs':settings['crs'],'cad_axes':settings['cad_axes'],'contour_step_m':cfg.contour_step,'dedup_m':cfg.dedup_xy,'z_min':settings['z_min'],'z_max':settings['z_max'],'exclude_zero':True,'prefer_recognized_layers':False},doc=doc)
+lo=np.array(g['origin'])-settings['max_edge'];hi=np.array(g['origin'])+[(g['nx']-1)*g['step'],(g['ny']-1)*g['step']]+settings['max_edge']
+points=points[np.all((points[:,:2]>=lo)&(points[:,:2]<=hi),axis=1)]
+dtm=DTM(points,settings['max_edge'])
+summary['source_counts']=counts
+summary['raw_height_layers']=dict(Counter(e.dxf.layer for e in doc.blocks.get(settings['terrain_block']) if e.dxftype() in ('POLYLINE','LWPOLYLINE')))
+# Counterfactual only: use the drawing's two contour layers, without replacing
+# the generated user model or changing source geometry.
+contours,contour_counts,_=read_survey_block(source,{'terrain_block':settings['terrain_block'],'terrain_crs':settings['terrain_crs'],'crs':settings['crs'],'cad_axes':settings['cad_axes'],'contour_step_m':cfg.contour_step,'dedup_m':cfg.dedup_xy,'z_min':settings['z_min'],'z_max':settings['z_max'],'exclude_zero':True,'layers':['CAA001','CAA002']},doc=doc)
+contours=contours[np.all((contours[:,:2]>=lo)&(contours[:,:2]<=hi),axis=1)]
+clean_dtm=DTM(contours,settings['max_edge'])
+gx,gy=np.meshgrid(g['origin'][0]+np.arange(g['nx'])*g['step'],g['origin'][1]+np.arange(g['ny'])*g['step'])
+clean=clean_dtm.z_many(np.column_stack([gx.ravel(),gy.ravel()])).reshape(z.shape)
+clean_edges=np.concatenate([np.abs(np.diff(clean,axis=0)).ravel(),np.abs(np.diff(clean,axis=1)).ravel()])
+clean_edges=clean_edges[np.isfinite(clean_edges)]
+summary['contour_only_diagnostic']={'layers':['CAA001','CAA002'],'counts':contour_counts,'valid_nodes':int(np.isfinite(clean).sum()),'z_range':[float(np.nanmin(clean)),float(np.nanmax(clean))],'max_adjacent_dz':float(clean_edges.max()),'edges_over_5m':int((clean_edges>5).sum()),'edges_over_8m':int((clean_edges>8).sum()),'not_applied':True}
+typed_cfg=GroundConfig(use_profile_points=False,z_min=settings['z_min'],z_max=settings['z_max'])
+typed_cfg.allowed_terrain_layers=typed_cfg.layers.spot_block+typed_cfg.layers.spot_text+typed_cfg.layers.contour
+typed=read_terrain_dxf(source,typed_cfg,document=doc)
+typed_local=typed.points[np.all((typed.points[:,:2]>=lo)&(typed.points[:,:2]<=hi),axis=1)]
+typed_dtm=DTM(typed_local,settings['max_edge'])
+typed_grid=typed_dtm.z_many(np.column_stack([gx.ravel(),gy.ravel()])).reshape(z.shape)
+typed_edges=np.concatenate([np.abs(np.diff(typed_grid,axis=0)).ravel(),np.abs(np.diff(typed_grid,axis=1)).ravel()])
+typed_edges=typed_edges[np.isfinite(typed_edges)]
+summary['typed_layer_diagnostic']={'layers':list(typed_cfg.allowed_terrain_layers),'counts':typed.stats['by_source'],'valid_nodes':int(np.isfinite(typed_grid).sum()),'z_range':[float(np.nanmin(typed_grid)),float(np.nanmax(typed_grid))],'max_adjacent_dz':float(typed_edges.max()),'edges_over_5m':int((typed_edges>5).sum()),'edges_over_8m':int((typed_edges>8).sum()),'covered_piers':sum(typed_dtm.z_at(p['X'],p['Y']) is not None for p in data['piers']),'not_applied':True}
+summary['pier_diagnostic']=[]
+for pier in data['piers']:
+    observed=clean_dtm.z_at(pier['X'],pier['Y'])
+    summary['pier_diagnostic'].append({'id':pier['id'],'generated':pier['ground']['center'],'contour_only':observed,'difference':None if observed is None else observed-pier['ground']['center']})
+for e in summary['worst_edges']:
+    triangles=[]
+    for row,col in [(e['row'],e['col']),(e['to_row'],e['to_col'])]:
+        xy=np.array(g['origin'])+np.array([col,row])*g['step']; k=int(dtm.tri.find_simplex(xy-dtm.origin))
+        triangles.append({'xy':xy.tolist(),'vertices':points[dtm.tri.simplices[k]].tolist() if k>=0 else [],'predicted':dtm.z_at(*xy)})
+    e['source_triangles']=triangles
+# Raw constant-height CAD entities retained in the selected TIN: look for
+# closed design outlines vs survey contours near the worst interpolation.
+raw=[]
+for entity in doc.blocks.get(settings['terrain_block']):
+    if entity.dxftype() not in ('LWPOLYLINE','POLYLINE'):continue
+    pts=np.array([tuple(p) for p in (entity.vertices_in_wcs() if entity.dxftype()=='LWPOLYLINE' else entity.points_in_wcs())])
+    if len(pts)<2 or np.ptp(pts[:,2])>1e-6:continue
+    raw.append({'handle':entity.dxf.handle,'layer':entity.dxf.layer,'z':float(pts[0,2]),'vertices':len(pts),'closed':bool(entity.closed if entity.dxftype()=='LWPOLYLINE' else entity.is_closed),'xy_min':pts[:,:2].min(0).tolist(),'xy_max':pts[:,:2].max(0).tolist(),'xyz':pts.tolist()})
+for e in summary['worst_edges'][:6]:
+    for tri in e['source_triangles']:
+        tri['entities']=[]
+        for north,east,height in tri['vertices']:
+            candidates=[]
+            for entity in raw:
+                if abs(entity['z']-height)>1e-3:continue
+                xy=np.array(entity['xyz'])[:,:2]; d=float(np.linalg.norm(xy-[east,north],axis=1).min())
+                candidates.append((d,entity))
+            candidates.sort(key=lambda r:r[0])
+            tri['entities'].append([{'distance_to_vertex':round(d,4),**{k:v for k,v in ent.items() if k!='xyz'}} for d,ent in candidates[:2]])
+(out/'report.json').write_text(json.dumps(summary,ensure_ascii=False,indent=2))
+print(json.dumps({k:v for k,v in summary.items() if k not in ('worst_edges','provenance')},ensure_ascii=False))
+for e in summary['worst_edges'][:3]:print(json.dumps(e,ensure_ascii=False))
+try:
+    import matplotlib
+    matplotlib.use('Agg')
+    import matplotlib.pyplot as plt
+    fig,axes=plt.subplots(1,3,figsize=(16,5),constrained_layout=True)
+    extent=[g['origin'][0],g['origin'][0]+(g['nx']-1)*g['step'],g['origin'][1],g['origin'][1]+(g['ny']-1)*g['step']]
+    for ax,values,title in [(axes[0],z,'Generated: all height-bearing layers'),(axes[1],typed_grid,'Diagnostic: survey spots + contour layers')]:
+        im=ax.imshow(values,origin='lower',extent=extent,vmin=20,vmax=47,aspect='auto',cmap='terrain');ax.set_title(title);ax.set_xlabel('Northing (m)');ax.set_ylabel('Easting (m)')
+    fig.colorbar(im,ax=axes[:2],label='EL (m)',shrink=.7)
+    row=summary['worst_edges'][0]['to_row'];axis=g['origin'][0]+np.arange(g['nx'])*g['step']
+    axes[2].plot(axis,z[row],label='Generated');axes[2].plot(axis,typed_grid[row],label='Survey spots + contours');axes[2].set_title(f'Transect E={g["origin"][1]+row*g["step"]:.3f} m');axes[2].set_xlabel('Northing (m)');axes[2].set_ylabel('EL (m)');axes[2].legend();axes[2].grid(alpha=.25)
+    fig.savefig(out/'comparison.png',dpi=150);plt.close(fig)
+except ImportError:
+    print('Plot dependency unavailable; numeric report preserved.')