krz/omaha-metro-blotter

 1"""Refresh raw_data/boundaries.geojson, the map's geographic context.
 2
 3City limits and the county line are static reference geometry, not archive
 4data, so this is a rare manual refresh rather than part of the daily pull.
 5Douglas County publishes its own city limits and boundary; Sarpy publishes
 6municipal boundaries under a different field name. Council Bluffs sits in
 7Pottawattamie County, Iowa, which publishes neither, so it is unoutlined.
 8"""
 9
10import json
11import ssl
12import urllib.parse
13import urllib.request
14from pathlib import Path
15
16import certifi
17
18OUT = Path(__file__).parent / "raw_data" / "boundaries.geojson"
19CTX = ssl.create_default_context(cafile=certifi.where())
20# ~55 m; enough to keep the shapes honest without shipping full survey detail
21TOLERANCE = 0.0005
22
23LAYERS = [
24    {"kind": "county", "name": "Douglas County", "name_field": None,
25     "url": "https://services.arcgis.com/pDAi2YK0L0QxVJHj/arcgis/rest/services"
26            "/Douglas_County_Boundary/FeatureServer/0"},
27    {"kind": "city", "name_field": "town",
28     "url": "https://dcgis.org/server/rest/services/Hosted"
29            "/City_Limits_(source)_view/FeatureServer/0"},
30    {"kind": "city", "name_field": "NAME",
31     "url": "https://geodata.sarpy.gov/arcgis/rest/services/Cadastral"
32            "/LandRecordsDynamic/MapServer/39"},
33]
34
35
36def fetch(layer):
37    params = {"where": "1=1", "outFields": "*", "outSR": "4326",
38              "maxAllowableOffset": TOLERANCE, "f": "geojson"}
39    url = layer["url"] + "/query?" + urllib.parse.urlencode(params)
40    # Sarpy's MapServer 403s the default Python-urllib agent
41    req = urllib.request.Request(url, headers={"User-Agent": "omaha-metro-blotter/1.0"})
42    with urllib.request.urlopen(req, timeout=120, context=CTX) as r:
43        return json.load(r)
44
45
46def main():
47    features = []
48    for layer in LAYERS:
49        gj = fetch(layer)
50        for f in gj.get("features", []):
51            field = layer["name_field"]
52            name = (f.get("properties") or {}).get(field) if field else layer["name"]
53            features.append({
54                "type": "Feature",
55                "properties": {"name": name, "kind": layer["kind"]},
56                "geometry": f["geometry"],
57            })
58        print(f"  {len(gj.get('features', []))} from {layer['url'].split('/')[2]}")
59
60    OUT.write_text(json.dumps(
61        {"type": "FeatureCollection", "features": features},
62        separators=(",", ":")))
63    verts = sum(len(r) for f in features
64                for poly in ([f["geometry"]["coordinates"]]
65                             if f["geometry"]["type"] == "Polygon"
66                             else f["geometry"]["coordinates"])
67                for r in poly)
68    print(f"  wrote {OUT.relative_to(OUT.parent.parent)}: "
69          f"{len(features)} features, {verts} vertices, "
70          f"{OUT.stat().st_size / 1024:.0f} KiB")
71
72
73if __name__ == "__main__":
74    main()