#!/usr/bin/env python3
import json, math, os, itertools

D='/sessions/upbeat-quirky-clarke/mnt/outputs/sfmap'
P=json.load(open(D+'/base_parks2.json'))['parks']
LAT0=37.77
MX=111320.0*math.cos(math.radians(LAT0)); MY=110574.0

def shoe(r):
    s=0.0
    for i in range(len(r)-1):
        s+=r[i][0]*r[i+1][1]-r[i+1][0]*r[i][1]
    return s/2.0
def area_m2(pg):
    a=abs(shoe(pg['outer']))
    for h in pg['holes']: a-=abs(shoe(h))
    return a*MX*MY
def parea(p): return sum(area_m2(pg) for pg in p['polys'])
def verts(p): return sum(len(r) for pg in p['polys'] for r in [pg['outer']]+pg['holes'])

print('=== TOTALS ===')
npoly=sum(len(p['polys']) for p in P)
nhole=sum(len(pg['holes']) for p in P for pg in p['polys'])
nvert=sum(verts(p) for p in P)
sz=os.path.getsize(D+'/base_parks2.json')
print('parks %d  polygons %d  holes %d  vertices %d  size %d B (%.1f KB)'%(len(P),npoly,nhole,nvert,sz,sz/1024))

print('\n=== WINDING ===')
bad_o=bad_h=0
for p in P:
    for pg in p['polys']:
        if shoe(pg['outer'])<=0: bad_o+=1; print('  outer not CCW:',p['n'])
        for h in pg['holes']:
            if shoe(h)>=0: bad_h+=1; print('  hole not CW:',p['n'])
print('outer rings not CCW: %d   holes not CW: %d'%(bad_o,bad_h))

print('\n=== HOLES ===')
for p in P:
    for i,pg in enumerate(p['polys']):
        for h in pg['holes']:
            print('  %-34s poly %d  hole %3d pts  %.0f m2 (%.2f ac)'%(p['n'],i,len(h),abs(shoe(h))*MX*MY,abs(shoe(h))*MX*MY/4046.86))

print('\n=== GOLDEN GATE PARK SECTIONS ===')
G=[p for p in P if p['n'].startswith('Golden Gate Park - Section')]
G.sort(key=lambda p:p['n'])
tot=0.0; tv=0
print('  %-32s %5s %6s %8s %10s'%('section','polys','verts','acres_db','area_km2'))
for p in G:
    a=parea(p)/1e6; tot+=a; tv+=verts(p)
    print('  %-32s %5d %6d %8.1f %10.4f'%(p['n'],len(p['polys']),verts(p),p['acres'],a))
print('  %-32s %5s %6d %8.1f %10.4f'%('TOTAL','',tv,sum(p['acres'] for p in G),tot))
print('  total acres from geometry: %.1f'%(tot*1e6/4046.86))
xs=[q[0] for p in G for pg in p['polys'] for q in pg['outer']]
ys=[q[1] for p in G for pg in p['polys'] for q in pg['outer']]
print('  lon extent %.6f .. %.6f  => %.0f m E-W'%(min(xs),max(xs),(max(xs)-min(xs))*MX))
print('  lat extent %.6f .. %.6f  => %.0f m N-S'%(min(ys),max(ys),(max(ys)-min(ys))*MY))

print('\n=== GGP SECTION 1 POLYGONS (Panhandle check) ===')
s1=[p for p in G if p['n'].endswith('Section 1')][0]
for i,pg in enumerate(s1['polys']):
    xs=[q[0] for q in pg['outer']]; ys=[q[1] for q in pg['outer']]
    print('  poly %d: %3d pts, %2d holes, %.3f ac, lon %.4f..%.4f lat %.4f..%.4f'%(
        i,len(pg['outer']),len(pg['holes']),area_m2(pg)/4046.86,min(xs),max(xs),min(ys),max(ys)))

print('\n=== GGP SEAM CHECK (min vertex-to-vertex distance, metres) ===')
def pts(p): return [q for pg in p['polys'] for q in pg['outer']]
names=[p['n'].replace('Golden Gate Park - Section','S') for p in G]
res=[]
for (i,a),(j,b) in itertools.combinations(list(enumerate(G)),2):
    A=pts(a); B=pts(b); best=1e18; ncoin=0
    for x1,y1 in A:
        for x2,y2 in B:
            d=math.hypot((x1-x2)*MX,(y1-y2)*MY)
            if d<best: best=d
            if d<0.5: ncoin+=1
    res.append((best,names[i],names[j],ncoin))
res.sort()
print('  adjacent pairs (min dist < 50 m):')
for d,a,b,nc in res:
    if d<50: print('    %-4s <-> %-4s  min = %6.2f m   coincident vertices: %d'%(a,b,d,nc))
print('  non-adjacent pairs (min dist >= 50 m):')
for d,a,b,nc in res:
    if d>=50: print('    %-4s <-> %-4s  min = %8.1f m'%(a,b,d))

print('\n=== SPOT CHECKS ===')
for want in ['Mission Dolores Park','Glen Canyon Park','San Francisco Zoo']:
    for p in P:
        if p['n']==want:
            print('  %-22s acres_db %7.2f   geom %7.2f ac (%.4f km2)  polys %d holes %d verts %d'%(
                p['n'],p['acres'],parea(p)/4046.86,parea(p)/1e6,len(p['polys']),
                sum(len(pg['holes']) for pg in p['polys']),verts(p)))

print('\n=== EXTREMES / SANITY ===')
sm=sorted(P,key=lambda p:p['acres'])[:3]
for p in sm: print('  smallest: %-32s %.3f ac  geom %.3f ac'%(p['n'],p['acres'],parea(p)/4046.86))
allx=[q[0] for p in P for pg in p['polys'] for r in [pg['outer']]+pg['holes'] for q in r]
ally=[q[1] for p in P for pg in p['polys'] for r in [pg['outer']]+pg['holes'] for q in r]
print('  dataset bbox: lon %.6f..%.6f  lat %.6f..%.6f'%(min(allx),max(allx),min(ally),max(ally)))
print('  max decimals lon:',max(len(str(q).split('.')[-1]) for q in allx))
names=[p['n'] for p in P]
print('  duplicate names:',[n for n in set(names) if names.count(n)>1])

# compare with old file
old=json.load(open(D+'/base_parks.json'))['parks']
print('\n=== OLD base_parks.json COMPARISON ===')
ov=sum(len(r) for p in old for r in p['rings'])
print('  old: %d parks, %d rings, %d vertices'%(len(old),sum(len(p['rings']) for p in old),ov))
print('  new: %d parks, %d rings, %d vertices'%(len(P),npoly+nhole,nvert))
og={p['n']:sum(len(r) for r in p['rings']) for p in old if p['n'].startswith('Golden Gate Park')}
for p in G:
    print('  %-32s old %4d verts -> new %4d verts'%(p['n'],og.get(p['n'],0),verts(p)))