Recipe 9 — Cascading-Impact Analysis

# Import required libraries
import os
from getpass import getpass
import pystac_client
import pystac
import pandas as pd
import geopandas as gpd
import folium
from datetime import datetime, timedelta
from typing import Optional, List, Dict, Any

# ============================================================================
# AUTHENTICATION & CONNECTION TO MONTANDON STAC API
# ============================================================================

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

# 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
try:
    catalog = pystac_client.Client.open(STAC_API_URL, headers=AUTH_HEADERS)
    print(f"Connected to: {catalog.title}")
    collections = list(catalog.get_collections())
    print(f"Available collections: {len(collections)}")
except Exception as e:
    print(f"Connection failed: {e}")
    catalog = None

# Define earthquake hazard codes
EARTHQUAKE_HAZARD_CODES = [
    "GH0101",           # UNDRR-ISC 2025: Earthquake (Seismic cluster)
    "EQ",               # GLIDE
    "nat-geo-ear-gro",  # EM-DAT: Earthquake > Ground movement
    "GH0001",           # UNDRR-ISC 2020: Earthquake
    "GH0002",           # UNDRR-ISC 2020: Ground Shaking
]

# Define search parameters
target_date = "2023-02-06"
start_datetime = f"{target_date}T00:00:00Z"
end_datetime = f"{target_date}T23:59:59Z"

# Bounding box for Turkey-Syria earthquake region - EXPANDED
# [west, south, east, north] - larger area to capture more aftershocks
TURKEY_SYRIA_BBOX = [34.0, 35.5, 40.0, 39.0]  # Expanded from [35.0, 36.0, 39.0, 38.5]

# Aftershock search duration (days after main event)
AFTERSHOCK_DAYS = 14  # Extended from 5 days to capture more aftershocks

print(f"Search Configuration:")
print(f"  Date: {target_date}")
print(f"  Bounding Box: {TURKEY_SYRIA_BBOX}")
print(f"  Aftershock window: {AFTERSHOCK_DAYS} days")


# ============================================================================
# HELPER FUNCTIONS
# ============================================================================

def get_magnitude(item: pystac.Item) -> Optional[float]:
    """
    Extract magnitude from eq:magnitude or monty:hazard_detail.severity_value
    """
    mag = item.properties.get("eq:magnitude")
    if mag is not None:
        return mag
    
    hazard_detail = item.properties.get("monty:hazard_detail")
    if hazard_detail and isinstance(hazard_detail, dict):
        mag = hazard_detail.get("severity_value")
        if mag is not None:
            return mag
    
    return None


def get_monty_correlation_id(item: pystac.Item) -> Optional[str]:
    """Extract the correlation ID from a STAC item."""
    return item.properties.get("monty:correlation_id")


def get_monty_hazard_detail(item: pystac.Item) -> Dict[str, Any]:
    """Extract hazard detail from a STAC Item."""
    return item.properties.get("monty:hazard_detail", {})


def get_monty_impact_detail(item: pystac.Item) -> Dict[str, Any]:
    """Extract impact detail from a STAC Item."""
    return item.properties.get("monty:impact_detail", {})

# Search for primary earthquake events using BBOX (more reliable than country codes)
# Uses authenticated catalog connection

if not catalog:
    print("No catalog connection. Please check authentication.")
    usgs_hazards = []
    gdacs_emdat_events = []
    primary_events = []
    events_with_mag = []
else:
    # Search USGS hazards with BBOX
    usgs_search = catalog.search(
        collections=["usgs-hazards"],
        bbox=TURKEY_SYRIA_BBOX,
        datetime=f"{start_datetime}/{end_datetime}",
        max_items=100
    )
    usgs_hazards = list(usgs_search.items())
    
    # Search GDACS and EMDAT with BBOX
    gdacs_emdat_search = catalog.search(
        collections=["gdacs-events", "emdat-events", "gdacs-hazards", "emdat-hazards"],
        bbox=TURKEY_SYRIA_BBOX,
        datetime=f"{start_datetime}/{end_datetime}",
        max_items=50
    )
    gdacs_emdat_events = list(gdacs_emdat_search.items())
    
    # Combine all events
    primary_events = usgs_hazards + gdacs_emdat_events
    print(f"Found {len(usgs_hazards)} USGS + {len(gdacs_emdat_events)} GDACS/EMDAT = {len(primary_events)} total events")
    
    # Filter events with magnitude
    events_with_mag = []
    for item in primary_events:
        mag = get_magnitude(item)
        if mag is not None and mag > 0:
            events_with_mag.append(item)
    
    if not events_with_mag:
        events_with_mag = primary_events

# Select the event with the highest magnitude
if events_with_mag:
    main_event = max(events_with_mag, key=lambda x: get_magnitude(x) or 0)
    main_magnitude = get_magnitude(main_event) or 0

    # Extract earthquake details
    eq_depth = main_event.properties.get("eq:depth")
    eq_alert = main_event.properties.get("eq:alert")
    eq_tsunami = main_event.properties.get("eq:tsunami")
    
    print(f"Primary Event: {main_event.properties.get('title')}")
    print(f"  Magnitude: {main_magnitude} | Depth: {eq_depth or 'N/A'} km | Alert: {eq_alert or 'N/A'}")
    
    # Set up aftershock search window - uses AFTERSHOCK_DAYS variable
    event_time = pd.to_datetime(main_event.datetime)
    aftershock_start = event_time + timedelta(minutes=10)
    aftershock_end = event_time + timedelta(days=AFTERSHOCK_DAYS)
    print(f"  Aftershock window: {aftershock_start.date()} to {aftershock_end.date()}")
else:
    print("No events found.")
    main_event = None
    main_magnitude = 0
    eq_depth = None
    eq_alert = None
    eq_tsunami = None

# Known epicenters (matching the helper function)
known_epicenters = {
    "us6000jllz": [37.014, 37.226],  # 7.8M Pazarcik
    "us6000jlqa": [37.203, 38.024],  # 7.5M Elbistan
}

# --- 2a. Aftershocks ---
# Increase max_items to capture more aftershocks
usgs_aftershock_search = catalog.search(
    collections=["usgs-hazards"],
    bbox=TURKEY_SYRIA_BBOX,
    datetime=f"{aftershock_start.isoformat()}/{aftershock_end.isoformat()}",
    max_items=2000  # Increased from 500 to capture more aftershocks
)
usgs_aftershocks = list(usgs_aftershock_search.items())
print(f"Found {len(usgs_aftershocks)} USGS aftershocks")

# Search GDACS hazards for additional aftershocks
gdacs_aftershock_search = catalog.search(
    collections=["gdacs-hazards"],
    bbox=TURKEY_SYRIA_BBOX,
    datetime=f"{aftershock_start.isoformat()}/{aftershock_end.isoformat()}",
    max_items=2000  # Increased from 500
)
gdacs_aftershocks = list(gdacs_aftershock_search.items())

# Combine all aftershocks
aftershocks = usgs_aftershocks + gdacs_aftershocks
print(f"Aftershocks: {len(usgs_aftershocks)} USGS + {len(gdacs_aftershocks)} GDACS = {len(aftershocks)} total")

# --- 2b. Landslides ---
LANDSLIDE_HAZARD_CODES = [
    "GH0007",           # UNDRR-ISC 2020: Landslide (Earthquake Trigger)
    "GH0301",           # UNDRR-ISC 2025: Landslide
    "LS",               # GLIDE
    "nat-hyd-mmw-lan",  # EM-DAT: Hydrological > Mass movement (wet) > Landslide
    "nat-geo-mmd-lan"   # EM-DAT: Geophysical > Mass movement (dry) > Landslide
]

landslide_code_filters = [{"op": "a_contains", "args": [{"property": "monty:hazard_codes"}, code]} for code in LANDSLIDE_HAZARD_CODES]

landslide_filter = {
    "op": "and",
    "args": [
        {
            "op": "or",
            "args": landslide_code_filters
        },
        {
            "op": "t_intersects",
            "args": [{"property": "datetime"}, {"interval": [aftershock_start.isoformat(), aftershock_end.isoformat()]}]
        }
    ]
}

search_landslides = catalog.search(
    collections=["usgs-hazards", "gdacs-hazards", "emdat-hazards"],
    bbox=TURKEY_SYRIA_BBOX,
    filter=landslide_filter,
    filter_lang="cql2-json",
    max_items=100
)
landslides = list(search_landslides.items())
print(f"Found {len(landslides)} landslides in the region")

# Impact search window
impact_end = event_time + timedelta(days=30)

impact_filter = {
    "op": "and",
    "args": [
        # Filter by Country Codes
        {
            "op": "or",
            "args": [
                {"op": "a_contains", "args": [{"property": "monty:country_codes"}, "TUR"]},
                {"op": "a_contains", "args": [{"property": "monty:country_codes"}, "SYR"]}
            ]
        },
        # Filter by Time
        {
            "op": "t_intersects",
            "args": [{"property": "datetime"}, {"interval": [start_datetime, impact_end.isoformat()]}]
        }
    ]
}

search_impacts = catalog.search(
    collections=["idmc-idu-impacts", "emdat-impacts", "gdacs-impacts"],
    bbox=TURKEY_SYRIA_BBOX,
    filter=impact_filter,
    filter_lang="cql2-json",
    max_items=100
)

impacts = list(search_impacts.items())
print(f"Found {len(impacts)} impact records")

# Helper to get coordinates from geometry (Point or Polygon)
def get_event_coordinates(item):
    """
    Extract coordinates from STAC item.
    For known major earthquakes, uses verified epicenter coordinates instead of polygon centroids.
    """
    if not item.geometry:
        return None
    
    # MANUAL OVERRIDE: Known epicenters for Turkey-Syria 2023 earthquakes
    # These are scientifically verified coordinates, more accurate than ShakeMap polygon centroids
    known_epicenters = {
        "us6000jllz": [37.014, 37.226],  # 7.8M Pazarcik/Kahramanmaras - [lon, lat]
        "us6000jlqa": [37.203, 38.024],  # 7.5M Elbistan - [lon, lat]
    }
    
    # Check if this is one of the known events
    for event_id, coords in known_epicenters.items():
        if event_id in item.id.lower():
            return coords
    
    # Fall back to geometry-based extraction
    geom_type = item.geometry.get("type")
    coords = item.geometry.get("coordinates")
    
    if not coords:
        return None
        
    if geom_type == "Point":
        return coords
    elif geom_type == "Polygon":
        # For polygons, calculate centroid (simple average of outer ring)
        # Note: This may not be accurate for large ShakeMap coverage areas
        try:
            outer_ring = coords[0]
            lons = [p[0] for p in outer_ring]
            lats = [p[1] for p in outer_ring]
            centroid_lon = sum(lons) / len(lons)
            centroid_lat = sum(lats) / len(lats)
            return [centroid_lon, centroid_lat]
        except (IndexError, TypeError):
            return None
    return None

# Helper function to add simulated intensity contours
def add_intensity_contours(map_obj, coords, magnitude):
    # Define zones (Color, Radius Multiplier in km)
    zones = [
        ("red", 40),      # Severe Shaking
        ("orange", 80),   # Strong Shaking
        ("yellow", 160)   # Moderate Shaking
    ]
    
    # Adjust radii based on magnitude relative to M7.0
    # Ensure a minimum scale to avoid tiny circles for smaller major events
    scale_factor = max(magnitude / 7.0, 0.5)
    
    for color, base_radius in zones:
        radius_km = base_radius * scale_factor
        folium.Circle(
            location=[coords[1], coords[0]],
            radius=radius_km * 1000,  # Meters
            color=color,
            weight=2,
            fill=True,
            fill_color=color,
            fill_opacity=0.15, # Slightly more transparent to handle overlap
            popup=f"Simulated Intensity Zone (~{int(radius_km)} km) for M{magnitude}"
        ).add_to(map_obj)

# Initialize Map centered on the main event
main_coords = get_event_coordinates(main_event)

if main_coords:
    # Folium uses [lat, lon]
    m = folium.Map(location=[main_coords[1], main_coords[0]], zoom_start=6)
    
    # Add contours for Main Event
    add_intensity_contours(m, main_coords, main_magnitude)
        
else:
    print("Warning: Could not determine main event coordinates. Defaulting to Turkey.")
    m = folium.Map(location=[37.0, 37.0], zoom_start=6)

# Create a detailed popup for the main event
main_popup_html = f"""<div style='width: 200px'>
<b>Main Earthquake Event</b><br>
<b>Magnitude:</b> {main_magnitude}<br>
<b>Date:</b> {main_event.datetime}<br>
<b>Depth:</b> {eq_depth if eq_depth else 'N/A'} km<br>
<b>Alert:</b> {eq_alert if eq_alert else 'N/A'}<br>
<b>Tsunami:</b> {'Yes' if eq_tsunami == 1 else 'No' if eq_tsunami == 0 else 'N/A'}<br>
<b>ID:</b> {main_event.id}
</div>"""

# Add Main Event Marker
if main_coords:
    folium.Marker(
        [main_coords[1], main_coords[0]],
        popup=folium.Popup(main_popup_html, max_width=250),
        tooltip=f"Main Event - Magnitude {main_magnitude}",
        icon=folium.Icon(color="red", icon="star")
    ).add_to(m)

# Initialize counters
major_event_count = 0
aftershock_count = 0

# Add Aftershocks and Major Events
print("Processing aftershocks...")
for item in aftershocks:
    coords = get_event_coordinates(item)
    if coords:
        mag = get_magnitude(item)
        if mag is None:
            continue
            
        # Check for Major Event (e.g., M >= 7.0)
        if mag >= 7.0:
            major_event_count += 1
            
            # Add Marker
            popup_html = f"""<div style='width: 150px'>
            <b>Major Event</b><br>
            <b>Magnitude:</b> {mag}<br>
            <b>Date:</b> {item.datetime}
            </div>"""
            
            folium.Marker(
                [coords[1], coords[0]],
                popup=folium.Popup(popup_html, max_width=200),
                tooltip=f"Major Event - Mag {mag}",
                icon=folium.Icon(color="red", icon="star")
            ).add_to(m)
            
            # Add Contours for Major Event
            add_intensity_contours(m, coords, mag)
            
        else:
            aftershock_count += 1
            popup_html = f"""<div style='width: 150px'>
            <b>Aftershock</b><br>
            <b>Magnitude:</b> {mag}<br>
            <b>Date:</b> {item.datetime}
            </div>"""
            
            folium.CircleMarker(
                [coords[1], coords[0]],
                radius=max(mag * 1.5, 3),  # Size by magnitude, min 3
                color="orange",
                fill=True,
                fill_opacity=0.6,
                popup=folium.Popup(popup_html, max_width=200),
                tooltip=f"Aftershock - Mag {mag}"
            ).add_to(m)

print(f"Map shows {major_event_count} major events (M>=7.0) and {aftershock_count} aftershocks")

# Add Landslides
for item in landslides:
    coords = get_event_coordinates(item)
    if coords:
        popup_html = f"""<div style='width: 150px'>
        <b>Landslide</b><br>
        <b>Date:</b> {item.datetime}<br>
        <b>Location:</b> {item.properties.get('title', 'N/A')}
        </div>"""
        
        folium.Marker(
            [coords[1], coords[0]],
            popup=folium.Popup(popup_html, max_width=200),
            tooltip="Landslide",
            icon=folium.Icon(color="green", icon="leaf")
        ).add_to(m)

# Add Impacts
for item in impacts:
    coords = get_event_coordinates(item)
    if coords:
        detail = item.properties.get("monty:impact_detail", {})
        
        impact_type = detail.get('type', 'Unknown')
        impact_value = detail.get('value', 'N/A')
        impact_unit = detail.get('unit', '')
        
        popup_html = f"""<div style='width: 150px'>
        <b>Impact</b><br>
        <b>Type:</b> {impact_type}<br>
        <b>Value:</b> {impact_value} {impact_unit}
        </div>"""
        
        folium.Marker(
            [coords[1], coords[0]],
            popup=folium.Popup(popup_html, max_width=200),
            tooltip=f"{impact_type}: {impact_value} {impact_unit}",
            icon=folium.Icon(color="blue", icon="info-sign")
        ).add_to(m)

# Add a legend
legend_html = '''<div style="position: fixed; 
                 bottom: 50px; right: 50px; width: 220px; height: auto; 
                 background-color: white; z-index:9999; font-size:14px;
                 border:2px solid grey; border-radius: 5px; padding: 10px">
                 <p style="margin-bottom: 5px;"><b>Legend</b></p>
                 <p style="margin: 3px;"><span style="color: red;">★</span> Main/Major Event (M7+)</p>
                 <p style="margin: 3px;"><span style="color: red; opacity: 0.5;">●</span> Severe Shaking Zone</p>
                 <p style="margin: 3px;"><span style="color: orange; opacity: 0.5;">●</span> Strong Shaking Zone</p>
                 <p style="margin: 3px;"><span style="color: yellow; opacity: 0.5;">●</span> Moderate Shaking Zone</p>
                 <p style="margin: 3px;"><span style="color: orange;">●</span> Aftershocks</p>
                 <p style="margin: 3px;"><span style="color: green;">🍃</span> Landslides</p>
                 <p style="margin: 3px;"><span style="color: blue;">ⓘ</span> Impacts</p>
                 </div>'''

m.get_root().html.add_child(folium.Element(legend_html))

# Display Map
m