Recipe 7 — Earthquake Cluster Maps

# Install necessary packages if not already installed
# !pip install pystac-client folium ipywidgets pandas pystac-monty

import os
import datetime
from getpass import getpass
import folium
from folium.plugins import MarkerCluster
import ipywidgets as widgets
from IPython.display import display, clear_output
import pandas as pd
from pystac_client import Client
import pystac

# Montandon STAC API URL (CORRECT URL with /stac suffix)
STAC_API_URL = "https://montandon-eoapi-stage.ifrc.org/stac"

# ============================================================================
# AUTHENTICATION
# ============================================================================
# First try to get token from environment variable
api_token = os.getenv('MONTANDON_API_TOKEN')

# If not set, prompt user to enter token
if api_token is None:
    print("=" * 70)
    print("AUTHENTICATION REQUIRED")
    print("=" * 70)
    print("\nThe Montandon STAC API requires a Bearer Token for authentication.")
    print("\nHow to get your token:")
    print("  1. Visit: https://goadmin-stage.ifrc.org/")
    print("  2. Log in with your IFRC credentials")
    print("  3. Generate an API token from your account settings")
    print("\nAlternatively, set the MONTANDON_API_TOKEN environment variable:")
    print("  PowerShell: $env:MONTANDON_API_TOKEN = 'your_token_here'")
    print("  Bash: export MONTANDON_API_TOKEN='your_token_here'")
    print("\n" + "=" * 70)
    api_token = getpass("Enter your Montandon API Token: ")

# Create authentication headers for pystac_client
AUTH_HEADERS = {
    "Authorization": f"Bearer {api_token}"
}

# Connect to the STAC API using pystac_client
print(f"\nConnecting to STAC API: {STAC_API_URL}")
try:
    catalog = Client.open(STAC_API_URL, headers=AUTH_HEADERS)
    print(f"✅ Connected successfully!")
    print(f"   Catalog ID: {catalog.id}")
    print(f"   Catalog Title: {catalog.title}")
    
    # List available collections (pystac_client handles pagination)
    collections = list(catalog.get_collections())
    print(f"   Available Collections: {len(collections)}")
except Exception as e:
    print(f"❌ Failed to connect: {e}")
    catalog = None

# Helper functions to extract Monty Extension data from STAC Items

def get_hazard_detail(item: pystac.Item):
    """Extracts hazard detail from a STAC Item."""
    hazard_dict = item.properties.get("monty:hazard_detail")
    if not hazard_dict:
        return None
    
    # Create a simple object to hold the hazard details
    class HazardDetail:
        def __init__(self, data):
            self.cluster = data.get("cluster")
            self.severity_value = data.get("severity_value")
            self.severity_unit = data.get("severity_unit")
            self.severity_label = data.get("severity_label")
            self.estimate_type = data.get("estimate_type")
    
    return HazardDetail(hazard_dict)


def get_impact_detail(item: pystac.Item):
    """Extracts impact detail from a STAC Item."""
    impact_dict = item.properties.get("monty:impact_detail")
    if not impact_dict:
        return None
    
    # Create a simple object to hold the impact details
    class ImpactDetail:
        def __init__(self, data):
            self.category = data.get("category")
            self.type = data.get("type")
            self.value = data.get("value")
            self.unit = data.get("unit")
            self.estimate_type = data.get("estimate_type")
    
    return ImpactDetail(impact_dict)


def get_monty_correlation_id(item: pystac.Item):
    """Extract the correlation ID from a STAC item."""
    return item.properties.get("monty:correlation_id")

# Define earthquake hazard codes across different classification systems
# These codes identify earthquakes in the STAC items
EARTHQUAKE_HAZARD_CODES = [
    # GLIDE classification
    "EQ",
    # EM-DAT CRED classification
    "nat-geo-ear-gro",  # Natural > Geophysical > Earthquake > Ground movement
    "nat-geo-ear-tsu",  # Natural > Geophysical > Earthquake > Tsunami
    # UNDRR-ISC 2025 Hazard Information Profiles
    "GH0101",           # Earthquake (Seismic cluster)
    # UNDRR-ISC 2020 (Historical)
    "GH0001",           # Earthquake
    "GH0002",           # Ground Shaking (Earthquake)
    "GH0003",           # Liquefaction (Earthquake Trigger)
    "GH0004",           # Earthquake Surface Rupture
    "GH0005",           # Subsidence and Uplift (Earthquake Trigger)
    "GH0006",           # Tsunami (Earthquake Trigger)
    "GH0007",           # Landslide (Earthquake Trigger)
]

# Define all event collections to search across (excluding USGS)
EVENT_COLLECTIONS = [
    "emdat-events",
    "desinventar-events",
    "gdacs-events",
    "gidd-events",
    "glide-events",
    "ifrc-events",
    "pdc-events"
]

# Define all hazard collections to search across (excluding USGS)
HAZARD_COLLECTIONS = [
    "gdacs-hazards",
    "pdc-hazards",
    "emdat-hazards",
]

# Define all impact collections to search across (excluding USGS)
IMPACT_COLLECTIONS = [
    "emdat-impacts", 
    "desinventar-impacts",
    "gdacs-impacts",
    "idmc-gidd-impacts",
    "idmc-idu-impacts",
    "ifrc-event-impacts",
    "pdc-impacts"
]


def is_earthquake(item: pystac.Item) -> bool:
    """
    Check if a STAC item is an earthquake based on its hazard codes.
    Returns True if any of the item's hazard codes match earthquake codes.
    """
    hazard_codes = item.properties.get("monty:hazard_codes", [])
    
    # Check if any hazard code matches our earthquake codes
    for code in hazard_codes:
        if code in EARTHQUAKE_HAZARD_CODES:
            return True
        # Also check for partial matches (some codes might have variations)
        code_upper = code.upper()
        if code_upper.startswith("EQ") or code_upper.startswith("GH01"):
            return True
    
    return False


def fetch_items(collection_ids, days=30, filter_earthquakes=True):
    """
    Fetch items from multiple collections for the last N days.
    Uses the authenticated catalog connection.
    pystac_client handles pagination automatically.
    
    Args:
        collection_ids: List of collection IDs to search across
        days: Number of days to look back
        filter_earthquakes: If True, only return earthquake-related items
    
    Returns:
        List of items from all collections
    """
    if not catalog:
        print("❌ No catalog connection available. Please check authentication.")
        return []
    
    # Calculate the date range
    end_date = datetime.datetime.now(datetime.timezone.utc)
    start_date = end_date - datetime.timedelta(days=days)
    datetime_range = f"{start_date.isoformat()}/{end_date.isoformat()}"
    
    print(f"Searching across collections: {collection_ids}")
    print(f"Date range: {start_date.date()} to {end_date.date()}")
    if filter_earthquakes:
        print("Filtering for EARTHQUAKE events only")
    
    all_items = []
    collection_counts = {}
    
    # Search each collection separately and track results
    for collection_id in collection_ids:
        try:
            # Use the authenticated catalog for search
            search = catalog.search(
                collections=[collection_id],
                datetime=datetime_range,
                max_items=500
            )
            
            # pystac_client handles pagination automatically
            items = list(search.items())
            
            # Filter for earthquakes if requested
            if filter_earthquakes:
                earthquake_items = [item for item in items if is_earthquake(item)]
                collection_counts[collection_id] = {"total": len(items), "earthquakes": len(earthquake_items)}
                all_items.extend(earthquake_items)
                print(f"  {collection_id}: {len(earthquake_items)} earthquakes out of {len(items)} total items")
            else:
                collection_counts[collection_id] = {"total": len(items), "earthquakes": len(items)}
                all_items.extend(items)
                print(f"  Found {len(items)} items in {collection_id}")
                
        except Exception as e:
            print(f"  ⚠️ Error searching {collection_id}: {e}")
            collection_counts[collection_id] = {"total": 0, "earthquakes": 0}
    
    total_earthquakes = sum(c.get("earthquakes", 0) for c in collection_counts.values())
    total_items = sum(c.get("total", 0) for c in collection_counts.values())
    
    print(f"\n✅ Total: Found {total_earthquakes} earthquakes out of {total_items} total items")
    
    return all_items


# Fetch all EARTHQUAKE data from multiple sources
if catalog:
    print("=" * 60)
    print("FETCHING EARTHQUAKE EVENTS FROM ALL SOURCES")
    print("=" * 60)
    event_items = fetch_items(EVENT_COLLECTIONS, days=30, filter_earthquakes=True)

    print("\n" + "=" * 60)
    print("FETCHING EARTHQUAKE HAZARDS FROM ALL SOURCES")
    print("=" * 60)
    hazard_items = fetch_items(HAZARD_COLLECTIONS, days=30, filter_earthquakes=True)

    print("\n" + "=" * 60)
    print("FETCHING EARTHQUAKE IMPACTS FROM ALL SOURCES")
    print("=" * 60)
    # Impacts might not have hazard codes, so we'll fetch all and filter by correlation later
    impact_items = fetch_items(IMPACT_COLLECTIONS, days=30, filter_earthquakes=False)
else:
    print("❌ Cannot fetch data - no catalog connection")
    event_items = []
    hazard_items = []
    impact_items = []

def get_geometry_centroid(geometry):
    """
    Extract centroid coordinates from any geometry type.
    Returns (longitude, latitude) tuple or (None, None) if no geometry.
    """
    if not geometry:
        return None, None
    
    geom_type = geometry.get("type")
    coords = geometry.get("coordinates")
    
    if not coords:
        return None, None
    
    if geom_type == "Point":
        # Point: [lon, lat] or [lon, lat, alt]
        return coords[0], coords[1]
    
    elif geom_type == "Polygon":
        # Polygon: [[[lon, lat], [lon, lat], ...]]
        # Calculate centroid by averaging all exterior ring coordinates
        ring = coords[0]  # Exterior ring
        if ring:
            avg_lon = sum(c[0] for c in ring) / len(ring)
            avg_lat = sum(c[1] for c in ring) / len(ring)
            return avg_lon, avg_lat
    
    elif geom_type == "MultiPolygon":
        # MultiPolygon: [[[[lon, lat], ...]], [[[lon, lat], ...]]]
        all_coords = []
        for polygon in coords:
            ring = polygon[0]  # Exterior ring of each polygon
            all_coords.extend(ring)
        if all_coords:
            avg_lon = sum(c[0] for c in all_coords) / len(all_coords)
            avg_lat = sum(c[1] for c in all_coords) / len(all_coords)
            return avg_lon, avg_lat
    
    elif geom_type == "LineString":
        # LineString: [[lon, lat], [lon, lat], ...]
        if coords:
            avg_lon = sum(c[0] for c in coords) / len(coords)
            avg_lat = sum(c[1] for c in coords) / len(coords)
            return avg_lon, avg_lat
    
    elif geom_type == "MultiPoint":
        # MultiPoint: [[lon, lat], [lon, lat], ...]
        if coords:
            avg_lon = sum(c[0] for c in coords) / len(coords)
            avg_lat = sum(c[1] for c in coords) / len(coords)
            return avg_lon, avg_lat
    
    return None, None


def safe_float(value):
    """Safely convert a value to float, returning None if conversion fails."""
    if value is None:
        return None
    try:
        val = float(value)
        return val if val > 0 else None
    except (ValueError, TypeError):
        return None


def extract_magnitude(props):
    """
    Extract earthquake magnitude from various possible property locations.
    Different data sources store magnitude in different fields.
    """
    # 1. Try eq:magnitude (USGS standard)
    mag = safe_float(props.get("eq:magnitude"))
    if mag is not None:
        return mag
    
    # 2. Try monty:hazard_detail.severity_value
    hazard_detail = props.get("monty:hazard_detail", {})
    if hazard_detail:
        mag = safe_float(hazard_detail.get("severity_value"))
        if mag is not None:
            return mag
    
    # 3. Try magnitude directly in properties (some sources)
    mag = safe_float(props.get("magnitude"))
    if mag is not None:
        return mag
    
    # 4. Try to extract from title (e.g., "M 5.2 - Location")
    title = props.get("title", "")
    if title:
        import re
        # Match patterns like "M5.2", "M 5.2", "Magnitude 5.2"
        match = re.search(r'[Mm](?:agnitude)?\s*(\d+\.?\d*)', title)
        if match:
            mag = safe_float(match.group(1))
            if mag is not None:
                return mag
    
    return None


def extract_depth(props):
    """
    Extract earthquake depth from various possible property locations.
    """
    # 1. Try eq:depth (USGS standard)
    depth = safe_float(props.get("eq:depth"))
    if depth is not None:
        return depth
    
    # 2. Try depth directly in properties
    depth = safe_float(props.get("depth"))
    if depth is not None:
        return depth
    
    # 3. Try monty:hazard_detail for depth info
    hazard_detail = props.get("monty:hazard_detail", {})
    if hazard_detail:
        depth = safe_float(hazard_detail.get("depth"))
        if depth is not None:
            return depth
    
    return None


def events_to_dataframe(events):
    """
    Convert STAC event items to a DataFrame.
    Handles events from multiple sources with different geometry types and properties.
    """
    data = []
    for item in events:
        props = item.properties
        
        # Get coordinates from any geometry type
        lon, lat = get_geometry_centroid(item.geometry)
        
        # Skip items with no valid geometry
        if lon is None or lat is None:
            print(f"Skipping item {item.id} - no valid geometry")
            continue
        
        # Get collection/source info
        source = item.collection_id if hasattr(item, 'collection_id') else "unknown"
        
        # Get hazard codes for classification
        hazard_codes = props.get("monty:hazard_codes", [])
        hazard_type = hazard_codes[0] if hazard_codes else "Unknown"
        
        # Extract magnitude using helper function
        magnitude = extract_magnitude(props)
        
        # Extract depth using helper function
        depth = extract_depth(props)
        
        # Get country codes
        country_codes = props.get("monty:country_codes", [])
        countries = ", ".join(country_codes) if country_codes else "Unknown"
        
        # Get hazard detail for additional info
        hazard_detail = props.get("monty:hazard_detail", {})
        severity_unit = hazard_detail.get("severity_unit", "") if hazard_detail else ""
        severity_label = hazard_detail.get("severity_label", "") if hazard_detail else ""
        
        data.append({
            "id": item.id,
            "title": props.get("title", "N/A"),
            "time": item.datetime,
            "source": source,
            "hazard_type": hazard_type,
            "magnitude": magnitude,
            "depth": depth,
            "severity_unit": severity_unit,
            "severity_label": severity_label,
            "countries": countries,
            "correlation_id": props.get("monty:correlation_id", "N/A"),
            "longitude": lon,
            "latitude": lat,
            "stac_item": item
        })
    
    df = pd.DataFrame(data)
    
    # Sort by time descending (most recent first)
    if not df.empty and "time" in df.columns:
        df = df.sort_values(by="time", ascending=False)
    
    # Print statistics about magnitude/depth availability
    if not df.empty:
        mag_count = df['magnitude'].notna().sum()
        depth_count = df['depth'].notna().sum()
        print(f"Created DataFrame with {len(df)} events from {df['source'].nunique()} sources")
        print(f"  - Events with magnitude: {mag_count} ({100*mag_count/len(df):.1f}%)")
        print(f"  - Events with depth: {depth_count} ({100*depth_count/len(df):.1f}%)")
        print(f"Sources: {df['source'].value_counts().to_dict()}")
    
    return df

events_df = events_to_dataframe(event_items)
events_df.head(10)

# Color mapping for different data sources (excluding USGS)
SOURCE_COLORS = {
    "emdat-events": "purple",
    "desinventar-events": "darkblue",
    "gdacs-events": "green",
    "gidd-events": "brown",
    "glide-events": "orange",
    "ifrc-events": "magenta",
    "pdc-events": "navy",
    "gfd-events": "teal"
}

def create_global_map(events_df):
    """
    Create a global map showing all earthquake events from multiple sources.
    Events are colored by source and sized by magnitude/severity.
    """
    if events_df.empty:
        print("No earthquake events to display")
        return folium.Map(location=[0, 0], zoom_start=2)
    
    # Filter out rows with invalid coordinates
    valid_df = events_df.dropna(subset=["latitude", "longitude"])
    
    if valid_df.empty:
        print("No events with valid coordinates")
        return folium.Map(location=[0, 0], zoom_start=2)
    
    center_lat = valid_df["latitude"].mean()
    center_lon = valid_df["longitude"].mean()
    
    m = folium.Map(location=[center_lat, center_lon], zoom_start=2, tiles="CartoDB positron")
    
    # Calculate statistics for title
    mag_available = valid_df['magnitude'].notna().sum()
    depth_available = valid_df['depth'].notna().sum()
    
    # Add title
    title_html = f"""
    <h3 align="center" style="font-size:18px; margin-top:10px;">
        <b>Earthquake Events from Multiple Sources</b>
    </h3>
    <p align="center" style="font-size:12px; color:gray;">
        Showing {len(valid_df)} earthquakes from {valid_df['source'].nunique()} data sources<br>
        Magnitude data: {mag_available}/{len(valid_df)} events | Depth data: {depth_available}/{len(valid_df)} events
    </p>
    """
    m.get_root().html.add_child(folium.Element(title_html))
    
    marker_cluster = MarkerCluster(name="Earthquakes").add_to(m)
    
    for _, row in valid_df.iterrows():
        mag = row.get("magnitude")
        depth = row.get("depth")
        source = row.get("source", "unknown")
        hazard_type = row.get("hazard_type", "Unknown")
        
        # Get color based on source
        color = SOURCE_COLORS.get(source, "gray")
        
        # Size based on magnitude/severity
        if mag is not None and mag > 0:
            if mag >= 8.0:
                radius = 14
            elif mag >= 7.0:
                radius = 11
            elif mag >= 6.0:
                radius = 9
            elif mag >= 5.0:
                radius = 7
            elif mag >= 4.0:
                radius = 5
            else:
                radius = 4
        else:
            radius = 4  # Default size for events without magnitude
            
        # Create popup content with magnitude and depth always shown
        popup_content = f"""
        <div style="font-size: 14px; min-width: 200px;">
            <b>{row['title']}</b><br>
            <hr style="margin: 5px 0;">
            <b>Magnitude:</b> {f'{mag:.1f}' if mag is not None else 'Not available'}<br>
            <b>Depth:</b> {f'{depth:.1f} km' if depth is not None else 'Not available'}<br>
            <hr style="margin: 5px 0;">
            <b>Source:</b> {source}<br>
            <b>Hazard Type:</b> {hazard_type}<br>
            <b>Time:</b> {row['time']}<br>
            <b>Countries:</b> {row.get('countries', 'N/A')}<br>
            <b>Correlation ID:</b> {row.get('correlation_id', 'N/A')}<br>
        </div>
        """
        
        # Tooltip with magnitude and source info
        if mag is not None and mag > 0:
            tooltip = f"M{mag:.1f}"
            if depth is not None:
                tooltip += f" | Depth: {depth:.1f}km"
            tooltip += f" | {source}"
        else:
            tooltip = f"{source}: {row['title'][:30]}..."
        
        folium.CircleMarker(
            location=[row["latitude"], row["longitude"]],
            radius=radius,
            color=color,
            fill=True,
            fill_opacity=0.7,
            popup=folium.Popup(popup_content, max_width=350),
            tooltip=tooltip
        ).add_to(marker_cluster)
    
    # Add a legend (excluding USGS)
    legend_html = """
    <div style="position: fixed; bottom: 50px; left: 50px; z-index: 1000; 
                background-color: white; padding: 10px; border: 2px solid gray;
                border-radius: 5px; font-size: 12px;">
        <b>Earthquake Data Sources</b><br>
    """
    for source, color in SOURCE_COLORS.items():
        source_name = source.replace("-events", "").upper()
        legend_html += f'<i style="background:{color}; width:12px; height:12px; display:inline-block; border-radius:50%;"></i> {source_name}<br>'
    
    legend_html += """<hr style="margin: 5px 0;">
        <b>Circle Size = Magnitude</b><br>
        <span style="font-size:10px;">Larger = Stronger earthquake</span>
    """
    legend_html += "</div>"
    
    m.get_root().html.add_child(folium.Element(legend_html))
    
    return m

global_map = create_global_map(events_df)
global_map

# Define hazard and impact collections to search (excluding USGS)
HAZARD_COLLECTIONS_DETAIL = [
    "gdacs-hazards",
    "pdc-hazards",
    "emdat-hazards"
]

IMPACT_COLLECTIONS_DETAIL = [
    "emdat-impacts", 
    "desinventar-impacts",
    "gdacs-impacts",
    "idmc-gidd-impacts",
    "idmc-idu-impacts",
    "ifrc-event-impacts",
    "pdc-impacts"
]

def find_related_items_by_correlation(event_item: pystac.Item):
    """
    Find related hazard and impact items using the correlation ID.
    Uses the authenticated catalog connection.
    pystac_client handles pagination automatically.
    """
    if not catalog:
        print("No catalog connection available")
        return [], [], []
    
    correlation_id = get_monty_correlation_id(event_item)
    
    if not correlation_id:
        print(f"No correlation ID found for event {event_item.id}")
        return [], [], []
    
    print(f"Searching for items with correlation ID: {correlation_id}")
    
    # Build CQL2 filter for correlation ID
    filter_dict = {
        "op": "=",
        "args": [
            {"property": "monty:correlation_id"},
            correlation_id
        ]
    }
    
    # Search for hazards across all hazard collections
    related_hazards = []
    hazard_sources = []
    try:
        hazard_search = catalog.search(
            collections=HAZARD_COLLECTIONS_DETAIL,
            filter=filter_dict,
            max_items=100
        )
        related_hazards = list(hazard_search.items())
        
        for hazard in related_hazards:
            source = hazard.collection_id if hasattr(hazard, 'collection_id') else "unknown"
            if source not in hazard_sources:
                hazard_sources.append(source)
        
        print(f"   Found {len(related_hazards)} hazards from sources: {hazard_sources}")
    except Exception as e:
        print(f"   Error searching hazards: {e}")
    
    # Search for impacts across ALL impact collections
    related_impacts = []
    impact_sources = []
    try:
        impact_search = catalog.search(
            collections=IMPACT_COLLECTIONS_DETAIL,
            filter=filter_dict,
            max_items=100
        )
        related_impacts = list(impact_search.items())
        
        for impact in related_impacts:
            source = impact.collection_id if hasattr(impact, 'collection_id') else "unknown"
            if source not in impact_sources:
                impact_sources.append(source)
        
        print(f"   Found {len(related_impacts)} impacts from sources: {impact_sources}")
    except Exception as e:
        print(f"   Error searching impacts: {e}")
    
    return related_hazards, related_impacts, impact_sources


def add_geometry_to_map(m, geometry, color, popup, tooltip, fill_opacity=0.3):
    """
    Add any geometry type to a folium map.
    Handles Point, Polygon, MultiPolygon, LineString, etc.
    """
    if not geometry:
        return
    
    geom_type = geometry.get("type") if isinstance(geometry, dict) else geometry["type"]
    coords = geometry.get("coordinates") if isinstance(geometry, dict) else geometry["coordinates"]
    
    if not coords:
        return
    
    if geom_type == "Point":
        folium.CircleMarker(
            location=[coords[1], coords[0]],
            radius=8,
            color=color,
            fill=True,
            fill_opacity=0.6,
            popup=popup,
            tooltip=tooltip,
        ).add_to(m)
    
    elif geom_type == "Polygon":
        ring = coords[0]
        folium_coords = [[c[1], c[0]] for c in ring]
        folium.Polygon(
            locations=folium_coords,
            color=color,
            weight=2,
            fill=True,
            fill_opacity=fill_opacity,
            popup=popup,
            tooltip=tooltip,
        ).add_to(m)
    
    elif geom_type == "MultiPolygon":
        for polygon in coords:
            ring = polygon[0]
            folium_coords = [[c[1], c[0]] for c in ring]
            folium.Polygon(
                locations=folium_coords,
                color=color,
                weight=2,
                fill=True,
                fill_opacity=fill_opacity,
                popup=popup,
                tooltip=tooltip,
            ).add_to(m)
    
    elif geom_type == "LineString":
        folium_coords = [[c[1], c[0]] for c in coords]
        folium.PolyLine(
            locations=folium_coords,
            color=color,
            weight=3,
            popup=popup,
            tooltip=tooltip,
        ).add_to(m)


# Function to create a map showing hazards and impacts for a selected earthquake event
def create_detail_map(event_id):
    """
    Create a detailed map showing hazards and impacts for a selected earthquake event.
    Uses the authenticated catalog connection for searches.
    """
    # Find the selected event
    selected_event = next((item for item in event_items if item.id == event_id), None)

    if not selected_event:
        print(f"Event not found: {event_id}")
        return folium.Map(location=[0, 0], zoom_start=2)
    
    # Find related hazards and impacts using CORRELATION ID
    related_hazards, related_impacts, impact_sources = find_related_items_by_correlation(selected_event)

    # Get event coordinates using the helper function
    event_lon, event_lat = get_geometry_centroid(selected_event.geometry)
    
    if event_lon is None or event_lat is None:
        print(f"No valid geometry for event {event_id}")
        return folium.Map(location=[0, 0], zoom_start=2)

    # Create the map centered on the event
    m = folium.Map(location=[event_lat, event_lon], zoom_start=5, tiles="CartoDB positron")

    # Get event details
    props = selected_event.properties
    title = props.get("title", "Earthquake Details")
    correlation_id = get_monty_correlation_id(selected_event) or "N/A"
    source = selected_event.collection_id if hasattr(selected_event, 'collection_id') else "unknown"
    hazard_codes = props.get("monty:hazard_codes", [])
    country_codes = props.get("monty:country_codes", [])
    
    # Extract magnitude and depth using helper functions
    magnitude = extract_magnitude(props)
    depth = extract_depth(props)
    
    # Format magnitude and depth strings
    mag_str = f"{magnitude:.1f}" if magnitude is not None else "Not available"
    depth_str = f"{depth:.1f} km" if depth is not None else "Not available"
    
    sources_str = ", ".join(impact_sources) if impact_sources else "None found"
    
    title_html = f"""
    <h3 align="center" style="font-size:20px">
        <b>{title}</b>
    </h3>
    <p align="center" style="font-size:12px; color:gray;">
        Magnitude: {mag_str} | Depth: {depth_str}<br>
        Source: {source} | Correlation ID: {correlation_id}<br>
        Hazards: {len(related_hazards)} | Impacts: {len(related_impacts)} (from: {sources_str})
    </p>
    """
    m.get_root().html.add_child(folium.Element(title_html))

    # Create popup content for the event
    date_time = selected_event.datetime.strftime("%Y-%m-%d %H:%M:%S UTC") if selected_event.datetime else "N/A"
    
    event_popup = f"""
    <div style="font-size: 14px; min-width: 250px;">
        <b>{title}</b><br>
        <hr style="margin: 5px 0;">
        <b>Magnitude:</b> {mag_str}<br>
        <b>Depth:</b> {depth_str}<br>
        <hr style="margin: 5px 0;">
        <b>Source:</b> {source}<br>
        <b>Time:</b> {date_time}<br>
        <b>ID:</b> {selected_event.id}<br>
        <b>Correlation ID:</b> {correlation_id}<br>
        <b>Hazard Codes:</b> {', '.join(hazard_codes) if hazard_codes else 'N/A'}<br>
        <b>Countries:</b> {', '.join(country_codes) if country_codes else 'N/A'}<br>
    </div>
    """

    # Add event marker/geometry with magnitude in tooltip
    tooltip_text = f"Earthquake: M{mag_str}" if magnitude else "Earthquake Event"
    
    add_geometry_to_map(
        m, 
        selected_event.geometry, 
        color="red", 
        popup=folium.Popup(event_popup, max_width=300),
        tooltip=tooltip_text,
        fill_opacity=0.5
    )

    # Color mapping for hazard sources
    HAZARD_COLORS = {
        "gdacs-hazards": "darkorange",
        "pdc-hazards": "coral",
        "emdat-hazards": "tomato"
    }

    # Add hazard geometries
    for hazard in related_hazards:
        hazard_detail = get_hazard_detail(hazard)
        hazard_source = hazard.collection_id if hasattr(hazard, 'collection_id') else "unknown"
        color = HAZARD_COLORS.get(hazard_source, "orange")
        
        # Get hazard magnitude/severity
        hazard_mag = None
        if hazard_detail and hazard_detail.severity_value:
            hazard_mag = hazard_detail.severity_value

        hazard_popup = f"""
        <b>{hazard.properties.get("title", "Hazard")}</b><br>
        <b>Source:</b> {hazard_source}<br>
        <b>ID:</b> {hazard.id}<br>
        <b>Magnitude/Severity:</b> {hazard_mag if hazard_mag else 'N/A'} {hazard_detail.severity_unit if hazard_detail else ''}<br>
        <b>Cluster:</b> {hazard_detail.cluster if hazard_detail else 'N/A'}<br>
        <b>Estimate Type:</b> {hazard_detail.estimate_type if hazard_detail else 'N/A'}<br>
        """

        add_geometry_to_map(
            m,
            hazard.geometry,
            color=color,
            popup=folium.Popup(hazard_popup, max_width=300),
            tooltip=f"Hazard: M{hazard_mag}" if hazard_mag else f"Hazard ({hazard_source})",
            fill_opacity=0.2
        )

    # Color mapping for impact sources
    IMPACT_COLORS = {
        "emdat-impacts": "purple",
        "desinventar-impacts": "darkblue",
        "gdacs-impacts": "green",
        "gfd-impacts": "cyan",
        "idmc-gidd-impacts": "brown",
        "idmc-idu-impacts": "olive",
        "ifrc-event-impacts": "magenta",
        "pdc-impacts": "navy"
    }

    # Add impact geometries
    for impact in related_impacts:
        impact_detail = get_impact_detail(impact)
        impact_source = impact.collection_id if hasattr(impact, 'collection_id') else "unknown"
        color = IMPACT_COLORS.get(impact_source, "blue")
        
        # Determine label based on impact type
        if impact_detail and impact_detail.type == "imptypdeat":
            label = f"Fatalities ({impact_source})"
        elif impact_detail and impact_detail.type == "imptypcost":
            label = f"Economic Loss ({impact_source})"
        else:
            label = f"Impact ({impact_source})"

        impact_popup = f"""
        <b>{impact.properties.get("title", "Impact")}</b><br>
        <b>Source:</b> {impact_source}<br>
        <b>ID:</b> {impact.id}<br>
        <b>Category:</b> {impact_detail.category if impact_detail else 'N/A'}<br>
        <b>Type:</b> {impact_detail.type if impact_detail else 'N/A'}<br>
        <b>Value:</b> {impact_detail.value if impact_detail else 'N/A'} {impact_detail.unit if impact_detail else ''}<br>
        <b>Estimate Type:</b> {impact_detail.estimate_type if impact_detail else 'N/A'}<br>
        """

        add_geometry_to_map(
            m,
            impact.geometry,
            color=color,
            popup=folium.Popup(impact_popup, max_width=300),
            tooltip=label,
            fill_opacity=0.3
        )

    # Add layer control
    folium.LayerControl().add_to(m)
    return m

# Step 4 — Interactive earthquake detail map
from IPython.display import HTML
import ipywidgets as widgets

# Build dropdown options from events fetched in Step 1
event_options = [
    (f"{item.properties.get('title', 'Unknown')} ({item.collection_id})", item.id)
    for item in event_items
]

event_dropdown = widgets.Dropdown(
    options=event_options,
    description="Select Earthquake:",
    style={"description_width": "initial"},
    layout=widgets.Layout(width="80%"),
)

render_button = widgets.Button(
    description="Render Map",
    button_style="primary",
    icon="globe",
)

map_output = widgets.Output()


def on_render_click(_):
    event_id = event_dropdown.value
    if not event_id:
        return

    # Disable immediately — prevents queuing up multiple renders while the
    # API call is in progress (this was causing the 9-maps problem)
    render_button.disabled = True
    render_button.description = "Loading..."

    try:
        # Show a loading message straight away
        map_output.clear_output(wait=False)
        with map_output:
            print(f"Loading map for: {event_id[:80]}...")

        # Fetch hazards/impacts and build the Folium map (involves API calls)
        detail_map = create_detail_map(event_id)

        # Swap loading message for the rendered map
        map_output.clear_output(wait=False)
        with map_output:
            display(HTML(detail_map._repr_html_()))
    except Exception as e:
        map_output.clear_output(wait=False)
        with map_output:
            print(f"Error loading map: {e}")
    finally:
        render_button.disabled = False
        render_button.description = "Render Map"


render_button.on_click(on_render_click)

if event_items:
    display(widgets.VBox([event_dropdown, render_button, map_output]))
else:
    print("No earthquake events available to display.")