Tutorial 03: GeoAnalytics with Open-Source Tools¶

# save as: scripts/geo/load_earthquakes.py
import json
import psycopg2
from psycopg2.extras import execute_values

# Connection parameters
conn = psycopg2.connect(
    host="csa-pg-geo-dev.postgres.database.azure.com",
    dbname="geodb",
    user="csaadmin",
    password="YourStr0ng_P@ssw0rd!",
    sslmode="require"
)

with open("data/geospatial/earthquakes_30day.geojson") as f:
    data = json.load(f)

rows = []
for feature in data["features"]:
    props = feature["properties"]
    coords = feature["geometry"]["coordinates"]  # [lon, lat, depth]
    rows.append((
        props["ids"].strip(",").split(",")[0] if props.get("ids") else props["code"],
        props.get("mag"),
        coords[2] if len(coords) > 2 else None,
        props.get("time"),
        props.get("place"),
        props.get("type"),
        coords[0],  # longitude
        coords[1],  # latitude
    ))

cur = conn.cursor()
execute_values(cur, """
    INSERT INTO geo.earthquake_events (event_id, magnitude, depth_km, event_time, place, event_type, geom)
    VALUES %s
    ON CONFLICT (event_id) DO NOTHING
""", rows, template="""(
    %s, %s, %s,
    TO_TIMESTAMP(%s / 1000.0),
    %s, %s,
    ST_SetSRID(ST_MakePoint(%s, %s), 4326)
)""")

conn.commit()
cur.execute("SELECT count(*) FROM geo.earthquake_events;")
print(f"Loaded {cur.fetchone()[0]} earthquake events")

cur.close()
conn.close()

# save as: scripts/geo/load_counties.py
import geopandas as gpd
from sqlalchemy import create_engine

engine = create_engine(
    "postgresql://csaadmin:YourStr0ng_P@ssw0rd!@"
    "csa-pg-geo-dev.postgres.database.azure.com:5432/geodb",
    connect_args={"sslmode": "require"}
)

counties = gpd.read_file("data/geospatial/counties/")
counties = counties.to_crs(epsg=4326)

# Select and rename columns
counties_clean = counties[["GEOID", "NAME", "STATEFP", "COUNTYFP", "ALAND", "AWATER", "geometry"]].copy()
counties_clean.columns = ["geoid", "name", "state_fips", "county_fips", "aland", "awater", "geom"]

# Add state names
import json
state_fips_map = {
    "01": "Alabama", "02": "Alaska", "04": "Arizona", "05": "Arkansas",
    "06": "California", "08": "Colorado", "09": "Connecticut", "10": "Delaware",
    "12": "Florida", "13": "Georgia", "15": "Hawaii", "16": "Idaho",
    "17": "Illinois", "18": "Indiana", "19": "Iowa", "20": "Kansas",
    # ... (abbreviated for tutorial, full map in the script)
}
counties_clean["state_name"] = counties_clean["state_fips"].map(state_fips_map).fillna("Other")

counties_clean.to_postgis("county_boundaries", engine, schema="geo", if_exists="replace", index=False)

print(f"Loaded {len(counties_clean)} county boundaries")

# save as: scripts/geo/postgis_to_geoparquet.py
import geopandas as gpd
from sqlalchemy import create_engine
import os

engine = create_engine(
    "postgresql://csaadmin:YourStr0ng_P@ssw0rd!@"
    "csa-pg-geo-dev.postgres.database.azure.com:5432/geodb",
    connect_args={"sslmode": "require"}
)

# Read the materialized view
print("Reading county earthquake summary from PostGIS...")
gdf = gpd.read_postgis(
    "SELECT * FROM geo.mv_county_earthquake_summary",
    engine,
    geom_col="geom"
)

print(f"Records: {len(gdf)}")
print(f"CRS: {gdf.crs}")
print(f"Columns: {list(gdf.columns)}")

# Write to GeoParquet (Silver layer)
os.makedirs("data/geospatial/silver", exist_ok=True)
output_path = "data/geospatial/silver/county_earthquake_summary.parquet"

gdf.to_parquet(output_path, engine="pyarrow", index=False)
print(f"\nWrote GeoParquet: {output_path}")
print(f"File size: {os.path.getsize(output_path) / 1024 / 1024:.2f} MB")

# Verify by reading back
gdf_verify = gpd.read_parquet(output_path)
print(f"Verified: {len(gdf_verify)} records, CRS={gdf_verify.crs}")

# save as: scripts/geo/geopandas_analysis.py
import geopandas as gpd
import pandas as pd

gdf = gpd.read_parquet("data/geospatial/silver/county_earthquake_summary.parquet")

# Top earthquake-prone states
state_summary = gdf.groupby("state").agg(
    total_earthquakes=("earthquake_count", "sum"),
    max_magnitude=("max_magnitude", "max"),
    counties_affected=("earthquake_count", lambda x: (x > 0).sum()),
    total_counties=("earthquake_count", "count"),
).sort_values("total_earthquakes", ascending=False)

print("Top 10 Earthquake-Prone States:")
print(state_summary.head(10).to_string())

# Calculate area-based density (earthquakes per 1000 sq km)
gdf["area_sqkm"] = gdf.to_crs(epsg=5070).geometry.area / 1e6
gdf["eq_density"] = gdf["earthquake_count"] / (gdf["area_sqkm"] / 1000)

# Counties with highest earthquake density
high_density = gdf[gdf["earthquake_count"] > 0].nlargest(10, "eq_density")
print("\nHighest Earthquake Density (per 1000 sq km):")
print(high_density[["county", "state", "earthquake_count", "area_sqkm", "eq_density"]].to_string())

# save as: scripts/geo/h3_indexing.py
import geopandas as gpd
import pandas as pd
import h3
import json

# Read earthquake data
with open("data/geospatial/earthquakes_30day.geojson") as f:
    eq_data = json.load(f)

records = []
for feature in eq_data["features"]:
    coords = feature["geometry"]["coordinates"]
    props = feature["properties"]
    lat, lon = coords[1], coords[0]

    # Generate H3 indexes at multiple resolutions
    h3_res4 = h3.latlng_to_cell(lat, lon, 4)   # ~1,770 km2
    h3_res6 = h3.latlng_to_cell(lat, lon, 6)   # ~36 km2
    h3_res8 = h3.latlng_to_cell(lat, lon, 8)   # ~0.74 km2

    records.append({
        "event_id": props.get("code", ""),
        "magnitude": props.get("mag"),
        "latitude": lat,
        "longitude": lon,
        "place": props.get("place"),
        "h3_res4": h3_res4,
        "h3_res6": h3_res6,
        "h3_res8": h3_res8,
    })

df = pd.DataFrame(records)
print(f"Total events: {len(df)}")
print(f"Unique H3 cells (res4): {df['h3_res4'].nunique()}")
print(f"Unique H3 cells (res6): {df['h3_res6'].nunique()}")
print(f"Unique H3 cells (res8): {df['h3_res8'].nunique()}")

# Aggregate by H3 resolution 4
h3_agg = df.groupby("h3_res4").agg(
    earthquake_count=("event_id", "count"),
    avg_magnitude=("magnitude", "mean"),
    max_magnitude=("magnitude", "max"),
).reset_index()

# Convert H3 cells to polygons for visualization
def h3_to_polygon(h3_index):
    from shapely.geometry import Polygon
    boundary = h3.cell_to_boundary(h3_index)
    # h3 returns (lat, lng) pairs, Shapely needs (lng, lat)
    return Polygon([(lng, lat) for lat, lng in boundary])

h3_agg["geometry"] = h3_agg["h3_res4"].apply(h3_to_polygon)
h3_gdf = gpd.GeoDataFrame(h3_agg, geometry="geometry", crs="EPSG:4326")

# Save H3 aggregation as GeoParquet
import os
os.makedirs("data/geospatial/gold", exist_ok=True)
h3_gdf.to_parquet("data/geospatial/gold/earthquake_h3_res4.parquet", index=False)

print(f"\nH3 Aggregation (Resolution 4):")
print(h3_gdf.nlargest(10, "earthquake_count")[
    ["h3_res4", "earthquake_count", "avg_magnitude", "max_magnitude"]
].to_string(index=False))

# Databricks notebook cell
%pip install apache-sedona geopandas h3

from sedona.spark import SedonaContext

sedona = SedonaContext.builder().getOrCreate()
print(f"Sedona version: {sedona.version}")

# Databricks notebook: notebooks/geoanalytics/01_spatial_analysis.py

from sedona.spark import SedonaContext
from pyspark.sql.functions import col, count, avg, max as spark_max

sedona = SedonaContext.builder().getOrCreate()

# Read county boundaries from GeoParquet
counties_df = sedona.read.format("geoparquet").load(
    "abfss://silver@csadlsdev.dfs.core.windows.net/geospatial/county_earthquake_summary.parquet"
)

counties_df.printSchema()
counties_df.show(5, truncate=False)
print(f"Total counties: {counties_df.count()}")

from sedona.sql.st_functions import ST_Buffer, ST_Intersects, ST_Area
from pyspark.sql import functions as F

# Read raw earthquake points
earthquakes_df = sedona.read.format("geojson").load(
    "abfss://bronze@csadlsdev.dfs.core.windows.net/geospatial/earthquakes_30day.geojson"
)

# Register as temp views
counties_df.createOrReplaceTempView("counties")
earthquakes_df.createOrReplaceTempView("earthquakes")

# Spatial join: earthquakes within 100km buffer of each county centroid
result = sedona.sql("""
    SELECT
        c.county,
        c.state,
        COUNT(e.*) AS nearby_earthquakes,
        AVG(e.properties.mag) AS avg_magnitude,
        MAX(e.properties.mag) AS max_magnitude,
        ST_Area(ST_Transform(c.geom, 'EPSG:4326', 'EPSG:5070')) / 1e6 AS area_sqkm
    FROM counties c
    JOIN earthquakes e
    ON ST_Intersects(
        ST_Buffer(ST_Transform(c.geom, 'EPSG:4326', 'EPSG:5070'), 100000),
        ST_Transform(e.geometry, 'EPSG:4326', 'EPSG:5070')
    )
    GROUP BY c.county, c.state, c.geom
    HAVING COUNT(e.*) >= 5
    ORDER BY nearby_earthquakes DESC
""")

result.show(20, truncate=False)

# Use Sedona's built-in H3 support (Sedona 1.6+)
h3_result = sedona.sql("""
    SELECT
        ST_H3CellIDs(geometry, 4, false)[0] AS h3_index,
        COUNT(*) AS event_count,
        AVG(properties.mag) AS avg_mag,
        MAX(properties.mag) AS max_mag
    FROM earthquakes
    GROUP BY ST_H3CellIDs(geometry, 4, false)[0]
    ORDER BY event_count DESC
    LIMIT 20
""")

h3_result.show(truncate=False)

# save as: scripts/geo/visualize_results.py
import geopandas as gpd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np

# Read Gold-layer data
counties = gpd.read_parquet("data/geospatial/silver/county_earthquake_summary.parquet")
h3_hexes = gpd.read_parquet("data/geospatial/gold/earthquake_h3_res4.parquet")

# Filter to continental US for cleaner visualization
continental = counties[
    ~counties["state"].isin(["Alaska", "Hawaii", "Puerto Rico", "Other"])
].copy()

fig, axes = plt.subplots(1, 2, figsize=(20, 8))

# Plot 1: Earthquake count by county
ax1 = axes[0]
continental.plot(
    column="earthquake_count",
    cmap="YlOrRd",
    legend=True,
    legend_kwds={"label": "Earthquake Count", "shrink": 0.6},
    edgecolor="gray",
    linewidth=0.1,
    ax=ax1,
    missing_kwds={"color": "lightgrey"},
)
ax1.set_title("Earthquakes by County (Last 30 Days)", fontsize=14, fontweight="bold")
ax1.set_xlim(-130, -65)
ax1.set_ylim(24, 50)
ax1.set_axis_off()

# Plot 2: H3 hexagonal aggregation
ax2 = axes[1]
continental.plot(ax=ax2, color="lightgrey", edgecolor="white", linewidth=0.1)
h3_us = h3_hexes[
    (h3_hexes.geometry.centroid.x > -130) & (h3_hexes.geometry.centroid.x < -65) &
    (h3_hexes.geometry.centroid.y > 24) & (h3_hexes.geometry.centroid.y < 50)
]
h3_us.plot(
    column="earthquake_count",
    cmap="plasma",
    legend=True,
    legend_kwds={"label": "Events per H3 Cell", "shrink": 0.6},
    alpha=0.7,
    edgecolor="black",
    linewidth=0.3,
    ax=ax2,
)
ax2.set_title("H3 Hexagonal Aggregation (Resolution 4)", fontsize=14, fontweight="bold")
ax2.set_xlim(-130, -65)
ax2.set_ylim(24, 50)
ax2.set_axis_off()

plt.tight_layout()
plt.savefig("data/geospatial/gold/earthquake_analysis.png", dpi=150, bbox_inches="tight")
plt.show()
print("Saved visualization to data/geospatial/gold/earthquake_analysis.png")

# save as: scripts/geo/publish_gold.py
import geopandas as gpd
import os

# Read Silver data
counties = gpd.read_parquet("data/geospatial/silver/county_earthquake_summary.parquet")
h3_hexes = gpd.read_parquet("data/geospatial/gold/earthquake_h3_res4.parquet")

gold_dir = "data/geospatial/gold"
os.makedirs(gold_dir, exist_ok=True)

# Gold dataset 1: County earthquake risk profile
risk_profile = counties[counties["earthquake_count"] > 0].copy()
risk_profile["risk_category"] = risk_profile["max_magnitude"].apply(
    lambda m: "High" if m >= 5.0 else ("Medium" if m >= 3.5 else "Low")
)
risk_profile.to_parquet(f"{gold_dir}/county_earthquake_risk.parquet", index=False)
print(f"Gold - County Risk: {len(risk_profile)} records")

# Gold dataset 2: H3 aggregation (already created)
print(f"Gold - H3 Hexagons: {len(h3_hexes)} records")

# Gold dataset 3: State-level summary (non-spatial)
import pandas as pd
state_summary = counties.groupby("state").agg(
    total_earthquakes=("earthquake_count", "sum"),
    max_magnitude=("max_magnitude", "max"),
    counties_affected=("earthquake_count", lambda x: (x > 0).sum()),
    total_counties=("earthquake_count", "count"),
).reset_index()
state_summary.to_parquet(f"{gold_dir}/state_earthquake_summary.parquet", index=False)
print(f"Gold - State Summary: {len(state_summary)} records")