"""Rebase NASA GISTEMP annual anomalies; stdlib only; not a new climate reconstruction.

Place this script beside nasa-gistemp.csv and run python analyze_baselines.py.
Source: https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv
NASA GISS/GISTEMP v4, snapshot retrieved 2026-10-26. Original reference 1951-1980.
Annual J-D values are already rounded to 0.01 degrees Celsius. Ignore incomplete years.
"""
from pathlib import Path
import csv,json,statistics,hashlib

HERE=Path(__file__).resolve().parent
SOURCE=HERE/'nasa-gistemp.csv'
lines=SOURCE.read_text(encoding='utf-8-sig').splitlines()
rows=list(csv.DictReader(lines[1:]))
annual={int(row['Year']):float(row['J-D']) for row in rows if row['J-D'].strip()!='***'}
periods={'1951-1980':(1951,1980),'1981-2010':(1981,2010),'1991-2020':(1991,2020)}
means={}
for label,(start,end) in periods.items():
    values=[annual[year] for year in range(start,end+1)]
    assert len(values)==30
    means[label]=statistics.mean(values)
series={label:{year:value-offset for year,value in annual.items()} for label,offset in means.items()}
def slope(values):
    xs=[year for year in values if 1970<=year<=2025]
    assert len(xs)==56
    xm=statistics.mean(xs);ym=statistics.mean(values[x] for x in xs)
    return sum((x-xm)*(values[x]-ym) for x in xs)/sum((x-xm)**2 for x in xs)*10
slopes={label:slope(values) for label,values in series.items()}
raw_slope=slope(annual)
assert max(abs(value-raw_slope) for value in slopes.values())<1e-12
for label,values in series.items():
    assert abs(statistics.mean(values[year] for year in range(periods[label][0],periods[label][1]+1)))<1e-12
    assert abs((values[2025]-values[1970])-(annual[2025]-annual[1970]))<1e-12
last=max(annual)
latest=rows[-1]
result={'sourceUrl':'https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv','sourceSha256':hashlib.sha256(SOURCE.read_bytes()).hexdigest(),'retrievedDate':'2026-10-26','unit':'degrees Celsius anomaly','completeAnnualRange':[min(annual),max(annual)],'annualCount':len(annual),'incompleteYearsExcluded':[int(r['Year']) for r in rows if r['J-D'].strip()=='***'],'baselineMeansRelativeToPublished1951_1980':means,'latestCompleteYear':last,'latestPublishedAnomaly':annual[last],'latestRebasedAnomalies':{label:values[last] for label,values in series.items()},'trendPeriod':[1970,2025],'trendUnit':'degrees Celsius per decade, descriptive OLS','rawTrend':raw_slope,'rebasedTrends':slopes,'endpointChange1970_2025':annual[2025]-annual[1970],'checks':{'eachBaselineHas30Years':True,'baselineMeansZeroWithinTolerance':True,'slopesInvariantTolerance':1e-12,'endpointChangeInvariant':True},'limitations':['Published annual values rounded to 0.01 C; rebased values rounded for presentation','No new station/ocean processing, spatial estimates, uncertainty ensemble or forecast','No reconstruction of an 1850-1900 baseline','No comparison equating annual anomalies with monthly anomalies or year-to-date means']}
(HERE/'results.json').write_text(json.dumps(result,indent=2),encoding='utf-8')
with (HERE/'rebased-annual.csv').open('w',newline='',encoding='utf-8') as f:
    writer=csv.writer(f);writer.writerow(['year','published_1951_1980_C']+[label+'_C' for label in periods])
    for year in annual:writer.writerow([year,annual[year]]+[round(series[label][year],8) for label in periods])
print(json.dumps(result,indent=2))