build_geojson.py

# -*- coding: utf-8 -*-
"""Build a web-ready GeoJSON of the 909 public REGA non-Saudi-ownership zones.
Source: rega_zones_api_raw.json (the REGA "Saudi Properties" portal API).
Rounds coordinates to 5 decimals (~1m). Area shown is REGA's OFFICIAL zoneArea
field, unit-normalized to km2 (see rega_area_km2) — never our own recomputation.

Paths resolve relative to this file, or override via env var:
  REGA_DATA_DIR  directory holding the raw REGA API dump (default: ./data)
"""
import json, math, os
HERE = os.path.dirname(os.path.abspath(__file__))
DATA_DIR = os.environ.get('REGA_DATA_DIR', os.path.join(HERE, 'data'))
SRC = os.path.join(DATA_DIR, 'rega_zones_api_raw.json')
OUT = os.path.join(HERE, 'public', 'zones.geojson')
raw=json.load(open(SRC,encoding='utf-8'))
def lst(k):
    v=raw[k]['data']; return (v.get('data') if isinstance(v,dict) else v) or []
regions={r['id']:r['name'] for r in lst('regions')}
cats={m['id']:m.get('mainZoneName') for m in lst('mainZonesGroups')}
sz=lst('subZones')

R=6371.0088
def ring_area(ring):
    if len(ring)<4: return 0.0
    s=0.0
    for i in range(len(ring)-1):
        lo1,la1=ring[i]; lo2,la2=ring[i+1]
        s+=math.radians(lo2-lo1)*(2+math.sin(math.radians(la1))+math.sin(math.radians(la2)))
    return abs(s*R*R/2.0)
def mp_area(coords):
    a=0.0
    for poly in coords:
        if poly:
            a+=ring_area(poly[0])-sum(ring_area(h) for h in poly[1:])
    return a
def rega_area_km2(zone_area, geom_km2):
    """Return REGA's official zoneArea expressed in km2.
    REGA stores zoneArea in MIXED units per record (m2 for some zones, km2 for
    others, even within the same category), so the raw value is ambiguous by a
    factor of 1e6. We pick the interpretation whose magnitude matches the polygon
    REGA also published. The geometry is used ONLY to disambiguate the unit; the
    value returned is REGA's own figure, not a recomputation. Unknown -> None."""
    if zone_area is None or geom_km2 <= 0:
        return None
    return min((zone_area, zone_area/1e6),
               key=lambda v: abs(math.log((v+1e-12)/(geom_km2+1e-12))))
def rnd(coords):
    return [[[ [round(x,5),round(y,5)] for x,y in ring] for ring in poly] for poly in coords]

# friendly short category labels + ids for styling
CATLABEL={1:'Urban boundaries',2:'Giga / mega-projects',3:'Cities & economic zones',
 4:'Riyadh City',5:'Jeddah Governorate',6:'Makkah City',7:'Al Madinah',8:'AlUla Governorate'}

feats=[]; bycat={}
for z in sz:
    g=z.get('geometry')
    if not g or g.get('type')!='MultiPolygon': continue
    cid=z.get('mainZoneCategoryId')
    # Show REGA's official area (unit-normalized), never our own recomputation.
    area=rega_area_km2(z.get('zoneArea'), mp_area(g['coordinates']))
    area=round(area,6) if area is not None else None  # 6 dp = ~1 m2, keeps small zones exact
    bycat[cid]=bycat.get(cid,0)+1
    feats.append({"type":"Feature",
      "properties":{"name":z.get('name'),"nameAr":z.get('nameAr'),"catId":cid,
        "cat":CATLABEL.get(cid,cats.get(cid)),"region":regions.get(z.get('regionId')),
        "code":z.get('zoneCode'),"areaKm2":area},
      "geometry":{"type":"MultiPolygon","coordinates":rnd(g['coordinates'])}})

fc={"type":"FeatureCollection","features":feats}
with open(OUT,'w',encoding='utf-8') as f:
    json.dump(fc,f,ensure_ascii=False,separators=(',',':'))
import os
print("WROTE",OUT,"| features:",len(feats),"| size MB:",round(os.path.getsize(OUT)/1e6,2))
print("by category:",{CATLABEL[k]:v for k,v in sorted(bycat.items())})