Montandon: Earthquakes

import os
import pandas as pd
from pystac_client import Client
import geopandas as gpd
from shapely.geometry import Point, Polygon, shape
from lonboard import viz

from datetime import datetime, timedelta, timezone

STAC_API_URL = "https://montandon-eoapi.ifrc.org/stac"
API_TOKEN = os.getenv('MONTANDON_API_TOKEN')
if not API_TOKEN:
    raise RuntimeError("MONTANDON_API_TOKEN is not set — required to authenticate with the Montandon STAC API.")

auth_headers = {"Authorization": f"Bearer {API_TOKEN}"}

# Connect to STAC API with authentication
try:
    client = Client.open(STAC_API_URL, headers=auth_headers)
    print(f"\n[OK] Connected to: {STAC_API_URL}")
    print(f"[OK] API Title: {client.title}")
    print(f"[OK] Authentication: Bearer Token (OpenID Connect)")
except Exception as e:
    print(f"\n[ERROR] Authentication failed: {e}")

client

collections_list = list(client.get_collections())

events_collection = [c for c in collections_list if '-events' in c.id]

for col in events_collection:
    print(f"{col.title} | {col.id}")

now = datetime.now(timezone.utc)
one_week_ago = two_days_ago = now - timedelta(hours=168)
magnitude = 4 # minimum magnitude to search for

search_hazards = client.search(
    collections=["usgs-hazards"],
    datetime=f"{one_week_ago.isoformat()}/{now.isoformat()}",
    max_items=1000,
)

earthquake_hazards = [
    item for item in search_hazards.items()
    if item.properties.get("monty:hazard_detail", {}).get("severity_value", 0) > magnitude
]

earthquake_hazards[0]

# Get the shape of an item
len(earthquake_hazards)

# ShakeMap polygons are only generated for significant events near populated areas.
# For everything else, fetch the epicenter Point from usgs-events as a fallback.
corr_ids = {item.properties['monty:corr_id'] for item in earthquake_hazards}

search_events = client.search(
    collections=["usgs-events"],
    datetime=f"{one_week_ago.isoformat()}/{now.isoformat()}",
    max_items=1000,
)

epicenters = {
    item.properties['monty:corr_id']: shape(item.geometry)
    for item in search_events.items()
    if item.properties.get('monty:corr_id') in corr_ids
}
print(f"Epicenter Points found for {len(epicenters)} / {len(corr_ids)} earthquakes")

def item_geometry(item):
    """Return ShakeMap polygon if valid, else epicenter Point from usgs-events."""
    geom = shape(item.geometry)
    if not geom.is_empty and geom.bounds != (0.0, 0.0, 0.0, 0.0):
        return geom  # ShakeMap polygon
    bbox = item.bbox
    if bbox and not all(v == 0 for v in bbox):
        return Point((bbox[0] + bbox[2]) / 2, (bbox[1] + bbox[3]) / 2)
    corr_id = item.properties.get('monty:corr_id')
    return epicenters.get(corr_id)  # epicenter Point, or None if missing from events

gdf = gpd.GeoDataFrame(
    [{
        'id': item.id,
        'monty_corr_id': item.properties['monty:corr_id'],
        'datetime': pd.to_datetime(item.properties['datetime']),
        'title': item.properties['title'],
        'magnitude': item.properties['eq:magnitude'],
        'geometry': item_geometry(item),
    } for item in earthquake_hazards],
    geometry='geometry',
    crs='EPSG:4326',
).sort_values('datetime', ascending=False)

null_geom = gdf['geometry'].isna().sum()
print(f"{len(gdf)} earthquakes total, {null_geom} with no usable geometry")

earthquakes = gdf
earthquakes

print(f"There have been {len(earthquakes)} recent earthquakes with a magnitude > 5, with {len(earthquakes['monty_corr_id'].unique())} Montandon correlation IDs.")

from lonboard import Map, PolygonLayer, ScatterplotLayer
from lonboard.layer_extension import DataFilterExtension
import numpy as np

# Split by geometry type: ShakeMap polygons for large events, epicenter Points for the rest
poly_gdf = earthquakes[earthquakes.geometry.geom_type == 'Polygon'].copy()
point_gdf = earthquakes[earthquakes.geometry.geom_type == 'Point'].copy()

# Scale radius by magnitude — each unit doubles the circle, matching the log scale of magnitude.
# M4 → ~20 km, M5 → ~40 km, M6 → ~80 km, M6.6 → ~122 km
radii = (2 ** (point_gdf['magnitude'] - 4) * 20_000).round().astype(int).values

layers = []
if len(poly_gdf):
    layers.append(PolygonLayer.from_geopandas(
        poly_gdf,
        get_fill_color=[220, 80, 0, 120],
        get_line_color=[180, 40, 0, 220],
        get_filter_value=poly_gdf['magnitude'].values,
        extensions=[DataFilterExtension(filter_size=1)],
        filter_range=(0, 10)
    ))
if len(point_gdf):
    layers.append(ScatterplotLayer.from_geopandas(
        point_gdf,
        get_radius=radii,
        get_fill_color=[220, 80, 0, 180],
        radius_min_pixels=3,
        get_filter_value=point_gdf['magnitude'].values,
        extensions=[DataFilterExtension(filter_size=1)],
        filter_range=(0, 10)
    ))

print(f"{len(poly_gdf)} ShakeMap polygons, {len(point_gdf)} epicenter points")
m = Map(layers=layers, show_tooltip=True)

from manywidgets import RangeSlider
from manywidgets.lonboard import FilterBinder
from IPython.display import display

slider = RangeSlider(label="Magnitude", min=0, max=10, low=2, high=10, step=1)
binders = [FilterBinder(slider, layer) for layer in layers]
display(*binders)

slider

m