# Dams from space, step 2: builds each DWS dam's full-supply outline from the OpenStreetMap water polygons chosen in
# outline-selection.json (the union when a dam is several polygons), checks it against DWS's capacity (the implied mean
# depth, capacity / area, must be 0.5 to 120 m) and writes outlines.geojson (ODbL: offered on /open with this method).
#   python3 -I data/build/dam-surface-from-space/outlines.py
# Needs the workshop (port 5433) with dams_tmp.osm_water loaded from the pinned extract (outline-match.sql says how).
import json, subprocess, sys
D = 'data/build/dam-surface-from-space'
sel = json.load(open(f'{D}/outline-selection.json'))
recs = json.load(open('data/records/dam-levels.json'))
# Each station's newest weekly record (the release version is no longer a report date, so matching on it found none;
# 13 Nov 2026).
dws = {}
for r in recs['records']:
    if r['station'] not in dws or r['report_date'] > dws[r['station']]['report_date']: dws[r['station']] = r
rows = [(s[0], s[1].replace('|', '/'), '{' + ','.join(str(f) for f in s[2]) + '}') for s in sel['sel']]
import tempfile, os
tmp = tempfile.NamedTemporaryFile('w', suffix='.csv', delete=False)
tmp.write(''.join(f'{a}|"{b}"|{c}\n' for a, b, c in rows)); tmp.close()
sql = f"create temp table s (station text, how text, fids int[]);\n\\copy s from '{tmp.name}' with (format csv, delimiter '|')\n"
sql += """
select json_agg(json_build_object('station', s.station, 'how', s.how,
  'osm', (select json_agg(coalesce(w.osm_id, 'w' || w.osm_way_id) order by w.ogc_fid) from dams_tmp.osm_water w where w.ogc_fid = any(s.fids)),
  'osm_name', (select string_agg(distinct w.name, '; ') from dams_tmp.osm_water w where w.ogc_fid = any(s.fids)),
  'km2', round((st_area(g::geography) / 1e6)::numeric, 4),
  'centre', json_build_array(round(st_x(st_pointonsurface(g))::numeric, 5), round(st_y(st_pointonsurface(g))::numeric, 5)),
  'bbox', json_build_array(round(st_xmin(g)::numeric, 5), round(st_ymin(g)::numeric, 5), round(st_xmax(g)::numeric, 5), round(st_ymax(g)::numeric, 5)),
  'geom', st_asgeojson(st_simplifypreservetopology(g, 0.00005), 5)::json) order by s.station)
from s, lateral (select st_union(w.geom) g from dams_tmp.osm_water w where w.ogc_fid = any(s.fids)) u;
"""
out = subprocess.run(['psql', '-p', '5433', '-d', 'unmapped', '-At', '-q', '-v', 'ON_ERROR_STOP=1'], input=sql, capture_output=True, text=True)
os.unlink(tmp.name)
if out.returncode: sys.exit(out.stderr)
items = json.loads(out.stdout.strip())
feats, bad = [], []
for it in items:
    r = dws[it['station']]
    depth = r['full_supply_capacity'] / it['km2'] if it['km2'] else None  # million m3 / km2 = m
    it['implied_depth_m'] = round(depth, 1) if depth else None
    if depth is None or not 0.5 <= depth <= 120: bad.append((it['station'], r['dam'], it['km2'], r['full_supply_capacity'], it['implied_depth_m']))
    g = it.pop('geom')
    feats.append({'type': 'Feature', 'id': it['station'], 'properties': {'dam': r['dam'], 'river': r['river'], 'capacity_mcm': r['full_supply_capacity'], **it}, 'geometry': g})
fc = {'type': 'FeatureCollection', 'name': 'dam-outlines',
      'attribution': '(c) OpenStreetMap contributors, ODbL 1.0 (https://www.openstreetmap.org/copyright). Derived by Unmapped Dataworks: one outline per DWS dam station, method in outlines.py.',
      'source': 'Geofabrik south-africa-261031.osm.pbf', 'unmatched': sel['unmatched'], 'features': feats}
open(f'{D}/outlines.geojson', 'w').write(json.dumps(fc, separators=(',', ':')))
print(f'{len(feats)} outlines, {len(sel["unmatched"])} unmatched; outside 0.5 to 120 m implied depth: {len(bad)}')
for b in bad: print('  ', b)