#!/usr/bin/env python3
"""weather.py - station weather for the Dyatlov Pass dates and wind-chill tables.

(a) GHCN-Daily (NOAA NCEI, https://www.ncei.noaa.gov/pub/data/ghcn/daily/by_station/<ID>.csv.gz,
    downloaded 2026-09-14 into raw/ghcnd/) for the stations nearest the pass:
    RSM00023921 IVDEL 60.6831N 60.4500E 93 m (WMO 23921)            ~125 km SSE of the tent
    RSM00023724 NJAKSIMVOL 62.4300N 60.8694E 51 m (WMO 23724)        ~105 km NE
    RSM00023711 TROICKO-PECHERSKOE 62.70N 56.20E 135 m (WMO 23711)   ~200 km NW
    RSM00023914 CHERDYN' 60.40N 56.52E 208 m (WMO 23914)             ~215 km SW
    RSM00023418 PECHORA 65.1167N 57.10E 53 m (WMO 23418)             ~390 km N
    RSM00028044 SEROV 59.60N 60.533E 132 m (TAVG only, partial)      ~245 km S
    (RSM00023813 UST-UNJA and RSM00028033 KARPINSK have precipitation only; Burmantovo (WMO 23904 is
     actually KOIGORODOK in GHCN) has no GHCN-Daily record.)
    Values are tenths of degC / tenths of mm in the files; converted here. GHCN day = the station's
    climatological day (Russian stations: 24 h ending at 18 UTC = 23:00 local on the given date),
    so TMIN of '1959-02-02' covers the night 1-2 Feb.
(b) Burmantovo 1 Feb 1959 hourly table read from the scan of the loose sheet in Maslennikov's 2nd notebook
    (dyatlovpass.com/maslennikov-notebook-2; the site's own converted table has one slip at 00:00).
(c) Wind chill: JAG/TI 2001 formula (Osczevski & Bluestein 2005), Twc = 13.12 + 0.6215 T - 11.37 V^0.16 + 0.3965 T V^0.16,
    T in degC, V = 10 m wind in km/h; valid for T <= 10 degC and V >= 4.8 km/h.
Output: out/ghcn_daily_1959.csv, out/weather_tables.md, out/windchill_table.csv
"""
import csv, gzip, math, os, collections
from geo_common import RAW, OUT, haversine, load_coordinates

C = load_coordinates(); tent = C["tent_dyatlovpass_map"]
STATIONS = [("RSM00023921", "Ivdel", 60.6831, 60.4500, 93), ("RSM00023724", "Nyaksimvol", 62.4300, 60.8694, 51),
            ("RSM00023711", "Troitsko-Pechorsk", 62.70, 56.20, 135), ("RSM00023914", "Cherdyn", 60.40, 56.52, 208),
            ("RSM00023418", "Pechora", 65.1167, 57.10, 53), ("RSM00028044", "Serov", 59.60, 60.533, 132)]
PERIODS = [("1959-01-28", "1959-02-05"), ("1959-02-24", "1959-03-10")]

def in_periods(d):
    ds = f"{d[:4]}-{d[4:6]}-{d[6:]}"
    return any(a <= ds <= b for a, b in PERIODS)

rows = []
tables = {}
for sid, name, lat, lon, elev in STATIONS:
    dist = haversine(tent["lat"], tent["lon"], lat, lon) / 1000
    tab = collections.defaultdict(dict)
    with gzip.open(os.path.join(RAW, "ghcnd", f"{sid}.csv.gz"), "rt") as f:
        for r in csv.reader(f):
            if r[1].startswith("1959") and in_periods(r[1]):
                tab[r[1]][r[2]] = (int(r[3]) / 10.0, r[4], r[5], r[6])
    tables[sid] = (name, dist, tab)
    for d in sorted(tab):
        v = tab[d]
        rows.append({"station": sid, "name": name, "dist_km": round(dist), "date": f"{d[:4]}-{d[4:6]}-{d[6:]}",
                     "tmax_c": v.get("TMAX", ("",))[0], "tmin_c": v.get("TMIN", ("",))[0], "tavg_c": v.get("TAVG", ("",))[0],
                     "prcp_mm": v.get("PRCP", ("",))[0], "snwd_mm": v.get("SNWD", ("",))[0],
                     "flags": ";".join(f"{e}:{''.join(v[e][1:3])}" for e in v if v[e][1] or v[e][2])})
with open(os.path.join(OUT, "ghcn_daily_1959.csv"), "w", newline="") as f:
    w = csv.DictWriter(f, fieldnames=list(rows[0].keys())); w.writeheader(); w.writerows(rows)

md = ["# Station weather (generated by weather.py)", "", "## GHCN-Daily, 28 Jan - 5 Feb 1959 (degC; PRCP mm; TMIN of a date = the preceding night)", ""]
hdr = "| date | " + " | ".join(f"{n} ({round(d)} km) max/min/avg" for sid, (n, d, _) in tables.items()) + " |"
md += [hdr, "|" + "---|" * (len(tables) + 1)]
def cell(v):
    def g(e): return f"{v[e][0]:.1f}" if e in v else "-"
    return f"{g('TMAX')} / {g('TMIN')} / {g('TAVG')}" + (f" (P {v['PRCP'][0]:.1f})" if "PRCP" in v and v["PRCP"][0] > 0 else "")
for a, b in PERIODS:
    if a == "1959-02-24": md += ["", "## GHCN-Daily, 24 Feb - 10 Mar 1959", "", hdr, "|" + "---|" * (len(tables) + 1)]
    days = sorted({d for _, (_, _, t) in tables.items() for d in t if a <= f"{d[:4]}-{d[4:6]}-{d[6:]}" <= b})
    for d in days:
        md.append(f"| {d[:4]}-{d[4:6]}-{d[6:]} | " + " | ".join(cell(t.get(d, {})) for _, (_, _, t) in tables.items()) + " |")

# (b) Burmantovo 1 Feb 1959 - read directly from the scan of the loose sheet in Maslennikov's 2nd notebook
#     (raw/maslennikov/Dyatlov-pass-Maslennikov-loose-pages-Burmantovo-weather-report.jpg; dyatlovpass.com/maslennikov-notebook-2).
#     Columns on the sheet: hour, sky, wind direction, wind (m/s), temperature (degC).
burm = [("15:00", "overcast (пасм)", "N", 5, -8), ("18:00", "cloudy (обл)", "N", 1, -10), ("19:00", "cloudy", "N", 3, -11), ("21:00", "cloudy", "N", 1, -13),
        ("23:00", "cloudy", "W", 1, -5), ("00:00 (2 Feb)", "cloudy", "NW", 3, -15), ("03:00 (2 Feb)", "clear (ясно)", "W", 3, -21)]
md += ["", "## Burmantovo (village on the Lozva, ~78 km SE of the tent, ~180 m a.s.l.) 1 Feb 1959 - loose sheet in Maslennikov's 2nd notebook", "",
       "Read from the scan (primary source). dyatlovpass.com's converted table (burmantovo-weather-february-1959) agrees except at 00:00, where it prints -5 degF (= -20.6 degC) for the sheet's -15 degC (= +5 degF).",
       "The -5 degC at 23:00 (wind veering N -> W) followed by -15 at 00:00 (NW) and -21 at 03:00 (clear, W) is the signature of a cold-front passage between 23:00 and 00:00 at Burmantovo.", "",
       "| local time | sky | wind dir | wind m/s | temp degC |", "|---|---|---|---|---|"]
for t, sky, wd, v, tc in burm:
    md.append(f"| {t} | {sky} | {wd} | {v} | {tc} |")

# (c) wind chill
def wct(T, v_ms):
    V = v_ms * 3.6
    if V < 4.8: return T
    return 13.12 + 0.6215 * T - 11.37 * V ** 0.16 + 0.3965 * T * V ** 0.16
temps = [-10, -15, -20, -25, -30]; winds = [2, 5, 8, 10, 15, 20, 25, 30]
md += ["", "## Wind chill (JAG/TI 2001), degC, for 10-m wind speeds", "", "| air T \\ wind | " + " | ".join(f"{w} m/s" for w in winds) + " |", "|" + "---|" * (len(winds) + 1)]
with open(os.path.join(OUT, "windchill_table.csv"), "w", newline="") as f:
    w = csv.writer(f); w.writerow(["air_temp_c"] + [f"wind_{x}_ms" for x in winds])
    for T in temps:
        vals = [wct(T, x) for x in winds]
        w.writerow([T] + [round(v, 1) for v in vals])
        md.append(f"| {T} | " + " | ".join(f"{v:.0f}" for v in vals) + " |")
md += ["", "Frostbite guidance (Environment Canada, Osczevski & Bluestein 2005): exposed skin freezes in <30 min below about -28 wind chill, <10 min below about -40, <5 min below about -48 (with 1-2 m/s local wind higher than 10-m wind not accounted for).",
       "", "Note: the JAG/TI formula uses the 10-m wind; face-level wind is ~2/3 of it. Wind chill is not an air temperature: it is the equivalent calm-air temperature giving the same facial heat loss for a walking adult; it does not describe hypothermia rate for a whole body in wet clothes.",
       "", "Also computed: WCT at the Pigoltsina 2019 values (see write-up)."]
open(os.path.join(OUT, "weather_tables.md"), "w").write("\n".join(md))
print("\n".join(md))