Which estates have the most EV chargers?

"""prep_ev.py: every number on the EV chargers page starts here.

Inputs (all in data/raw/):
  ev.json          LTA DataMall, EV Charging Points Batch file
  pa.geojson       URA Master Plan 2019 Planning Area Boundary (No Sea), data.gov.sg
  pop.json         SingStat table C020123, General Household Survey 2025 (residents by planning area)

    python prep/prep_ev.py --no-download    use the three files already in data/raw/
    python prep/prep_ev.py                  download them first. LTA needs an AccountKey: set LTA_ACCOUNT_KEY, or set LTA_FILE_URL to the
                                            presigned "Link" that the EVCBatch call returns (it works for only a few minutes).

Writes data/ev_numbers.json. page.py counts everything again with plain Python (no shapely) and stops on any difference."""
import json, os, re, sys, urllib.request
from pathlib import Path
from shapely.geometry import shape, Point, box
from shapely.strtree import STRtree

HERE = Path(__file__).resolve().parent.parent
RAW = HERE / "data" / "raw"

MIN_RESIDENTS = 20000      # an area is ranked only if at least this many residents live there
EV_READY_EIGHT = ["Ang Mo Kio", "Bedok", "Choa Chu Kang", "Jurong West", "Punggol", "Queenstown", "Sembawang", "Tengah"]   # LTA news release, March 2021
HDB_NAME = re.compile(r"^BLK |MULTI STOREY CAR PARK")                                                                       # names that point to an HDB block or carpark
LON0, LON1, LAT0, LAT1 = 103.60, 104.09, 1.205, 1.475                                                                       # the map frame
W = 800.0


def fetch(url, dest, headers=None):
    req = urllib.request.Request(url, headers=headers or {"User-Agent": "Mozilla/5.0"})
    with urllib.request.urlopen(req, timeout=120) as r:
        dest.write_bytes(r.read())


def download():
    RAW.mkdir(parents=True, exist_ok=True)
    link = os.environ.get("LTA_FILE_URL")
    if not link:
        key = os.environ["LTA_ACCOUNT_KEY"]
        req = urllib.request.Request("https://datamall2.mytransport.sg/ltaodataservice/EVCBatch", headers={"AccountKey": key, "accept": "application/json"})
        link = json.load(urllib.request.urlopen(req, timeout=60))["value"][0]["Link"]
    fetch(link, RAW / "ev.json")
    u = json.load(urllib.request.urlopen("https://api-open.data.gov.sg/v1/public/api/datasets/d_4765db0e87b9c86336792efe8a1f7a66/poll-download", timeout=60))["data"]["url"]
    fetch(u, RAW / "pa.geojson")
    fetch("https://tablebuilder.singstat.gov.sg/api/table/tabledata/C020123?limit=5000", RAW / "pop.json")


def compute():
    ev = json.load(open(RAW / "ev.json"))
    feats = sorted(json.load(open(RAW / "pa.geojson"))["features"], key=lambda f: f["properties"]["PLN_AREA_N"])
    areas = [(f["properties"]["PLN_AREA_N"].title(), shape(f["geometry"])) for f in feats]
    tree = STRtree([g for _, g in areas])

    def area_of(lon, lat):
        pt = Point(lon, lat)
        for i in tree.query(pt):
            if areas[i][1].contains(pt):
                return areas[i][0], False
        return areas[int(tree.nearest(pt))][0], True      # a point just off the coast: nearest area

    rows = json.load(open(RAW / "pop.json"))["Data"]["row"]
    residents = {}
    for r in rows:
        if r["rowText"].endswith(" - Total"):
            v = [c for c in r["columns"][0]["columns"] if c["key"] == "Total"][0]["value"]
            residents[r["rowText"][:-8].strip().title()] = int(v) if v.isdigit() else 0

    stats = {n: dict(n=n, pts=0, locs=0, hdb=0, dc=0, res=residents.get(n, 0)) for n, _ in areas}
    sites, nearest, unnamed, max_site = [], 0, 0, 0
    for loc in ev["evLocationsData"]:
        a, near = area_of(loc["longtitude"], loc["latitude"])     # (sic) LTA spells the field longtitude
        nearest += near
        pts = dc = 0
        for ch in loc["chargingPoints"]:
            for pt in ch["plugTypes"]:
                k = len(pt["evIds"])                              # one evCpId = one charging point (connector)
                pts += k
                dc += k if pt["current"] == "DC" else 0
        s = stats[a]
        s["pts"] += pts; s["locs"] += 1; s["dc"] += dc
        named_hdb = bool(HDB_NAME.search(loc["name"].upper()))
        s["hdb"] += pts if named_hdb else 0
        unnamed += pts if not loc["name"].strip() else 0
        max_site = max(max_site, pts)
        sites.append((loc["longtitude"], loc["latitude"], pts))

    # ---- map geometry: planning areas simplified, and the locations ----
    kx = W / (LON1 - LON0)
    H = (LAT1 - LAT0) * kx
    proj = lambda lon, lat: ((lon - LON0) * kx, (LAT1 - lat) * kx)
    clip = box(LON0, LAT0, LON1, LAT1)
    out_areas = []
    for n, geom in areas:
        g = geom.intersection(clip).simplify(0.0005, preserve_topology=True)
        polys = [g] if g.geom_type == "Polygon" else list(getattr(g, "geoms", []))
        d = []
        for p in polys:
            if p.geom_type != "Polygon" or p.is_empty:
                continue
            for ring in [p.exterior] + list(p.interiors):
                d.append("M" + "L".join(f"{x:.1f},{y:.1f}" for x, y in (proj(a, b) for a, b in ring.coords)) + "Z")
        if not d:
            continue
        c = g.representative_point(); cx, cy = proj(c.x, c.y)
        s = stats[n]
        out_areas.append(dict(n=n, d="".join(d), cx=round(cx, 1), cy=round(cy, 1), res=s["res"], pts=s["pts"], locs=s["locs"], hdb=s["hdb"], dc=s["dc"]))
    # areas with no drawn outline would break the page's totals
    assert {a["n"] for a in out_areas} == set(stats), "an area has no outline in the map frame"
    locs = []
    for lon, lat, p in sites:
        if LON0 <= lon <= LON1 and LAT0 <= lat <= LAT1:
            x, y = proj(lon, lat); locs.append([round(x, 1), round(y, 1), p])
    return dict(
        updated=ev["LastUpdatedTime"], n_locs=len(sites), n_points=sum(s["pts"] for s in stats.values()), dc_total=sum(s["dc"] for s in stats.values()),
        nearest=nearest, unnamed_pts=unnamed, max_site=max_site, min_residents=MIN_RESIDENTS, ev_ready_eight=EV_READY_EIGHT,
        map=dict(w=W, h=round(H, 1), lon=[LON0, LON1], lat=[LAT0, LAT1]), off_map=len(sites) - len(locs),
        areas=out_areas, locs=locs)


if __name__ == "__main__":
    if "--no-download" not in sys.argv:
        download()
    d = compute()
    (HERE / "data" / "ev_numbers.json").write_text(json.dumps(d, ensure_ascii=False, separators=(",", ":")), encoding="utf-8")
    print("wrote data/ev_numbers.json:", d["n_points"], "points,", d["n_locs"], "locations, list time stamp", d["updated"])