import json, math, sys
sys.path.insert(0,'.')
import fed_src as S

LO,HI,LA,HA = -122.52,-122.35,37.70,37.84

def sarea(r):
    s=0.0
    for i in range(len(r)-1):
        x1,y1=r[i]; x2,y2=r[i+1]
        s+=x1*y2-x2*y1
    return s/2.0

def close(r):
    r=[list(map(float,p)) for p in r]
    if r[0]!=r[-1]: r.append(r[0][:])
    return r

def rnd(r): return [[round(p[0],6),round(p[1],6)] for p in r]

# Sutherland-Hodgman clip to bbox
def clip(ring):
    def cl(pts, inside, inter):
        out=[]
        n=len(pts)
        for i in range(n):
            a=pts[i]; b=pts[(i+1)%n]
            ia,ib=inside(a),inside(b)
            if ia: out.append(a)
            if ia!=ib: out.append(inter(a,b))
        return out
    pts=ring[:-1]
    edges=[(lambda p: p[0]>=LO, lambda a,b:[LO,a[1]+(b[1]-a[1])*(LO-a[0])/(b[0]-a[0])]),
           (lambda p: p[0]<=HI, lambda a,b:[HI,a[1]+(b[1]-a[1])*(HI-a[0])/(b[0]-a[0])]),
           (lambda p: p[1]>=LA, lambda a,b:[a[0]+(b[0]-a[0])*(LA-a[1])/(b[1]-a[1]),LA]),
           (lambda p: p[1]<=HA, lambda a,b:[a[0]+(b[0]-a[0])*(HA-a[1])/(b[1]-a[1]),HA])]
    for ins,itr in edges:
        if not pts: return None
        pts=cl(pts,ins,itr)
    if len(pts)<3: return None
    pts.append(pts[0][:])
    return pts

def mkpolys(rings):
    """rings: list of rings (first is an outer). classify by sign vs first."""
    rings=[close(r) for r in rings]
    s0=sarea(rings[0])
    polys=[]; cur=None
    for r in rings:
        if cur is None or (sarea(r)*s0)>0:
            cur={"outer":r,"holes":[]}; polys.append(cur)
        else:
            cur["holes"].append(r)
    return polys

def norm(polys):
    out=[]
    for p in polys:
        o=clip(p["outer"])
        if o is None: continue
        if sarea(o)<0: o=o[::-1]          # outer CCW
        hs=[]
        for h in p["holes"]:
            hc=clip(h)
            if hc is None: continue
            if sarea(hc)>0: hc=hc[::-1]   # hole CW
            hs.append(rnd(hc))
        out.append({"outer":rnd(o),"holes":hs})
    return out

parks=[]
for oid,(n,ag,ac,acres) in S.META.items():
    if oid in S.R:
        polys=mkpolys(S.R[oid])
    else:
        polys=[]
        for poly in S.MP[oid]:
            polys+=mkpolys(poly)
    np_=norm(polys)
    if not np_: continue
    parks.append({"n":n,"agency":ag,"access":ac,"acres":acres,"polys":np_})

parks.sort(key=lambda p:-p["acres"])
src="https://gis.cnra.ca.gov/arcgis/rest/services/Boundaries/CPAD_AgencyLevel/MapServer/0 (CPAD Holdings by Agency Level, California Protected Areas Database / GreenInfo Network)"
json.dump({"source":src,"parks":parks},open("base_parks_fed.json","w"),separators=(",",":"))

# ---- verification ----
def pip(pt,poly):
    x,y=pt
    def inr(r):
        c=False; n=len(r)-1
        for i in range(n):
            x1,y1=r[i]; x2,y2=r[i+1]
            if (y1>y)!=(y2>y):
                xi=x1+(y-y1)*(x2-x1)/(y2-y1)
                if x<xi: c=not c
        return c
    for p in poly:
        if inr(p["outer"]) and not any(inr(h) for h in p["holes"]): return True
    return False

def find(n):
    for p in parks:
        if p["n"]==n: return p
    return None

tests=[("Presidio contains 37.7980,-122.4660","Presidio",(-122.4660,37.7980),True),
 ("Presidio contains 37.7930,-122.4750","Presidio",(-122.4750,37.7930),True),
 ("Presidio EXCLUDES Presidio Hts 37.7885,-122.4520","Presidio",(-122.4520,37.7885),False),
 ("Crissy Field contains 37.8035,-122.4650","Crissy Field",(-122.4650,37.8035),True),
 ("Lands End contains 37.7825,-122.5050","Golden Gate NRA - Lands End",(-122.5050,37.7825),True)]
for label,name,pt,exp in tests:
    p=find(name); got=pip(pt,p["polys"]) if p else None
    print(("PASS" if got==exp else "FAIL"),"|",label,"| got",got)

# Fort Funston: either half
ff=[find("Golden Gate NRA - Fort Funston (north)"),find("Golden Gate NRA - Fort Funston (south)")]
got=any(pip((-122.5030,37.7150),p["polys"]) for p in ff if p)
print(("PASS" if got else "FAIL"),"| Fort Funston contains 37.7150,-122.5030 | got",got)

# area of Presidio in km2 (equal-area approx via local scaling)
def km2(poly):
    tot=0.0
    for p in poly:
        for r,sgn in [(p["outer"],1)]+[(h,-1) for h in p["holes"]]:
            lat0=sum(q[1] for q in r)/len(r)
            k=111.32*math.cos(math.radians(lat0))
            s=0.0
            for i in range(len(r)-1):
                x1,y1=r[i]; x2,y2=r[i+1]
                s+=(x1*k)*(y2*110.574)-(x2*k)*(y1*110.574)
            tot+=sgn*abs(s/2.0)
    return tot
print("Presidio area: %.2f km2 (%.0f acres)"%(km2(find("Presidio")["polys"]),km2(find("Presidio")["polys"])*247.105))

bad=0; nh=0; no=0
for p in parks:
    for q in p["polys"]:
        no+=1
        if sarea(q["outer"])<=0: bad+=1
        for h in q["holes"]:
            nh+=1
            if sarea(h)>=0: bad+=1
print("rings: %d outers, %d holes, winding violations: %d"%(no,nh,bad))
print("units:",len(parks))
agset=sorted(set(p["agency"] for p in parks))
print("agencies:",agset)
print("recpark units:",[p["n"] for p in parks if "Recreation and Park" in p["agency"]])