troggle/core/position_utils.py

import bisect
from pathlib import Path
import xml.etree.ElementTree as ET

def load_and_clean_gpx(filename):
    if not Path(filename).exists:
        print("No file")
    else:
        print(f"Loading '{filename}'")
    # Parse GPX
    tree = ET.parse(filename)
    root = tree.getroot()
    
    # Extract points (handling namespaces)
    ns = {'gpx': 'http://www.topografix.com/GPX/1/0'}
    raw_points = []
    for pt in root.findall('.//gpx:trkpt', ns):
        lat = float(pt.get('lat'))
        lon = float(pt.get('lon'))
        raw_points.append((lon, lat)) # Using (x, y) order for easier logic

    n_raw = len(raw_points)
    if not raw_points:
        return []

    # 1. Remove exactly duplicate consecutive points
    cleaned_points = [raw_points[0]]
    for i in range(1, len(raw_points)):
        if raw_points[i] != raw_points[i-1]:
            cleaned_points.append(raw_points[i])
            
    n_clean =len(cleaned_points)
    print(f"read {n_raw} points and returned {n_clean} cleaned points")
    return cleaned_points



def save_cleaned_gpx(points, output_filename="cleaned_border.gpx"):
    """
    Saves a list of (lon, lat) tuples as a GPX track.
    """
    # Create the root element with necessary namespaces
    root = ET.Element("gpx", 
                      version="1.1", 
                      creator="CaveAreaProcessor",
                      xmlns="http://www.topografix.com/GPX/1/1")
    
    # Create the track structure
    trk = ET.SubElement(root, "trk")
    name = ET.SubElement(trk, "name")
    name.text = "Cleaned Border 1626-1623"
    
    trkseg = ET.SubElement(trk, "trkseg")
    
    # Add each point as a track point
    for lon, lat in points:
        # Note: GPX uses lat/lon attributes, we provide them as strings
        ET.SubElement(trkseg, "trkpt", lat=str(lat), lon=str(lon))
    
    # Write to file
    tree = ET.ElementTree(root)
    # Use indenting for readability if available (Python 3.9+)
    if hasattr(ET, "indent"):
        ET.indent(tree, space="\t", level=0)
        
    tree.write(output_filename, encoding="windows-1252", xml_declaration=True)
    print(f"Cleaned track exported to '{output_filename}'")


def split_into_monotonic_segments(points):
    """
    Splits the track into segments that are monotonic in Longitude (E/W).
    This allows us to use binary search on each segment.
    """
    if not points: return []
    
    segments = []
    current_segment = [points[0]]
    
    for i in range(1, len(points)):
        p_prev = points[i-1]
        p_curr = points[i]
        
        # Check if we should start a new segment based on direction change in Lon
        if len(current_segment) > 1:
            prev_dir = current_segment[-1][0] - current_segment[-2][0]
            curr_dir = p_curr[0] - p_prev[0]
            
            # If direction changed (and wasn't zero before), split
            if (prev_dir > 0 and curr_dir < 0) or (prev_dir < 0 and curr_dir > 0):
                segments.append(current_segment)
                current_segment = [p_prev]
        
        current_segment.append(p_curr)
    
    segments.append(current_segment)
    print(f"{len(segments)} segments in border.")
    return segments


def generate_boundary_segments():
    points = load_and_clean_gpx('1623-6_border.gpx')
    # save_cleaned_gpx(points, "cleaned_border_output.gpx") # done once, not needed
    mono_segments = split_into_monotonic_segments(points)
    
    # Pre-calculate global bounds and segment metadata
    ALL_LONS = [p[0] for p in points]
    MIN_LON, MAX_LON = min(ALL_LONS), max(ALL_LONS)

    # Prepare segments for binary search
    # We store each segment as (sorted_lons, corresponding_lats)
    PREPARED_SEGMENTS = []
    for seg in mono_segments:
        lons = [p[0] for p in seg]
        lats = [p[1] for p in seg]
        
        # Ensure lons are strictly increasing for bisect
        is_increasing = lons[-1] > lons[0]
        if not is_increasing:
            lons.reverse()
            lats.reverse()
            
        PREPARED_SEGMENTS.append({
            'lons': lons,
            'lats': lats,
            'min_lon': min(lons),
            'max_lon': max(lons)
        })
    return PREPARED_SEGMENTS, MIN_LON, MAX_LON 

def which_area(lat, lon):
    # Initialize state on the function object itself if it doesn't exist
    # Yes, this is an attrubte *on a function*, an unusual python capbility
    if not hasattr(which_area, "data"):
        prepared_segments, min_l, max_l = generate_boundary_segments()
        which_area.data = {
            "prepared_segments": prepared_segments,
            "min_lon": min_l,
            "max_lon": max_l
        }
    
    if lon < which_area.data["min_lon"] or lon > which_area.data["max_lon"]:
        return False, "None"
    
    # Area 1626 is North, Area 1623 is South
    # We find the boundary latitude at this longitude
    boundary_lat = None
    
    for seg in which_area.data["prepared_segments"]:
        # Check if lon is within this monotonic segment
        if seg['min_lon'] <= lon <= seg['max_lon']:
            lons = seg['lons']
            lats = seg['lats']
            
            # Binary search to find the segment indices O(log n)
            idx = bisect.bisect_right(lons, lon)
            
            if idx == 0:
                boundary_lat = lats[0]
            elif idx == len(lons):
                boundary_lat = lats[-1]
            else:
                # Simple Linear Interpolation (no square roots/Pythagoras)
                # lat = y1 + (y2 - y1) * (x - x1) / (x2 - x1)
                x1, x2 = lons[idx-1], lons[idx]
                y1, y2 = lats[idx-1], lats[idx]
                
                # Basic ratio calculation
                t = (lon - x1) / (x2 - x1)
                boundary_lat = y1 + t * (y2 - y1)
            
            # Once we find the boundary at this lon, we compare
            # Note: If the line doubles back, this logic uses the first segment found.
            break 

    if boundary_lat is None:
        return False, "None"
    
    # Compare input lat to boundary lat
    # North of boundary = 1626, South = 1623
    area = "1626" if lat >= boundary_lat else "1623"
    return True, area