City-Lifeline-Cockpit/tools/city-massing/audit-user-sample-quality.py

import json, math, hashlib
from pathlib import Path
from PIL import Image, ImageDraw, ImageFont
ROOT=Path('tools/city-massing/city/cad-bim'); OUT=Path('.scratch/cockpit-shell-panel-layout/verify/artifacts/t35/user-samples'); OUT.mkdir(parents=True,exist_ok=True)
SOURCE=Path('D:/GIS/GIS/yuanlingPic3'); Z=19; N=2**Z; SIZE=1280; WINDOW_TILES=5
TARGETS=[('positive-fine',[110.40957689672715,28.4430329383677]),('negative-estate',[110.40570023361066,28.434788752462875])]
def xy(lon,lat):return ((lon+180)/360*N*256,(1-math.asinh(math.tan(math.radians(lat)))/math.pi)/2*N*256)
def dist(a,b):return math.hypot((a[0]-b[0])*math.cos(math.radians((a[1]+b[1])/2)),a[1]-b[1])*111320
def coords(g):
 out=[]
 def rec(v):
  if isinstance(v,list) and len(v)>=2 and isinstance(v[0],(int,float)) and isinstance(v[1],(int,float)):out.append(v);return
  if isinstance(v,list):
   for x in v:rec(x)
 if g and 'coordinates' in g:rec(g['coordinates'])
 return out
def center(g):
 p=coords(g)
 return [sum(x[0] for x in p)/len(p),sum(x[1] for x in p)/len(p)] if p else None
def line(draw,pts,color,width=2):draw.line([(round(x),round(y)) for x,y in pts],fill=color,width=width,joint='curve')
def pix(c,x0,y0):p=xy(*c);return[p[0]-x0*256,p[1]-y0*256]
features=json.loads(Path('tools/city-massing/city/cad-buildings.geojson').read_text())['features'];plan=json.loads((ROOT/'plan.json').read_text());planby={b['id']:b for b in plan['buildings']};seated=json.loads((ROOT/'seated-plan.json').read_text());seatby={b['id']:b for b in seated['buildings']};exclude={b['id']:b for b in seated['excluded']};roots=['pilot','batch-urban8000-dedup','batch-urban8000b','batch-urban8000c','batch-urban-final','batch-urban-tail','batch-unfitted-fallback'];active=set()
for n in roots:
 p=ROOT/f'lod-{n}-roofed2/manifest.json'
 if p.exists():
  for c in json.loads(p.read_text())['cells']:active.update(c['buildingIds'])
font=ImageFont.truetype('C:/Windows/Fonts/arial.ttf',13);reports=[]
for name,target in TARGETS:
 near=[]
 for i,f in enumerate(features):
  c=center(f['geometry'])
  if c is None:continue
  d=dist(c,target)
  if d>190:continue
  id=f'cad-{i}';b=planby.get(id,{});s=seatby.get(id,{});ex=exclude.get(id,{});tier=b.get('tier',ex.get('tier','not-planned'))
  near.append({'id':id,'distanceM':round(d,1),'center':c,'surveyAreaSqm':b.get('areaSqm'),'rectFillRatio':b.get('rectFillRatio'),'activeSourceId':id in active,'planned':id in planby,'seated':id in seatby,'tier':tier,'kind':b.get('kind'),'archetype':b.get('archetype'),'estimatedFloors':b.get('floors'),'estimatedHeightM':s.get('heightM'),'terrainReliefM':s.get('terrainRelief',ex.get('relief')),'exclusionReason':ex.get('reason')})
 near.sort(key=lambda x:x['distanceM'])
 for b in near:b['class']=('active' if b['activeSourceId'] else 'seated-not-active' if b['seated'] else 'planned-not-seated' if b['planned'] else 'not-in-seated-plan')
 counts={k:sum(b['class']==k for b in near) for k in ['active','seated-not-active','planned-not-seated','not-in-seated-plan']}
 p=xy(*target);cx=math.floor(p[0]/256);cy=math.floor(p[1]/256);x0=cx-WINDOW_TILES//2;y0=cy-WINDOW_TILES//2
 image=Image.new('RGB',(SIZE,SIZE),'#fff');tiles=[]
 for ty in range(WINDOW_TILES):
  for tx in range(WINDOW_TILES):
   path=SOURCE/str(Z)/str(x0+tx)/f'{N-1-(y0+ty)}.png'
   if not path.exists():continue
   raw=path.read_bytes()
   with Image.open(path) as tile:image.paste(tile.convert('RGB').resize((SIZE//WINDOW_TILES,SIZE//WINDOW_TILES)),(tx*SIZE//WINDOW_TILES,ty*SIZE//WINDOW_TILES))
   tiles.append({'path':str(path),'sha256':hashlib.sha256(raw).hexdigest(),'bytes':len(raw)})
 scale=(SIZE/WINDOW_TILES)/256
 def point(c):q=pix(c,x0,y0);return[q[0]*scale,q[1]*scale]
 draw=ImageDraw.Draw(image); colors={'active':'#24e777','seated-not-active':'#24bfff','planned-not-seated':'#ff4c56','not-in-seated-plan':'#ffc24a'}
 for b in near:
  f=features[int(b['id'][4:])];pts=[point(c) for c in coords(f['geometry'])]
  if len(pts)<3:continue
  line(draw,pts+[pts[0]],colors[b['class']],2 if b['class']=='active' else 1)
  if b['class']=='seated-not-active' and (b['estimatedFloors'] or 0)>=4 and (b['surveyAreaSqm'] or 0)>=150:
   q=point(b['center']);draw.text((q[0],q[1]),b['id'],fill='#fff4aa',font=font,stroke_width=1,stroke_fill='#111')
 q=point(target);draw.ellipse((q[0]-8,q[1]-8,q[0]+8,q[1]+8),outline='#f000ff',width=3);draw.line((q[0]-14,q[1],q[0]+14,q[1]),fill='#f000ff',width=2);draw.line((q[0],q[1]-14,q[0],q[1]+14),fill='#f000ff',width=2)
 out=OUT/f'{name}-coverage.jpg';image.save(out,quality=96)
 reports.append({'name':name,'target':target,'radiusM':190,'image':str(out),'tileSources':tiles,'counts':counts,'buildings':near,'highRiseNotActive':[b for b in near if b['class']=='seated-not-active' and (b['estimatedFloors'] or 0)>=4 and (b['surveyAreaSqm'] or 0)>=150]})
result={'targets':reports,'legend':{'active':'green','seated-not-active':'cyan','planned-not-seated':'red','not-in-seated-plan':'amber','userCoordinate':'magenta cross'},'areaSource':'plan.json areaSqm; not source GeoJSON properties','heightSource':'estimated from footprint area; not measured','modelChanged':False,'browserValidated':False}
(OUT/'report-v2.json').write_text(json.dumps(result,indent=2));print(json.dumps({'passed':True,'targets':[{'name':r['name'],'counts':r['counts'],'highRiseNotActive':[(b['id'],b['surveyAreaSqm'],b['estimatedFloors'],b['terrainReliefM']) for b in r['highRiseNotActive']]} for r in reports],'modelChanged':False}))