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

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.

### The Problems

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.

### Solution

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.

Travel Time Estimation Using Quadkeys

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.

### Implementation

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

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.

Road Network Edge Matching With TrianglesMore on Road Network Matching

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.

https://medium.com/media/9a750dcacd2d7e78ca022b489c7dc0de/href

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):

https://medium.com/media/ec9ae3b2afb260052f2ae3ff5455b0e4/href

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.

https://medium.com/media/4df0d93389a435d4fc399a8ae673b770/href

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.

https://medium.com/media/fcf7b5d14d9280664123c0bf9dd2ec8f/href

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

https://medium.com/media/ce42fa5f00a0db6fc711a0eaecacfdf5/href

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.

https://medium.com/media/e188387243071d0156ded58aef80ba14/href

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

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.

https://medium.com/media/52f96118744423e30f414f865493187a/href

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.

https://medium.com/media/22333a655f5407f43af56410394d24fc/href

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.

https://medium.com/media/34bc93041cfc75e919d7841aa40dbf26/href

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.

https://medium.com/media/28fdb7bea1bf3c5740eaa357ae3d365a/href

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.

https://medium.com/media/d6f2ab4a125597a3ac5c10d6330cb9e1/href

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.

https://medium.com/media/8707a2ac9399790682104269b6cd05de/href

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.

### Probabilities and Predictions

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.

https://medium.com/media/aea28c175ef7fe7a53777a4971f22f57/href

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).

https://medium.com/media/1bb65aa7b34d8430baa10b1f7d86ae7b/href

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.

https://medium.com/media/df6e7390b0e2317e40b4eed53ac5a00f/href

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**.

https://medium.com/media/f5ca47415a9ec49299074f3282e54680/href

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**.

https://medium.com/media/0266e4c91cc7b5093f625d6da91e33d3/href

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

https://medium.com/media/6cf9aafd70ad5615896ab32d26ecb819/href

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

https://medium.com/media/519adc0488aebd49769027d291aef4fa/href

### The Application

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

**Figure 28** — The Streamlit application shows the two described features in action. The input trajectory is in blue, and you can draw it using the tool menu on the left-hand side of the map. Once drawn, the code computes its probability and displays it at the bottom. The three red trajectories are the three most-likely fifty-edge predictions for where the source trajectory may evolve. You get a popup with the computed probability by clicking on each trajectory. (Image source: Author)

### Conclusion

In this article, I presented a means to predict a vehicle’s future trajectory when driving through a digitally mapped road network. Using a historical trajectory database, this method assigns a probability to any trip and also predicts the most likely directions for the near future. Consequently, this method can detect unlikely or even novel trajectories never seen before.

We start with an extensive database of vehicle trajectories from the area of interest. Each path is a chronological sequence of geospatial coordinates (latitude and longitude) and other relevant properties such as speed. We typically collect these trajectories from onboard GPS receivers and centrally compile them into a database.

GPS samples are noisy due to unavoidable errors that occur during signal measurement. Natural and artificial obstacles, such as urban canyons, can significantly decrease the signal’s reception accuracy and increase geolocation errors. Fortunately, workable solutions solve this issue by probabilistically matching the GPS samples to a digital map. This is what map matching is all about.

By matching the noisy GPS samples to a known digital map, we not only correct the accuracy problem by projecting each instance to the map’s most likely road segment, but we also get a discrete sequence of existing map-defined locations that the vehicle most likely went through. This last result is instrumental for our trajectory prediction because it essentially converts a set of noisy GPS coordinates into a clean and well-known collection of points in the digital map. These digital markers are fixed and never change, and by projecting the GPS sample sequence into them, we get a string of well-known tokens that we can later use for prediction.

We calculate all probabilities using the known token sequence frequencies for arbitrary trajectories and their future evolution. The result is a couple of Python scripts, one for data preparation and another for data input and visualization using the Streamlit platform.

### Notes

The original paper authors licensed the dataset under the Apache 2.0 License (see the VED and EVED GitHub repositories). Note that this also applies to derivative work.

### References

**[1]** Zhang, S., Fatih, D., Abdulqadir, F., Schwarz, T., & Ma, X. (2022). Extended vehicle energy dataset (eVED): an enhanced large-scale dataset for deep learning on vehicle trip energy consumption. *arXiv*. https://doi.org/10.48550/arXiv.2203.08630

**[2]** Efficient and fast map matching with Valhalla — Sandeep Pandey (ikespand.github.io)

**[3]** Map Matching done right using Valhalla’s Meili | by Serge Zotov | Towards Data Science

João Paulo Figueira is a Data Scientist at tb.lx by Daimler Truck in Lisbon, Portugal.

Map-Matching for Trajectory Prediction was originally published in Towards Data Science on Medium, where people are continuing the conversation by highlighting and responding to this story.