道路折线上经纬度点间行驶距离计算方案及API需求

Alright, let's tackle this problem head-on—you need to calculate the actual road-traveled distance between two GPS points along a predefined road polyline, not straight-line distance, and you want options that either skip spatial queries (for high throughput) or use BigQuery/PostGIS if that's the only feasible path. Here's how to pull this off:

Since you mentioned high data ingestion rates, this pure in-memory method will avoid database round-trips and keep things fast. The core idea is to preprocess your road polyline to build a cumulative distance map, then use that to calculate the distance between any two projected points on the road.

Step 1: Preprocess the Road Polyline #

First, we'll precompute the total distance from the start of the road to each vertex in the polyline. This only needs to be done once per road (or when the road geometry changes):

• Store your road's GPS points in order (as a list of (lat, lon) tuples, matching the direction of travel).
• Calculate the cumulative distance array: each entry represents the total miles from the road's starting point to that vertex.

Step 2: Real-Time Distance Calculation #

For each pair of reported GPS points:

1. Project each point onto the closest segment of the road polyline.
2. Calculate the total distance from the road's start to each projected point.
3. The difference between these two totals is your traveled distance.

Here's a Python implementation to illustrate:

import math

def haversine(lat1, lon1, lat2, lon2):
    """Calculate straight-line distance between two GPS points (in miles)"""
    earth_radius_miles = 3956
    d_lat = math.radians(lat2 - lat1)
    d_lon = math.radians(lon2 - lon1)
    a = (math.sin(d_lat / 2) ** 2 +
         math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *
         math.sin(d_lon / 2) ** 2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    return earth_radius_miles * c

def preprocess_road(road_points):
    """Build cumulative distance array for the road polyline"""
    cumulative_distances = [0.0]
    for i in range(1, len(road_points)):
        lat1, lon1 = road_points[i-1]
        lat2, lon2 = road_points[i]
        segment_distance = haversine(lat1, lon1, lat2, lon2)
        cumulative_distances.append(cumulative_distances[-1] + segment_distance)
    return cumulative_distances

def get_projected_total_distance(point, road_points, cumulative_distances):
    """Find the total distance from road start to the point's projection on the road"""
    min_point_distance = float('inf')
    projected_total = 0.0

    for i in range(1, len(road_points)):
        p1 = road_points[i-1]
        p2 = road_points[i]
        seg_length = haversine(p1[0], p1[1], p2[0], p2[1])
        
        # Calculate distance from point to each segment endpoint
        dist_to_p1 = haversine(point[0], point[1], p1[0], p1[1])
        dist_to_p2 = haversine(point[0], point[1], p2[0], p2[1])
        
        # Use Heron's formula to find projection distance from p1 to the point's projection
        s = (dist_to_p1 + dist_to_p2 + seg_length) / 2
        area = math.sqrt(max(s * (s - dist_to_p1) * (s - dist_to_p2) * (s - seg_length), 0))
        proj_dist_from_p1 = (2 * area) / seg_length if seg_length != 0 else 0

        # Determine if projection is within the segment, or use closest endpoint
        if 0 <= proj_dist_from_p1 <= seg_length:
            current_total = cumulative_distances[i-1] + proj_dist_from_p1
            current_point_dist = haversine(point[0], point[1], 
                                         p1[0] + (p2[0]-p1[0])*(proj_dist_from_p1/seg_length),
                                         p1[1] + (p2[1]-p1[1])*(proj_dist_from_p1/seg_length))
        else:
            if dist_to_p1 < dist_to_p2:
                current_total = cumulative_distances[i-1]
                current_point_dist = dist_to_p1
            else:
                current_total = cumulative_distances[i]
                current_point_dist = dist_to_p2

        # Track the closest projection
        if current_point_dist < min_point_distance:
            min_point_distance = current_point_dist
            projected_total = current_total

    return projected_total

def calculateMilesBetweenPointsAlongRoad(latlon1, latlon2, road_polyline):
    """Main function to compute road distance between two points"""
    cumulative_dist = preprocess_road(road_polyline)
    dist1 = get_projected_total_distance(latlon1, road_polyline, cumulative_dist)
    dist2 = get_projected_total_distance(latlon2, road_polyline, cumulative_dist)
    return abs(dist2 - dist1)

Pros & Cons #

• Pros: Blazing fast for high-volume data, no database dependency, easy to scale in memory.
• Cons: Requires preprocessing roads (only a one-time cost), works best if each GPS point maps to a single known road (no dynamic road matching).

If you need to handle dynamic road data or match points to multiple roads automatically, using spatial databases is a solid fallback.

PostGIS Implementation #

1. Store your road polyline as a LINESTRING (with SRID 4326 for GPS coordinates) in a table (e.g., roads with columns id and geom).
2. For each GPS point pair, find the closest road segment, project the points onto it, and calculate the distance difference.

Example SQL:

WITH road_segments AS (
    SELECT 
        id,
        geom,
        ST_Length(geom::geography) / 1609.34 AS seg_length_miles, -- Convert meters to miles
        ST_StartPoint(geom) AS start_point,
        ST_EndPoint(geom) AS end_point
    FROM roads
),
point_a_proj AS (
    SELECT 
        rs.id,
        ST_LineLocatePoint(rs.geom, ST_SetSRID(ST_MakePoint(lonA, latA), 4326)) AS fraction_a,
        (rs.seg_length_miles * ST_LineLocatePoint(rs.geom, ST_SetSRID(ST_MakePoint(lonA, latA), 4326))) AS dist_from_start_a
    FROM road_segments rs
    ORDER BY ST_Distance(ST_SetSRID(ST_MakePoint(lonA, latA), 4326), rs.geom)
    LIMIT 1
),
point_b_proj AS (
    SELECT 
        rs.id,
        ST_LineLocatePoint(rs.geom, ST_SetSRID(ST_MakePoint(lonB, latB), 4326)) AS fraction_b,
        (rs.seg_length_miles * ST_LineLocatePoint(rs.geom, ST_SetSRID(ST_MakePoint(lonB, latB), 4326))) AS dist_from_start_b
    FROM road_segments rs
    ORDER BY ST_Distance(ST_SetSRID(ST_MakePoint(lonB, latB), 4326), rs.geom)
    LIMIT 1
)
SELECT ABS(pb.dist_from_start_b - pa.dist_from_start_a) AS distance_miles
FROM point_a_proj pa
JOIN point_b_proj pb ON pa.id = pb.id;

BigQuery Implementation #

BigQuery uses GEOGRAPHY types for spatial data. The approach mirrors PostGIS, with slight syntax differences:

Example SQL:

DECLARE point_a GEOGRAPHY DEFAULT ST_GEOGPOINT(lonA, latA);
DECLARE point_b GEOGRAPHY DEFAULT ST_GEOGPOINT(lonB, latB);

WITH closest_road AS (
    SELECT 
        path,
        ST_Length(path) / 1609.34 AS total_road_miles -- Convert meters to miles
    FROM roads
    ORDER BY ST_Distance(point_a, path) + ST_Distance(point_b, path)
    LIMIT 1
),
proj_a_fraction AS (
    SELECT ST_LineLocatePoint(path, point_a) AS frac_a
    FROM closest_road
),
proj_b_fraction AS (
    SELECT ST_LineLocatePoint(path, point_b) AS frac_b
    FROM closest_road
)
SELECT ABS((cr.total_road_miles * pf.frac_b) - (cr.total_road_miles * pa.frac_a)) AS distance_miles
FROM closest_road cr, proj_a_fraction pa, proj_b_fraction pf;

Pros & Cons #

• Pros: Handles dynamic road data, automatically matches points to roads, great for batch processing.
• Cons: Adds database query overhead, less ideal for ultra-high throughput real-time data.

内容的提问来源于stack exchange，提问作者Jason