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(): filename = '1623-6_border.gpx' # Check if we are running inside Django try: from django.conf import settings # If settings is configured, use the troggle path if hasattr(settings, 'TROGGLE_PATH'): gpx_path = Path(settings.TROGGLE_PATH) / 'core' / filename else: gpx_path = Path(filename) except (ImportError, Exception): # Fallback for standalone testing outside of Django gpx_path = Path(filename) points = load_and_clean_gpx(gpx_path) # 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