#!/usr/bin/env python3
"""Estonia's cereal harvests, 1980-2026, from Statistics Estonia's table PM0281.

Our summary for https://factchecked.science/estonia-growing-season/ (Fact-Checked Science, October 26, 2026):
yearly cereal area, production and yield for the whole country; means by decade and period; least-squares and
Theil-Sen trends; the trend since 2015 with a Mann-Kendall test, and five-year means of production; record years;
and other crops, 2000-2009 against 2016-2025.

Input: estonia-crops.csv, next to this script (Statistics Estonia, PM0281: Agricultural land and crops by county,
whole country, latest update 30 September 2026; the 2026 values are preliminary, production from expert
assessments). Statistics Estonia's open data are licensed CC BY-SA 4.0; so is that file.

Usage:
    python3 crops_summary.py               # summarise estonia-crops.csv
    python3 crops_summary.py --fetch       # first download a fresh copy from Statistics Estonia's API
Standard library only (Python 3.8 or later).
"""
import argparse
import csv
import json
import math
import os
import statistics
import sys
import urllib.request

HERE = os.path.dirname(os.path.abspath(__file__))
CSV_FILE = os.path.join(HERE, "estonia-crops.csv")
API = ("https://andmed.stat.ee/api/v1/en/stat/majandus/pellumajandus/pellumajandussaaduste-tootmine/"
       "taimekasvatussaaduste-tootmine/PM0281.px")
CROPS = ["CRC0000", "CRC1100", "CRC1210", "CRC1300", "CRC1410", "CRC1600", "CRR1000", "CRI1110", "CRP0000", "CRARA"]


def fetch():
    """Download PM0281 for the whole country and rewrite estonia-crops.csv."""
    query = {"query": [
        {"code": "Maakondade ja üksuste jaotus", "selection": {"filter": "item", "values": ["EE_T"]}},
        {"code": "Põllumajandusmaa/põllukultuur", "selection": {"filter": "item", "values": CROPS}},
        {"code": "Näitaja", "selection": {"filter": "item", "values": ["AR_C", "AR_HA_PR", "PR", "YI"]}}],
        "response": {"format": "json-stat2"}}
    req = urllib.request.Request(API, data=json.dumps(query).encode("utf-8"), method="POST",
                                 headers={"Content-Type": "application/json", "User-Agent": "crops_summary.py (factchecked.science)"})
    with urllib.request.urlopen(req, timeout=120) as r:
        js = json.loads(r.read())
    ids, sizes, vals = js["id"], js["size"], js["value"]
    cats = {}
    for d in ids:
        c = js["dimension"][d]["category"]
        cats[d] = [(k, c.get("label", {}).get(k, k)) for k in sorted(c["index"], key=lambda k: c["index"][k])]
    stride = {d: 1 for d in ids}
    for i, d in enumerate(ids):
        for n in sizes[i + 1:]:
            stride[d] *= n
    crop, ind, year = "Põllumajandusmaa/põllukultuur", "Näitaja", "Vaatlusperiood"
    rows = []
    for ci, (cc, cn) in enumerate(cats[crop]):
        for ii, (ic, iname) in enumerate(cats[ind]):
            for yi, (yc, _) in enumerate(cats[year]):
                v = vals[ci * stride[crop] + ii * stride[ind] + yi * stride[year]]
                if v is not None:
                    rows.append((int(yc), cc, cn, ic, iname, v))
    rows.sort(key=lambda r: (r[1], r[3], r[0]))
    with open(CSV_FILE, "w", newline="", encoding="utf-8") as f:
        w = csv.writer(f, lineterminator="\n")
        w.writerow(["year", "crop_code", "crop", "indicator_code", "indicator", "value"])
        for r in rows:
            w.writerow([r[0], r[1], r[2], r[3], r[4], int(r[5]) if float(r[5]).is_integer() else r[5]])
    print(f"downloaded {len(rows)} values into {CSV_FILE}")


def load():
    data = {}
    with open(CSV_FILE, encoding="utf-8") as f:
        for row in csv.DictReader(f):
            data.setdefault(row["crop_code"], {"name": row["crop"]}).setdefault(row["indicator_code"], {})[int(row["year"])] = float(row["value"])
    return data


def ols(xs, ys):
    mx, my = statistics.fmean(xs), statistics.fmean(ys)
    return sum((x - mx) * (y - my) for x, y in zip(xs, ys)) / sum((x - mx) ** 2 for x in xs)


def theil_sen(xs, ys):
    return statistics.median((ys[j] - ys[i]) / (xs[j] - xs[i]) for i in range(len(xs)) for j in range(i + 1, len(xs)))


def mean_of(s, y0, y1):
    v = [s[y] for y in range(y0, y1 + 1) if y in s]
    return statistics.fmean(v) if v else None


def mann_kendall_p(ys):
    """Two-sided Mann-Kendall p, normal approximation with tie correction (as in gsl_rerun.py)."""
    n = len(ys)
    s = sum((ys[j] > ys[i]) - (ys[j] < ys[i]) for i in range(n - 1) for j in range(i + 1, n))
    ties = {}
    for v in ys:
        ties[v] = ties.get(v, 0) + 1
    var = (n * (n - 1) * (2 * n + 5) - sum(t * (t - 1) * (2 * t + 5) for t in ties.values())) / 18.0
    z = 0.0 if s == 0 else (s - (1 if s > 0 else -1)) / math.sqrt(var)
    return math.erfc(abs(z) / math.sqrt(2))


def main():
    ap = argparse.ArgumentParser(description="Summarise Estonia's cereal statistics (Statistics Estonia, PM0281).")
    ap.add_argument("--fetch", action="store_true", help="download a fresh copy of the table first")
    a = ap.parse_args()
    if hasattr(sys.stdout, "reconfigure"):
        sys.stdout.reconfigure(encoding="utf-8")
    if a.fetch:
        fetch()
    d = load()
    c = d["CRC0000"]
    print("Cereals for grain, whole country (Statistics Estonia, PM0281; 2026 preliminary)")
    print("year   area (ha)   production (t)   yield (kg/ha)")
    for y in sorted(set(c.get("PR", {})) | set(c.get("YI", {}))):
        print(f"{y}  {c.get('AR_C', {}).get(y, float('nan')):10,.0f}  {c.get('PR', {}).get(y, float('nan')):14,.0f}  {c.get('YI', {}).get(y, float('nan')):12,.0f}")
    print("\nPeriod means: production (t), yield (kg/ha), area (ha)")
    for y0, y1 in [(1980, 1989), (1990, 1999), (2000, 2009), (2010, 2019), (2016, 2025), (2026, 2026)]:
        print(f"  {y0}-{y1}: {mean_of(c['PR'], y0, y1):12,.0f}  {mean_of(c['YI'], y0, y1):7,.0f}  {mean_of(c['AR_C'], y0, y1):10,.0f}")
    print("\nTrends (per year)")
    for y0, y1 in [(1992, 2025), (2000, 2025), (2010, 2025)]:
        for ind, unit in (("PR", "t"), ("YI", "kg/ha")):
            xs = [y for y in range(y0, y1 + 1) if y in c[ind]]
            ys = [c[ind][y] for y in xs]
            print(f"  {ind} {y0}-{y1}: least squares {ols(xs, ys):+,.1f} {unit}, Theil-Sen {theil_sen(xs, ys):+,.1f} {unit} ({len(xs)} years)")
    print("\nSince 2015 (2026 left out, being preliminary): Theil-Sen trend and Mann-Kendall p")
    for ind, unit in (("PR", "t"), ("YI", "kg/ha")):
        xs = [y for y in range(2015, 2026) if y in c[ind]]
        ys = [c[ind][y] for y in xs]
        print(f"  {ind} 2015-2025: {theil_sen(xs, ys):+,.1f} {unit} a year, p = {mann_kendall_p(ys):.2f} ({len(xs)} years)")
    print("Five-year means of production (million t): " + "; ".join(
        f"{y0}-{y0 + 4} {mean_of(c['PR'], y0, y0 + 4) / 1e6:.2f}" for y0 in (2001, 2006, 2011, 2016, 2021)))
    for ind in ("PR", "YI"):
        top = sorted(((y, v) for y, v in c[ind].items()), key=lambda kv: -kv[1])[:5]
        low = sorted(((y, v) for y, v in c[ind].items()), key=lambda kv: kv[1])[:3]
        print(f"  highest {ind}: {top}; lowest: {low}")
    print("\nOther crops: mean 2000-2009 -> mean 2016-2025 (production t; yield kg/ha; area ha)")
    for code in ["CRC1100", "CRC1210", "CRC1300", "CRC1410", "CRC1600", "CRR1000", "CRI1110", "CRP0000"]:
        e = d[code]
        parts = []
        for ind in ("PR", "YI", "AR_C"):
            s = e.get(ind, {})
            parts.append(f"{mean_of(s, 2000, 2009):,.0f} -> {mean_of(s, 2016, 2025):,.0f}")
        print(f"  {e['name']}: " + "; ".join(parts))


if __name__ == "__main__":
    main()