find_longest_view_systematic.py

"""
Systematic Algorithm to Find Longest Unobstructed View from England Coast
- Sample coastline at configurable intervals (default: 100km)
- Two-phase checking for speed:
  * Phase 1: Coarse check every 100km to filter out paths < 7000km
  * Phase 2: Detailed check every 50m for promising paths only
- Parallel processing using all available CPU cores
- Comprehensive timing information
"""

import geopandas as gpd
import numpy as np
from pyproj import Geod
from shapely.geometry import Point, LineString, MultiLineString
from pathlib import Path
from tqdm import tqdm
import json
import time
from datetime import datetime
from multiprocessing import Pool, cpu_count
from functools import partial
import argparse


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...")
    import requests
    import zipfile
    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'

    with open(zip_path, 'wb') as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)

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

    zip_path.unlink()
    return str(shapefile_path)


def load_geodata():
    """Load England and world countries data from Natural Earth."""
    print("Loading geographical data...")
    shapefile_path = download_natural_earth_data()
    world = gpd.read_file(shapefile_path)

    # Load detailed UK regions
    data_dir = Path('data')
    admin1_path = data_dir / 'ne_10m_admin_1_states_provinces.shp'

    if not admin1_path.exists():
        print("Downloading detailed UK regions data...")
        import requests
        import zipfile
        url = 'https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_1_states_provinces.zip'
        response = requests.get(url, stream=True)
        zip_path = data_dir / 'ne_admin1.zip'

        with open(zip_path, 'wb') as f:
            for chunk in response.iter_content(chunk_size=8192):
                f.write(chunk)

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

        zip_path.unlink()

    admin1 = gpd.read_file(str(admin1_path))
    england = admin1[admin1['name'] == 'England'].copy()

    if len(england) == 0:
        print("Warning: England not found, using UK")
        england = world[world['NAME'] == 'United Kingdom'].copy()

    return england, world


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

    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 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)

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

    return points, coastline, total_length_km


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

    ahead_dist = min(point_distance + delta, total_length)
    point_ahead = coastline.interpolate(ahead_dist)
    behind_dist = max(point_distance - delta, 0)
    point_behind = coastline.interpolate(behind_dist)

    azimuth_tangent, _, _ = geod.inv(point_behind.x, point_behind.y,
                                      point_ahead.x, point_ahead.y)

    azimuth_perp1 = (azimuth_tangent + 90) % 360
    azimuth_perp2 = (azimuth_tangent - 90) % 360

    test_dist = 1000
    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)

    in_uk1 = uk_geom.contains(test_point1)
    in_uk2 = uk_geom.contains(test_point2)

    if not in_uk1:
        return azimuth_perp1
    elif not in_uk2:
        return azimuth_perp2
    else:
        return azimuth_perp1


def coarse_ray_check(start_point, azimuth, min_distance_km, all_countries):
    """
    Quick coarse check every 100km to filter out unpromising paths.
    Returns True if path looks good (no land hit before min_distance_km), False otherwise.
    """
    geod = Geod(ellps='WGS84')

    # Check every 100km up to the minimum distance threshold
    coarse_interval_km = 100
    num_checks = int(min_distance_km / coarse_interval_km)

    for i in range(1, num_checks + 1):
        distance_km = i * coarse_interval_km

        lon, lat, _ = geod.fwd(start_point.x, start_point.y, azimuth, distance_km * 1000)
        test_point = Point(lon, lat)

        # Check which country this point is in
        possible_matches_idx = list(all_countries.sindex.intersection(test_point.bounds))
        possible_matches = all_countries.iloc[possible_matches_idx]

        for idx, country in possible_matches.iterrows():
            if country.geometry.contains(test_point):
                # Hit ANY land before minimum distance - reject this path
                return False

    # No land hits in coarse check - looks promising
    return True


def check_ray_intersections(start_point, azimuth, max_distance_km, all_countries, check_interval_m=50):
    """
    Cast a ray and check for land intersections at regular intervals.
    Returns (target_country, distance_km) if clear path found, else (None, None)
    """
    geod = Geod(ellps='WGS84')

    # Check every 50m from start to max distance
    num_checks = int((max_distance_km * 1000) / check_interval_m)

    target_country = None
    target_distance_km = None

    for i in range(1, num_checks):
        distance_m = i * check_interval_m
        distance_km = distance_m / 1000

        lon, lat, _ = geod.fwd(start_point.x, start_point.y, azimuth, distance_m)
        test_point = Point(lon, lat)

        # Check which country this point is in
        possible_matches_idx = list(all_countries.sindex.intersection(test_point.bounds))
        possible_matches = all_countries.iloc[possible_matches_idx]

        for idx, country in possible_matches.iterrows():
            if country.geometry.contains(test_point):
                country_name = country['NAME']

                # Skip UK
                if country_name == 'United Kingdom':
                    return None, None, None

                # If we haven't found a target yet, this is it
                if target_country is None:
                    target_country = country_name
                    target_distance_km = distance_km
                    return target_country, target_distance_km, test_point

                # If we found a different country, there's an intersection
                if country_name != target_country:
                    return None, None, None

    # If we found a target and no intersections, return it
    if target_country:
        return target_country, target_distance_km, test_point

    return None, None, None


def check_single_point(point, coastline, uk_geom, all_countries, max_distance_km=20000, min_distance_km=7000.0, check_interval_m=50):
    """
    Check the perpendicular outward direction from a single coastline point.
    Uses two-phase checking: coarse filter first, then detailed check.
    Returns view info if valid, None otherwise.
    """
    azimuth = calculate_outward_direction(point, coastline, uk_geom)

    # Phase 1: Coarse check - quickly filter out paths that won't beat current best
    if not coarse_ray_check(point, azimuth, min_distance_km, all_countries):
        # Path hits land before minimum distance - skip detailed check
        return None

    # Phase 2: Detailed check - only for promising paths
    country, distance, hit_point = check_ray_intersections(
        point, azimuth, max_distance_km, all_countries, check_interval_m=check_interval_m
    )

    if country and distance:
        return {
            'point': point,
            'azimuth': azimuth,
            'country': country,
            'distance': distance,
            'hit_point': hit_point
        }

    return None


def process_point_worker(point, coastline, uk_geom, all_countries, max_distance_km=20000, min_distance_km=7000.0, check_interval_m=50):
    """
    Worker function for parallel processing.
    Wraps check_single_point for use with multiprocessing.
    """
    try:
        return check_single_point(point, coastline, uk_geom, all_countries, max_distance_km, min_distance_km, check_interval_m)
    except Exception as e:
        print(f"Error processing point {point}: {e}")
        return None


def main(spacing_km=100.0, min_threshold_km=7000.0, check_interval_m=50):
    """Main execution with systematic checking."""
    start_time = time.time()

    print("=" * 70)
    print("SYSTEMATIC ALGORITHM: Find Longest Unobstructed View")
    print("=" * 70)
    print(f"Start time: {datetime.now().strftime('%H:%M:%S')}")
    print(f"Coastline spacing: {spacing_km} km")
    print(f"Coarse check threshold: {min_threshold_km} km")
    print(f"Fine-grain check interval: {check_interval_m} m")
    print()

    # Load data
    t0 = time.time()
    print("STEP 1: Loading geographical data...")
    england, all_countries = load_geodata()
    print(f"  ✓ Loaded in {time.time()-t0:.1f}s")

    t0 = time.time()
    print("STEP 2: Building spatial index...")
    _ = all_countries.sindex
    print(f"  ✓ Built in {time.time()-t0:.1f}s")

    uk_geom = england.geometry.iloc[0]

    # Sample coastline
    t0 = time.time()
    print(f"STEP 3: Sampling coastline at {spacing_km}km intervals...")
    coastline_points, coastline, total_length_km = get_coastline_points(uk_geom, spacing_km=spacing_km)
    print(f"  ✓ Coastline length: {total_length_km:.0f} km")
    print(f"  ✓ Generated {len(coastline_points)} sample points")
    print(f"  ✓ Sampled in {time.time()-t0:.1f}s")

    # Main scan with parallel processing
    print()
    print("=" * 70)
    print("STEP 4: SYSTEMATIC SCAN (PARALLEL + TWO-PHASE)")
    print("  - Checking perpendicular outward direction from each point")
    print("  - Phase 1: Coarse check every 100km (quick filter)")
    print(f"  - Phase 2: Detailed check every {check_interval_m}m (only if promising)")
    num_cores = cpu_count()
    print(f"  - Using {num_cores} CPU cores")
    print("=" * 70)
    print()

    all_views = []
    scan_start = time.time()

    # Create partial function with fixed parameters
    worker_func = partial(
        process_point_worker,
        coastline=coastline,
        uk_geom=uk_geom,
        all_countries=all_countries,
        max_distance_km=20000,
        min_distance_km=min_threshold_km,
        check_interval_m=check_interval_m
    )

    # Process points in parallel
    with Pool(processes=num_cores) as pool:
        # Use imap for progress tracking
        results = []
        for i, result in enumerate(tqdm(
            pool.imap(worker_func, coastline_points),
            total=len(coastline_points),
            desc="Scanning coastline points",
            unit="point"
        )):
            if result:
                all_views.append(result)
            results.append(result)

            # Progress update every 10 points
            if (i + 1) % 10 == 0:
                elapsed = time.time() - scan_start
                rate = (i + 1) / elapsed
                remaining = (len(coastline_points) - i - 1) / rate
                print(f"  [{i+1}/{len(coastline_points)}] "
                      f"Rate: {rate:.2f} pts/s, "
                      f"ETA: {remaining:.1f}s, "
                      f"Found: {len(all_views)} views")

    scan_time = time.time() - scan_start
    print(f"\n  ✓ Scan completed in {scan_time/60:.1f} minutes")
    print(f"  ✓ Found {len(all_views)} valid views")

    if len(all_views) == 0:
        print("\n❌ No valid views found!")
        return

    # Find longest
    all_views.sort(key=lambda x: x['distance'], reverse=True)
    longest_view = all_views[0]

    # Results
    print()
    print("=" * 70)
    print("RESULTS: LONGEST UNOBSTRUCTED VIEW FROM ENGLAND")
    print("=" * 70)
    print(f"Country: {longest_view['country']}")
    print(f"Distance: {longest_view['distance']:.2f} km")
    print(f"From: ({longest_view['point'].y:.4f}°N, {longest_view['point'].x:.4f}°E)")
    print(f"Azimuth: {longest_view['azimuth']:.1f}°")
    print("=" * 70)

    # Save results
    geod = Geod(ellps='WGS84')
    end_lon, end_lat, _ = geod.fwd(
        longest_view['point'].x, longest_view['point'].y,
        longest_view['azimuth'], longest_view['distance'] * 1000
    )

    json_data = {
        'longest_view': {
            'country': longest_view['country'],
            'distance_km': float(longest_view['distance']),
            'start': {
                'lon': float(longest_view['point'].x),
                'lat': float(longest_view['point'].y)
            },
            'end': {
                'lon': float(end_lon),
                'lat': float(end_lat)
            },
            'azimuth': float(longest_view['azimuth'])
        },
        'top_10_for_reference': [
            {
                'country': v['country'],
                'distance_km': float(v['distance']),
                'start': {'lon': float(v['point'].x), 'lat': float(v['point'].y)},
                'azimuth': float(v['azimuth'])
            }
            for v in all_views[:10]
        ],
        'metadata': {
            'coastline_spacing_km': float(spacing_km),
            'coarse_check_interval_km': 100,
            'fine_check_interval_m': int(check_interval_m),
            'min_distance_threshold_km': float(min_threshold_km),
            'method': 'two_phase_perpendicular_outward',
            'parallel_processing': True,
            'num_cores': num_cores,
            'total_scan_time_seconds': scan_time,
            'total_views_found': len(all_views)
        }
    }

    # Generate unique filename with timestamp and parameters
    timestamp = datetime.now().strftime('%Y%m%d_%H%M%S')
    output_file = f'longest_view_result_{timestamp}_spacing{spacing_km}km_check{check_interval_m}m_thresh{min_threshold_km}km.json'
    with open(output_file, 'w') as f:
        json.dump(json_data, f, indent=2)

    total_time = time.time() - start_time
    print(f"\n✓ Results saved to {output_file}")
    print(f"✓ Total execution time: {total_time/60:.1f} minutes")
    print(f"✓ Finished at: {datetime.now().strftime('%H:%M:%S')}")


if __name__ == '__main__':
    parser = argparse.ArgumentParser(
        description='Find longest unobstructed view from England coast',
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument(
        '--spacing',
        type=float,
        default=100.0,
        help='Coastline sampling interval in km (smaller = more accurate but slower)'
    )
    parser.add_argument(
        '--threshold',
        type=float,
        default=5000.0,
        help='Minimum distance threshold for coarse check in km (paths hitting land before this are filtered out)'
    )
    parser.add_argument(
        '--check-interval',
        type=int,
        default=50,
        help='Fine-grain checking interval in meters along the ray path'
    )

    args = parser.parse_args()
    main(spacing_km=args.spacing, min_threshold_km=args.threshold, check_interval_m=args.check_interval)