Continuous Map Matching to Paths
Under Travel Time Constraints

A naive attempt for map matching is to assign each measurement $M_i$ to the closest node in $G$ or to the closest point on an edge in $G$ [30]. However, this method does not produce sensible matching paths reliably as it neglects the necessary continuity of the movement in the graph [22]. A more sophisticated and wide-spread method for map matching is the following one [23, 11, 19], which we refer to as DAG algorithm:

• $\mathrm{■}$

  For each $M_i$, select a finite set of candidate locations $L_i$ that exist in $G$.

• $\mathrm{■}$

  Construct a layered DAG where layer $i$ contains a node for each element in $L_i$ and where there are edges $(a,b)$ for each $a\in L_i,b\in L_{i+1}$. Edge costs correspond to the minimum travel time between the respective nodes in the graph.

• $\mathrm{■}$

  Prune edges that do not obey given travel time constraints between consecutive locations.

• $\mathrm{■}$

  Find the shortest path $\pi$ from a node in layer $1$ to a node in layer $k$ and return it.

The advantage of this approach is its simplicity. However, the drawback of restricting candidate locations to a finite candidate location set in the underlying graph is that the quality of the returned path and the efficiency of the method crucially depend on the density and distribution of the candidate locations. Most commonly, $L_i$ is either the set of the $c$-closest nodes to $p_i$ in $V$ for some $c\in \mathbb{N}$ or the set of nodes within a disk $D_i$ of radius $r$ centered at $p_i$ (where $r$ models the location measurement uncertainty).

Potential issues that arise when using these methods in the DAG approach are illustrated in Figure 2.

If the candidate set is formed by the $c$ nodes that are closest to the measured location, selecting a sensible value for $c$ is challenging. If there are densely sampled paths in the graph (which often is the case for curvy routes in real networks), then even with a large value of $c$ all candidate locations might be on one road; and good alternatives, which are slightly further away, could be completely missed. Setting $c$ to a large value might result in a mapping to locations that are far away from the measured position. It also increases the number of nodes in the constructed DAG and with that the running time of the approach. When using disks $D_i$ of radius $r$ to determine the candidate locations as $D_i\cap V$, feasible mapping locations can be missed because $D_i$ might only intersect an edge but contains none of its endpoints. Furthermore, the running time for computing the DAG might also become huge if many nodes are contained in the disks, as shortest paths have to be computed between all pairs of candidate locations for consecutive layers. In general, using only a fixed set of candidate locations becomes particularly problematic when travel time constraints are taken into account. It might be the case that there exists a path ${\pi }_{i,i+1}$ in $G$ that connects a location close to $p_i$ to a location close to $p_{i+1}$ with a travel time of at most $t_{i+1}-t_i$. However, if the respective start and end locations are on edges and not on nodes, it might be that only a superpath of ${\pi }_{i,i+1}$ is represented in the DAG with a travel time that exceeds $t_{i+1}-t_i$. Thus, the algorithm cannot identify a feasible matching path even though there is one. Clearly, all of the discussed issues can already arise if $G$ is a path graph. And even for this special case, the described DAG algorithm is the prevailing standard [13].

In this paper, we propose a continuous map matching approach that does not require to fix a finite set of candidate locations for each measurement a priori. Instead, we allow $M_i$ to be mapped to any location in $D_i$ that exists in $G$. This prohibits the usage of the DAG algorithm, since there can now be an infinite number of candidate locations in each layer. We prove that our novel algorithm always identifies a feasible matching path (if one exists) and that it does so even faster in practice than the DAG algorithm which might fail to produce a correct result due to its restricted candidate location set.

Many approaches for map matching have been proposed in the literature, using different input models, objective functions, and assumptions about the underlying embedded graph.

In [23, 9, 19], the DAG algorithm is used with a more refined candidate selection strategy. Here, first the $c$ closest edges to the measured position or the edges within a certain radius are selected and then a suitable point per edge is chosen. However, the potential issues described above for finite candidate location sets also apply to these strategies. In [11], engineering methods are described to accelerate the baseline algorithm with candidate locations being selected as nodes in disks of radius $r$. In particular, methods for fast one-to-many shortest path computations are shown to reduce the running time significantly. We will show that for path graphs, even faster solutions are possible.

A completely different approach for map matching is to interpret $M$ as a polyline by linear interpolation between the measured positions and to find the path in $G$ with minimum Fréchet distance to that polyline [16, 4]. With this objective, there is no need to specify any restricted candidate location sets and one gets a continuous mapping sequence. However, the linear interpolation between position measurements might not be realistic, especially if measurements are sparse. Moreover, the computation of the exact Fréchet distance is computationally expensive. Thus, for fast query times, one has to resort to approximations [6]. Finally, it is unclear how to incorporate travel time constraints into this framework. Finally, as the Fréchet distance is a bottleneck measure, there is usually a large set of matching paths with the same objective function value. Thus, oftentimes there needs to be a second map matching step with another objective function or a combined objective to select the final matching path [29]. Other polyline distance measures, as the weak Fréchet distance or the vertex-monotone Fréchet distance [3, 29] as well as dynamic time warping [28], have also been applied in the map matching context. But they all suffer from similar issues and, in particular, the respective methods are not designed to identify a feasible matching path with respect to time constraints.

There is also a plethora of HMM or machine learning based methods for map matching [22, 24, 18, 20, 31]. However, those also do not come with quality guarantees and they either require a fine tuned model or high quality training data to work well in practice.

In [1, 27], map matching of bus stop sequences to a transport network are discussed. They also leverage the DAG algorithm but incorporate additional requirements as e.g turn restrictions or consistency among different bus lines. However, travel time constraints are also not taken into account.

The following algorithm identifies for each $M_i$ its possible mapping locations that are part of some feasible matching path. The idea is to start with all intervals in $D_i\cap P$ and to then exclude points that do not obey the travel time constraints with respect to $M_{i-1}$ or $M_{i+1}$.

1. 1.

   For each $i\in [k]$, compute the intervals $I_i$ of $P$ inside the disk $D_i$ and sort them increasingly.

2. 2.

   Forward sweep: For $i=2,\mathrm{\dots },k$, consider the intervals in $I_i$. For each interval $[a,b]$, find the latest endpoint $b^{\prime }$ of an interval in $I_{i-1}$ which is prior to $a$. If no such $b^{\prime }$ exists, delete the interval. Set $b=\min (b,b^{\prime }+t_i-t_{i-1})$. If $[a,b]=\mathrm{\varnothing }$, delete the interval.

3. 3.

   Backward sweep: For $i=k-1,\mathrm{\dots },1$, consider the intervals in $I_i$. For each interval $[a,b]$, find the earliest starting point $a^{\prime }$ of an interval in $I_{i+1}$ which starts after $b$. If no such $a^{\prime }$ exists, delete the interval. Set $a=\max (a,a^{\prime }-t_{i+1}+t_i)$. If $[a,b]=\mathrm{\varnothing }$, delete the interval.

The three steps are illustrated in Figure 3. The correctness of the algorithm is shown below.

If we get rid of the assumption of non-overlapping disks, only intervals within one disk remain non-overlapping while intervals in $I_i$ and $I_{i+1}$ might now indeed share points. This impacts the interval update procedure during the sweeps. In the forward sweep, for an interval $[a,b]\in I_i$ we can no longer only consider the interval in $I_{i-1}$ with the latest endpoint prior to $a$. Instead, we still find this interval via binary search but then proceed to check the successive intervals in the sorted list for intersection with $[a,b]$ until we reach the end of the list or the first interval that starts after $b$. We then subdivide $[a,b]$ at the starting points of its overlapping intervals. For each subinterval $[a_s,b_s]$ created in this way that starts with an overlap with some interval $[a^{\prime },b^{\prime }]$, we set $b_s$ to $\min (b_s,b^{\prime }+t_i-t_{i-1})$. Figure 4 illustrates the process. Only the first subinterval might not have an overlapping portion with an interval from $I_{i-1}$. In this case, we can simply handle it with the procedure described above for non-overlapping intervals. If intervals in $I_i$ overlap after this process, they can be merged into one. The backward sweep works analogue.

While the update procedure for non-overlapping disks can only decrease the total number of intervals, the one that takes care of overlaps might also increase it. However, as the starting points of the new intervals are inherited from starting points of intervals in $I_{i-1}$, there can be at most $kn$ interval starting points per disk during the course of the algorithm. While multiple intervals from $I_{i-1}$ might be considered for a single interval $[a,b]$ from $I_i$, the set of intervals from $I_i$ strictly contained $[a,b]$ it is only considered for $[a,b]$ and not for any other interval from $I_i$ as those are all pairwise non-overlapping. Therefore, apart from the binary searches, the number of interval access operations is bounded by $|I_{i-1}|+|I_i|\in 𝒪(kn)$. Thus, the total running time per sweep is in $𝒪(k^{2}n\log nk)$.

In [17], it was discussed how to preprocess a given set of $n$ line segments in ${ℝ}^2$ such that for a query circle $C$ of radius $r$ the subset of segments that intersect $C$ can be reported efficiently. There are two types of segment-circle intersections (both visible in Figure 3):

1. (i)

   The segment intersects $C$ once. Here, one segment endpoint is inside $C$ and the other one is outside of $C$.

2. (ii)

   The segment intersects $C$ twice. Here, both segment endpoints are outside of $C$ but the segment point closest to the center of $C$ is inside $C$.

It is suggested in [17] to build a separate data structure for each type. We review and analyze these data structures in more detail below. For the first intersection type, we propose a significantly faster and simpler data structure for our use case in which the segments are known to form a path. For the second intersection type, we show that the method proposed in [17] only works for segments that are longer than the circle diameter. We then describe how to handle the missing case without increasing the preprocessing or query time. This result applies to arbitrary segment sets and might be of independent interest.

Subsequently, we describe in more detail how to leverage the resulting data structures for improving the COMMA running time.

The goal is now to leverage our new segment-circle intersection data structure for accelerating the interval computation phase of COMMA. For each disk $D_i$, we query the data structure with the disk boundary circle $C_i$ to retrieve the set of segments that intersect $C_i$. Then, we compute the explicit intersection points for each segment. To get the final intervals, we also need the travel time from the start point of $P$ to the respective intersection point. To have access to this value in constant time, we also compute a prefix sum array for the edge costs along $P$. Querying this data structure with an edge segment index and a point on that segment only demands to interpolate the cost assigned to that segment in the array up to the given point. The resulting intersection point timestamps are sorted to construct the respective intervals. Using this procedure, we can compute the set of intervals for each $D_i$ without having to traverse any edges that are fully contained in $D_i$. Note that the DAG based algorithm always has to consider the endpoints of all of these segments to identify the candidate location sets. Thus, the running time of COMMA is significantly reduced while the $𝒪(k{n}^2)$ running time of the baseline remains unchanged in our realistic model.

In the evaluation, we use two types of data sets: generated paths and measurement sequences as well as GTFS data. The former is created as follows. To get a path $P$ of length $n$, we add vectors with coordinates that are uniformly sampled at random. The time values for the points are accumulated in the same way. The corresponding measurement sequence $M$ is generated as follows. We first randomly sample $k$ values in $[0,1]$ and sort them ascendingly. These values are used to find positions along the path with an associated time value. Perturbing said points gives us our generated input. Figure 6 shows a small example and Figure 7 a large one.

For real data, we used GTFS bus timetables. We extract the coordinates present in shapes.txt of a bus line, and its corresponding stops in stops.txt (defining a trip of the vehicle). The time values for the trip are present in stop_times.txt. Using said times together with the shape_dist_traveled data in stop_times.txt, we can induce valid time values for the shape of the bus line. An example is shown in Figure 1.

For the DAG algorithm as well as for COMMA, we need to retrieve the intervals of $P$ inside the disks $D_1,\mathrm{\dots },D_k$, and for the DAG algorithm additionally the contained segment endpoints as those serve as candidate locations. However, it might happen that the DAG algorithm does not find a feasible solution when restricted to segment endpoints inside the disks. Thus, we also implement a sampling variant, which additionally places possible locations along the edges inside $D_i$ equidistantly (see again Figure 6). This allows us to study what sampling rates are needed to find feasible mapping sequences.

Naively, the intervals and segment endpoints can be retrieved in time $𝒪(kn)$ for both algorithms. But as discussed in our theoretical analysis, interval computation might easily dominate the running time of COMMA and thus using a suitable segment-circle intersection data structure is advised. However, the data structures we used in the theoretical analysis are quite involved, using transformations to three dimensions as well as several nested data structures (especially for intersections of type (ii)) [17]. To the best of our knowledge, no library implementation is available. Engineering these data structures would certainly be interesting, but a naive implementation is expected to be rather slow. Thus, we decided to use CGAL’s AABB Tree in both algorithms, DAG and COMMA. It is a hierarchical data structure where each node represents an axis-aligned bounding box (AABB) that encloses the set of geometric primitives that are represented by the respective leaf nodes. In our case, the leaf nodes store the path segments. Thus, the search algorithm is similar to the one we described in Section 4 only that here the oracle that decides whether to explore a branch of the tree might be false positive in case only the bounding box of the segment and the disk intersect but not the elements themselves. But it allows us to handle both possible types of segment-circle intersection in one data structure as well as segment endpoint extraction for the DAG baseline.

In Figure 8, we evaluate the running time for interval and segment extraction on generated paths with varying values of $n$ and $k$. We observe that the running time of the naive approach increases linear with both $n$ and $k$ and is in the order of minutes already for small parameter values, while the AABB Tree provides results very quickly for all choices of $n$ and $k$.

Next, we compare COMMA to the DAG baseline on real and generated input data. In Figure 9, a comparison of the algorithms on GTFS bus stop mapping instances is provided. For the baseline, the running time difference between using or not using the AABB tree is rather small, as most of the effort is required to construct the DAG. We remark that the unsuccessful runs of the baseline due to a too low sampling rate to identify a feasible match are not even included in the timings. For COMMA, the running time is already better on the larger instances without the AABB tree, and with the help of the AABB tree it achieves speed-ups of up to four orders of magnitude. This allows for much faster GTFS data processing especially for data sets with thousands of routes.

For further scalability studies, we use the above described generator. In Figure 10, left, we observe that the DAG baseline without additional samples fails to produce a feasible solution on the example instance. With a higher sampling rate and lower sample distance, the failure rate decreases but at the same time the running time increases significantly, as the DAG gets much larger. COMMA produces a feasible mapping path and outperforms the DAG algorithm by four orders of magnitude. A similar behavior was observed across all generated instances. The main reason for this huge speed-up is illustrated in Figure 10, right. The number of intervals that COMMA has to consider is significantly smaller than the number of locations that the DAG algorithm needs to process except for very small radii. As the DAG algorithm additionally scales worse with the number of locations to consider than COMMA, the running time gap gets more pronounced the more measurements need to processed and the larger the disks are.

Executing the DAG baseline on the instance shown in Figure 7 with $n=10000$, $k=1000$ does not yield a feasible result even with a sampling distance as small as $0.004$, although the DAG already contains $3.5$ million candidate nodes and $796$ million edges (that respect the time bound for consecutive locations). COMMA, however, scales very well and can easily identify feasible mapping paths for such instances. Indeed, for the generated instances depicted in Figure 8 with up to $n=100000$ segments and up to $k=10000$ measurements, COMMA produces the solution in less than a second. Thus, COMMA is also suitable to be used as an efficient subroutine in map matching algorithms for general graphs that select a set of candidate paths and then check whether they constitute feasible matching paths for GPS sequences with potentially thousands of measurements.

To further illustrate the scalability of COMMA, Table 1 shows results for very large $n$ and $k$ and also provides a breakdown of the running time. Even for millions of path segments and measurements, the running time stays within a few minutes. We see that the interval computation dominates the overall running time even when using the AABB tree, while the sweeps and the path extraction are very efficient.

Figure 11 provides a detailed analysis of the number of segment bounding boxes checked by the AABB tree traversal and the actual number of intervals that were identified for varying $n$ and $k$.

Clearly, there is quite a significant gap between those two values. Thus, it might indeed be worthwhile to implement and engineer a segment-circle intersection data structure with better theoretical guarantees in future work to cater for huge instances. However, for typical instances derived from GTFS data or from GPS measurements with up to a few thousand segments and measurements, leveraging the AABB tree is sufficient to achieve very low query times.