"""Before the seal: how noisy is a 5-day window against its flanks? Reads ONLY rows from
August to October of the stations named on the command line (seven of the frame, every fifth), so no spring number is seen before the rule is sealed."""
import io, sys, zipfile, urllib.request, re, datetime as dt
import numpy as np
U = "https://opendata.dwd.de/climate_environment/CDC/observations_germany/climate/daily/kl/historical/"
idx = urllib.request.urlopen(U, timeout=60).read().decode()
want = sys.argv[1:]
out = {}
for sid in want:
    name = re.search(rf"tageswerte_KL_{sid}_\d+_\d+_hist\.zip", idx).group(0)
    z = zipfile.ZipFile(io.BytesIO(urllib.request.urlopen(U + name, timeout=120).read()))
    f = [n for n in z.namelist() if n.startswith("produkt_klima_tag")][0]
    rows = {}
    for line in io.TextIOWrapper(z.open(f), "latin1"):
        p = [x.strip() for x in line.split(";")]
        if p[0] == "STATIONS_ID" or len(p) < 15:
            continue
        d = p[1]
        if d[4:6] not in ("08", "09", "10"):
            continue  # nothing outside August to October is even parsed
        t = float(p[13])
        if t > -900:
            rows[dt.date(int(d[:4]), int(d[4:6]), int(d[6:8]))] = t
    out[sid] = rows
    print(sid, len(rows), flush=True)
def contrast(rows, y, m, d):
    s = dt.date(y, m, d)
    w = [rows.get(s + dt.timedelta(i)) for i in range(5)]
    fl = [rows.get(s + dt.timedelta(i)) for i in list(range(-10, 0)) + list(range(5, 15))]
    if None in w or sum(x is None for x in fl) > 2:
        return None
    a = np.mean([x for x in fl[:10] if x is not None])
    b = np.mean([x for x in fl[10:] if x is not None])
    return (a + b) / 2 - np.mean(w)
for m, d in ((9, 1), (9, 11), (9, 21), (10, 1)):
    ys = []
    for y in range(1925, 2026):
        v = [contrast(r, y, m, d) for r in out.values()]
        v = [x for x in v if x is not None]
        if len(v) >= len(out) / 2:
            ys.append(np.mean(v))
    ys = np.array(ys)
    print(m, d, "n", len(ys), "mean %.3f sd %.3f se %.3f" % (ys.mean(), ys.std(ddof=1), ys.std(ddof=1) / len(ys) ** .5))