If you have collected GPS tracks, you know that the results can have varying accuracy. The track points collected along a route are not always on the road and can be jittery.

If you are a logistics, delivery or a cab company – this poses a big problem. The distance you compute using these points will not be accurate – especially if the points are spaced apart. Also, you cannot compare tracks collected at different devices or people since they will have different geometries even if they were on the same route.

A solution to this problem is to snap each point to the nearest road segment. Though this may sound easy in principle, but doing it accurately is challenging. You cannot pick the nearest road segment for a point – because the nearest point maybe on an intersecting street. You need to consider the the route between previous and the next points to find the most plausible snapping location.

Fortunately, an open-source project called Open Source Routing Machine (OSRM) has solve this problem with fast and scalable algorithms. We can use OSRM’s match service to snap the gps points to the most appropriate road segment. OSRM engine uses data from OpenStreetMap (OSM) project. OSM has pretty good street network coverage in most parts of the world and is constantly improving. By leveraging both open data from OSM and open routing algorithms from OSRM – we can implement a snapping service.

OSRM works by taking the input via a HTTP API, computing the results and returning them via a JSON object.

Running OSRM Service

OSRM provides a demo server and a demo HTTP service. But I have found that the the demo server is often overloaded and not suitable for uses other than occasional testing.

If you want to use OSRM engine in your project, the best option is to run your own service on your computer or server. Running your own instance of a service may sound scary, but it is quite straightforward to set it up using Docker. The documentation has pretty good instructions. Here are the steps I followed to run a local instance with data for the city of Bengaluru, India.

Get the Data

An easy way to get OpenStreetMap extracts at a city-level is Interline. If you need country and continent level data, they can be downloaded from GeoFabrik.

I signed-up for a free API key and got the extract downloaded for Bangalore as bengaluru_india.osm.pbf file. I created a new folder on my system, copied the data file there, started Docker and ran the following commands in a terminal. The only change from the documentation is the –max-matching-size parameter which I increased to 5000 so we can match large GPS tracks.

docker run -t -v "${PWD}:/data" osrm/osrm-backend osrm-extract -p /opt/car.lua /data/bengaluru_india.osm.pbf

docker run -t -v "${PWD}:/data" osrm/osrm-backend osrm-partition /data/bengaluru_india.osrm

docker run -t -v "${PWD}:/data" osrm/osrm-backend osrm-customize /data/bengaluru_india.osrm

docker run -t -i -p 5000:5000 -v "${PWD}:/data" osrm/osrm-backend osrm-routed --algorithm mld --max-matching-size 5000 /data/bengaluru_india.osrm

After running the last command, a server will start on your machine and it can take requests for matching at URL http://127.0.0.1:5000

The format of match request is as follows, where the key part being the {coordinates} parameter which are the coordinates of each point on the track as a strong of the format longitude1, latitude1;longitude2, latitude2.

/match/v1/{profile}/{coordinates}?steps={true|false}&geometries={polyline|polyline6|geojson}&overview={simplified|full|false}&annotations={true|false}

We need to compile this URL programmatically by reading the GPS tracks and send it to the local match service we started in the previous step. The result also needs to be processed and converted to a track line for visualizing. This is where QGIS comes in. Using PyQGIS, we can write a processing script that makes this interaction easy and intuitive.

Matching the GPS Track

Open QGIS. Go to Processing → Toolbox → Create New Script

Copy/Paste the following code in the script editor and save it as snap_to_road.py

import requests
from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing, QgsProcessingAlgorithm, 
    QgsProcessingParameterFeatureSource, QgsProcessingParameterFeatureSink,
    QgsProcessingParameterString, QgsProcessingParameterNumber, QgsWkbTypes,
    QgsGeometry, QgsFeatureSink, QgsFields, QgsPoint, QgsFeature)
from PyQt5.QtXml import QDomDocument
class ExportLayoutAlgorithm(QgsProcessingAlgorithm):
    """Exports the current map view to PDF"""
    INPUT = 'INPUT'
    OUTPUT = 'OUTPUT'
    SERVICE = 'SERVICE'
    TOLERANCE = 'TOLERANCE'
    
    def flags(self):
          return super().flags() | QgsProcessingAlgorithm.FlagNoThreading
    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                'INPUT',
                self.tr('Input vector layer'),
                types=[QgsProcessing.TypeVectorPoint]
            )
        )
        
        self.addParameter(
            QgsProcessingParameterString(
                self.SERVICE,
                self.tr('OSRM Service URL'),
                'http://127.0.0.1:5000'
            )
        )
        
        self.addParameter(
            QgsProcessingParameterNumber(
                self.TOLERANCE,
                self.tr('Snapping Tolerance (meters)'),
                QgsProcessingParameterNumber.Integer,
                10
            )
        )
        
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                'Snapped Line',
                QgsProcessing.TypeVectorLine
            )
        )
    def processAlgorithm(self, parameters, context, feedback):
        source = self.parameterAsSource(parameters, self.INPUT, context)
        service = self.parameterAsString(parameters, self.SERVICE, context)
        tolerance = self.parameterAsInt(parameters, self.TOLERANCE, context)
        
        sink, dest_id = self.parameterAsSink(
            parameters,
            self.OUTPUT,
            context,
            QgsFields(),
            QgsWkbTypes.LineString,
            source.sourceCrs()
            )
        
        # Compute the number of steps to display within the progress bar and
        # get features from source
        total = 100.0 / source.featureCount() if source.featureCount() else 0
        features = source.getFeatures()
        
        coordinate_list = []
        for current, f in enumerate(features):
            # Stop the algorithm if cancel button has been clicked
            if feedback.isCanceled():
                break
            geom = f.geometry().asPoint()
            coordinates = '{},{}'.format(geom.x(), geom.y())
            coordinate_list.append(coordinates)
            feedback.setProgress(int(current * total))
        coordinate_str = ';'.join(coordinate_list)
        radius = ['{}'.format(tolerance)]
        radius_str = ';'.join(radius*len(coordinate_list))
        service_url = '/match/v1/driving/{}'.format(coordinate_str)
        request_url = service + service_url
        payload = {'geometries': 'geojson', 'steps': 'false', 'radiuses': radius_str}
        r = requests.get(request_url, params=payload)
        results = r.json()
        
        for match in results['matchings']:
            coords = match['geometry']['coordinates']
            point_list = [QgsPoint(coord[0], coord[1]) for coord in coords]
            out_f = QgsFeature()
            out_f.setGeometry(QgsGeometry.fromPolyline(point_list))
            sink.addFeature(out_f, QgsFeatureSink.FastInsert)
            
        return {self.OUTPUT: sink} 
    def name(self):
        return 'snap_to_roads'
    def displayName(self):
        return self.tr('Snap to Roads')
        
    def shortHelpString(self):
        return self.tr('Snaps GPS Trackpoints to OSM roads using OSRM service')
    def group(self):
        return self.tr(self.groupId())
    def groupId(self):
        return ''
    def tr(self, string):
        return QCoreApplication.translate('Processing', string)
    def createInstance(self):
        return ExportLayoutAlgorithm()

Once saved, a new algorithm will appear in Processing → Toolbox → Scripts → Snap To Roads. Load your GPS track points in QGIS and double-click the script to run it.

The resulting snapped road line will be added to the QGIS Layers panel. You can see that OSRM worked like a charm and results are exactly as one would expect.

If you want to try out the algorithm, you can download sample_gps_track.gpx. Get the Bengaluru OSM extract from Interline. Do leave a comment and let me know if you run into issues.