#!/usr/bin/env python3
"""A rough digitization of Fig. 3a of Yuan et al. (2026), "The Qianlong XXXV event", Front. Astron. Space Sci.
13:1857812, https://doi.org/10.3389/fspas.2026.1857812 (CC BY 4.0), to estimate how far single storms scatter about
the paper's fitted curve. For https://factchecked.science/qianlong-1770-storm/ (Fact-Checked Science, 27 October
2026). Standard library only: it decodes the PNG itself and runs nothing from it.

It fetches the figure from the publisher's image service, as a PNG 1,772 pixels wide (about 0.3 MB, verified
HTTPS), or reads your own copy of that PNG:
  python3 -I digitize_fig3.py [fig3.png] [out.json]
URL: https://www.frontiersin.org/api/ipx/w=1772&f=png/https://www.frontiersin.org/files/Articles/1857812/xml-images/fspas-13-1857812-g003.webp
The copy we used has SHA-256 1bc45e3a2719ea17bf1beeded5eb7d86fb94fff369937e5ed2580683b4a6ac03; the script says
whether yours matches (if the service re-renders the figure differently, the pixel counts can shift a little).

Rules, fixed before the first run:
  - Axes, from the image: pixel columns 51 and 851 are -Dst = 0 and 3500 nT; pixel rows 483 and 24 are 0 and 90
    degrees MLAT (Fig. 3a's frame and tick labels).
  - A "blue" pixel has B - R > 50 and G - R > 30: the light-blue fill of the markers of the fitting subset. The
    fitted line and its confidence band are violet, with G about equal to R, and are left out. Only pixels inside
    the frame of panel a, above 2 degrees, are used.
  - Blue pixels are grouped into 8-connected blobs of at least 6 pixels; a blob counts as round(area / median blob
    area) markers at its centroid (markers that overlap cannot be separated, so their spread is lost).
  - Markers below -Dst = 200 nT are dropped (the fitting subset is -Dst > 250 nT; this allows for marker size).
  - The curve is the paper's, MLAT = 71.53 - 0.4252 (-Dst)^(2/3), with the coefficient solved from its two results
    (see reconstruct.py; it prints 0.43).
  - Reported: the markers recovered (the paper fitted 37 storms), R^2 of the curve on them (the paper: 0.563), the
    residual standard deviation s in degrees, and the -Dst range that a band of +/- 1.96 s about the curve gives at
    20.40 and 28.30 degrees: a rough 95% band for single storms, before the curve's own uncertainty is added.
"""
import hashlib
import json
import math
import statistics
import struct
import sys
import urllib.request
import zlib

URL = ("https://www.frontiersin.org/api/ipx/w=1772&f=png/"
       "https://www.frontiersin.org/files/Articles/1857812/xml-images/fspas-13-1857812-g003.webp")
OURS = "1bc45e3a2719ea17bf1beeded5eb7d86fb94fff369937e5ed2580683b4a6ac03"
X0, X1, Y0, Y1 = 51, 851, 483, 24
A, C = 0.4252, 71.53


def png_rgb(data):
    """Decode an 8-bit, non-interlaced RGB or RGBA PNG into rows of (r, g, b) tuples."""
    if data[:8] != b"\x89PNG\r\n\x1a\n":
        sys.exit("not a PNG")
    i, idat, w = 8, [], None
    while i < len(data):
        n = struct.unpack(">I", data[i:i + 4])[0]
        kind, body = data[i + 4:i + 8], data[i + 8:i + 8 + n]
        if kind == b"IHDR":
            w, h, depth, ctype, _, _, inter = struct.unpack(">IIBBBBB", body)
            if depth != 8 or ctype not in (2, 6) or inter != 0:
                sys.exit("unsupported PNG (need 8-bit RGB or RGBA, not interlaced)")
            bpp = 3 if ctype == 2 else 4
        elif kind == b"IDAT":
            idat.append(body)
        elif kind == b"IEND":
            break
        i += 12 + n
    raw = zlib.decompress(b"".join(idat))
    stride = w * bpp
    prev = bytearray(stride)
    rows = []
    for y in range(h):
        f = raw[y * (stride + 1)]
        line = bytearray(raw[y * (stride + 1) + 1:(y + 1) * (stride + 1)])
        if f == 1:
            for x in range(bpp, stride):
                line[x] = (line[x] + line[x - bpp]) & 255
        elif f == 2:
            for x in range(stride):
                line[x] = (line[x] + prev[x]) & 255
        elif f == 3:
            for x in range(stride):
                left = line[x - bpp] if x >= bpp else 0
                line[x] = (line[x] + ((left + prev[x]) >> 1)) & 255
        elif f == 4:
            for x in range(stride):
                a = line[x - bpp] if x >= bpp else 0
                b = prev[x]
                c = prev[x - bpp] if x >= bpp else 0
                p = a + b - c
                pa, pb, pc = abs(p - a), abs(p - b), abs(p - c)
                line[x] = (line[x] + (a if pa <= pb and pa <= pc else b if pb <= pc else c)) & 255
        elif f != 0:
            sys.exit(f"bad PNG filter {f}")
        rows.append([(line[x], line[x + 1], line[x + 2]) for x in range(0, stride, bpp)])
        prev = line
    return rows


def to_data(px, py):
    return (px - X0) / (X1 - X0) * 3500.0, (Y0 - py) / (Y0 - Y1) * 90.0


def curve(d):
    return C - A * d ** (2 / 3)


def inv(lat):
    return ((C - lat) / A) ** 1.5


def main():
    sys.stdout.reconfigure(encoding="utf-8")
    args = sys.argv[1:]
    src = args[0] if args and args[0].lower().endswith(".png") else URL
    outp = args[-1] if args and args[-1].lower().endswith(".json") else None
    if src == URL:
        req = urllib.request.Request(URL, headers={"User-Agent": "factchecked.science digitize_fig3.py"})
        with urllib.request.urlopen(req, timeout=120) as r:
            data = r.read()
    else:
        with open(src, "rb") as f:
            data = f.read()
    digest = hashlib.sha256(data).hexdigest()
    print(f"read {src}: {len(data):,} bytes, SHA-256 {digest} ({'the same file as ours' if digest == OURS else 'NOT the file we used'})")
    img = png_rgb(data)
    h, w = len(img), min(X1, len(img[0]))
    mask = [[False] * w for _ in range(h)]
    for y in range(Y1 + 2, Y0 - 9):
        row = img[y]
        for x in range(X0 + 1, w):
            r, g, b = row[x]
            mask[y][x] = (b - r > 50) and (g - r > 30)
    seen = [[False] * w for _ in range(h)]
    blobs = []
    for y0 in range(h):
        for x0 in range(w):
            if mask[y0][x0] and not seen[y0][x0]:
                stack, pts = [(y0, x0)], []
                seen[y0][x0] = True
                while stack:
                    y, x = stack.pop()
                    pts.append((y, x))
                    for dy in (-1, 0, 1):
                        for dx in (-1, 0, 1):
                            yy, xx = y + dy, x + dx
                            if 0 <= yy < h and 0 <= xx < w and mask[yy][xx] and not seen[yy][xx]:
                                seen[yy][xx] = True
                                stack.append((yy, xx))
                if len(pts) >= 6:
                    blobs.append((len(pts), sum(p[1] for p in pts) / len(pts), sum(p[0] for p in pts) / len(pts)))
    single = statistics.median(a for a, _, _ in blobs)
    pts = []
    for a, x, y in blobs:
        pts += [to_data(x, y)] * max(1, int(round(a / single)))
    pts = [(d, lat) for d, lat in pts if d >= 200]
    res = [lat - curve(d) for d, lat in pts]
    mean_lat = sum(lat for _, lat in pts) / len(pts)
    ss_res = sum(r * r for r in res)
    s = math.sqrt(ss_res / (len(res) - 2))
    r2 = 1 - ss_res / sum((lat - mean_lat) ** 2 for _, lat in pts)
    out = {"blobs": len(blobs), "single_marker_area_px": single, "markers": len(pts),
           "dst_range": [min(d for d, _ in pts), max(d for d, _ in pts)], "r2_of_paper_curve": r2, "residual_sd_deg": s,
           "points": [[round(d), round(lat, 1)] for d, lat in sorted(pts)]}
    print(f"{len(blobs)} blobs, a single marker about {single} pixels: {len(pts)} markers recovered (the paper fitted 37),"
          f" -Dst {out['dst_range'][0]:.0f} to {out['dst_range'][1]:.0f} nT")
    print(f"R^2 of the paper's curve on them {r2:.2f} (the paper: 0.563); scatter about the curve s = {s:.1f} deg")
    for name, lat in (("1770 (20.40 deg)", 20.40), ("1859 (28.30 deg)", 28.30)):
        lo, hi = lat + 1.96 * s, lat - 1.96 * s
        out[name] = {"central": inv(lat), "pi95": [inv(lo), inv(hi) if hi < C else None]}
        print(f"  {name}: {inv(lat):,.0f} nT on the curve; +/- 1.96 s gives {inv(lo):,.0f} to {inv(hi):,.0f} nT")
    if outp:
        with open(outp, "w", encoding="utf-8", newline="\n") as f:
            json.dump(out, f, indent=1)
        print("wrote", outp)


if __name__ == "__main__":
    main()