import re, math, json

def parse_multipolygon(wkt):
    """Return list of rings (each list of (lon,lat)) for a MULTIPOLYGON."""
    body = wkt[wkt.index('('):]
    rings = []
    for m in re.finditer(r'\(([-0-9eE., ]+)\)', body):
        pts = []
        for pair in m.group(1).split(','):
            pair = pair.strip()
            if not pair: continue
            a, b = pair.split()
            pts.append((float(a), float(b)))
        if len(pts) >= 4:
            rings.append(pts)
    return rings

def ring_area_m2(ring):
    """Planar shoelace area in m^2 using local equirectangular scaling."""
    if len(ring) < 3: return 0.0
    lat0 = sum(p[1] for p in ring)/len(ring)
    kx = 111320.0*math.cos(math.radians(lat0)); ky = 110540.0
    s = 0.0
    for i in range(len(ring)):
        x1,y1 = ring[i][0]*kx, ring[i][1]*ky
        x2,y2 = ring[(i+1)%len(ring)][0]*kx, ring[(i+1)%len(ring)][1]*ky
        s += x1*y2 - x2*y1
    return abs(s)/2.0

def _perp(p, a, b):
    if a == b:
        return math.hypot(p[0]-a[0], p[1]-a[1])
    dx, dy = b[0]-a[0], b[1]-a[1]
    t = ((p[0]-a[0])*dx + (p[1]-a[1])*dy)/(dx*dx+dy*dy)
    t = max(0.0, min(1.0, t))
    return math.hypot(p[0]-(a[0]+t*dx), p[1]-(a[1]+t*dy))

def dp(pts, tol):
    """Douglas-Peucker on an open polyline (degrees)."""
    if len(pts) < 3: return pts[:]
    keep = [False]*len(pts); keep[0] = keep[-1] = True
    stack = [(0, len(pts)-1)]
    while stack:
        i, j = stack.pop()
        if j <= i+1: continue
        dmax, idx = -1.0, -1
        for k in range(i+1, j):
            d = _perp(pts[k], pts[i], pts[j])
            if d > dmax: dmax, idx = d, k
        if dmax > tol:
            keep[idx] = True
            stack.append((i, idx)); stack.append((idx, j))
    return [p for p, k in zip(pts, keep) if k]

def dp_ring(ring, tol):
    """Simplify a closed ring, keeping it closed."""
    r = ring[:]
    if r[0] == r[-1]: r = r[:-1]
    if len(r) < 4: return r
    # anchor at the two most distant-ish points to avoid collapsing the ring
    out = dp(r + [r[0]], tol)
    if out[0] == out[-1]: out = out[:-1]
    return out

def despike(ring, max_spike_m=200.0, cos_thresh=-0.85):
    """Drop vertices that form a near-retracing needle (simplification artifact)."""
    r = ring[:]
    if r[0] == r[-1]: r = r[:-1]
    changed = True
    while changed and len(r) > 4:
        changed = False
        lat0 = sum(p[1] for p in r)/len(r)
        kx = 111320.0*math.cos(math.radians(lat0)); ky = 110540.0
        out, n, i = [], len(r), 0
        skip = set()
        for i in range(n):
            if i in skip: continue
            p0, p1, p2 = r[(i-1) % n], r[i], r[(i+1) % n]
            v1 = ((p0[0]-p1[0])*kx, (p0[1]-p1[1])*ky)
            v2 = ((p2[0]-p1[0])*kx, (p2[1]-p1[1])*ky)
            l1 = math.hypot(*v1); l2 = math.hypot(*v2)
            if l1 > 0 and l2 > 0:
                c = (v1[0]*v2[0] + v1[1]*v2[1])/(l1*l2)
                if c < cos_thresh and min(l1, l2) < max_spike_m:
                    skip.add(i); changed = True
        r = [p for i, p in enumerate(r) if i not in skip]
    return r

def dissolve_two(a, b):
    """Merge two rings that share a contiguous chain of identical vertices."""
    A = a[:-1] if a[0] == a[-1] else a[:]
    B = b[:-1] if b[0] == b[-1] else b[:]
    sa, sb = set(A), set(B)
    shared = sa & sb
    if len(shared) < 2: return None

    def split(R):
        n = len(R)
        inS = [p in shared for p in R]
        if all(inS) or not any(inS): return None
        start = next(i for i in range(n) if not inS[i] and inS[(i-1) % n])
        rot = [R[(start+k) % n] for k in range(n)]
        inR = [p in shared for p in rot]
        # rot = [non-shared block ..., shared block ...]
        try:
            first_shared = next(i for i in range(n) if inR[i])
        except StopIteration:
            return None
        if not all(inR[first_shared:]): return None   # shared run not contiguous
        nonshared = rot[:first_shared]
        sh = rot[first_shared:]
        return [sh[-1]] + nonshared + [sh[0]]

    pa, pb = split(A), split(B)
    if pa is None or pb is None: return None
    if pb[0] != pa[-1]:
        pb = pb[::-1]
    if pb[0] != pa[-1] or pb[-1] != pa[0]:
        return None
    return pa + pb[1:-1]