"""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.')