Hastebin

#!/usr/bin/env python3
"""
Route Buffer Coverage Generator

This script takes an existing GPX route file and generates a return path with a zig-zag pattern
that ensures coverage of a specified buffer zone around the original route. This is useful for
forcing navigation apps to download all map tiles in the vicinity of a planned route.

Usage:
    python route_buffer_coverage.py input_gpx buffer_width spacing output_gpx

Arguments:
    input_gpx:     Path to the input GPX file containing the original route
    buffer_width:  Width of the buffer zone to cover on each side of the route (in meters)
    spacing:       Distance between zig-zag lines in meters (default: 200)
    output_gpx:    Path to save the output GPX file
"""

import argparse
import math
import xml.etree.ElementTree as ET
import numpy as np
from datetime import datetime

# Define GPX namespaces
NAMESPACES = {
    'gpx': 'http://www.topografix.com/GPX/1/1',
    'xsi': 'http://www.w3.org/2001/XMLSchema-instance'
}

def haversine_distance(lat1, lon1, lat2, lon2):
    """Calculate the great-circle distance between two points in meters."""
    R = 6371000  # Earth radius in meters
    
    # Convert decimal degrees to radians
    lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])
    
    # Haversine formula
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    c = 2 * math.asin(math.sqrt(a))
    
    return R * c

def calculate_bearing(lat1, lon1, lat2, lon2):
    """Calculate the bearing from point 1 to point 2 in degrees."""
    # Convert to radians
    lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])
    
    # Calculate bearing
    y = math.sin(lon2 - lon1) * math.cos(lat2)
    x = math.cos(lat1) * math.sin(lat2) - math.sin(lat1) * math.cos(lat2) * math.cos(lon2 - lon1)
    bearing = math.atan2(y, x)
    
    # Convert to degrees
    bearing = math.degrees(bearing)
    bearing = (bearing + 360) % 360
    
    return bearing

def calculate_perpendicular_bearing(bearing):
    """Calculate perpendicular bearings (90 degrees to left and right)."""
    left_bearing = (bearing - 90) % 360
    right_bearing = (bearing + 90) % 360
    return left_bearing, right_bearing

def calculate_new_point(lat, lon, bearing, distance):
    """Calculate new GPS coordinates given a starting point, bearing (in degrees), and distance (in meters)."""
    R = 6371000  # Earth radius in meters
    
    # Convert to radians
    lat1 = math.radians(lat)
    lon1 = math.radians(lon)
    bearing = math.radians(bearing)
    
    # Calculate angular distance
    angular_distance = distance / R
    
    # Calculate new latitude
    lat2 = math.asin(math.sin(lat1) * math.cos(angular_distance) +
                     math.cos(lat1) * math.sin(angular_distance) * math.cos(bearing))
    
    # Calculate new longitude
    lon2 = lon1 + math.atan2(math.sin(bearing) * math.sin(angular_distance) * math.cos(lat1),
                           math.cos(angular_distance) - math.sin(lat1) * math.sin(lat2))
    
    # Convert back to degrees
    lat2 = math.degrees(lat2)
    lon2 = math.degrees(lon2)
    
    return lat2, lon2

def parse_gpx_route(gpx_file):
    """Parse the input GPX file and extract route points."""
    try:
        tree = ET.parse(gpx_file)
        root = tree.getroot()
        
        waypoints = []
        
        # Try to find trackpoints
        for trkpt in root.findall('.//gpx:trkpt', NAMESPACES):
            lat = float(trkpt.get('lat'))
            lon = float(trkpt.get('lon'))
            waypoints.append((lat, lon))
        
        # If no trackpoints found, try routepoints
        if not waypoints:
            for rtept in root.findall('.//gpx:rtept', NAMESPACES):
                lat = float(rtept.get('lat'))
                lon = float(rtept.get('lon'))
                waypoints.append((lat, lon))
        
        return waypoints
    except Exception as e:
        print(f"Error parsing GPX file: {e}")
        return []

def generate_buffer_zigzag(route_points, buffer_width, spacing):
    """Generate a zig-zag pattern to cover the buffer zone around the route."""
    if len(route_points) < 2:
        print("Error: Route must contain at least 2 points.")
        return []
    
    # Start with the original route
    result_points = route_points.copy()
    zigzag_points = []
    
    # Calculate accumulated distance along the route
    accumulated_distances = [0]
    for i in range(1, len(route_points)):
        prev_lat, prev_lon = route_points[i-1]
        curr_lat, curr_lon = route_points[i]
        segment_distance = haversine_distance(prev_lat, prev_lon, curr_lat, curr_lon)
        accumulated_distances.append(accumulated_distances[-1] + segment_distance)
    
    total_route_length = accumulated_distances[-1]
    
    # Calculate the interval for perpendicular crossings
    # Use buffer_width as a basis for the interval to ensure good coverage
    interval = buffer_width
    
    # Number of perpendicular crossings
    num_crossings = max(3, int(total_route_length / interval))
    
    # Create perpendicular lines at regular intervals
    for i in range(num_crossings + 1):
        # Calculate distance along the route for this crossing
        target_distance = (i * total_route_length) / num_crossings
        
        # Find the segment containing this distance
        segment_idx = 0
        while segment_idx < len(accumulated_distances) - 1 and accumulated_distances[segment_idx + 1] < target_distance:
            segment_idx += 1
        
        # If we're at the end of the route
        if segment_idx >= len(route_points) - 1:
            segment_idx = len(route_points) - 2
        
        # Calculate interpolation factor within this segment
        segment_start_dist = accumulated_distances[segment_idx]
        segment_end_dist = accumulated_distances[segment_idx + 1]
        segment_length = segment_end_dist - segment_start_dist
        
        # Avoid division by zero for very short segments
        if segment_length < 0.1:
            interpolation_factor = 0
        else:
            interpolation_factor = (target_distance - segment_start_dist) / segment_length
        
        # Interpolate to find the exact crossing point
        start_lat, start_lon = route_points[segment_idx]
        end_lat, end_lon = route_points[segment_idx + 1]
        
        crossing_lat = start_lat + interpolation_factor * (end_lat - start_lat)
        crossing_lon = start_lon + interpolation_factor * (end_lon - start_lon)
        
        # Calculate bearing of the route at this point
        bearing = calculate_bearing(start_lat, start_lon, end_lat, end_lon)
        
        # Calculate perpendicular bearings
        left_bearing, right_bearing = calculate_perpendicular_bearing(bearing)
        
        # Create perpendicular line crossing the route at this point
        crossing_points = []
        
        # Center point (on the route)
        crossing_points.append((crossing_lat, crossing_lon))
        
        # Left side points (multiple points to ensure full coverage)
        for dist in range(0, int(buffer_width) + 1, int(spacing)):
            new_lat, new_lon = calculate_new_point(crossing_lat, crossing_lon, left_bearing, dist)
            crossing_points.append((new_lat, new_lon))
        
        # Right side points (multiple points to ensure full coverage)
        for dist in range(0, int(buffer_width) + 1, int(spacing)):
            new_lat, new_lon = calculate_new_point(crossing_lat, crossing_lon, right_bearing, dist)
            crossing_points.append((new_lat, new_lon))
        
        # Add this set of crossing points to our zigzag collection
        zigzag_points.extend(crossing_points)
    
    # Create a route back to the start using the zigzag points
    # First, get the first and last points of the original route
    start_point = route_points[0]
    
    # Add a path back to start from the end of the zigzags
    result_points.extend(zigzag_points)
    result_points.append(start_point)  # Close the loop
    
    return result_points

def write_gpx(points, output_file):
    """Write the points to a GPX file."""
    root = ET.Element('gpx')
    root.set('version', '1.1')
    root.set('creator', 'Route Buffer Coverage Generator')
    root.set('xmlns', 'http://www.topografix.com/GPX/1/1')
    
    # Add metadata
    metadata = ET.SubElement(root, 'metadata')
    name = ET.SubElement(metadata, 'name')
    name.text = 'Route with Buffer Coverage'
    time = ET.SubElement(metadata, 'time')
    time.text = datetime.utcnow().strftime('%Y-%m-%dT%H:%M:%SZ')
    
    # Create track
    trk = ET.SubElement(root, 'trk')
    trk_name = ET.SubElement(trk, 'name')
    trk_name.text = 'Combined Route'
    trkseg = ET.SubElement(trk, 'trkseg')
    
    # Add track points
    for lat, lon in points:
        trkpt = ET.SubElement(trkseg, 'trkpt')
        trkpt.set('lat', str(lat))
        trkpt.set('lon', str(lon))
    
    # Write to file
    tree = ET.ElementTree(root)
    ET.indent(tree, space="  ")  # Pretty print XML (Python 3.9+)
    tree.write(output_file, encoding='utf-8', xml_declaration=True)

def main():
    parser = argparse.ArgumentParser(description="Generate a return path with zig-zags to cover a buffer zone around an existing route.")
    parser.add_argument("input_gpx", help="Path to the input GPX file containing the original route")
    parser.add_argument("buffer_width", type=float, help="Width of the buffer zone to cover on each side of the route (in meters)")
    parser.add_argument("-s", "--spacing", type=float, default=200.0, help="Spacing between points along perpendicular lines (default: 200)")
    parser.add_argument("-i", "--interval", type=float, default=None, help="Distance between perpendicular crossings in meters (default: equal to buffer_width)")
    parser.add_argument("output_gpx", help="Path to save the output GPX file")
    
    args = parser.parse_args()
    
    # Parse input GPX
    route_points = parse_gpx_route(args.input_gpx)
    if not route_points:
        print("Error: No route points found in the input GPX file.")
        return
    
    print(f"Parsed {len(route_points)} points from input route.")
    
    # Set interval if not specified
    interval = args.interval if args.interval else args.buffer_width
    
    # Generate buffer coverage
    result_points = generate_buffer_zigzag(route_points, args.buffer_width, args.spacing)
    
    # Write output GPX
    write_gpx(result_points, args.output_gpx)
    
    print(f"Generated route with buffer coverage ({len(result_points)} points).")
    print(f"Buffer width: {args.buffer_width}m")
    print(f"Crossing interval: {interval}m")
    print(f"Output written to: {args.output_gpx}")

if __name__ == "__main__":
    main()