Recipe 1 — Connecting to the Montandon STAC API

# Import libraries
import pandas as pd
from pystac_client import Client
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
import seaborn as sns
from typing import List, Dict, Any
import warnings
import os
from getpass import getpass
import json
warnings.filterwarnings('ignore')

# Set display options for better readability
pd.set_option('display.max_rows', 20)
pd.set_option('display.max_columns', None)
pd.set_option('display.max_colwidth', 50)
pd.set_option('display.width', None)

# Plotting style
sns.set_style('whitegrid')
plt.rcParams['figure.figsize'] = (12, 6)

print("All libraries imported successfully!")

# Connect to Montandon STAC API with Authentication
STAC_API_URL = "https://montandon-eoapi-stage.ifrc.org/stac"

# Get authentication token
# Option 1: From environment variable (recommended for automation)
api_token = os.getenv('MONTANDON_API_TOKEN')

# Option 2: Prompt user for token if not in environment
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)
    
    # Prompt for token (hidden input)
    api_token = getpass("Enter your Montandon API Token: ")
    
    if not api_token or api_token.strip() == "":
        raise ValueError("API token is required to access the Montandon STAC API")

# Create authentication headers
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)")
    print(f"[OK] Auth Provider: https://goadmin-stage.ifrc.org/o/.well-known/openid-configuration")
except Exception as e:
    print(f"\n[ERROR] Authentication failed: {e}")
    print("\nPlease check:")
    print("  1. Your token is valid and not expired")
    print("  2. You have the correct permissions")
    print("  3. The API endpoint is accessible")
    raise

# Check API conformance and capabilities
# Note: The API supports STAC specifications - conformance details available via API endpoints
print("API Connection Verified")
print("=" * 60)
print(f"Connected to: {STAC_API_URL}")
print(f"API Title: {client.title}")
print("\nThe API is ready for data queries.")
print("You can now explore collections and search for disaster event data.")

# Define helper function using pystac_client's client.search() method
from typing import Optional

def search_stac(
    collections: Optional[List[str]] = None,
    max_items: int = 100,
    bbox: Optional[List[float]] = None,
    datetime_range: Optional[str] = None,
    query: Optional[Dict[str, Any]] = None,
    sortby: Optional[List[str]] = None,
    filter_body: Optional[Dict[str, Any]] = None,
    filter_lang: Optional[str] = None,
) -> list:
    """
    Search STAC API using pystac_client's client.search() method.
    
    This uses the standard pystac_client search interface which handles
    pagination automatically.
    
    Parameters:
    -----------
    collections : list of str
        Collection IDs to search
    max_items : int
        Maximum number of results to return
    bbox : list of float
        Bounding box [min_lon, min_lat, max_lon, max_lat]
    datetime_range : str
        ISO 8601 datetime range string (e.g., '2024-01-01/2024-12-31')
    query : dict
        Query parameters for filtering
    sortby : list of str
        Sort order (e.g., [{'field': 'properties.datetime', 'direction': 'desc'}])
    filter_body : dict
        CQL2 filter expression
    filter_lang : str
        Filter language (e.g., 'cql2-json')
    
    Returns:
    --------
    list : List of pystac Item objects matching the search criteria
    """
    # Build search parameters
    search_params = {
        "max_items": max_items
    }
    
    if collections:
        search_params["collections"] = collections
    if bbox:
        search_params["bbox"] = bbox
    if datetime_range:
        search_params["datetime"] = datetime_range
    if query:
        search_params["query"] = query
    if sortby:
        search_params["sortby"] = sortby
    if filter_body:
        search_params["filter"] = filter_body
    if filter_lang:
        search_params["filter_lang"] = filter_lang
    
    # Use client.search() - pagination is handled automatically
    search = client.search(**search_params)
    
    # Return list of items
    return list(search.items())

print("✅ PySTAC helper function initialized (using client.search())")

# Get all collections using pystac_client
# The client handles authentication and pagination automatically

print("Fetching all collections...")

# Get all collections using the client's get_collections() method
collections_list = list(client.get_collections())

print(f"\n{'='*60}")
print(f"Total Collections Available: {len(collections_list)}")
print(f"{'='*60}\n")

# Create a DataFrame for better visualization
collection_data = []
for coll in collections_list:
    desc = coll.description or 'N/A'
    if desc and desc != 'N/A' and len(desc) > 80:
        desc = desc[:80] + '...'
    
    collection_data.append({
        'Collection ID': coll.id,
        'Title': coll.title or 'N/A',
        'Description': desc,
        'License': coll.license or 'N/A'
    })

collections_df = pd.DataFrame(collection_data)
collections_df

# Categorize collections
events_collections = [c for c in collections_list if '-events' in c.id]
hazards_collections = [c for c in collections_list if '-hazards' in c.id]
impacts_collections = [c for c in collections_list if '-impacts' in c.id]

print("Collections by Type:\n")
print(f"Events Collections ({len(events_collections)}):")
for col in events_collections:
    print(f"   • {col.id}")

print(f"Hazards Collections ({len(hazards_collections)}):")
for col in hazards_collections:
    print(f"   • {col.id}")

print(f"Impacts Collections ({len(impacts_collections)}):")
for col in impacts_collections:
    print(f"   • {col.id}")

# Extract data sources
sources = {}
for col in collections_list:
    source = col.id.split('-')[0]
    if source not in sources:
        sources[source] = []
    sources[source].append(col.id)

print("Collections by Data Source:\n")
for source, cols in sorted(sources.items()):
    print(f"{source.upper()} ({len(cols)} collections):")
    for col in sorted(cols):
        print(f"   • {col}")
    print()

# Get a specific collection using client.get_collection()
collection_id = 'gdacs-hazards'

# Fetch collection details using pystac_client
gdacs_hazards = client.get_collection(collection_id)

print(f"Collection Details: {collection_id}\n")
print(f"Title: {gdacs_hazards.title or 'N/A'}")
print(f"Description: {gdacs_hazards.description or 'N/A'}")
print(f"License: {gdacs_hazards.license or 'N/A'}")

# Spatial extent
if gdacs_hazards.extent:
    if gdacs_hazards.extent.spatial:
        print(f"\nSpatial Extent:")
        print(f"   Bounding Box: {gdacs_hazards.extent.spatial.bboxes}")
    
    if gdacs_hazards.extent.temporal:
        print(f"\nTemporal Extent:")
        print(f"   Intervals: {gdacs_hazards.extent.temporal.intervals}")

# Additional properties (keywords, providers, etc.)
if gdacs_hazards.keywords:
    print(f"\nKeywords: {', '.join(gdacs_hazards.keywords)}")

if gdacs_hazards.providers:
    print(f"\nProviders:")
    for provider in gdacs_hazards.providers:
        print(f"   • {provider.name} ({provider.roles})")

# Simple search using client.search()
try:
    search = client.search(
        collections=["gdacs-hazards"],
        max_items=10
    )
    items = list(search.items())
    
    print(f"✅ Retrieved {len(items)} items from gdacs-hazards collection\n")
    
    # Display first item details
    if items:
        first_item = items[0]
        print("📄 First Item Details:")
        print(f"ID: {first_item.id}")
        print(f"Collection: {first_item.collection_id}")
        print(f"Geometry Type: {first_item.geometry.get('type', 'None')}")
        print(f"\nProperties:")
        for key, value in list(first_item.properties.items())[:10]:  # Show first 10 properties
            print(f"  • {key}: {value}")
    else:
        print("No items found in the search results")
        
except Exception as e:
    print(f"❌ Search failed: {e}")

# Direct search using pystac_client's client.search() method

# Perform search using client.search()
search = client.search(
    collections=["gdacs-hazards"],
    max_items=5
)
items_list = list(search.items())

print(f"📡 PySTAC Client Response")
print(f"\n📊 Results:")
print(f"  • Number of Items: {len(items_list)}")
print(f"  • Collection: gdacs-hazards")

# Display feature IDs
if items_list:
    print(f"\n Item IDs:")
    for item in items_list:
        print(f"  • {item.id}")

# Define date range - last 90 days
end_date = datetime.now()
start_date = end_date - timedelta(days=90)

# Format as ISO 8601
datetime_range = f"{start_date.strftime('%Y-%m-%dT%H:%M:%SZ')}/{end_date.strftime('%Y-%m-%dT%H:%M:%SZ')}"

print(f"📅 Searching for events between:")
print(f"   Start: {start_date.strftime('%Y-%m-%d')}")
print(f"   End: {end_date.strftime('%Y-%m-%d')}\n")

# Search with datetime filter using client.search()
search = client.search(
    collections=['gdacs-events'],
    datetime=datetime_range,
    max_items=20
)
items = list(search.items())

print(f"✅ Found {len(items)} events in the last 90 days")

# Display event details
if items:
    print("\n📋 Recent Events:")
    for i, item in enumerate(items[:5], 1):
        event_date = item.properties.get('datetime', 'N/A')
        event_name = item.properties.get('title', item.id)
        print(f"  {i}. {event_name} - {event_date}")

# Define bounding box for Southeast Asia [min_lon, min_lat, max_lon, max_lat]
southeast_asia_bbox = [95.0, -10.0, 141.0, 28.0]

print("🗺️ Searching for events in Southeast Asia")
print(f"   Bounding Box: {southeast_asia_bbox}\n")

# Search with bbox using client.search()
search = client.search(
    collections=['gdacs-events', 'emdat-events', 'usgs-events'],
    bbox=southeast_asia_bbox,
    max_items=30
)
items = list(search.items())

print(f"✅ Found {len(items)} events in Southeast Asia")

# Analyze by collection
if items:
    collection_counts = {}
    for item in items:
        col = item.collection_id
        collection_counts[col] = collection_counts.get(col, 0) + 1
    
    print("\n📊 Events by Collection:")
    for col, count in sorted(collection_counts.items()):
        print(f"  • {col}: {count} events")

# Search for recent events in Europe
europe_bbox = [-10.0, 35.0, 40.0, 71.0]
recent_datetime = "2024-01-01T00:00:00Z/2024-12-31T23:59:59Z"

print("🔍 Searching for 2024 events in Europe\n")

search = client.search(
    collections=['gdacs-events', 'emdat-events'],
    bbox=europe_bbox,
    datetime=recent_datetime,
    max_items=50
)
items = list(search.items())

print(f"✅ Found {len(items)} events in Europe during 2024")

# Display sample
if items:
    print("\n📋 Sample Events:")
    for i, item in enumerate(items[:3], 1):
        print(f"\n  {i}. ID: {item.id}")
        print(f"     Collection: {item.collection_id}")
        print(f"     Title: {item.properties.get('title', 'N/A')}")
        print(f"     Date: {item.properties.get('datetime', 'N/A')}")
        print(f"     Country: {item.properties.get('monty:country_codes', 'N/A')}")

# CQL2 filter for Spain flood events
spain_flood_late_oct_filter = {
    "op": "and",
    "args": [
        {
            "op": "a_contains",
            "args": [
                {"property": "monty:country_codes"},
                "ESP"
            ]
        },
        {
            "op": "a_overlaps",
            "args": [
                {"property": "monty:hazard_codes"},
                ["nat-hyd-flo-flo", "FL"]
            ]
        },
        {
            "op": "t_intersects",
            "args": [
                {"property": "datetime"},
                {"interval": ["2024-10-20T00:00:00Z", "2024-11-05T23:59:59Z"]}
            ]
        }
    ]
}

# Search using client.search() with CQL2 filter
search = client.search(
    collections=["gdacs-events", "emdat-events", "glide-events", "reference-events"],
    filter=spain_flood_late_oct_filter,
    filter_lang="cql2-json",
    max_items=100
)
spain_flood_late_oct_events = list(search.items())

print(f"✅ Found {len(spain_flood_late_oct_events)} flood events in Spain (Oct 20 - Nov 5, 2024)\n")

if spain_flood_late_oct_events:
    print(f"📋 All {len(spain_flood_late_oct_events)} Spain Late October/Early November Flood Events:\n")
    for i, event in enumerate(spain_flood_late_oct_events, 1):
        print(f"\n{i}. Event Details:")
        print(f"   ID: {event.id}")
        print(f"   Title: {event.properties.get('title', 'N/A')}")
        print(f"   Date: {event.properties.get('datetime', 'N/A')}")
        print(f"   Collection: {event.collection_id}")
        print(f"   Countries: {event.properties.get('monty:country_codes', 'N/A')}")
        print(f"   Hazard Codes: {event.properties.get('monty:hazard_codes', 'N/A')}")
        print(f"   Correlation ID: {event.properties.get('monty:corr_id', 'N/A')}")

else:    print("⚠️ No flood events found in Spain during Oct 20 - Nov 5, 2024")

# Get correlation IDs from Spain flood events and find related hazards
if spain_flood_late_oct_events:
    # Collect unique correlation IDs
    correlation_ids = set()
    for event in spain_flood_late_oct_events:
        corr_id = event.properties.get('monty:corr_id')
        if corr_id:
            correlation_ids.add(corr_id)
    
    print(f"🔗 Found {len(correlation_ids)} unique correlation IDs from Spain flood events\n")
    print("="*80)
    
    # Search for hazards for each correlation ID
    all_hazards = []
    hazard_collections_used = ["glide-hazards", "gdacs-hazards", "usgs-hazards", "reference-hazards"]
    
    for idx, corr_id in enumerate(sorted(correlation_ids), 1):
        print(f"\n📍 CORRELATION ID {idx}/{len(correlation_ids)}: {corr_id}")
        print("="*80)
        
        # Search for related hazards using CQL2 filter
        hazard_filter = {
            "op": "=",
            "args": [
                {"property": "monty:corr_id"},
                corr_id
            ]
        }
        
        hazard_search = client.search(
            collections=hazard_collections_used,
            filter=hazard_filter,
            filter_lang="cql2-json",
            max_items=50
        )
        hazards = list(hazard_search.items())
        
        print(f"\n✅ Found {len(hazards)} hazard(s) across collections: {', '.join(hazard_collections_used)}\n")
        
        if hazards:
            all_hazards.extend(hazards)
            
            # Group hazards by collection for clearer display
            hazards_by_collection = {}
            for hazard in hazards:
                coll = hazard.collection_id
                if coll not in hazards_by_collection:
                    hazards_by_collection[coll] = []
                hazards_by_collection[coll].append(hazard)
            
            # Display hazards grouped by collection
            for collection_name, coll_hazards in sorted(hazards_by_collection.items()):
                print(f"📦 Collection: {collection_name} ({len(coll_hazards)} hazard{'s' if len(coll_hazards) > 1 else ''})")
                print("-" * 80)
                
                for i, hazard in enumerate(coll_hazards, 1):
                    print(f"\n  {i}. Hazard ID: {hazard.id}")
                    print(f"     Date: {hazard.properties.get('datetime', 'N/A')}")
                    print(f"     Title: {hazard.properties.get('title', 'N/A')}")
                    
                    # Display hazard detail if available
                    hazard_detail = hazard.properties.get('monty:hazard_detail', {})
                    if hazard_detail:
                        severity = hazard_detail.get('severity_value', 'N/A')
                        severity_unit = hazard_detail.get('severity_unit', '')
                        cluster = hazard_detail.get('cluster', 'N/A')
                        print(f"     Severity: {severity} {severity_unit}")
                        print(f"     Cluster: {cluster}")
                    else:
                        print(f"     Severity: N/A")
                        print(f"     Cluster: N/A")
                
                print()
        else:
            print("⚠️ No hazards found for this correlation ID")
    
    print("\n" + "="*80)
    print(f"\n🎯 FINAL SUMMARY")
    print("="*80)
    print(f"Total Correlation IDs: {len(correlation_ids)}")
    print(f"Total Hazards Found: {len(all_hazards)}")
    
    # Summary by collection
    if all_hazards:
        print(f"\n📊 Hazards by Collection:")
        collection_summary = {}
        for hazard in all_hazards:
            coll = hazard.collection_id
            collection_summary[coll] = collection_summary.get(coll, 0) + 1
        
        for coll, count in sorted(collection_summary.items()):
            print(f"  • {coll}: {count} hazard{'s' if count > 1 else ''}")
    
else:
    print("ℹ️ No Spain flood events found to search for hazards")

# Search for impacts related to Spain flood events
if spain_flood_late_oct_events:
    print("🔍 Searching for impacts related to Spain flood events\n")
    print("="*80)
    
    # Define impact filter for Spain floods
    impact_filter = {
        "op": "and",
        "args": [
            {
                "op": "a_contains",
                "args": [
                    {"property": "monty:country_codes"},
                    "ESP"
                ]
            },
            {
                "op": "a_overlaps",
                "args": [
                    {"property": "monty:hazard_codes"},
                    ["nat-hyd-flo-flo", "FL"]
                ]
            },
            {
                "op": "t_intersects",
                "args": [
                    {"property": "datetime"},
                    {"interval": ["2024-10-20T00:00:00Z", "2024-11-15T23:59:59Z"]}
                ]
            }
        ]
    }
    
    # Search across impact collections
    impact_collections = ["gdacs-impacts", "emdat-impacts", "idmc-gidd-impacts", "idmc-idu-impacts"]
    
    impact_search = client.search(
        collections=impact_collections,
        filter=impact_filter,
        filter_lang="cql2-json",
        max_items=100
    )
    impacts = list(impact_search.items())
    
    print(f"✅ Found {len(impacts)} impact(s) across collections: {', '.join(impact_collections)}\n")

# Diagnostic: Examine the actual structure of impact items
if impacts:
    print("🔍 DIAGNOSTIC: Examining Impact Data Structure\n")
    print("="*80)
    
    # Sample a few impacts from different collections
    sample_impacts = {}
    for impact in impacts[:10]:
        coll = impact.collection_id
        if coll not in sample_impacts:
            sample_impacts[coll] = impact
    
    for collection_name, impact in sample_impacts.items():
        print(f"\n📦 Collection: {collection_name}")
        print(f"Impact ID: {impact.id}")
        print("-" * 80)
        
        # Show all properties
        print("\n📋 All Properties:")
        for key, value in impact.properties.items():
            if isinstance(value, dict):
                print(f"  {key}:")
                for sub_key, sub_value in value.items():
                    print(f"    {sub_key}: {sub_value}")
            else:
                print(f"  {key}: {value}")
        
        # Show monty:impact_detail specifically
        impact_detail = impact.properties.get('monty:impact_detail', {})
        if impact_detail:
            print(f"\n⚡ monty:impact_detail structure:")
            for key, val in impact_detail.items():
                print(f"    {key}: {val}")
        else:
            print(f"\n⚠️ monty:impact_detail is empty or missing")
        
        print("\n" + "="*80)
else:
    print("No impacts to examine")

# Create a DataFrame with all impact data
if impacts:
    impact_data = []
    
    for impact in impacts:
        impact_detail = impact.properties.get('monty:impact_detail', {})
        
        # Extract all relevant fields
        impact_record = {
            'Impact ID': impact.id,
            'Collection': impact.collection_id,
            'Correlation ID': impact.properties.get('monty:corr_id', 'N/A'),
            'Date': impact.properties.get('datetime', 'N/A'),
            'Title': impact.properties.get('title', 'N/A'),
            'Category': impact_detail.get('category', 'N/A') if impact_detail else 'N/A',
            'Type': impact_detail.get('type', 'N/A') if impact_detail else 'N/A',
            'Value': impact_detail.get('value', None) if impact_detail else None,
            'Unit': impact_detail.get('unit', 'N/A') if impact_detail else 'N/A',
            'Estimate Type': impact_detail.get('estimate_type', 'N/A') if impact_detail else 'N/A'
        }
        
        impact_data.append(impact_record)
    
    # Create DataFrame
    impacts_df = pd.DataFrame(impact_data)
    
    # Display summary
    print("Impact Data Summary")
    print("="*80)
    print(f"Total Impacts: {len(impacts_df)}")
    print(f"\nDataFrame Shape: {impacts_df.shape}")
    print(f"Columns: {list(impacts_df.columns)}\n")
    
    # Display the DataFrame
    print("="*80)
    print("\nAll Impacts Table:\n")
    display(impacts_df)
    
    # Summary statistics by category and type
    print("\n" + "="*80)
    print("\nSummary by Category and Type:\n")
    
    summary = impacts_df.groupby(['Category', 'Type']).agg({
        'Value': ['count', 'sum', 'mean', 'min', 'max']
    }).round(0)
    
    # Reset index to make Category and Type regular columns
    summary = summary.reset_index()
    
    # Flatten the multi-level column names
    summary.columns = ['Category', 'Type', 'Count', 'Sum', 'Mean', 'Min', 'Max']
    
    display(summary)
    
else:
    print("No impacts available to create DataFrame")

# Your practice code here
# Example: Search for floods in your region