gps-wizard

Line Intersection

The following code was generated by Google’s AI mode.

“calculate intersection lat lon pairs python library pypy”

Pre-requisites

pip install shapely pyproj

Python Code

Slightly different Google searches produced different examples, although very similar in nature.

Example 1

from shapely.geometry import LineString
from shapely.ops import transform
import pyproj

def calculate_line_intersection_lat_lon(line1_coords, line2_coords):
    """
    Calculates the intersection point of two lines defined by lat/lon pairs.

    Args:
        line1_coords (list): A list of two (lon, lat) tuples for the first line.
        line2_coords (list): A list of two (lon, lat) tuples for the second line.

    Returns:
        tuple or None: The (longitude, latitude) of the intersection point,
                       or None if they do not intersect.
    """
    # Define the line segments using Shapely's LineString
    line1 = LineString(line1_coords)
    line2 = LineString(line2_coords)

    # Set up the coordinate transformation
    # EPSG:4326 is WGS 84 (lat/lon)
    # A local UTM zone is better for accurate planar calculations
    project = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:32632", always_xy=True).transform

    # Project the line strings for accurate intersection calculation
    projected_line1 = transform(project, line1)
    projected_line2 = transform(project, line2)

    # Find the intersection using Shapely's `intersection` method
    intersection = projected_line1.intersection(projected_line2)

    # Check if an intersection exists and is a Point
    if intersection.is_empty or intersection.geom_type != 'Point':
        return None

    # Project the intersection point back to lat/lon
    reverse_project = pyproj.Transformer.from_crs("EPSG:32632", "EPSG:4326", always_xy=True).transform
    lon, lat = reverse_project(intersection.x, intersection.y)

    return (lon, lat)

# Define two line segments in (longitude, latitude) format
# These lines are chosen to intersect
line_a = [(-122.4194, 37.7749), (-122.0, 37.0)] # Roughly San Francisco to Santa Cruz
line_b = [(-122.5, 37.3), (-122.1, 37.8)] # A line crossing the first

# Calculate the intersection
intersection_point = calculate_line_intersection_lat_lon(line_a, line_b)

if intersection_point:
    print(f"Intersection found at (Lon, Lat): {intersection_point}")
else:
    print("No single intersection point found.")

The original code used EPSG:32632 (WGS 84 / UTM zone 32N) but has been modified to use EPSG:2154 (RGF93 v1 / Lambert-93).

Example 2

from shapely.geometry import LineString
from shapely.ops import transform
import pyproj

# Define the two lines using latitude and longitude coordinates
# Line 1: From Miami to London
line1 = LineString([(-80.1918, 25.7617), (-0.1278, 51.5074)])
# Line 2: From Bogota to Dublin
line2 = LineString([(-74.0721, 4.7110), (-6.2603, 53.3498)])

# Set up the coordinate transformation.
# This projects the WGS84 coordinates (EPSG:4326) to a planar system (e.g., Mercator).
# The projection chosen will influence accuracy, so a local one is sometimes better for small areas.
# For a global scale, a generic projection works for calculating intersections.
project_to_planar = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:2154", always_xy=True).transform

# Project the lines to a planar coordinate system
planar_line1 = transform(project_to_planar, line1)
planar_line2 = transform(project_to_planar, line2)

# Calculate the intersection in the planar system
intersection_planar = planar_line1.intersection(planar_line2)

# If an intersection exists, transform the coordinates back to WGS84
if not intersection_planar.is_empty:
    project_to_latlon = pyproj.Transformer.from_crs("EPSG:2154", "EPSG:4326", always_xy=True).transform
    intersection_latlon = transform(project_to_latlon, intersection_planar)
    
    print(f"Intersection point(s): {intersection_latlon.wkt}")
else:
    print("Lines do not intersect.")

The original code used EPSG:3395 (WGS 84 / World Mercator) but has been modified to use EPSG:2154 (RGF93 v1 / Lambert-93).

Notes

The use of EPSG:4326 is correct as it corresponds to WGS84 - World Geodetic System 1984.

Two planar coordinate systems are suitable for distance calculations in La Palme: