Introduction

Map matching involves the projection of raw geographic coordinates to map data. For instance, we can match recorded GPS data to the actual course of the road, or in real time navigation, we can match our raw position with the route we are following. Let’s see how we can match our location to a track on the map.

In this example we have a GPX track trough the town center of Cleves (Germany) (which is already map matched in GIS). The red pin indicates our raw position. We want to match it to the nearest point on the track, indicated by the green pin.

 

 

 

Finding the closest point on a line segment

To find the closest point we can use some simple vector operations. Let’s see how this works in a simple Cartesian plane, without geo-coordinates.

A and B are the closest track points. P is the raw position that we want to match.

First, we represent the line segment AB as well as AP as vectors.

Vector \(\vec{AB}\) is \(B-A=(4-1, 6-2)=(3, 4)\) and vector \(\vec{AP}\) is \(P-A=(2-1, 5-2)=(1, 3)\).

To find the closest point D on AB we should project the vector \(\vec{AP}\) onto \(\vec{AB}\).

This is done by calculating the dot product of \(\vec{AP}\) and \(\vec{AB}\), and then dividing by the square of the magnitude of \(\vec{AB}\). Mathematically, the projection scalar is given by:
\[\text{proj}_{\vec{AB}}(\vec{AP}) = \left( \frac{\vec{AP} \cdot \vec{AB}}{\| \vec{AB} \|^2} \right) \vec{AB}\]
The dot product of \(\vec{AP}\) and \(\vec{AB}\) is \(\vec{AP} \cdot \vec{AB} = (1 \times 3) + (3 \times 4) = 15\). The magnitude squared of \(\vec{AB}\) is \(\| \vec{AB} \|^2=3^2+4^2=25\). The projection scalar is then \( \frac{15}{25}=0.6\).

To find the closest point we scale \(\vec{AB}\) by the projection scalar: \(0.6 \times \vec{AB} = (1.8, 2.4)\). We add this scaled vector to A to find point D: \(D = A + (1.8, 2.4) = (2.8, 4.4)\).

Finally, we can calculate the perpendicular distance from point \(P\) to the line segment with the Euclidean distance formula: \(\sqrt{(2 – 2.8^2) + (5 – 4.4^2)} = 1\).

Since only short distances apply for our purpose, we can ignore the curvature of the earth, so we can use this linear algebra to find the closest point on the track. (Use the Haversine formula for long distances.)

 

Python implementation

Let’s see how we can implement our methodology in Python code. You can find the colab here: track_matching.ipynb.

We start by parsing a GPX file. Make sure that gpxpy is installed. Upload the GPX if you’re using colab. In this example the file is Kleve.gpx.

!pip install gpxpy

Import it, as well as numpy and folium. We’ll use Folium to display our data on the map.

import gpxpy
import folium
import numpy as np

We’ll parse the GPX to a list of route points. Then we’ll display a map, centered to the first route point.

gpx_file = open('Kleve.gpx', 'r')
gpx = gpxpy.parse(gpx_file)

route_points = []
for track in gpx.tracks:
    for segment in track.segments:
        for point in segment.points:
            route_points.append((point.latitude, point.longitude))

map_center = route_points[0]
map = folium.Map(location=map_center, zoom_start=14)

folium.PolyLine(route_points).add_to(map)

map.save('map_gpx.html')

Define a function to calculate the simple Euclidean, as discussed above. You don’t need to be a flat-earther to apply this approach to points that are close to each other.

def calculate_euclidean_distance(lat1, lon1, lat2, lon2):
    """
    Calculate the approximate distance in meters between two points (lat1, lon1) and (lat2, lon2)
    on the earth, assuming a flat earth.
    """
    METERS_PER_DEGREE_LATITUDE = 111320
    METERS_PER_DEGREE_LONGITUDE = METERS_PER_DEGREE_LATITUDE * np.cos(np.radians(lat1))

    delta_lat_meters = (lat2 - lat1) * METERS_PER_DEGREE_LATITUDE
    delta_lon_meters = (lon2 - lon1) * METERS_PER_DEGREE_LONGITUDE

    return np.sqrt(delta_lat_meters**2 + delta_lon_meters**2)

Now we can a implement function to perform our vector operations. It returns the closest point on a segment and the distance to that point from our raw location. In a real-life situation (like a navigation application) you’ll want to discard the result if the distance exceeds a certain value (like 20 m) and consider our location a deviation from the route that we don’t want to match.

def point_line_distance(point, line_start, line_end):
    """
    Calculate the perpendicular distance from `point` to the line segment `line_start`-`line_end`
    and check if the perpendicular intersection point lies within the segment.
    """
    line_vec = np.array(line_end) - np.array(line_start)
    point_vec = np.array(point) - np.array(line_start)
    line_len_squared = np.dot(line_vec, line_vec)

    projection_scale = np.dot(point_vec, line_vec) / line_len_squared
    projection_point = np.array(line_start) + projection_scale * line_vec

    if 0 <= projection_scale <= 1:
        # The projection point lies within the segment
        distance = calculate_euclidean_distance(point[0], point[1], projection_point[0], projection_point[1])
        return distance, projection_point.tolist()
    else:
        # Projection point lies outside the segment, use the nearest endpoint
        dist_to_start = calculate_euclidean_distance(point[0], point[1], line_start[0], line_start[1])
        dist_to_end = calculate_euclidean_distance(point[0], point[1], line_end[0], line_end[1])

        if dist_to_start < dist_to_end:
            return dist_to_start, line_start
        else:
            return dist_to_end, line_end

Finally, we can implement a top function that only takes the track and the raw location. It returns the matched point on the track.

def find_closest_point_on_route(route_points, my_location):
    """
    Find the closest point on the route (a series of line segments) to the given location.
    """
    closest_point = None
    min_distance = float('inf')

    for i in range(len(route_points) - 1):
        segment_start, segment_end = route_points[i], route_points[i + 1]
        distance, point_on_segment = point_line_distance(my_location, segment_start, segment_end)

        if distance < min_distance:
            min_distance = distance
            closest_point = point_on_segment

    return closest_point

To see it in action, feed it with a raw location and place a marker at the matched location it returns.

location = (51.78962, 6.14120)

folium.Marker(
    location=(location),
    icon=folium.Icon(icon="remove", color="red")
).add_to(map)

closest_point = find_closest_point_on_route(route_points, location)

folium.Marker(
    location=(closest_point),
    icon=folium.Icon(icon="ok", color="green")
).add_to(map)

map