import json, math

# Approximate San Francisco land boundary (lon, lat), traced counterclockwise.
# Covers the city proper incl. Presidio, excludes Farallones. Treasure Island added separately.
SF = [
 (-122.5115,37.7790),(-122.5093,37.7860),(-122.5060,37.7900),(-122.4990,37.7905),
 (-122.4880,37.7975),(-122.4790,37.8055),(-122.4720,37.8100),(-122.4650,37.8110),
 (-122.4560,37.8065),(-122.4480,37.8055),(-122.4400,37.8085),(-122.4310,37.8085),
 (-122.4240,37.8070),(-122.4160,37.8085),(-122.4060,37.8080),(-122.4000,37.8055),
 (-122.3935,37.7975),(-122.3880,37.7930),(-122.3830,37.7860),(-122.3865,37.7800),
 (-122.3860,37.7730),(-122.3890,37.7680),(-122.3860,37.7600),(-122.3785,37.7550),
 (-122.3760,37.7480),(-122.3660,37.7395),(-122.3555,37.7240),(-122.3760,37.7130),
 (-122.3830,37.7080),(-122.3960,37.7080),(-122.4080,37.7085),(-122.4240,37.7090),
 (-122.4400,37.7080),(-122.4560,37.7085),(-122.4680,37.7060),(-122.4790,37.7080),
 (-122.4880,37.7130),(-122.4930,37.7180),(-122.5040,37.7250),(-122.5090,37.7400),
 (-122.5105,37.7600),(-122.5115,37.7790),
]
TI = [(-122.3760,37.8010),(-122.3630,37.8010),(-122.3610,37.8280),(-122.3730,37.8290),(-122.3760,37.8010)]

def inside(pt, poly):
    x,y = pt; n=len(poly); c=False; j=n-1
    for i in range(n):
        xi,yi = poly[i]; xj,yj = poly[j]
        if ((yi>y)!=(yj>y)) and (x < (xj-xi)*(y-yi)/(yj-yi)+xi): c = not c
        j=i
    return c

# Hex grid, ~400m spacing (flat-top hexes, centers on a staggered lattice)
LAT0 = 37.758
M_PER_DEG_LAT = 110574.0
M_PER_DEG_LON = 111320.0*math.cos(math.radians(LAT0))
SPACING = 200.0                      # center-to-center, m
dx = SPACING/M_PER_DEG_LON           # horizontal step (deg lon)
dy = SPACING*math.sqrt(3)/2/M_PER_DEG_LAT  # row step (deg lat)

lon_min,lon_max = -122.5135,-122.3540
lat_min,lat_max =  37.7050, 37.8320

cells=[]; row=0; lat=lat_min
while lat <= lat_max:
    off = (dx/2) if row%2 else 0.0
    lon = lon_min+off
    while lon <= lon_max:
        if inside((lon,lat),SF) or inside((lon,lat),TI):
            cells.append({"id":len(cells),"lon":round(lon,6),"lat":round(lat,6),"row":row})
        lon += dx
    lat += dy; row += 1

json.dump({"spacing_m":SPACING,"dx":dx,"dy":dy,"cells":cells}, open("grid200.json","w"))
print("cells:", len(cells))

# Emit batched Open-Elevation query strings (~250 pts per URL to stay under URL length limits)
B=250
batches=[cells[i:i+B] for i in range(0,len(cells),B)]
with open("elev_urls200.txt","w") as f:
    for b in batches:
        f.write("https://api.open-elevation.com/api/v1/lookup?locations="+"|".join(f'{c["lat"]},{c["lon"]}' for c in b)+"\n")
print("elevation batches:", len(batches))