#!/usr/bin/env python3
"""Read the seven temperature curves of Ocko et al. 2017, Fig. 2B, out of the PDF's own vectors.

Ocko, King, Andreen, Bardunias, Turner, Soar & Mahadevan (2017), "Solar-powered ventilation
of African termite mounds", J Exp Biol 220:3260-3269, doi:10.1242/jeb.160895. No data were
deposited with the paper. But Fig. 2B is drawn as vector paths, so the curves can be read
back exactly as plotted, with no pixel tracing and no guessing at colours.

Input:  the authors' lab copy of the published PDF,
        https://softmath.seas.harvard.edu/wp-content/uploads/2019/10/2017-10.pdf
        sha256 30427f10954efda92b6d49ef0cd214bec3b6daf6cce9e34d6b642f0c02b9b46c
        (not committed: it is the journal's copyright; fetch it yourself)
Needs:  pdftocairo (poppler-utils)
Output: data/ocko2017-fig2b.json

Axis calibration is read from the figure's own grid: the vertical grid lines at 0 h and 24 h
and the tick marks at 15 and 40 degC. Each curve is identified by stroke colour (the four
sides) or dash pattern (the three black curves), matching the figure's legend.

    python3 extract-fig2b.py path/to/2017-10.pdf
"""
import hashlib, json, re, subprocess, sys, tempfile, os

SHA = "30427f10954efda92b6d49ef0cd214bec3b6daf6cce9e34d6b642f0c02b9b46c"
pdf = sys.argv[1]
if hashlib.sha256(open(pdf, "rb").read()).hexdigest() != SHA:
    sys.exit("not the PDF this extraction was written against (sha256 differs)")

svgpath = os.path.join(tempfile.mkdtemp(), "p4.svg")
subprocess.run(["pdftocairo", "-svg", "-f", "4", "-l", "4", pdf, svgpath], check=True)
svg = open(svgpath).read()

def paths():
    for m in re.finditer(r"<path ([^>]*)/>", svg):
        a = m.group(1)
        d = re.search(r' d="([^"]*)"', a)
        if not d:
            continue
        st = re.search(r'stroke="([^"]*)"', a)
        da = re.search(r'stroke-dasharray="([^"]*)"', a)
        tr = re.search(r'transform="matrix\(([^)]*)\)"', a)
        M = [float(x) for x in tr.group(1).split(",")] if tr else [1, 0, 0, 1, 0, 0]
        toks = re.findall(r"[MLCZ]|-?[\d.]+(?:e-?\d+)?", d.group(1))
        pts, cur, cmd = [], [], None
        for t in toks:
            if t in "MLCZ":
                cmd = t
                continue
            cur.append(float(t))
            if len(cur) == {"M": 2, "L": 2, "C": 6}[cmd]:
                x, y = cur[-2], cur[-1]
                a_, b_, c_, d_, e_, f_ = M
                pts.append((a_ * x + c_ * y + e_, b_ * x + d_ * y + f_))
                cur = []
        yield (st.group(1) if st else None, da.group(1) if da else None, pts)

# The panel: grid verticals span y 236.51..347.11 (B's frame). 0 h and 24 h are the frame's
# left and right edges; 15 and 40 degC are the bottom and top ticks. Found by listing the
# short tick paths and the frame (see README); asserted below rather than trusted.
X0, X24, Y15, Y40 = 76.95, 295.45, 347.11, 236.51
COLOURS = {"92.98": "north", "14.95": "east", "23.22": "south", "63.59": "west"}
DASH = {None: "center", "0.92152 0.92152": "top", "0.15359 0.46076": "nest"}

curves = {}
for stroke, dash, pts in paths():
    if len(pts) < 30 or not stroke:
        continue
    xs = [p[0] for p in pts]
    if abs(min(xs) - X0) > 0.5 or abs(max(xs) - 295.14) > 0.5:
        continue  # only the full-day curves of panel B start at its left edge and end at its right
    ys = [p[1] for p in pts]
    if not (Y40 <= min(ys) and max(ys) <= Y15):
        continue
    name = None
    for k, v in COLOURS.items():
        if stroke.startswith("rgb(" + k):
            name = v
    if stroke.startswith("rgb(6.28"):
        name = DASH.get(dash)
    if not name:
        continue
    assert name not in curves, name
    curves[name] = [[round(24 * (x - X0) / (X24 - X0), 4), round(15 + 25 * (Y15 - y) / (Y15 - Y40), 4)] for x, y in pts]

assert sorted(curves) == sorted(["north", "east", "south", "west", "center", "top", "nest"]), sorted(curves)
out = {
    "source": "Ocko et al. 2017, J Exp Biol 220:3260, Fig. 2B, read from the PDF's vector paths",
    "pdf": "https://softmath.seas.harvard.edu/wp-content/uploads/2019/10/2017-10.pdf",
    "pdf_sha256": SHA,
    "units": {"x": "hour of day", "y": "degC"},
    "sensors": "iButton DS1922L, 0.5 degC accuracy (caption). Sides: ~1 m up, 5-10 cm below the surface. "
               "Center: central axis ~1 m up. Top: central axis ~15 cm from the top. Nest: ~0.5 m below ground.",
    "caveat": "One mound, one day, April-May 2014 or 2015 (the paper does not say which day). "
              "These are the curves as the authors plotted them, not their raw logger readings.",
    "curves": curves,
}
json.dump(out, open(os.path.join(os.path.dirname(os.path.abspath(__file__)), "data/ocko2017-fig2b.json"), "w"), indent=1)
print({k: len(v) for k, v in curves.items()})