Haze in Singapore by year and part of the island

"""combine_psi.py - joins the two NEA 24-hour PSI sources into one hourly file for the haze page.
Run (inside the build folder):  python3 -I prep/combine_psi.py <out.csv>
  reads data/raw/historical_24hr_psi.csv and data/raw/realtime/*.json
1. Historical file (data.gov.sg, NEA 24-hour PSI, one row an hour, 1 Apr 2014 01:00 to 1 Sep 2026 00:00, 108,864 rows).
   Four rows have a timestamp that does not fit their place in the file (see FIXES). Each is set to the one hour that is
   missing between its neighbours. All four readings are below 70, so no count of hours at 101 or higher changes.
2. Real-time files (api-open.data.gov.sg /v2/real-time/api/psi?date=YYYY-MM-DD, one file for each day, fetched 13 Oct 2026).
   Used from 1 Sep 2026 00:00 on. For 25 to 31 Aug, where the two sources overlap, the differences are counted and printed.
Standard library only."""
import csv, glob, json, os, sys
from datetime import datetime, timedelta
HERE = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
R = ['north', 'south', 'east', 'west', 'central']
out = sys.argv[1]
# _id of the row -> the hour it belongs to (the file is in time order, so the hour is the one missing at that place)
FIXES = {17616: '2016-04-04 00:00',   # written "2016-04-04 010:00:00"
         16956: '2016-03-07 12:00',   # written 7/3/2016 0:00 twice; 12 rows after the first
         35640: '2018-04-25 00:00',   # written 24/4/2018 0:00 twice; 24 rows after the first
         38355: '2018-08-16 03:00'}   # written 16/8/2018 2:00 twice
FMT = '%Y-%m-%d %H:%M'
hist = {}
rows_n = 0
for row in csv.DictReader(open(os.path.join(HERE, 'data/raw/historical_24hr_psi.csv'), newline='')):
    rows_n += 1
    i = int(row['_id'])
    if i in FIXES:
        ts = datetime.strptime(FIXES[i], FMT)
    else:
        ts = datetime.strptime(row['timestamp'], '%d/%m/%Y %H:%M')
    assert ts not in hist, ('two rows for one hour', ts)
    hist[ts] = [int(row[r]) for r in R]
ts_sorted = sorted(hist)
assert all(b - a == timedelta(hours=1) for a, b in zip(ts_sorted, ts_sorted[1:])), 'history is not one row per hour'
rt = {}
for fn in sorted(glob.glob(os.path.join(HERE, 'data/raw/realtime/*.json'))):
    for it in json.load(open(fn))['data']['items']:
        p = it['readings']['psi_twenty_four_hourly']
        rt[datetime.strptime(it['timestamp'][:16].replace('T', ' '), FMT)] = [p[r] for r in R]
overlap = [t for t in rt if t in hist]
diff = {r: sum(1 for t in overlap if hist[t][k] != rt[t][k]) for k, r in enumerate(R)}
print('overlap rows', len(overlap), 'differences', diff)
cut = datetime(2026, 9, 1, 0, 0)
rows = [(t, hist[t], 'historical') for t in ts_sorted if t < cut] + [(t, rt[t], 'realtime') for t in sorted(rt) if t >= cut]
assert all(b[0] - a[0] == timedelta(hours=1) for a, b in zip(rows, rows[1:])), 'a gap between the two sources'
with open(out, 'w', newline='') as f:
    f.write('ts,' + ','.join(R) + ',src\n')
    for t, v, s in rows:
        f.write(t.strftime(FMT) + ':00,' + ','.join(str(x) for x in v) + ',' + s + '\n')
print('rows', len(rows), rows[0][0], rows[-1][0], 'rows from the history', sum(1 for r in rows if r[2] == 'historical'), 'rows from the real-time feed', sum(1 for r in rows if r[2] == 'realtime'))

"""prep_haze.py - builds data/data.json for the CompareSG haze page.
Run (inside the build folder):  python3 -I prep/prep_haze.py <hourly.csv> <out.json>
Input : the hourly file made by prep/combine_psi.py (NEA 24-hour PSI, one row an hour, five regions: the data.gov.sg
        historical file to 31 Aug 2026 23:00, then NEA's real-time feed).
A 'view' is one series with one value an hour: a region, 'any' (the highest of the five readings that hour) or
'all5' (the lowest of the five, so 101 or more means every part was at 101 or more).
For each year, window (the whole year, or 1 Jan to the last day's date) and view it counts the hours at 101 or higher
(and at 201 or higher), the days with such an hour, the highest reading (and the first hour it was reached) and the longest
unbroken run of hours at 101 or higher. Standard library only."""
import csv, json, sys
from datetime import datetime
src, out = sys.argv[1], sys.argv[2]
R = ['north', 'south', 'east', 'west', 'central']
VIEWS = ['any'] + R + ['all5']
rows = []
for r in csv.DictReader(open(src, newline='')):
    t = datetime.strptime(r['ts'], '%Y-%m-%d %H:%M:%S')
    v = {k: int(r[k]) for k in R}
    v['any'] = max(v[k] for k in R)
    v['all5'] = min(v[k] for k in R)
    rows.append({'ts': t, 'year': t.year, 'date': t.strftime('%Y-%m-%d'), 'md': t.month * 100 + t.day, **v})
last = rows[-1]['ts']
assert all(a['ts'] < b['ts'] for a, b in zip(rows, rows[1:]))
CUT = last.month * 100 + last.day
stamp = lambda t: t.strftime('%Y-%m-%d %H:%M')

def one(g, v):
    hi = [x[v] >= 101 for x in g]
    pk = max(x[v] for x in g)
    first = next(x for x in g if x[v] == pk)
    o = {'h': sum(hi), 'v': sum(1 for x in g if x[v] >= 201), 'd': len({x['date'] for x in g if x[v] >= 101}),
         'p': pk, 't': stamp(first['ts'])}
    if v == 'any':
        o['pr'] = [r for r in R if first[r] == pk][0]
    longest = cur = 0
    for f in hi:
        cur = cur + 1 if f else 0
        longest = max(longest, cur)
    o['r'] = longest
    o['open'] = bool(cur > 0 and cur == longest and g[-1]['ts'] == last)
    return o

res = {'last': stamp(last), 'cut': '%02d-%02d' % (CUT // 100, CUT % 100), 'views': VIEWS, 'Y': {}}
for y in sorted({x['year'] for x in rows}):
    g = [x for x in rows if x['year'] == y]
    res['Y'][y] = {'hrs': len(g), 'first': stamp(g[0]['ts']), 'full': {v: one(g, v) for v in VIEWS}}
    gw = [x for x in g if x['md'] <= CUT]
    res['Y'][y]['same'] = {v: one(gw, v) for v in VIEWS}
    res['Y'][y]['hrs_same'] = len(gw)

# the current year day by day, 1 Aug to the last day: the highest reading each day in every view
cal = []
dates = sorted({x['date'] for x in rows if x['year'] == last.year and x['ts'] >= datetime(last.year, 8, 1)})
for dt in dates:
    gg = [x for x in rows if x['date'] == dt]
    cal.append({'d': dt, 'n': len(gg), **{v: max(x[v] for x in gg) for v in VIEWS}})
res['cal'] = cal

# facts for the sentences
cur_y = last.year
a = [x for x in rows if x['year'] == cur_y and x['all5'] >= 101]
b = [x for x in rows if x['year'] < cur_y and x['all5'] >= 101]
res['all5_first_2026'] = stamp(a[0]['ts'])
res['all5_prev_last'] = stamp(b[-1]['ts'])
c = [x for x in rows if x['year'] == cur_y and x['central'] == 192]
res['c192'] = [stamp(c[0]['ts']), stamp(c[-1]['ts']), len(c)]
res['gaps'] = sum(1 for x, y2 in zip(rows, rows[1:]) if (y2['ts'] - x['ts']).total_seconds() != 3600)
json.dump(res, open(out, 'w'), separators=(',', ':'))
print('ok', last, 'gaps in hourly steps:', res['gaps'])