12 Essential Python Libraries for Geospatial Data Analysis (with Hands-On Examples)

import os
import requests
import geopandas as gpd

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

# 1. Fetch restaurants within 2km of Berlin center
params = {
    "categories": "catering.restaurant",
    "filter": "circle:13.4050,52.5200,2000",
    "limit": 200,
    "apiKey": API_KEY,
}
resp = requests.get("https://api.geoapify.com/v2/places", params=params)
data = resp.json()

# 2. Convert to GeoDataFrame (WGS84, EPSG:4326)
gdf = gpd.GeoDataFrame.from_features(data["features"], crs="EPSG:4326")

# 3. Reproject to a projected CRS (UTM Zone 33N, meters)
gdf_utm = gdf.to_crs(epsg=32633)

# 4. Build a 1 km buffer around each restaurant
buffers = gdf_utm.copy()
buffers["geometry"] = buffers.buffer(1000)  # 1000 meters

# 5. Spatial join: how many restaurants fall within each buffer
joined = gpd.sjoin(gdf_utm, buffers, predicate="within")
density = joined.groupby(joined.index_right).size()

# 6. Attach density counts back to restaurants
gdf_utm["nearby_restaurants"] = density.values

# 7. Inspect results
print(gdf_utm[["name", "nearby_restaurants"]].head())

# (Optional) Plot restaurants, sized by density
ax = gdf_utm.plot(
    column="nearby_restaurants",
    cmap="OrRd",
    legend=True,
    figsize=(8, 8),
    markersize=50,
    alpha=0.7
)
ax.set_title("Restaurant density within 1 km (Berlin center)")

import requests, fiona
from fiona.crs import CRS  # optional, but nice for explicit CRS

API_KEY = "YOUR_API_KEY" # Sign up for free at https://www.geoapify.com/ and replace with your Geoapify API key.

# Fetch a 15-minute driving isochrone around Berlin
params = {"lat": 52.52, "lon": 13.405, "type": "time", "mode": "drive", "range": 900, "apiKey": API_KEY}
isojson = requests.get("https://api.geoapify.com/v1/isoline", params=params, timeout=30).json()

# Assume polygon isoline; extract properties/geometry from the API response
features = isojson["features"]
# Build a schema that matches the GeoJSON properties and geometry type
# (You can introspect the first feature)
geom_type = features[0]["geometry"]["type"]          # e.g., 'Polygon' or 'MultiPolygon'
props = features[0].get("properties", {})            # dict of properties
schema = {
    "geometry": geom_type,
    "properties": {k: "str" if isinstance(v, str) else "float" if isinstance(v, float) else "int"
                   for k, v in props.items()}
}

# Write GeoJSON with Fiona
with fiona.open(
    "berlin_drive_15min_fiona.geojson",
    mode="w",
    driver="GeoJSON",
    schema=schema,
    crs=CRS.from_epsg(4326)  # Geoapify isolines are in WGS84
) as sink:
    sink.writerecords(features)

# Read it back
with fiona.open("berlin_drive_15min_fiona.geojson", "r") as src:
    print("Driver:", src.driver)     # 'GeoJSON'
    print("Features:", len(list(src)))

import os
import gzip
import json
import requests
from pathlib import Path

import rasterio
from rasterio.mask import mask
from shapely.geometry import shape, mapping

API_KEY = "YOUR_API_KEY" # Sign up for free at https://www.geoapify.com/ and replace with your Geoapify API key.

# 1) Download SRTM tile (HGT)
data_dir = Path("data")
data_dir.mkdir(exist_ok=True)
zip_path = data_dir / "N52E013.hgt.gz"
if not zip_path.exists():
    resp = requests.get(
        "https://s3.amazonaws.com/elevation-tiles-prod/skadi/N52/N52E013.hgt.gz",
        timeout=60,
    )
    resp.raise_for_status()
    zip_path.write_bytes(resp.content)

# Decompress .hgt
hgt_path = data_dir / "N52E013.hgt"
if not hgt_path.exists():
    with gzip.open(zip_path, "rb") as src, hgt_path.open("wb") as dst:
        dst.write(src.read())

# 2) Get 15-minute drive-time polygon
params = {"lat": 52.52, "lon": 13.405, "type": "time", "mode": "drive", "range": 900, "apiKey": API_KEY}
isojson = requests.get("https://api.geoapify.com/v1/isoline", params=params, timeout=30).json()
polygon = shape(isojson["features"][0]["geometry"])  # shapely geometry

# 3) Clip DEM to isochrone
with rasterio.open(hgt_path) as src:
    # Use GeoJSON-like dicts for mask:
    out_image, out_transform = mask(src, [mapping(polygon)], crop=True)
    out_meta = src.meta.copy()

    # IMPORTANT: write as GeoTIFF, not HGT; update dimensions/transform
    out_meta.update({
        "driver": "GTiff",
        "height": out_image.shape[1],
        "width":  out_image.shape[2],
        "transform": out_transform,
        # Optional niceties:
        "compress": "deflate",
        "tiled": True
    })

    # Mean elevation ignoring masked areas
    # (mask() returns a masked array; use .mask to ignore outside region)
    mean_elev = float(out_image[out_image != src.nodata].mean())
    print(f"Mean elevation within 15-minute drive: {mean_elev:.1f} m")

# 4) Save clipped DEM as GeoTIFF
clipped_path = data_dir / "berlin_isochrone_elevation.tif"
with rasterio.open(clipped_path, "w", **out_meta) as dest:
    dest.write(out_image)

with rasterio.open(clipped_path) as src:
    array = src.read(1)  # first band

    plt.imshow(array, cmap="terrain", vmin=0, vmax=100)  # fix scale 0–100
    cbar = plt.colorbar(label="Elevation (m)")
    cbar.set_ticks([0, 20, 40, 60, 80, 100])  # optional: custom ticks
    plt.title("Elevation within 15-min drive area (Berlin)")
    plt.show()

print(f"Clipped DEM written to {clipped_path}")

import os
import requests
import pandas as pd
from pyproj import CRS, Transformer

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

addresses = [
    "Brandenburg Gate, Berlin",
    "Museum Island, Berlin",
    "Berlin Central Station",
]

# 1. Geocode addresses with Geoapify (returns WGS84 coordinates)
results = []
for query in addresses:
    params = {"text": query, "format": "json", "apiKey": API_KEY}
    resp = requests.get("https://api.geoapify.com/v1/geocode/search", params=params, timeout=15)
    resp.raise_for_status()
    hit = resp.json()["results"][0]
    results.append({
        "query": query,
        "lat": hit["lat"],
        "lon": hit["lon"],
    })

df = pd.DataFrame(results)
print(df)

# 2. Prepare CRS objects
crs_wgs84 = CRS.from_epsg(4326)
crs_web_mercator = CRS.from_epsg(3857)
crs_utm33 = CRS.from_epsg(32633)

# 3. Build transformers (always use `always_xy=True` for lon/lat order)
transformer_web = Transformer.from_crs(crs_wgs84, crs_web_mercator, always_xy=True)
transformer_utm = Transformer.from_crs(crs_wgs84, crs_utm33, always_xy=True)

# 4. Apply transformations
x3857, y3857 = transformer_web.transform(df["lon"].values, df["lat"].values)
xutm, yutm = transformer_utm.transform(df["lon"].values, df["lat"].values)

df["x_3857"] = x3857
df["y_3857"] = y3857
df["x_utm"] = xutm
df["y_utm"] = yutm

print(df[["query", "x_3857", "y_3857", "x_utm", "y_utm"]])

import os
import requests
from rtree import index
from shapely.geometry import shape, Point

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

# 1. Fetch cafés within 1.5 km of Union Square (New York City)
params = {
    "categories": "catering.cafe",
    "filter": "circle:-73.9903,40.7359,1500",
    "limit": 500,
    "apiKey": API_KEY,
}
resp = requests.get("https://api.geoapify.com/v2/places", params=params, timeout=30)
resp.raise_for_status()
features = resp.json()["features"]

# 2. Build an Rtree spatial index
rtree_idx = index.Index()
point_geoms = []
properties = []
for idx_id, feature in enumerate(features):
    geom = shape(feature["geometry"])
    point_geoms.append(geom)
    properties.append(feature.get("properties", {}))
    rtree_idx.insert(idx_id, geom.bounds)

print(f"Loaded {len(point_geoms)} cafés into the index")

# 3. Query the 5 nearest cafés to a user location (14th St–Union Sq Station as example)
user_location = Point(-73.9901, 40.7359)
nearest_ids = list(rtree_idx.nearest(user_location.bounds, 5))

for rank, idx_id in enumerate(nearest_ids, start=1):
    prop = properties[idx_id]
    geom = point_geoms[idx_id]
    distance = user_location.distance(geom) * 111_139  # rough degrees-to-meters factor near NYC
    print(f"#{rank}: {prop.get('name', 'Unnamed cafe')} — ~{distance:.0f} m away")

import os
import requests
import folium

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

workspaces = [
    "Brooklyn Navy Yard, NY",
    "Industry City, Brooklyn",
    "New Lab, Brooklyn",
]

# 1. Geocode locations
geocoded = []
for place in workspaces:
    params = {"text": place, "format": "json", "apiKey": API_KEY}
    resp = requests.get("https://api.geoapify.com/v1/geocode/search", params=params, timeout=15)
    resp.raise_for_status()
    result = resp.json()["results"][0]
    geocoded.append({"name": place, "lat": result["lat"], "lon": result["lon"]})

# 2. Request a 20-minute walking isochrone around Brooklyn Navy Yard
iso_params = {
    "lat": geocoded[0]["lat"],
    "lon": geocoded[0]["lon"],
    "type": "time",
    "mode": "walk",
    "range": 1200,
    "apiKey": API_KEY,
}
iso_resp = requests.get("https://api.geoapify.com/v1/isoline", params=iso_params, timeout=30)
iso_geojson = iso_resp.json()

# 3. Build Folium map
center = [geocoded[0]["lat"], geocoded[0]["lon"]]
folium_map = folium.Map(
    location=center,
    zoom_start=13,
    tiles=f"https://maps.geoapify.com/v1/tile/maptiler-3d/{{z}}/{{x}}/{{y}}.png?apiKey={API_KEY}",
    attr="© Geoapify, © OpenMapTiles, © OpenStreetMap contributors"
)

# Add isochrone GeoJSON with styling
folium.GeoJson(
    iso_geojson,
    name="10-min walk",
    style_function=lambda feature: {
        "fillColor": "#1f78b4",
        "color": "#1f78b4",
        "weight": 2,
        "fillOpacity": 0.2,
    },
    tooltip=folium.GeoJsonTooltip(fields=["range"], aliases=["ISO-Time (s)"])
).add_to(folium_map)

# Add markers for each workspace
for place in geocoded:
    # Build a Geoapify marker URL (for example, material icon "work")
    marker_url = (
        f"https://api.geoapify.com/v2/icon/?type=awesome&color=%232ed375&size=50&icon=work&iconType=material&contentSize=15&noShadow&scaleFactor=2&apiKey={API_KEY}"
    )

    folium.Marker(
        [place["lat"], place["lon"]],
        popup=folium.Popup(place["name"], offset=(0, -15)),
        icon=folium.CustomIcon(marker_url, icon_size=(39, 50)),  # choose size
    ).add_to(folium_map)

folium.LayerControl().add_to(folium_map)
folium_map.save("brooklyn-workspaces.html")
print("Map saved to brooklyn-workspaces.html")

import os
import requests
import pandas as pd
import plotly.graph_objects as go
from shapely.geometry import shape

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

# 1. Fetch art galleries within 1 km of the High Line (Chelsea, NYC)
params = {
    "categories": "entertainment.culture.gallery",
    "filter": "circle:-74.0047,40.7480,1000",
    "limit": 200,
    "apiKey": API_KEY,
}
resp = requests.get("https://api.geoapify.com/v2/places", params=params, timeout=30)
resp.raise_for_status()
features = resp.json()["features"]

data = [
    {
        "name": feat["properties"].get("name", "Unnamed gallery"),
        "lat": feat["geometry"]["coordinates"][1],
        "lon": feat["geometry"]["coordinates"][0],
        "rating": feat["properties"].get("rating", "n/a"),
        "address": feat["properties"].get("formatted", "")
    }
    for feat in features
]
galleries = pd.DataFrame(data)

# 2. Request a 10-minute biking isochrone around the same point
iso_params = {
    "lat": 40.7480,
    "lon": -74.0047,
    "type": "time",
    "mode": "bicycle",
    "range": 600,
    "apiKey": API_KEY,
}
isojson = requests.get("https://api.geoapify.com/v1/isoline", params=iso_params, timeout=30).json()

geom = shape(isojson["features"][0]["geometry"])

if isinstance(geom, MultiPolygon):
    geom = max(geom.geoms, key=lambda g: g.area)

x, y = geom.exterior.coords.xy
lat, lon = list(y), list(x)

# 3. Build Plotly figure
fig = go.Figure()

# Isochrone polygon
fig.add_trace(
    go.Scattermap(
        lat=lat,
        lon=lon,
        mode="lines",
        line=dict(color="#1f78b4", width=2),
        fill="toself",
        fillcolor="rgba(31,120,180,0.2)",
        name="10-min bike reach"
    )
)

# Gallery points
fig.add_trace(
    go.Scattermap(
        lat=galleries["lat"],
        lon=galleries["lon"],
        mode="markers",
        marker=dict(size=9, color="#e31a1c"),
        name="Art galleries",
        text=galleries["name"],
        hovertemplate="<b>%{text}</b><br>Rating: %{customdata[0]}<br>%{customdata[1]}<extra></extra>",
        customdata=list(zip(galleries["rating"], galleries["address"]))
    )
)

# Layout (note: use map=, not mapbox=)
fig.update_layout(
    map=dict(
        style="open-street-map",  # can also use Geoapify tiles here
        center=dict(lat=40.7480, lon=-74.0047),
        zoom=14
    ),
    margin=dict(l=0, r=0, t=40, b=0),
    title="Chelsea art galleries within a 10-minute bike ride"
)

fig.write_html("chelsea-galleries.html", auto_open=False)
print("Interactive map saved to chelsea-galleries.html")

import os
import requests
from sqlalchemy import create_engine, Column, Integer, String
from sqlalchemy.orm import declarative_base, sessionmaker
from geoalchemy2 import Geometry
from geoalchemy2.shape import from_shape
from shapely.geometry import Point

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
DATABASE_URL = os.environ.get("DATABASE_URL")  # e.g., postgresql+psycopg://user:pass@localhost:5432/geodb

if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")
if not DATABASE_URL:
    raise RuntimeError("Set DATABASE_URL to your PostGIS connection string")

engine = create_engine(DATABASE_URL, echo=False, future=True)
SessionLocal = sessionmaker(bind=engine)
Base = declarative_base()

class Workspace(Base):
    __tablename__ = "workspaces"
    id = Column(Integer, primary_key=True)
    name = Column(String, nullable=False)
    address = Column(String, nullable=False)
    geom = Column(Geometry("POINT", srid=4326), nullable=False)

Base.metadata.create_all(engine)

workspace_names = [
    "WeWork Times Square, New York",
    "Industrious Bryant Park, New York",
    "The Yard: Herald Square, New York",
]

with SessionLocal() as session:
    for name in workspace_names:
        params = {"text": name, "format": "json", "apiKey": API_KEY}
        resp = requests.get("https://api.geoapify.com/v1/geocode/search", params=params, timeout=15)
        resp.raise_for_status()
        result = resp.json()["results"][0]
        point = Point(result["lon"], result["lat"])

        workspace = Workspace(
            name=name,
            address=result.get("formatted", name),
            geom=from_shape(point, srid=4326),
        )
        session.add(workspace)

    session.commit()

with SessionLocal() as session:
    times_square = Point(-73.9855, 40.7580)
    query = (
        session.query(Workspace)
        .filter(Workspace.geom.ST_DWithin(from_shape(times_square, srid=4326), 500))  # meters
    )

    print("Workspaces within 500 m of Times Square:")
    for ws in query:
        print(f"- {ws.name} @ {ws.address}")

import os
import requests
import geopandas as gpd
import libpysal
from esda import Moran_Local

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

# 1) Fetch POIs (police + community centers) around Chicago (10 km radius)
params = {
    "categories": "service.police,activity.community_center",
    "filter": "circle:-87.6298,41.8781,10000",
    "limit": 500,
    "apiKey": API_KEY,
}
resp = requests.get("https://api.geoapify.com/v2/places", params=params, timeout=30)
resp.raise_for_status()
features = resp.json()["features"]

pois = gpd.GeoDataFrame.from_features(features, crs="EPSG:4326")

# 2) Project to a planar CRS in meters (Illinois East / EPSG:26971)
pois = pois.to_crs(epsg=26971)

# 3) Create a fishnet grid (1 km x 1 km cells) covering the POI extent
cell_size = 1000  # meters
minx, miny, maxx, maxy = pois.total_bounds

# snap bounds to the grid
minx = np.floor(minx / cell_size) * cell_size
miny = np.floor(miny / cell_size) * cell_size
maxx = np.ceil(maxx / cell_size) * cell_size
maxy = np.ceil(maxy / cell_size) * cell_size

xs = np.arange(minx, maxx, cell_size)
ys = np.arange(miny, maxy, cell_size)

cells = [box(x, y, x + cell_size, y + cell_size) for x in xs for y in ys]
fishnet = gpd.GeoDataFrame(geometry=cells, crs=pois.crs)

# 4) Spatial join: count POIs per grid cell
joined = gpd.sjoin(fishnet, pois, how="left", predicate="contains")
counts = joined.groupby(level=0).size()              # group by fishnet index
fishnet["poi_count"] = counts.reindex(fishnet.index, fill_value=0).astype(int)

# 5) Build spatial weights (queen contiguity) and run Local Moran’s I
w = libpysal.weights.Queen.from_dataframe(fishnet)
w.transform = "r"

y = fishnet["poi_count"].values
moran_local = Moran_Local(y, w)

fishnet["local_I"] = moran_local.Is
fishnet["significant"] = moran_local.p_sim < 0.05
fishnet["cluster_type"] = moran_local.q  # 1=HH, 2=LH, 3=LL, 4=HL

print(fishnet[["poi_count", "local_I", "significant", "cluster_type"]].head())

import os
import requests
import osmnx as ox
import networkx as nx
from shapely.geometry import shape, MultiPolygon

API_KEY = os.getenv("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

# 1) Download pedestrian network (use a less ambiguous name)
place = "South Boston, Massachusetts, USA"
G = ox.graph_from_place(place, network_type="walk", simplify=True)

# 2) 15-minute walk isochrone from Geoapify
center_lat, center_lon = 42.3510, -71.0410
iso_params = {"lat": center_lat, "lon": center_lon, "type": "time", "mode": "walk", "range": 900, "apiKey": API_KEY}
isojson = requests.get("https://api.geoapify.com/v1/isoline", params=iso_params, timeout=30).json()

geom = shape(isojson["features"][0]["geometry"]).buffer(0)
if isinstance(geom, MultiPolygon):  # pick the largest piece if multipolygon
    geom = max(geom.geoms, key=lambda g: g.area)

# 3) Clip the network by the isochrone (OSMnx 2.x API)
subgraph = ox.truncate.truncate_graph_polygon(G, geom, truncate_by_edge=True)
subgraph = ox.truncate.largest_component(subgraph)

# 4) Edge betweenness centrality (OSMnx recommends using NetworkX for this)
#    Collapse parallel edges by minimizing 'length', then compute on an undirected view.
D = ox.convert.to_digraph(subgraph, weight="length")
Gu = D.to_undirected()
edge_bc = nx.edge_betweenness_centrality(Gu, weight="length", normalized=True)

# 5) Map results back to the edges GeoDataFrame
edges = ox.convert.graph_to_gdfs(subgraph, nodes=False).reset_index()
edges["pair"] = edges.apply(lambda r: tuple(sorted((r["u"], r["v"]))), axis=1)
edges["betweenness"] = edges["pair"].map(edge_bc).fillna(0.0)

print(edges[["u", "v", "length", "betweenness"]].head())

import os
import math
import gzip
import requests
from pathlib import Path

import numpy as np
import xarray as xr
import rioxarray as rxr
from shapely.geometry import shape, MultiPolygon, mapping
from rioxarray.merge import merge_arrays
import matplotlib.pyplot as plt

API_KEY = os.getenv("GEOAPIFY_API_KEY")

# --- 1) Get Boston boundary from Geoapify Place Details by OSM id/type ---
# OSM relation 2315704 = City of Boston
params = {
    "osm_id": 2315704,
    "osm_type": "r",
    # Ask for original polygon geometry (adds a feature alongside 'details')
    "features": "details,details.full_geometry",
    "apiKey": API_KEY,
}
resp = requests.get("https://api.geoapify.com/v2/place-details", params=params, timeout=30)
resp.raise_for_status()
fc = resp.json()["features"]

# Prefer the full/original geometry if present, otherwise fallback to details geometry
feat_full = next((f for f in fc if f["properties"].get("feature_type") == "details.full_geometry"), None)
feat_detl = next((f for f in fc if f["properties"].get("feature_type") == "details"), None)
if not (feat_full or feat_detl):
    raise RuntimeError("No boundary geometry returned by Place Details API.")

geom = shape((feat_full or feat_detl)["geometry"])
if isinstance(geom, MultiPolygon):
    geom = max(geom.geoms, key=lambda g: g.area)  # take largest part

# --- 2) Figure out which SRTM tiles we need (1°x1° HGT) and download them ---
def tile_id(sw_lat, sw_lon):
    lat_prefix = "N" if sw_lat >= 0 else "S"
    lon_prefix = "E" if sw_lon >= 0 else "W"
    return f"{lat_prefix}{abs(int(sw_lat)):02d}{lon_prefix}{abs(int(sw_lon)):03d}"

minx, miny, maxx, maxy = geom.bounds  # lon/lat bbox

lat_floors = range(math.floor(miny), math.floor(maxy) + 1)
lon_floors = range(math.floor(minx), math.floor(maxx) + 1)

data_dir = Path("data"); data_dir.mkdir(exist_ok=True)
hgt_paths = []

for lat0 in lat_floors:
    for lon0 in lon_floors:
        tid = tile_id(lat0, lon0)              # e.g., N42W071
        folder = tid[:3]                        # e.g., N42
        url = f"https://s3.amazonaws.com/elevation-tiles-prod/skadi/{folder}/{tid}.hgt.gz"
        gz_path = data_dir / f"{tid}.hgt.gz"
        hgt_path = data_dir / f"{tid}.hgt"
        
        if not gz_path.exists():
            r = requests.get(url, timeout=60)
            r.raise_for_status()
            gz_path.write_bytes(r.content)
        if not hgt_path.exists():
            with gzip.open(gz_path, "rb") as src, open(hgt_path, "wb") as dst:
                dst.write(src.read())
        hgt_paths.append(hgt_path)

if not hgt_paths:
    raise RuntimeError("No HGT tiles determined from boundary bbox.")

# --- 3) Open tiles as DataArrays, mosaic, and clip to Boston boundary ---
void_vals = (-32768, 32767, -9999, 65535)  # common SRTM/HGT sentinels
arrays = []
for p in hgt_paths:
    da = rxr.open_rasterio(p).squeeze(drop=True).astype("float32")  # -> dims ('y','x'), allow NaNs
    # Ensure CRS + nodata are set
    if da.rio.crs is None:
        da = da.rio.write_crs("EPSG:4326")
    # mask voids and clamp to a plausible elevation range
    da = da.where(~da.isin(void_vals))
    da = da.where((da > -500) & (da < 9000))
    
    arrays.append(da)

mosaic = merge_arrays(arrays)  # single DataArray in EPSG:4326
mosaic = mosaic.where(~mosaic.isin(void_vals))  # belt & suspenders
clipped = mosaic.rio.clip([mapping(geom)], crs="EPSG:4326")

# Mask NODATA explicitly (common SRTM sentinel)
clipped = clipped.where(clipped != clipped.rio.nodata)

# --- 4) Quick stats and preview ---
stats = dict(
    count=int(np.isfinite(clipped.values).sum()),
    min=float(np.nanmin(clipped.values)),
    mean=float(np.nanmean(clipped.values)),
    max=float(np.nanmax(clipped.values)),
)
print("Boston DEM stats (after masking):", stats)

# --- 5) Save outputs ---
# 5a) Single-band (grayscale) GeoTIFF with proper georeferencing
singleband_path = data_dir / "boston_dem_singleband.tif"
clipped.rio.to_raster(singleband_path, compress="deflate")
print("Wrote single-band DEM:", singleband_path)

# 5b) Colorized RGB GeoTIFF so external apps show the same colors

# Plot (color) in notebook
plt.figure(figsize=(8, 6))
# robust=True trims color limits to 2–98 percentile so extremes don't wash out the color
ax = clipped.plot(cmap="terrain", robust=True)
plt.title("Boston SRTM DEM (colorized)")
plt.tight_layout()
plt.show()

vmin, vmax = np.nanpercentile(clipped.values, [2, 98])
norm = plt.Normalize(vmin=vmin, vmax=vmax)
cmap = plt.get_cmap("terrain")

rgba = cmap(norm(clipped.values))            # H x W x 4
rgb = (rgba[..., :3] * 255).astype("uint8") # 0..255 uint8

rgb_path = data_dir / "boston_dem_rgb.tif"
with rasterio.open(
    rgb_path, "w",
    driver="GTiff",
    height=rgb.shape[0], width=rgb.shape[1],
    count=3, dtype="uint8",
    transform=clipped.rio.transform(),
    crs=clipped.rio.crs,
    compress="deflate"
) as dst:
    dst.write(rgb[..., 0], 1)
    dst.write(rgb[..., 1], 2)
    dst.write(rgb[..., 2], 3)

print("Wrote colorized RGB DEM:", rgb_path)

import os
import requests
import geopandas as gpd
import geoplot as gplt
import contextily as ctx
import matplotlib.pyplot as plt

API_KEY = os.environ.get("GEOAPIFY_API_KEY")
if not API_KEY:
    raise RuntimeError("Set the GEOAPIFY_API_KEY environment variable")

# ---- 1) Get Portland neighborhoods via Boundaries "consists-of" ----
# This ID corresponds to Portland; geometry=geometry_1000 returns simplified polygons.
PORTLAND_ID = "51ba64411daca75ec059c0092a6673c54640f00101f901d3d8020000000000c00206920308506f72746c616e64"
bnd_url = (
    "https://api.geoapify.com/v1/boundaries/consists-of"
    f"?id={PORTLAND_ID}&geometry=geometry_1000&apiKey={API_KEY}"
)
bnd_resp = requests.get(bnd_url, timeout=30)
bnd_resp.raise_for_status()
neighborhoods = gpd.GeoDataFrame.from_features(bnd_resp.json()["features"], crs="EPSG:4326")

# ---- 2) Fetch supermarkets inside Portland boundary ----
# Use the same ID with the Places API’s place filter.
places_params = {
    "categories": "commercial.supermarket",
    "filter": f"place:{PORTLAND_ID}",
    "limit": 500,
    "apiKey": API_KEY,
}
p_resp = requests.get("https://api.geoapify.com/v2/places", params=places_params, timeout=30)
p_resp.raise_for_status()
charging_gdf = gpd.GeoDataFrame.from_features(p_resp.json()["features"], crs="EPSG:4326")

# ---- 3) Spatial join: count chargers per neighborhood ----
joined = gpd.sjoin(neighborhoods, charging_gdf, predicate="contains", how="left")
counts = joined.groupby(joined.index).size()
neighborhoods["charger_count"] = neighborhoods.index.map(counts).fillna(0).astype(int)

# ---- 4) Plot a choropleth + points over a Geoapify basemap ----
neighborhoods_3857 = neighborhoods.to_crs(epsg=3857).explode(index_parts=False).reset_index(drop=True)
charging_3857 = charging_gdf.to_crs(epsg=3857)

fig, ax = plt.subplots(figsize=(10, 10))

gplt.choropleth(
    neighborhoods_3857,
    hue="charger_count",
    cmap="Blues",
    linewidth=0.3,
    edgecolor="white",
    legend=True,
    ax=ax,
)
gplt.pointplot(charging_3857, color="#f03b20", alpha=0.7, s=10, ax=ax)

# 1) set a solid extent (pad a bit so at least one tile is guaranteed)
minx, miny, maxx, maxy = neighborhoods_3857.total_bounds
pad = max((maxx - minx), (maxy - miny)) * 0.05  # 5% padding
ax.set_xlim(minx - pad, maxx + pad)
ax.set_ylim(miny - pad, maxy + pad)

# 2) pass source as URL string + correct CRS + explicit zoom
tile_url = f"https://maps.geoapify.com/v1/tile/dark-matter/{{z}}/{{x}}/{{y}}.png?apiKey={API_KEY}"
ctx.add_basemap(
    ax,
    source=tile_url,     # <- string, not dict
    crs="EPSG:3857",
    zoom=12,
    reset_extent=False,
    attribution="© Geoapify, © OpenMapTiles, © OpenStreetMap contributors",
)

ax.set_title("Supermarkets by Portland neighborhood")
ax.axis("off")
plt.tight_layout()
plt.savefig("portland-chargers.png", dpi=200)
print("Static map saved to portland-chargers.png")