"""Heat lost by the Atlantic between 26.5°N and 70°N in the public DEEP-C data, as used on
factchecked.science/amoc-heat-loss/: a partial check of Bell et al. (2026).

Usage:
  python3 -I deepc_region.py FMASS_FILE_OR_FOLDER [results.json]

  If the argument is a folder, the script downloads DEEP-C_Fmass_v5_198501-201711.nc (207 MB) into it from the
  University of Reading Research Data Archive and checks its SHA-256; if it is a file, it checks and reads that.

The data: Liu, C. and Allan, R. P. (2022), Reconstructions of the radiation fluxes at the top of atmosphere and net
surface energy flux: DEEP-C Version 5.0, University of Reading, doi:10.17864/1947.000347, CC BY 4.0. The file holds
the monthly net downward surface energy flux (W/m2), January 1985 to November 2017, built with ERA5. The preprint
uses version 5.1, 1985-2023, with ERA5, JRA55 and MERRA2, which is not public; so this is a partial check, not a
reproduction.

What it does (our choices, not the authors'):
- the region: the blue polygon of the preprint's Fig. 1(a), vertices read from the figure, in degrees east and north:
  (-80,26.5) (-5,26.5) (-5,40) (0,40) (0,50) (-5,50) (-5,60) (15,60) (15,70) (-70,70) (-70,60) (-75,60) (-75,40)
  (-80,40); a grid cell counts if its centre is inside;
- ocean cells only: the file has no land mask, so a cell counts as ocean if its monthly values over 1985-2016 vary
  with a standard deviation above 40 W/m2 (most land cells in this product vary by less than about 25, most ocean
  cells in the region by more than about 85, and some cells fall between); 60 and 80 are tried
  too, and change the trend by less than 0.002 PW a decade;
- cell areas on the file's N128 Gaussian grid, from latitude bounds halfway between grid latitudes;
- heat loss = minus the downward flux, integrated over the region, in PW; annual means of complete years, 1985-2016;
- the trend and its 95% interval with the preprint's own recipe (lags up to 19, as for its Fig. 1) and with AR(1).

Needs Python 3 with NumPy and SciPy in the interpreter's own site-packages (python3 -I ignores packages installed
with pip install --user). Reads the file as data only. Statistics: trendstats.py (same folder).
"""
import hashlib
import json
import os
import sys
import urllib.request

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import numpy as np  # noqa: E402
from scipy.io import netcdf_file  # noqa: E402

from trendstats import fit, summary  # noqa: E402

sys.stdout.reconfigure(encoding="utf-8")
URL = "https://researchdata.reading.ac.uk/347/5/DEEP-C_Fmass_v5_198501-201711.nc"
NAME = "DEEP-C_Fmass_v5_198501-201711.nc"
SHA256 = "d69b54767781ee0a7880e022ba2dc59c8bbb263494871000278be585af4bc1f3"
POLY = [(-80, 26.5), (-5, 26.5), (-5, 40), (0, 40), (0, 50), (-5, 50), (-5, 60), (15, 60), (15, 70),
        (-70, 70), (-70, 60), (-75, 60), (-75, 40), (-80, 40)]
R_EARTH = 6.371e6

if len(sys.argv) < 2:
    sys.exit(__doc__)
target = sys.argv[1]
out_json = sys.argv[2] if len(sys.argv) > 2 else None
path = os.path.join(target, NAME) if os.path.isdir(target) else target
if not os.path.isfile(path):
    print(f"downloading {URL} (207 MB) to {path}")
    req = urllib.request.Request(URL, headers={"User-Agent": "factchecked.science re-run (deepc_region.py)"})
    with urllib.request.urlopen(req, timeout=120) as r, open(path, "wb") as f:
        while True:
            chunk = r.read(1 << 20)
            if not chunk:
                break
            f.write(chunk)
h = hashlib.sha256()
with open(path, "rb") as f:
    for chunk in iter(lambda: f.read(1 << 22), b""):
        h.update(chunk)
if h.hexdigest() != SHA256:
    sys.exit(f"{path}: SHA-256 {h.hexdigest()} is not the {SHA256} of the file we used; stopping")
print(f"{NAME}: SHA-256 matches")


def inside(x, y, poly):
    """Even-odd rule, vectorised over arrays of points."""
    res = np.zeros(x.shape, bool)
    for i in range(len(poly)):
        x1, y1 = poly[i]
        x2, y2 = poly[(i + 1) % len(poly)]
        cross = (y1 > y) != (y2 > y)
        xint = x1 + (y - y1) * (x2 - x1) / ((y2 - y1) if y2 != y1 else 1e-12)
        res ^= cross & (x < xint)
    return res


with netcdf_file(path, "r", mmap=False) as f:
    F = f.variables["Fmass"].data.astype("f8")
    lat = f.variables["lat"].data.astype("f8").copy()
    lon = f.variables["lon"].data.astype("f8").copy()
    days = f.variables["time"].data.astype("f8").copy()
when = np.datetime64("1985-01-01") + days.astype("timedelta64[D]")
years = when.astype("datetime64[Y]").astype(int) + 1970
full = sorted(y for y in set(years.tolist()) if (years == y).sum() == 12)
print(f"{F.shape[0]} months, {str(when[0])} to {str(when[-1])}; complete years {full[0]}-{full[-1]}")

lat_b = np.concatenate([[90.0], 0.5 * (lat[1:] + lat[:-1]), [-90.0]])
band = np.abs(np.sin(np.radians(lat_b[:-1])) - np.sin(np.radians(lat_b[1:])))
area = (R_EARTH ** 2 * np.radians(360.0 / len(lon)) * band)[:, None] * np.ones((1, len(lon)))
lon180 = np.where(lon > 180, lon - 360, lon)
LON, LAT = np.meshgrid(lon180, lat)
region = inside(LON, LAT, POLY)
sd = F[np.isin(years, full)].std(0)

res = {"data": {"file": NAME, "sha256": SHA256, "doi": "10.17864/1947.000347", "licence": "CC BY 4.0"},
       "polygon_lon_lat": POLY, "masks": {}}
for sd_min in (40, 60, 80):
    m = region & (sd > sd_min)
    series = -(F[:, m] * area[m]).sum(1) / 1e15
    ann = np.array([series[years == y].mean() for y in full])
    r = fit(np.array(full) + 0.5, ann, kmax=19)
    s = summary(r, digits=4)
    s["area_million_km2"] = round(float(area[m].sum()) / 1e12, 2)
    s["cells"] = int(m.sum())
    s["per_unit_area_W_m2_per_decade"] = round(10 * r["slope_per_yr"] * 1e15 / float(area[m].sum()), 2)
    s["decade_mean_first"] = round(float(ann[:10].mean()), 3)
    s["decade_mean_last"] = round(float(ann[-10:].mean()), 3)
    s["decade_change_pct"] = round(100 * (ann[-10:].mean() - ann[:10].mean()) / ann[:10].mean(), 1)
    subs = {}
    for y0, y1 in ((1985, 2004), (2005, 2016)):
        sel = (np.array(full) >= y0) & (np.array(full) <= y1)
        subs[f"{y0}-{y1}"] = summary(fit(np.array(full)[sel] + 0.5, ann[sel], kmax=9), digits=4)
    s["sub_periods"] = subs
    if sd_min == 40:
        s["annual_heat_loss_PW"] = {str(y): round(float(v), 4) for y, v in zip(full, ann)}
    res["masks"][f"sd_over_{sd_min}"] = s
    print(f"\nocean mask: sd > {sd_min} W/m2: {s['cells']} cells, {s['area_million_km2']} million km2;"
          f" mean heat loss {s['mean']:.3f} PW")
    print(f"  trend {full[0]}-{full[-1]}: {s['slope_per_decade']:+.4f} PW/decade, 95% +/-{s['ci95_preprint_method']:.4f}"
          f" (preprint's recipe), +/-{s['ci95_ar1']:.4f} (AR(1)); {s['per_unit_area_W_m2_per_decade']:+.2f} W/m2 per decade")
    print(f"  fitted {s['fit_start']:.3f} -> {s['fit_end']:.3f} PW: {s['change_pct_of_fit_start']:+.1f}%;"
          f" decade means {s['decade_mean_first']:.3f} -> {s['decade_mean_last']:.3f} PW ({s['decade_change_pct']:+.1f}%)")
    for k, v in subs.items():
        print(f"  {k}: {v['slope_per_decade']:+.4f} PW/decade, 95% +/-{v['ci95_preprint_method']:.4f}")

if out_json:
    with open(out_json, "w", encoding="utf-8", newline="\n") as fh:
        json.dump(res, fh, indent=1, ensure_ascii=False)
        fh.write("\n")
    print("\nwrote", out_json)