Map-Matching for Trajectory Prediction | by João Paulo Figueira | Jul, 2023

Where are you going? Should you be going that way?

João Paulo Figueira
Towards Data Science
Photo by Mateusz Wacławek on Unsplash

This article presents a method to predict vehicle trajectories on a digital road network using a database of past trips sampled from noisy GPS sensors. Besides predicting future directions, this method also assigns a probability to an arbitrary sequence of locations.

Central to this idea is using a digital map unto which we project all sampled locations by aggregating them into individual trajectories and matching them to the map. This matching process discretizes the continuous GPS space into predetermined locations and sequences. After encoding these locations into unique geospatial tokens, we can more easily predict sequences, evaluate the probability of current observations and estimate future directions. This is the gist of this article.

What problems am I trying to solve here? If you need to analyze vehicle path data, you might need to answer questions like those in the article’s sub-heading.

Where are you going? Should you be going that way?

How do you evaluate the probability that the path under observation follows frequently traveled directions? This is an important question as, by answering it, you could program an automated system to classify trips according to their observed frequency. A new trajectory with a low score would cause concern and prompt immediate flagging.

How do you predict which maneuvers the vehicle will do next? Will it keep going straight ahead, or will it turn right at the next intersection? Where do you expect to see the vehicle in the next ten minutes or ten miles? Quick answers to these questions will assist an online tracking software solution in providing answers and insights to delivery planners, online route optimizers, or even opportunity charging systems.

The solution I am presenting here uses a database of historical trajectories, each consisting of a timed sequence of positions generated by the motion of a specific vehicle. Each positional record must contain time, position information, a reference to the vehicle identifier, and the trajectory identifier. A vehicle has many trajectories, and each trajectory has many positional records. A sample of our input data is depicted in Figure 1 below.

Figure 1 — The table above shows a small sample of a trajectory from the Extended Vehicle Energy Dataset. Although each row contains more properties than the ones displayed, we only need the implicit order and the locations. Note that there are many duplicated locations due to the sampling strategy. We will handle this issue later on. (Image source: Author)

I drew the data above from the Extended Vehicle Energy Dataset (EVED) [1] article. You can build the corresponding database by following the code in one of my previous articles.

Our first job is to match these trajectories to a supporting digital map. The purpose of this step is not only to eliminate the GPS data sampling errors but, most importantly, to coerce the acquired trip data to an existing road network where each node and edge are known and fixed. Each recorded trajectory is thus converted from a sequence of geospatial locations into another sequence of numeric tokens coinciding with the existing digital map nodes. Here, we will use open-sourced data and software, with map data sourced from OpenStreetMap (compiled by Geofabrik), the Valhalla map-matching package, and H3 as the geospatial tokenizer.

Edge Versus Node Matching

Map-matching is more nuanced than it might look at first sight. To illustrate what this concept entails, let us look at Figure 2 below.

Figure 2 — The map above shows a noisy trajectory taken from the EVED in blue. As you can see, it does not match the nearest road and needs matching to the map. Once we project the blue line vertices to the map, we obtain a sequence of projections of the original points to the inferred map edges, and you can see the resulting trajectory in red. Still, this path misses the underlying map in some places, most notably in the image’s center, where the red line jumps between roads. We aim to reconstruct the trip’s path on the map, as represented by the green line, that follows the underlying map nodes. (Image source: Author using Folium and OpenStreetMap imagery)

Figure 2 above shows that we can derive two trajectories from an original GPS sequence. We obtain the first trajectory by projecting the original GPS locations into the nearest (and most likely) road network segments. As you can see, the resulting polyline will only sometimes follow the road because the map uses graph nodes to define its basic shapes. By projecting the original locations to the map edges, we get new points that belong to the map but may stray from the map’s geometry when connected to the next ones by a straight line.

By projecting the GPS trajectory to the map nodes, we get a path that perfectly overlays the map, as shown by the green line in Figure 2. Although this path better represents the initially driven trajectory, it does not necessarily have a one-to-one location correspondence with the original. Fortunately, this will be fine for us as we will always map-match any trajectory to the map nodes, so we will always get coherent data, with one exception. The Valhalla map-matching code always edge-projects the initial and final trajectory points, so we will systematically discard them as they do not correspond to map nodes.

H3 Tokenization

Unfortunately, Valhalla does not report the unique road network node identifiers, so we must convert the node coordinates to unique integer tokens for later sequence frequency calculation. This is where H3 enters the picture by allowing us to encode the node coordinates into a sixty-four-bit integer uniquely. We pick the Valhalla-generated polyline, strip the initial and final points (these points are not nodes but edge projections), and map all remaining coordinates to level 15 H3 indices.

The Dual Graph

Using the process above, we convert each historical trajectory into a sequence of H3 tokens. The next step is to convert each trajectory to a sequence of token triplets. Three values in a sequence represent two consecutive edges of the prediction graph, and we want to know the frequencies of these, as they will be the core data for both the prediction and the probability assessment. Figure 3 below depicts this process visually.

Figure 3 — The list of geospatial tokens on the left is expanded to another list of triplets, representing a dual vision of the implicit graph. Each token is a node on the geospatial graph, and its sequence represents the edges. The transformed list considers each edge a node in the dual graph, and the middle token is the new edge, as shown in the right column. (Image source: Author)

The transformation above computes the dual of the road graph, reversing the roles of the original nodes and edges.

We can now start to answer the proposed questions.

Should you be going that way?

We need to know the vehicle trajectory up to a given moment to answer this question. We map-match and tokenize the trajectory using the same process as above and then compute each trajectory triplet frequency using the known historic frequencies. The final result is the product of all individual frequencies. If the input trajectory has an unknown triplet, its frequency will be zero as the final path probability.

A triplet probability is the ratio of counts of a specific sequence (A, B, C) to the count of all (A, B, *) triplets, as depicted in Figure 4 below.

Figure 4 — The triplet probability is the ratio of its frequency to the frequency of all triplets with the same two initial tokens. (Image source: Author)

The trip probability is just the product of individual trip triplets, as depicted in Figure 5 below.

Figure 5 — The trip probability is the simple product of all matched triplets. (Image source: Author)

Where are you going?

We use the same principles to answer this question but start with the last known triplet only. We can predict the k most likely successors using this triplet as input by enumerating all triplets that have as their first two tokens the last two of the input. Figure 6 below illustrates the process for triplet sequence generation and evaluation.

Figure 6 — In this fictitious case, the next most likely triplet is the one with the highest observed frequency: (B, C, D). (Image source: Author)

We can extract the top k successor triplets and repeat the process to predict the most likely trip.

We are ready to discuss the implementation details, starting with map-matching and some associated concepts. Next, we will see how to use the Valhalla toolset from Python, extract the matched paths and generate the token sequences. The data preprocessing step will be over once we store the result in the database.

Finally, I illustrate a simple user interface using Streamlit that calculates the probability of any hand-drawn trajectory and then projects it into the future.


Map-matching converts GPS coordinates sampled from a moving object’s path into an existing road graph. A road graph is a discrete model of the underlying physical road network consisting of nodes and connecting edges. Each node corresponds to a known geospatial location along the road, encoded as a latitude, longitude, and altitude tuple. Each directed edge connects adjacent nodes following the underlying road and contains many properties such as the heading, maximum speed, road type, and more. Figure 7 below illustrates the concept with a straightforward example.

Figure 7 — The picture above shows a tiny digital road network highlighting an intersection. Each red dot represents a known geospatial location along the existing road. The blue lines represent the connecting edges between the nodes. Note that these edges are usually directed and might also be multiple. (Image source: Author)

When successful, the map-matching process produces relevant and valuable information on the sampled trajectory. On the one hand, the process projects the sampled GPS points to locations along the most likely road graph edges. The map-matching process “corrects” the observed spots by squarely placing them over the inferred road graph edges. On the other hand, the method also reconstructs the sequence of graph nodes by providing the most likely path through the road graph corresponding to the sampled GPS locations. Note that, as previously explained, these outputs are different. The first output contains coordinates along the edges of the most likely path, while the second output consists of the reconstructed sequence of graph nodes. Figure 8 below illustrates the process.

Figure 8 — The diagram above illustrates the map-matching process, where the green dots represent the observed GPS coordinates, and the orange diamonds are the projected locations along the known edges. Note that, for the simplified example above, we can only safely reconstruct the path between nodes 2 and 3. This predicament is not as dire as it looks because, in actual maps, trajectories match many more edges than just one, so the information loss is minimal. (Image source: Author)

A byproduct of the map-matching process is the standardization of the input locations using a shared road network representation, especially when considering the second output type: the most likely sequence of nodes. When converting sampled GPS trajectories to a series of nodes, we make them comparable by reducing the inferred path to a series of node identifiers. We can think of these node sequences as phrases of a known language, where each inferred node identifier is a word, and their arrangement conveys behavioral information.

This is the fifth article where I explore the Extended Vehicle Energy Dataset¹ (EVED) [1]. This dataset is an enhancement and review of prior work and provides the map-matched versions of the original GPS-sampled locations (the orange diamonds in Figure 8 above).

Unfortunately, the EVED only contains the projected GPS locations and misses the reconstructed road network node sequences. In my previous two articles, I addressed the issue of rebuilding the road segment sequences from the transformed GPS locations without map-matching. I found the result somewhat disappointing, as I expected less than the observed 16% of defective reconstructions. You can follow this discussion from the articles below.

Now I am looking at the source map-matching tool to see how far it can go in correcting the defective reconstructions. So let’s put Valhalla through its paces. Below are the steps, references, and code I used to run Valhalla on a Docker container.

Valhalla Setup

Here I closely follow the instructions provided by Sandeep Pandey [2] on his blog.

First, make sure that you have Docker installed on your machine. To install the Docker engine, please follow the online instructions. If you work on a Mac, a great alternative is Colima.

Once installed, you must pull a Valhalla image from GitHub by issuing the following commands at your command line, as the shell code in Figure 9 below depicts.

Figure 9 — Pulling Valhalla’s docker image from the command line. (Image source: Author)

While executing the above commands, you may have to enter your GitHub credentials. Also, ensure you have cloned this article’s GitHub repository, as some files and folder structures refer to it.

Once done, you should open a new terminal window and issue the following command to start the Valhalla API server (MacOS, Linux, WSL):

Figure 10 — The above command runs the pulled Valhalla image in a Docker container. During first-time execution, the command also downloads and prepares the latest Geofabrik Michigan data file before starting. (Image source: Author)

The command line above explicitly states which OSM file to download from the Geofabrik service, the latest Michigan file. This specification means that when executed the first time, the server will download and process the file and generate an optimized database. In subsequent calls, the server omits these steps. When needed, delete everything under the target directory to refresh the downloaded data and spin up Docker again.

We can now call the Valhalla API with a specialized client.

Enter PyValhalla

This spin-off project simply offers packaged Python bindings to the fantastic Valhalla project.

Using the PyValhalla Python package is quite simple. We start with a neat install procedure using the following command line.

Figure 11 — You can install PyValhalla using PIP. (Image source: Author)

In your Python code, you must import the required references, instantiate a configuration from the processed GeoFabrik files and finally create an Actor object, your gateway to the Valhalla API.

Figure 12 — The code above shows how easy it is to set up PyValhalla on a Python application or notebook. (Image source: Author)

Before we call the Meili map-matching service, we must get the trajectory GPS locations using the function listed below in Figure 13.

Figure 13 — The function above loads the unique positions of a vehicle’s trajectory, returning a Pandas DataFrame with latitude, longitude, and timestamp. (Image source: Author)

We can now set up the parameter dictionary to pass into the PyValhalla call to trace the route. Please refer to the Valhalla documentation for more details on these parameters. The function below calls the map-matching feature in Valhalla (Meili) and is included in the data preparation script. It illustrates how to determine the inferred route from a Pandas data frame containing the observed GPS locations encoded as latitude, longitude, and time tuples.

Figure 14 — The function above accepts a PyValhalla Actor object and a Pandas DataFrame containing the source path and returns a map-matched string-encoded polyline. This string is later decoded into a list of geospatial locations corresponding to the digital map network nodes, except for the extremities, which are edge-projected. (Image source: Author)

The above function returns the matched path as a string-encoded polyline. As illustrated in the data preparation code below, we can easily decode the returned string using a PyValhalla library call. Note that this function returns a polyline whose first and last locations are projected to edges, not graph nodes. You will see these extremities removed by code later in the article.

Let us now look at the data preparation phase, where we convert all the trajectories in the EVED database into a set of map edge sequences, from where we can derive pattern frequencies.

Data preparation aims at converting the noisy GPS-acquired trajectories into sequences of geospatial tokens corresponding to known map locations. The main code iterates through the existing trips, processing one at a time.

In this article, I use an SQLite database to store all the data processing results. We start by filling the matched trajectory path. You can follow the description using the code in Figure 15 below.

Figure 15 — The code above contains the preprocessing data loop. This loop iterates through the known trajectories, computes their map-matched paths (if any), tokenizes the nodes, and expands them into triplets. The code stores all intermediary and final results in the database. (Image source: Author)

For each trajectory, we instantiate an object of the Actor type (line 9). This is an unstated requirement, as each call to the map-matching service requires a new instance. Next, we load the trajectory points (line 13) acquired by the vehicles’ GPS receivers with the added noise, as stated in the original VED article. In line 14, we make the map-matching call to Valhalla, retrieve the string-encoded matched path, and save it to the database. Next, we decode the string into a list of geospatial coordinates, remove the extremities (line 17) and then convert them to a list of H3 indices computed at level 15 (line 19). On line 23, we save the converted H3 indices and the original coordinates to the database for later reverse mapping. Finally, on lines 25 to 27, we generate a sequence of 3-tuples based on the H3 indices list and save them for later inference calculations.

Let’s go through each of these steps and explain them in detail.

Trajectory Loading

We have seen how to load each trajectory from the database (see Figure 13). A trajectory is a time-ordered sequence of sampled GPS locations encoded as a latitude and longitude pair. Note that we are not using the matched versions of these locations as provided by the EVED data. Here, we use the noisy and original coordinates as they existed in the initial VED database.

Map Matching

The code that calls the map-matching service is already presented in Figure 14 above. Its central issue is the configuration settings; other than that; it is a pretty straightforward call. Saving the resulting encoded string to the database is also simple.

Figure 16 — The code above saves the encoded polyline string to the new database. (Image source: Author)

On line 17 of the main loop (Figure 15), we decode the geometry string into a list of latitude and longitude tuples. Note that this is where we strip out the initial and final locations, as they are not projected to nodes. Next, we convert this list to its corresponding H3 token list on line 19. We use the maximum detail level to try and avoid overlaps and ensure a one-to-one relationship between H3 tokens and map graph nodes. We insert the tokens in the database in the following two lines. First, we save the whole token list associating it to the trajectory.

Figure 17 — The function above inserts the trajectory H3 token list in the database. (Image source: Author)

Next, we insert the mapping of node coordinates to H3 tokens to enable drawing polylines from a given list of tokens. This feature will be helpful later on when inferring future trip directions.

Figure 18 — We insert a mapping between H3 tokens and node coordinates to enable the reconstruction of a trajectory from given inferred tokens. (Image source: Author)

We can now generate and save the corresponding token triples. The function below uses the newly generated list of H3 tokens and expands it to another list of triples, as detailed in Figure 3 above. The expansion code is depicted in Figure 19 below.

Figure 19 — The code above converts a list of H3 tokens into a list of the corresponding triples. (Image source: Author)

After triplet expansion, we can finally save the final product to the database, as shown by the code in Figure 20 below. Through clever querying of this table, we will infer current trip probabilities and future most-likely trajectories.

Figure 20 — The function above saves the H3 triples to the database. This is the final step of the data preparation phase. We can now move on to exploring the information we collected. (Image source: Author)

We are now done with one cycle of the data preparation loop. Once the outer loop is completed, we have a new database with all the trajectories converted to token sequences that we can explore at will.

You can find the whole data preparation code in the GitHub repository.

We now turn to the problem of estimating existing trip probabilities and predicting future directions. Let’s start by defining what I mean by “existing trip probabilities.”

Trip Probabilities

We start with an arbitrary path projected into the road network nodes through map-matching. Thus, we have a sequence of nodes from the map and want to assess how probable that sequence is, using as a frequency reference the known trip database. We use the formula in Figure 5 above. In a nutshell, we compute the product of the probabilities of all individual token triplets.

To illustrate this feature, I implemented a simple Streamlit application that allows the user to draw an arbitrary trip over the covered Ann Arbor area and immediately compute its probability.

Once the user draws points on the map representing the trip or the hypothetical GPS samples, the code map matches them to retrieve the underlying H3 tokens. From then on, it’s a simple matter of computing the individual triplet frequencies and multiplying them to compute the total probability. The function in Figure 21 below computes the probability of an arbitrary trip.

Figure 21 — The function above computes an arbitrary path probability from the triplet frequency database. (Image source: Author)

The code gets support from another function that retrieves the successors of any existing pair of H3 tokens. The function listed below in Figure 22 queries the frequency database and returns a Python Counter object with the counts of all successors of the input token pair. When the query finds no successors, the function returns the None constant. Note how the function uses a cache to improve database access performance (code not listed here).

Figure 22 — The function above queries the frequency database for the known successors of any pair of H3 tokens and returns a Counter object with the counts of all successors. (Image source: Author)

I designed both functions such that the computed probability is zero when no known successors exist for any given node.

Let us look at how we can predict a trajectory’s most probable future path.

Predicting Directions

We only need the last two tokens from a given running trip to predict its most likely future directions. The idea involves expanding all the successors of that token pair and selecting the most frequent ones. The code below shows the function as the entry point to the directions prediction service.

Figure 23 — The function above populates a FeatureGroup object from Folium with the predicted paths of the existing user-provided trip. (Image source: Author)

The above function starts by retrieving the user-drawn trajectory as a list of map-matched H3 tokens and extracting the last pair. We call this token pair the seed and will expand it further in the code. At line 9, we call the seed-expansion function that returns a list of polylines corresponding to the input expansion criteria: the maximum branching per iteration and the total number of iterations.

Let us see how the seed expansion function works by following the code listed below in Figure 24.

Figure 24 — The seed expansion function uses the PredictedPath class to manage each iteration. Please see below for more details on this class. (Image source: Author)

By calling a path expansion function that generates the best successor paths, the seed expansion function iteratively expands paths, starting with the initial one. Path expansion operates by picking a path and generating the most probable expansions, as shown below in Figure 25.

Figure 25 — The path expansion function above iterates the most frequent successors to the current path. It creates a new path for each of the most frequent successors using a specialized function (see below). (Image source: Author)

The code generates new paths by appending the successor nodes to the source path, as shown in Figure 26 below.

Figure 26 — To generate a “child” path, we only need to append the successor node to an existing path, as shown below. Note that the code creates a copy of the original path before appending the new node. (Image source: Author)

The code implements predicted paths using a specialized class, as shown in Figure 27.

Figure 27 — The class above implements a predicted path with probability sorting support, creation from a seed token pair, and map polyline generation. (Image source: Author)

We can now see the resulting Streamlit application in Figure 28 below.

Source link

This post originally appeared on TechToday.