import json, glob, math, os

BBOX = (-122.52, 37.70, -122.35, 37.84)
LAT0 = 37.77
MX = 111320.0 * math.cos(math.radians(LAT0))   # m per deg lon
MY = 110574.0                                   # m per deg lat
TOL_M = 40.0
MIN_AREA_M2 = 10000.0

def to_m(p):  return (p[0]*MX, p[1]*MY)
def clip_rect(ring, b):
    xmin,ymin,xmax,ymax = b
    pts = list(ring)
    if len(pts)>1 and pts[0]==pts[-1]: pts = pts[:-1]
    edges = [
        (lambda p: p[0]>=xmin, 0, xmin),
        (lambda p: p[0]<=xmax, 0, xmax),
        (lambda p: p[1]>=ymin, 1, ymin),
        (lambda p: p[1]<=ymax, 1, ymax),
    ]
    for inside, axis, val in edges:
        if not pts: return []
        out = []
        n = len(pts)
        for i in range(n):
            cur, prv = pts[i], pts[i-1]
            ci, pi = inside(cur), inside(prv)
            if ci != pi:
                d = cur[axis]-prv[axis]
                t = 0.0 if d==0 else (val-prv[axis])/d
                ip = [prv[0]+(cur[0]-prv[0])*t, prv[1]+(cur[1]-prv[1])*t]
                ip[axis] = val
                if ci: out.append(ip); out.append(list(cur))
                else:  out.append(ip)
            elif ci:
                out.append(list(cur))
        pts = out
    return pts

def dp(pts, tol):
    if len(pts) < 3: return pts
    a, b = to_m(pts[0]), to_m(pts[-1])
    dx, dy = b[0]-a[0], b[1]-a[1]
    L = math.hypot(dx, dy)
    imax, dmax = 0, -1.0
    for i in range(1, len(pts)-1):
        p = to_m(pts[i])
        if L == 0:
            d = math.hypot(p[0]-a[0], p[1]-a[1])
        else:
            d = abs(dy*(p[0]-a[0]) - dx*(p[1]-a[1]))/L
        if d > dmax: imax, dmax = i, d
    if dmax <= tol: return [pts[0], pts[-1]]
    return dp(pts[:imax+1], tol)[:-1] + dp(pts[imax:], tol)

def simplify_ring(ring, tol):
    r = list(ring)
    if len(r)>1 and r[0]==r[-1]: r = r[:-1]
    if len(r) < 3: return []
    s = dp(r + [r[0]], tol)
    if len(s)>1 and s[0]==s[-1]: s = s[:-1]
    return s if len(s) >= 3 else []

def area_m2(ring):
    if len(ring) < 3: return 0.0
    s = 0.0
    for i in range(len(ring)):
        x1,y1 = to_m(ring[i]); x2,y2 = to_m(ring[(i+1)%len(ring)])
        s += x1*y2 - x2*y1
    return abs(s)/2.0

def process(rings_with_meta, tol=TOL_M):
    out, dropped_small, dropped_clip = [], 0, 0
    for ring, meta in rings_with_meta:
        c = clip_rect(ring, BBOX)
        if len(c) < 3: dropped_clip += 1; continue
        s = simplify_ring(c, tol)
        if len(s) < 3: dropped_clip += 1; continue
        if area_m2(s) < MIN_AREA_M2: dropped_small += 1; continue
        out.append({"sev": meta["sev"], "class": meta["class"],
                    "ring": [[round(p[0],5), round(p[1],5)] for p in s]})
    return out, dropped_small, dropped_clip

# ---------- LIQUEFACTION ----------
lq_rings, lq_holes = [], 0
for f in sorted(glob.glob('haz_raw/lq*.json')):
    for feat in json.load(open(f))['features']:
        coords = feat['geometry']['coordinates']
        lq_holes += len(coords) - 1
        lq_rings.append((coords[0], {"sev": 1.0,
            "class": "Liquefaction Zone of Required Investigation"}))
lq_polys, lq_s, lq_c = process(lq_rings)
lq = {
 "source": "https://services2.arcgis.com/zr3KAIbsRSUyARHG/ArcGIS/rest/services/CGS_Liquefaction_Zones/FeatureServer/0 - layer 'CGS Liquefaction Zones' (California Geological Survey, Seismic Hazards Mapping Act regulatory Seismic Hazard Zones)",
 "layer_type": "REGULATORY Seismic Hazard Zone for liquefaction (statutory 'Zone of Required Investigation' under the Seismic Hazards Mapping Act). Binary in/out zone - NOT a graded liquefaction-susceptibility surface, so no low/moderate/high classes exist.",
 "scenario": "n/a (regulatory zone, not a scenario). Quads: San Francisco North, San Francisco South, Hunters Point, Oakland West",
 "crs": "wgs84", "polygons": lq_polys}
json.dump(lq, open('haz_liquefaction.json','w'), separators=(',',':'))

# ---------- TSUNAMI ----------
ts_rings, ts_holes = [], 0
for feat in json.load(open('haz_raw/ts38.json'))['features']:
    coords = feat['geometry']['coordinates']
    ts_holes += len(coords) - 1
    ts_rings.append((coords[0], {"sev": 1.0, "class": feat['properties']['Evacuate']}))
ts_polys, ts_s, ts_c = process(ts_rings)
ts = {
 "source": "https://gis.conservation.ca.gov/server/rest/services/CGS/IW_Tsunami_Hazard_Area/FeatureServer/0 - layer 'CA_Tsunami_Hazard_Area' (California Geological Survey / Cal OES), OBJECTID 38 = San Francisco County",
 "layer_type": "Tsunami Hazard Area (maximum-considered tsunami inundation for evacuation planning; model results represent inundation exceeding a ~975-year average return period event). Binary hazard-area polygon.",
 "scenario": "Maximum considered tsunami, ~975-year average return period; San Francisco County sheet",
 "crs": "wgs84", "polygons": ts_polys}
json.dump(ts, open('haz_tsunami.json','w'), separators=(',',':'))

def stats(name, polys, holes, ds, dc):
    v = sum(len(p['ring']) for p in polys)
    print(f"{name}: {len(polys)} polygons, {v} vertices, "
          f"{os.path.getsize(name)} bytes, dropped_small={ds} dropped_by_clip={dc} "
          f"interior_holes_in_source_dropped={holes}")
stats('haz_liquefaction.json', lq_polys, lq_holes, lq_s, lq_c)
stats('haz_tsunami.json', ts_polys, ts_holes, ts_s, ts_c)