longest_views_map.py

"""
UK Longest Views Map
Finds the 5 longest direct views from the UK coastline to distant countries.
"""

import geopandas as gpd
import folium
from shapely.geometry import Point, LineString, MultiLineString
import numpy as np
from pyproj import Geod
import requests
import zipfile
from pathlib import Path
from tqdm import tqdm


def download_natural_earth_data():
    """Download Natural Earth data if not already cached."""
    data_dir = Path('data')
    data_dir.mkdir(exist_ok=True)

    shapefile_path = data_dir / 'ne_10m_admin_0_countries.shp'

    if shapefile_path.exists():
        print("Using cached Natural Earth data...")
        return str(shapefile_path)

    print("Downloading Natural Earth data...")
    url = 'https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries.zip'

    response = requests.get(url, stream=True)
    zip_path = data_dir / 'ne_countries.zip'

    total_size = int(response.headers.get('content-length', 0))
    with open(zip_path, 'wb') as f:
        with tqdm(total=total_size, unit='B', unit_scale=True, desc="Downloading") as pbar:
            for chunk in response.iter_content(chunk_size=8192):
                f.write(chunk)
                pbar.update(len(chunk))

    print("Extracting...")
    with zipfile.ZipFile(zip_path, 'r') as zip_ref:
        zip_ref.extractall(data_dir)

    zip_path.unlink()
    print("Data downloaded successfully!")
    return str(shapefile_path)


def load_geodata():
    """Load UK and world countries data from Natural Earth."""
    print("Loading geographical data...")
    shapefile_path = download_natural_earth_data()
    world = gpd.read_file(shapefile_path)
    uk = world[world['NAME'] == 'United Kingdom'].copy()
    all_countries = world.copy()
    return uk, all_countries, world


def get_coastline_points(uk_geom, spacing_km=1.0):
    """Extract points along the UK coastline at regular intervals."""
    geod = Geod(ellps='WGS84')

    # Get the exterior boundary
    if uk_geom.geom_type == 'MultiPolygon':
        lines = []
        for poly in uk_geom.geoms:
            lines.append(LineString(poly.exterior.coords))
        coastline = MultiLineString(lines)
    else:
        coastline = LineString(uk_geom.exterior.coords)

    # Calculate total coastline length
    if coastline.geom_type == 'MultiLineString':
        total_length_m = 0
        for line in coastline.geoms:
            coords_list = list(line.coords)
            for i in range(len(coords_list) - 1):
                lon1, lat1 = coords_list[i]
                lon2, lat2 = coords_list[i + 1]
                _, _, dist = geod.inv(lon1, lat1, lon2, lat2)
                total_length_m += dist
    else:
        coords_list = list(coastline.coords)
        total_length_m = 0
        for i in range(len(coords_list) - 1):
            lon1, lat1 = coords_list[i]
            lon2, lat2 = coords_list[i + 1]
            _, _, dist = geod.inv(lon1, lat1, lon2, lat2)
            total_length_m += dist

    total_length_km = total_length_m / 1000
    num_points = int(total_length_km / spacing_km)

    print(f"UK coastline length: ~{total_length_km:.0f} km")
    print(f"Sampling {num_points} points at {spacing_km} km intervals...")

    # Sample points along the coastline
    distances = np.linspace(0, coastline.length, num_points)
    points = [coastline.interpolate(distance) for distance in distances]

    return points, coastline


def calculate_outward_direction(point, coastline, uk_geom):
    """Calculate the outward-facing direction from a coastline point."""
    geod = Geod(ellps='WGS84')

    # Find the position of this point along the coastline
    point_distance = coastline.project(point)

    # Get a small distance ahead and behind on the coastline to estimate tangent
    delta = 0.001
    total_length = coastline.length

    # Get point slightly ahead
    ahead_dist = min(point_distance + delta, total_length)
    point_ahead = coastline.interpolate(ahead_dist)

    # Get point slightly behind
    behind_dist = max(point_distance - delta, 0)
    point_behind = coastline.interpolate(behind_dist)

    # Calculate azimuth of the tangent
    azimuth_tangent, _, _ = geod.inv(point_behind.x, point_behind.y,
                                      point_ahead.x, point_ahead.y)

    # Perpendicular directions are +90 and -90 degrees from tangent
    azimuth_perp1 = (azimuth_tangent + 90) % 360
    azimuth_perp2 = (azimuth_tangent - 90) % 360

    # Project a short distance in both perpendicular directions
    test_dist = 1000  # 1km for testing
    lon1, lat1, _ = geod.fwd(point.x, point.y, azimuth_perp1, test_dist)
    lon2, lat2, _ = geod.fwd(point.x, point.y, azimuth_perp2, test_dist)

    test_point1 = Point(lon1, lat1)
    test_point2 = Point(lon2, lat2)

    # Check which direction is outward (not contained in UK)
    in_uk1 = uk_geom.contains(test_point1)
    in_uk2 = uk_geom.contains(test_point2)

    # Pick the direction that's NOT in the UK (outward direction)
    if not in_uk1:
        return azimuth_perp1
    elif not in_uk2:
        return azimuth_perp2
    else:
        return azimuth_perp1


def find_facing_country_and_distance(point, azimuth, countries_gdf, max_distance_km=40075):
    """
    Cast a ray from the coastline point and find which country it hits and the distance.
    Excludes UK to avoid counting rays that curve back and hit the UK.
    Returns (country_name, distance_km, hit_point)
    """
    geod = Geod(ellps='WGS84')

    # Cast a ray up to Earth's circumference
    distances = np.linspace(10, max_distance_km, 400)

    for dist_km in distances:
        lon, lat, _ = geod.fwd(point.x, point.y, azimuth, dist_km * 1000)
        ray_point = Point(lon, lat)

        # Check which country this point falls in
        for idx, country in countries_gdf.iterrows():
            if country.geometry.contains(ray_point):
                country_name = country['NAME']
                # Skip UK - we only want views to other countries
                if country_name == 'United Kingdom':
                    continue
                return country_name, dist_km, ray_point

    return None, None, None


def create_great_circle_line(start_point, azimuth, distance_km, num_points=100):
    """Create a series of points along a great circle route."""
    geod = Geod(ellps='WGS84')
    points = []

    # Create points at regular intervals along the great circle
    distances = np.linspace(0, distance_km * 1000, num_points)
    for dist in distances:
        lon, lat, _ = geod.fwd(start_point.x, start_point.y, azimuth, dist)
        points.append([lat, lon])

    return points


def path_intersects_uk(start_point, azimuth, distance_km, uk_geom, num_check_points=200):
    """
    Check if a great circle path intersects the UK at any point.
    Uses dense sampling to detect intersections.
    Start checking from 10km out to avoid detecting the starting point on the coast.
    """
    geod = Geod(ellps='WGS84')

    # Check points from 10km to the destination
    # Use dense sampling to catch intersections
    distances = np.linspace(10, distance_km, num_check_points)

    for dist_km in distances:
        lon, lat, _ = geod.fwd(start_point.x, start_point.y, azimuth, dist_km * 1000)
        test_point = Point(lon, lat)

        # If any point along the path is in the UK, the path intersects
        if uk_geom.contains(test_point):
            return True

    return False


def create_map(uk, longest_views):
    """Create an interactive folium map showing the 5 longest views."""
    print("Creating map visualization...")

    # Center on UK
    m = folium.Map(
        location=[54.5, -3.5],
        zoom_start=4,
        tiles='OpenStreetMap'
    )

    # Add UK outline
    folium.GeoJson(
        uk,
        style_function=lambda x: {
            'fillColor': '#E8E8E8',
            'color': '#666666',
            'weight': 2,
            'fillOpacity': 0.3
        }
    ).add_to(m)

    # Colors for the 5 longest views
    colors = ['#E74C3C', '#3498DB', '#2ECC71', '#F39C12', '#9B59B6']

    # Add markers and great circle routes for each of the 5 longest views
    for i, view in enumerate(longest_views):
        color = colors[i]
        point = view['point']
        country = view['country']
        distance = view['distance']
        azimuth = view['azimuth']

        # Add marker at the UK coastline point
        folium.Marker(
            location=[point.y, point.x],
            popup=f"<b>View #{i+1}</b><br>To: {country}<br>Distance: {distance:.0f} km",
            tooltip=f"View #{i+1}: {country} ({distance:.0f} km)",
            icon=folium.Icon(color='red' if i == 0 else 'blue', icon='eye', prefix='fa')
        ).add_to(m)

        # Add circle marker for emphasis
        folium.CircleMarker(
            location=[point.y, point.x],
            radius=8,
            color=color,
            fill=True,
            fillColor=color,
            fillOpacity=0.7,
            weight=3
        ).add_to(m)

        # Draw full great circle route
        gc_points = create_great_circle_line(point, azimuth, distance, num_points=100)

        folium.PolyLine(
            gc_points,
            color=color,
            weight=2,
            opacity=0.6,
            popup=f"Great circle to {country} ({distance:.0f} km)"
        ).add_to(m)

        # Add marker at the destination point
        if len(gc_points) > 1:
            end_point = gc_points[-1]
            folium.CircleMarker(
                location=end_point,
                radius=6,
                color=color,
                fill=True,
                fillColor=color,
                fillOpacity=0.8,
                popup=f"Destination: {country}",
                tooltip=country
            ).add_to(m)

    # Add legend
    legend_html = '''
    <div style="position: fixed;
                top: 50px; right: 50px; width: 300px;
                background-color: white; z-index:9999; font-size:14px;
                border:2px solid grey; border-radius: 5px; padding: 15px">
    <h4 style="margin-top: 0; margin-bottom: 10px;">5 Longest Views from UK Coast</h4>
    '''

    for i, view in enumerate(longest_views):
        color = colors[i]
        country = view['country']
        distance = view['distance']
        legend_html += f'''
        <p style="margin: 8px 0;">
            <span style="background-color:{color}; padding: 5px 10px; margin-right: 8px;
                         border-radius: 3px; color: white; font-weight: bold;">#{i+1}</span>
            <b>{country}</b><br>
            <span style="margin-left: 45px; color: #666;">{distance:.0f} km</span>
        </p>
        '''

    legend_html += '</div>'
    m.get_root().html.add_child(folium.Element(legend_html))

    return m


def main():
    """Main execution function."""
    print("UK Longest Views Finder")
    print("=" * 50)

    # Load data
    uk, all_countries, world = load_geodata()

    # Get UK coastline points at 10km intervals
    uk_geom = uk.geometry.iloc[0]
    points, coastline = get_coastline_points(uk_geom, spacing_km=10.0)

    # Analyze each point
    print("Analyzing coastline points for longest views...")
    views_data = []

    for point in tqdm(points, desc="Processing coastline", unit="points"):
        # Calculate outward direction
        azimuth = calculate_outward_direction(point, coastline, uk_geom)

        # Find which country we're facing and the distance
        country, distance, hit_point = find_facing_country_and_distance(
            point, azimuth, all_countries
        )

        if country and country != 'United Kingdom':  # Exclude UK itself
            views_data.append({
                'point': point,
                'azimuth': azimuth,
                'country': country,
                'distance': distance,
                'hit_point': hit_point
            })

    # Filter out views where the great circle path intersects the UK
    print(f"\nFound {len(views_data)} potential views. Filtering out paths that intersect UK...")
    valid_views = []
    for view in tqdm(views_data, desc="Validating paths", unit="views"):
        if not path_intersects_uk(view['point'], view['azimuth'], view['distance'], uk_geom):
            valid_views.append(view)

    print(f"After filtering: {len(valid_views)} valid views remain")

    # Sort by distance and get top 5
    valid_views.sort(key=lambda x: x['distance'], reverse=True)
    top_5_views = valid_views[:5]

    # Print results
    print("\n" + "=" * 50)
    print("Top 5 Longest Views from UK Coastline:")
    print("=" * 50)
    for i, view in enumerate(top_5_views, 1):
        print(f"{i}. {view['country']}: {view['distance']:.0f} km")
        print(f"   From: ({view['point'].y:.4f}, {view['point'].x:.4f})")
        print(f"   Azimuth: {view['azimuth']:.1f}°")
        print()

    # Create and save map
    m = create_map(uk, top_5_views)
    output_file = 'longest_views_map.html'
    m.save(output_file)

    print(f"Map saved to {output_file}")
    print("Open this file in a web browser to view the interactive map!")


if __name__ == '__main__':
    main()