bottle_trajectory_demo.py

#!/usr/bin/env python3
"""
Mediterranean Bottle Trajectory Demonstration
===========================================

Creates conceptual trajectory maps showing expected bottle drift patterns
based on known Mediterranean circulation and bottle physics.
"""

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from datetime import datetime, timedelta

def mediterranean_current_field(lon, lat):
    """
    Simplified Mediterranean circulation model
    Based on known patterns: Atlantic Water eastward flow, cyclonic gyres
    """
    
    # Base eastward flow (Atlantic Water)
    u_base = 0.12  # m/s eastward
    v_base = 0.02  # m/s slight northward
    
    # Add gyre structures
    # Western Mediterranean gyre (around Balearic Sea)
    west_gyre_lon, west_gyre_lat = 5, 40
    west_gyre_radius = 5  # degrees
    west_gyre_strength = 0.08
    
    # Eastern Mediterranean circulation
    east_gyre_lon, east_gyre_lat = 32, 35
    east_gyre_radius = 4
    east_gyre_strength = 0.06
    
    # Calculate distance-based modifications
    west_dist = np.sqrt((lon - west_gyre_lon)**2 + (lat - west_gyre_lat)**2)
    east_dist = np.sqrt((lon - east_gyre_lon)**2 + (lat - east_gyre_lat)**2)
    
    # Apply gyre effects
    u_current = u_base + west_gyre_strength * np.exp(-west_dist/west_gyre_radius) * np.sin(2*np.pi*lat/10)
    v_current = v_base + east_gyre_strength * np.exp(-east_dist/east_gyre_radius) * np.cos(2*np.pi*lon/15)
    
    # Reduce flow near coasts (simplified)
    if lat < 32:  # Near North African coast
        u_current *= 1.2  # Enhance eastward flow
    if lat > 42:  # Near European coast
        u_current *= 0.7  # Reduce flow
    
    return u_current, v_current

def wind_field(lon, lat, season='summer'):
    """Mediterranean wind patterns by season"""
    
    if season == 'summer':
        # Etesian winds - northerly in eastern Med
        u_wind = 2.0 + 1.5 * np.sin(2*np.pi*lon/30)
        v_wind = -1.0 if lon > 25 else 1.0
    elif season == 'winter':
        # Stronger westerlies
        u_wind = 4.0 + 2.0 * np.sin(2*np.pi*lat/40)
        v_wind = 2.0 + np.cos(2*np.pi*lon/35)
    else:  # spring/autumn
        u_wind = 3.0
        v_wind = 1.5
    
    return u_wind, v_wind

def simulate_bottle_drift(start_lon, start_lat, season='summer', days=30):
    """
    Simulate bottle drift using simplified physics
    """
    
    # Time setup
    dt = 3600  # 1 hour time step
    total_steps = days * 24
    
    # Bottle properties
    wind_factor = 0.02  # 2% of wind speed (partially submerged bottle)
    current_factor = 1.0  # Full current drift
    
    # Storage arrays
    times = []
    lons = []
    lats = []
    
    # Initial conditions
    lon, lat = start_lon, start_lat
    
    for step in range(total_steps):
        # Current time
        current_time = datetime(2024, 6, 15) + timedelta(seconds=step*dt)
        
        # Get environmental conditions
        u_current, v_current = mediterranean_current_field(lon, lat)
        u_wind, v_wind = wind_field(lon, lat, season)
        
        # Add some random variability
        u_current += np.random.normal(0, 0.02)
        v_current += np.random.normal(0, 0.02)
        u_wind += np.random.normal(0, 0.5)
        v_wind += np.random.normal(0, 0.5)
        
        # Calculate total drift velocity
        u_total = current_factor * u_current + wind_factor * u_wind
        v_total = current_factor * v_current + wind_factor * v_wind
        
        # Convert to lat/lon change (rough approximation)
        delta_lon = u_total * dt / (111000 * np.cos(np.radians(lat)))  # m to degrees
        delta_lat = v_total * dt / 111000  # m to degrees
        
        # Update position
        lon += delta_lon
        lat += delta_lat
        
        # Simple boundary conditions (stay in Mediterranean)
        lon = np.clip(lon, 5, 37)
        lat = np.clip(lat, 30, 46)
        
        # Store position (every 6 hours for efficiency)
        if step % 6 == 0:
            times.append(current_time)
            lons.append(lon)
            lats.append(lat)
        
        # Simple stranding check (very rough)
        if (lon < 6 or lon > 36.5 or lat < 30.5 or lat > 45.5):
            break
    
    return np.array(times), np.array(lons), np.array(lats)

def create_mediterranean_overview():
    """Create overview map of the Mediterranean with release points"""
    
    fig = plt.figure(figsize=(16, 12))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    # Set Mediterranean extent
    ax.set_extent([5, 37, 30, 46], crs=ccrs.PlateCarree())
    
    # Add map features
    ax.add_feature(cfeature.LAND, color='lightgray', alpha=0.8)
    ax.add_feature(cfeature.OCEAN, color='lightblue', alpha=0.3)
    ax.add_feature(cfeature.COASTLINE, linewidth=1, color='black')
    ax.add_feature(cfeature.BORDERS, linewidth=0.5, color='gray')
    
    # Add gridlines
    gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
    gl.top_labels = False
    gl.right_labels = False
    
    # Release points
    release_points = {
        'Alexandria': (29.9187, 31.2001),
        'Port Said': (32.3019, 31.2653),
        'Damietta': (31.8160, 31.4165),
        'Marsa Matruh': (27.2373, 31.3543)
    }
    
    colors = ['red', 'blue', 'green', 'orange']
    
    # Simulate and plot trajectories for each release point
    for i, (location, (start_lon, start_lat)) in enumerate(release_points.items()):
        color = colors[i]
        
        # Plot release point
        ax.plot(start_lon, start_lat, 'o', color=color, markersize=12,
                markeredgecolor='black', markeredgewidth=2,
                transform=ccrs.PlateCarree())
        
        # Simulate multiple trajectories for uncertainty
        for traj in range(5):
            times, lons, lats = simulate_bottle_drift(start_lon, start_lat, 'summer', 30)
            
            if len(lons) > 1:
                ax.plot(lons, lats, color=color, alpha=0.6, linewidth=1.5,
                        transform=ccrs.PlateCarree())
        
        # Add location label
        ax.text(start_lon, start_lat - 0.5, location, ha='center', fontsize=10,
                fontweight='bold', transform=ccrs.PlateCarree(),
                bbox=dict(boxstyle='round,pad=0.3', facecolor='white', alpha=0.9))
    
    # Add destination markers
    destinations = {
        'Cyprus': (33.4299, 35.1264),
        'Crete': (24.8093, 35.2401),
        'Sicily': (14.0154, 37.2967),
        'Sardinia': (8.5, 40.0),
        'Turkey': (32.0, 36.5)
    }
    
    for dest, (lon, lat) in destinations.items():
        ax.plot(lon, lat, 's', color='darkred', markersize=8,
                transform=ccrs.PlateCarree())
        ax.text(lon, lat + 0.3, dest, ha='center', fontsize=9,
                transform=ccrs.PlateCarree(),
                bbox=dict(boxstyle='round,pad=0.2', facecolor='yellow', alpha=0.7))
    
    ax.set_title('Mediterranean Bottle Drift Trajectories\n2L Bottle Half-Filled with Rice - Conceptual Model',
                 fontsize=16, fontweight='bold', pad=20)
    
    # Add legend
    legend_elements = [plt.Line2D([0], [0], marker='o', color='w', 
                                  markerfacecolor=colors[i], markersize=10,
                                  label=location, markeredgecolor='black')
                       for i, location in enumerate(release_points.keys())]
    ax.legend(handles=legend_elements, loc='upper left', fontsize=12)
    
    plt.tight_layout()
    plt.savefig('mediterranean_bottle_trajectories_conceptual.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    print("✓ Saved mediterranean_bottle_trajectories_conceptual.png")

def create_seasonal_comparison():
    """Compare trajectories across seasons"""
    
    fig, axes = plt.subplots(2, 2, figsize=(20, 16), 
                            subplot_kw={'projection': ccrs.PlateCarree()})
    
    seasons = ['spring', 'summer', 'autumn', 'winter']
    season_names = ['Spring (March)', 'Summer (June)', 'Autumn (September)', 'Winter (December)']
    
    # Alexandria as test location
    start_lon, start_lat = 29.9187, 31.2001
    
    for i, (season, season_name) in enumerate(zip(seasons, season_names)):
        ax = axes[i//2, i%2]
        
        # Set extent
        ax.set_extent([25, 37, 30, 42], crs=ccrs.PlateCarree())
        
        # Add map features
        ax.add_feature(cfeature.LAND, color='lightgray', alpha=0.8)
        ax.add_feature(cfeature.OCEAN, color='lightblue', alpha=0.3)
        ax.add_feature(cfeature.COASTLINE, linewidth=0.8, color='black')
        
        # Plot multiple trajectories
        for traj in range(8):
            times, lons, lats = simulate_bottle_drift(start_lon, start_lat, season, 30)
            
            if len(lons) > 1:
                ax.plot(lons, lats, color='red', alpha=0.7, linewidth=2,
                        transform=ccrs.PlateCarree())
        
        # Plot release point
        ax.plot(start_lon, start_lat, 'ko', markersize=10, 
                markeredgecolor='white', markeredgewidth=2,
                transform=ccrs.PlateCarree())
        
        ax.set_title(f'{season_name}\nFrom Alexandria', fontsize=14, fontweight='bold')
        
        # Add grid
        gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
        gl.top_labels = False
        gl.right_labels = False
    
    fig.suptitle('Seasonal Variation in Bottle Drift\nConceptual Mediterranean Circulation Model',
                 fontsize=18, fontweight='bold')
    
    plt.tight_layout()
    plt.savefig('seasonal_bottle_drift_comparison.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    print("✓ Saved seasonal_bottle_drift_comparison.png")

def create_probability_heatmap():
    """Create landfall probability heatmap"""
    
    fig = plt.figure(figsize=(16, 10))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    
    # Set extent
    ax.set_extent([5, 37, 30, 46], crs=ccrs.PlateCarree())
    
    # Add map features
    ax.add_feature(cfeature.LAND, color='lightgray', alpha=0.9)
    ax.add_feature(cfeature.OCEAN, color='white')
    ax.add_feature(cfeature.COASTLINE, linewidth=1, color='black')
    
    # Collect landfall positions from multiple simulations
    all_final_lons = []
    all_final_lats = []
    
    release_points = {
        'Alexandria': (29.9187, 31.2001),
        'Port Said': (32.3019, 31.2653),
        'Damietta': (31.8160, 31.4165)
    }
    
    for location, (start_lon, start_lat) in release_points.items():
        for simulation in range(50):  # Many simulations for statistics
            times, lons, lats = simulate_bottle_drift(start_lon, start_lat, 'summer', 30)
            
            if len(lons) > 0:
                all_final_lons.append(lons[-1])
                all_final_lats.append(lats[-1])
    
    # Create heatmap
    if len(all_final_lons) > 0:
        hb = ax.hexbin(all_final_lons, all_final_lats, gridsize=20,
                       cmap='Reds', alpha=0.8, transform=ccrs.PlateCarree())
        
        # Add colorbar
        cbar = plt.colorbar(hb, ax=ax, shrink=0.8, aspect=30)
        cbar.set_label('Landing Probability Density', fontsize=12)
    
    # Plot release points
    for location, (lon, lat) in release_points.items():
        ax.plot(lon, lat, 'bo', markersize=12, markeredgecolor='white',
                markeredgewidth=2, transform=ccrs.PlateCarree())
    
    ax.set_title('Bottle Landing Probability Map\nBased on Conceptual Mediterranean Circulation',
                 fontsize=16, fontweight='bold')
    
    # Add gridlines
    gl = ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False)
    gl.top_labels = False
    gl.right_labels = False
    
    plt.tight_layout()
    plt.savefig('bottle_landing_probability_heatmap.png', dpi=300, bbox_inches='tight')
    plt.show()
    
    print("✓ Saved bottle_landing_probability_heatmap.png")

def analyze_drift_statistics():
    """Analyze drift statistics and create summary"""
    
    print("\nBOTTLE DRIFT ANALYSIS")
    print("=" * 50)
    
    release_points = {
        'Alexandria': (29.9187, 31.2001),
        'Port Said': (32.3019, 31.2653),
        'Damietta': (31.8160, 31.4165),
        'Marsa Matruh': (27.2373, 31.3543)
    }
    
    for location, (start_lon, start_lat) in release_points.items():
        print(f"\n{location} (Start: {start_lon:.2f}°E, {start_lat:.2f}°N):")
        
        distances = []
        final_positions = []
        
        # Run multiple simulations
        for sim in range(30):
            times, lons, lats = simulate_bottle_drift(start_lon, start_lat, 'summer', 30)
            
            if len(lons) > 1:
                final_lon, final_lat = lons[-1], lats[-1]
                final_positions.append((final_lon, final_lat))
                
                # Calculate distance
                delta_lon = final_lon - start_lon
                delta_lat = final_lat - start_lat
                distance_km = np.sqrt((delta_lon * 111 * np.cos(np.radians(start_lat)))**2 + 
                                    (delta_lat * 111)**2)
                distances.append(distance_km)
        
        if distances:
            mean_distance = np.mean(distances)
            std_distance = np.std(distances)
            
            # Calculate mean final position
            final_lons = [pos[0] for pos in final_positions]
            final_lats = [pos[1] for pos in final_positions]
            mean_final_lon = np.mean(final_lons)
            mean_final_lat = np.mean(final_lats)
            
            print(f"  Average drift distance: {mean_distance:.0f} ± {std_distance:.0f} km")
            print(f"  Average final position: {mean_final_lon:.2f}°E, {mean_final_lat:.2f}°N")
            
            # Determine likely destination region
            if mean_final_lon > 32.5:
                if mean_final_lat > 35:
                    region = "Turkish coast"
                else:
                    region = "Cyprus region"
            elif mean_final_lon > 28 and mean_final_lat > 35:
                region = "Greek islands"
            elif mean_final_lon < 20:
                region = "Italian waters"
            else:
                region = "Central Mediterranean"
            
            print(f"  Most likely destination: {region}")
            
            # Success probability (rough estimate)
            northeast_count = sum(1 for lon, lat in final_positions 
                                if lon > start_lon + 2 and lat > start_lat + 1)
            success_rate = northeast_count / len(final_positions) * 100
            print(f"  Northeast drift success: {success_rate:.0f}%")

def main():
    """Main demonstration function"""
    
    print("Mediterranean Bottle Trajectory Demonstration")
    print("=" * 50)
    print("Creating conceptual drift maps based on:")
    print("• Known Mediterranean circulation patterns")
    print("• 2L bottle half-filled with rice physics")
    print("• Seasonal wind and current variations")
    print("• Multiple trajectory uncertainty analysis")
    print()
    
    print("1. Creating overview trajectory map...")
    create_mediterranean_overview()
    
    print("2. Creating seasonal comparison...")
    create_seasonal_comparison()
    
    print("3. Creating probability heatmap...")
    create_probability_heatmap()
    
    print("4. Analyzing drift statistics...")
    analyze_drift_statistics()
    
    print("\n" + "=" * 50)
    print("DEMONSTRATION COMPLETE!")
    print("=" * 50)
    print("\nGenerated files:")
    print("• mediterranean_bottle_trajectories_conceptual.png")
    print("• seasonal_bottle_drift_comparison.png") 
    print("• bottle_landing_probability_heatmap.png")
    print("\nNote: These are conceptual demonstrations based on")
    print("simplified Mediterranean circulation models.")
    print("For accurate predictions, use real oceanographic data.")

if __name__ == "__main__":
    main()