[GIS] Polyline matching python

gpsNetworkpythonroad

I have two polylines, the first consists of individual GPS traces and the second is a series of road network poylines.

I wish to match the GPS trace polyline to that of a road network polyline or series of polylines in order to pass the time attribute from the GPS traces to the road network poyline(s). Of course, for a range of accuracy related reasons the GPS trace does not perfectly match the road network.

This may be visually inspected. The blue line shows the GPS trace polyline and the gray represents the road network polylines.

Due to the complexity of road junctions and directionality restrictions, a nearest neighbour match for the range of coordinates within the polyline won't work.

An example GPS trace polyline looks like this:

GPS_trace = [   
"polyline_decoded":[
          [51.5201035,-0.0965073],[51.5201,-0.09651],[51.52004,-0.09693],[51.5201406,-0.096996]],
"duration":15]

An example road polyline looks like this:

road_network_link = [
"polyline" : [
[51.5223,-0.0965],[51.521,-0.097],[51.52,-0.0979],[51.5212,-0.0970]]
"id" : "10200A"]

This appears to be the problem a good solution – https://developers.google.com/maps/documentation/roads/snap but I'd prefer to try myself first.

I have a large amount of GPS polylines and I am thus looking for a matching process that scales well. I'm aware this functionality exists in ArcGIS and QGIS.

My first inclination is to perform a shortest path algorithm for the GPS polyline points via the road link network network. Can anyone point me in the right direction?

Best Answer

I created route, placed points on it 250 m apart, and randomly spread them from original position using normal distribution with standard deviation of 30 m:

Iterations:

Compute shortest path along network between 2 network nodes nearest to first and last “GPS” points
Count number of links in that path, that are more than 100 m away from GPS points
Break if none (assume path found), otherwise increase the cost of travel through selected links by factor of 100 and go to 1

After 3 iterations:

Algorithm converged to this solution:

Related Solutions

[GIS] Automatic matching of two street networks

A robust method to match networks is described in Mustière, S., Devogele, T., Dec. 2008. Matching networks with different levels of detail. GeoInformatica 12 (4), 435-453.. It has been used at the French national mapping agency to match 2 geographical databases with different levels of detail (see image below). The purpose was to do exactly what you need: Transfer of attributes.

network matching algorithm

This matching process compares the two network elements taking into account geometrical and also topological criteria: Network elements are not matched only if they are closed to each other (using Hausfdorff distance) and with comparable shapes, but also if they are connected to other network elements that are matched together. One-to-many relations are used.

A good news: This process is implemented in the opensource GéOxygene library. This document describes how to use it. Bad news: You have to speak both Java and French to use it...

[GIS] Matching GPS tracks

As @Loxodromes said above, I too am not sure that an open source library exists to do this. But it's simple enough to do in Python if you're happy enough with the scripting. For example, if you have access to numpy and scipy you can use a KDTree to easily calculate points from trail A that are within some tolerance of points from trail B.

With a bit of work you can take this a bit further by stacking the points into a single array and playing with labelled groups. This has the bonus of coping with more than two base data sets for comparison, though note this is not memory friendly - if you've got a lot of points you might need to do some work to make this more memory efficient. This also assumes everything is in the same projection.

import numpy as np
import scipy.spatial

For this example I'll dummy up some data, but take a look at numpy.loadtxt to read in your CSVs.

np.random.seed(20140201)
num_pts = 50
points_a = np.vstack([
    np.linspace(0., 10., num=num_pts),
    np.linspace(10., 0., num=num_pts)
    ]).T

points_b = points_a + np.random.random([num_pts, 2]) - 0.5
points_c = points_a + np.random.random([num_pts, 2]) - 0.5
points_d = points_a + np.vstack([
    np.sin(np.linspace(0., 2 * np.pi, num_pts)),
    np.sin(np.linspace(0., 2 * np.pi, num_pts)),
    ]).T

all_trails = [points_a, points_b, points_c, points_d]

You'll also need to specify a tolerance

tolerance = 0.1

Then, so you can process all the points in bulk but still know what group they're in, stack the arrays.

labelled_pts = np.vstack([
    np.hstack([a, np.ones((a.shape[0], 1)) * i])
    for i, a in enumerate(all_trails)
])

You can now build a KDTree from the labelled points. Remember that you don't want the labels themselves in the tree - they're used later on to classify results

tree = scipy.spatial.KDTree(labelled_pts[:, :2])

You use the ball point algorithm to get all the points within tolerance of another set of points (which is conveniently also our input points).

points_within_tolerance = tree.query_ball_point(labelled_pts[:, :2], tolerance)

This returns an array of the same length as the incoming points, with each value in the array being a tuple of indexes of the found points in the tree. Because you put in our original set there will always be at least one match. However you can then build a simple vectorisation function to test whether each item in the tree matches a point from a different group.

vfunc = np.vectorize(lambda a: np.any(labelled_pts[a, 2] != labelled_pts[a[0], 2]))

matches = vfunc(points_within_tolerance)
matching_points = labelled_pts[matches, :2]

The vfunc simply returns a numpy array of the results of this function, in this case True or False which we can use to index out our points.

So now you have points on the GPS trails which cross, but you want to group points into contiguous segments of track that overlap. For that you can use the scipy hierarchical clustering methods to group the data into groups which are linked by at most the tolerance distance.

import scipy.cluster.hierarchy

clusters = scipy.cluster.hierarchy.fclusterdata(matching_points, tolerance, 'distance')

clusters is an array of the same length of your matched points containing cluster indexes for each point. This means it's easy to get back a table of x, y, original_trail, segment by stacking the output together.

print np.hstack([
    matching_points,              #x, y
    np.vstack([
        labelled_pts[matches, 2], #original_trail
        clusters                  #segment
    ]).T
])

Or you can draw up the clusters.

from itertools import cycle, izip
import matplotlib.pyplot as plt

for pts, colour in izip(all_trails, cycle(['blue', 'red', 'orange', 'green', 'pink'])):
    plt.scatter(pts[:, 0], pts[:, 1], c=colour)

for clust_idx, shape, size in izip(set(clusters), cycle(['o', 'v', '^', '<', '>', 's', 'p', '*', '8', 'd']), cycle([40, 50, 60])):
    plt.scatter(matching_points[clusters == clust_idx, 0], matching_points[clusters == clust_idx, 1], c='yellow', marker=shape, s=size)

plt.show()

Sample output of the clustering. The yellow shapes are the different clusters.

Hopefully this all makes sense!

Best Answer

Related Solutions

[GIS] Automatic matching of two street networks

[GIS] Matching GPS tracks

Related Question