Geospatial Analytics using Python and SQL

Let’s look at a ride-hailing example (e.g., Uber). When generating routes, Uber wants to minimize the distance and time to reach a destination:

There are two map matching approaches covered in this section:

Let’s use HMMs for map matching. This is the approach Uber takes for its map matching. Let’s visualize and explain the problem.

In essence, we have GPS points (O) and a given road network with road segments (R). We are trying to go from ambiguous GPS data on the left, to an inferred route, as depicted on the right. Thus, we are attempting to predict the probability of a GPS point belonging to a specific road segment at a certain point in time, in sequence of the driver’s route.

Let’s further set up our problem.

Our goal is to figure out for every GPS observation, which road segment, in theory, could have produced that GPS point.

To do this, candidate road segments are generated via kNN lookup (i.e., look up k closest road segments within search radius). Once candidate segments are identified, draw the projection of that GPS signal on road signal, which is a line perpendicular from signal to candidate road segment.

These projections are depicted as the transparent circles in the image on the right. Note that observations can have multiple segment candidates (e.g., O2 and O4 can be projected to R1/R2 or R2/R3, respectively).

The emission probability represents the likelihood of a vehicle present on certain road segments at certain moments.

For example, let’s take R3 and O5. We are trying to calculate the probability of seeing O5 if the car was traversing on R3. This is to be calculated by function of great circle distance between GPS signal and its projection for a given candidate road segment.

For one GPS point with m number of road segments nearby, there will be m emission probabilities representing the likelihood of this GPS trace on each road segment.

For GPS points G1, which have m nearby segments, and G2, which has n nearby segments, there are m * n  transition probabilities. These probabilities are in HMM, in which they pick up a sequence of states with maximum probabilities that are most likely to represent road segments on which the vehicle was moving.

The haversine distance between the GPS point and its snap point on the road segment. The emission probability indicates how likely the GPS would be observed if the vehicle is actually on the road segment.

We then calculating transition probability regarding one GPS point on a certain segment to another GPS point on a certain segment, calculated using the following formula. This formula is called DistanceDiff, and is defined as the following: the absolute value of the difference between the haversine distance of two GPS points and the routable distance between two snap points associated with the GPS points.

It is less likely that the vehicle was traversing through these two segments when emitting GPS positions if DistanceDiff is greater than others.

Now, we have to look at the transition probability of our HMM. In this example, the transition probability represents the likelihood of a vehicle moving from one road segment to another road segment over a certain duration.

An example transition probability to calculate is for R1, R2 – the probability the car traversed from R1 to R2.

An assumption we can make is that certain road transitions should be more likely than others. Any transition that provides a feasible route has a higher transition probability. A feasible route is the difference between the haversine distance and the shortest path distance between the two projections is closest to 0.

By calculating the transition probability, we can show the driver’s trip unfolding, following where they traversed. This concludes with the following trip recorded for our driver:

*In[54]:*

#Visualization from leuvenmapmatching import visualization as mmviz

#Examples from leuvenmapmatching.matcher.distance import DistanceMatcher from leuvenmapmatching.map.inmem import InMemMap from leuvenmapmatching.util.gpx import gpx_to_path

leuvenmapmatching.readthedocs.io/en/latest/usage/introduction.html#example-1-simple
A first, simple example. Some parameters are given to tune the algorithm. The max_dist and obs_noise are distances that indicate the maximal distance between observation and road segment and the expected noise in the measurements, respectively. The min_prob_norm prunes the lattice in that it drops paths that drop below 0.5 normalized probability. The probability is normalized to allow for easier reasoning about the probability of a path. It is computed as the exponential smoothed log probability components instead of the sum as would be the case for log likelihood.

*In[3]:*

#Define GPS Points to Map Match path = [(0.8, 0.7), (0.9, 0.7), (1.1, 1.0), (1.2, 1.5), (1.2, 1.6), (1.1, 2.0), (1.1, 2.3), (1.3, 2.9), (1.2, 3.1), (1.5, 3.2), (1.8, 3.5), (2.0, 3.7), (2.3, 3.5), (2.4, 3.2), (2.6, 3.1), (2.9, 3.1), (3.0, 3.2), (3.1, 3.8), (3.0, 4.0), (3.1, 4.3), (3.1, 4.6), (3.0, 4.9)]

matcher = DistanceMatcher(map_con, max_dist=2, obs_noise=1, min_prob_norm=0.5) states, _ = matcher.match(path) nodes = matcher.path_pred_onlynodes

print("States\n--") print(states) print("Nodes\n--") print(nodes) print("") matcher.print_lattice_stats()

*Out[3]:*

States

Stats lattice - nbr levels : 22 nbr lattice : 1002 avg lattice[level] : 45.54545454545455 min lattice[level] : 7 max lattice[level] : 97 avg obs distance : 0.15514927458475236 last logprob : -0.5464565099511667 last length : 22 last norm logprob : -0.024838932270507576

Now, we will visualize the above map, via leuvenmapmatching.readthedocs.io/en/latest/usage/visualisation.html

*In[7]:*

*Out[7]:* (None, None)

Let’s try visualizing an overlay of OpenStreetMaps:

*In[12]:*

*Out[12]:*

leuvenmapmatching.readthedocs.io/en/latest/usage/introduction.html#example-2-non-emitting-states
In case there are less observations that states (an assumption of HMMs), non-emittings states allow you to deal with this. States will be inserted that are not associated with any of the given observations if this improves the probability of the path.

It is possible to also associate a distribtion over the distance between observations and the non-emitting states (obs_noise_ne). This allows the algorithm to prefer nearby road segments. This value should be larger than obs_noise as it is mapped to the line between the previous and next observation, which does not necessarily run over the relevant segment. Setting this to infinity is the same as using pure non-emitting states that ignore observations completely.

*In[14]:*

#Define a Map to which Points are Mapped mapdb = InMemMap("mymap", graph={ "A": ((1, 0.00), ["B"]), "B": ((3, 0.00), ["A", "C"]), "C": ((4, 0.70), ["B", "D"]), "D": ((5, 1.00), ["C", "E"]), "E": ((6, 1.00), ["D", "F"]), "F": ((7, 0.70), ["E", "G"]), "G": ((8, 0.00), ["F", "H"]), "H": ((10, 0.0), ["G", "I"]), "I": ((10, 2.0), ["H"]) }, use_latlon=False)

matcher = DistanceMatcher(mapdb, max_dist_init=0.2, obs_noise=1, obs_noise_ne=10, non_emitting_states=True, only_edges=True) states, _ = matcher.match(path) nodes = matcher.path_pred_onlynodes

print("States\n--") print(states) print("Nodes\n--") print(nodes) print("") matcher.print_lattice_stats()

mmviz.plot_map(mapdb, matcher=matcher, show_labels=True, show_matching=True, filename="example_2_plot.png")

*Out[14]:*

States

Stats lattice - nbr levels : 3 nbr lattice : 40 avg lattice[level] : 13.333333333333334 min lattice[level] : 8 max lattice[level] : 16 avg obs distance : 0.26790850746762634 last logprob : -2.373678241605297 last length : 3 last norm logprob : -0.791226080535099 (None, None)

leuvenmapmatching.readthedocs.io/en/latest/usage/openstreetmap.html
The next section requires OSM data. OSM is open source maps and map data. Free!

*In[37]:*

Once we have a file containing the region we are interested in, we can select the roads we want to use to create a graph from. In this case we focus on `ways' with a `highway' tag. Those represent a variety of roads. For a more detailed filtering look at the possible values of the highway tag.

*In[48]:*

for entity in osmread.parse_file(str(xml_file)): if isinstance(entity, osmread.Way) and 'highway' in entity.tags: for node_a, node_b in zip(entity.nodes, entity.nodes[1:]): map_con.add_edge(node_a, node_b) # Some roads are one-way. We’ll add both directions. map_con.add_edge(node_b, node_a) if isinstance(entity, osmread.Node): map_con.add_node(entity.id, (entity.lat, entity.lon)) map_con.purge()

*Out[48]:*

*In[49]:*

#Remove Extra None from track track = [(lambda x: x[0:2])(x) for x in track]

matcher = DistanceMatcher(map_con, max_dist=100, max_dist_init=25, # meter min_prob_norm=0.001, non_emitting_length_factor=0.75, obs_noise=50, obs_noise_ne=75, # meter dist_noise=50, # meter non_emitting_states=True) states, lastidx = matcher.match(track)

*Out[49]:*

CPU times: user 2min 10s, sys: 1.61 s, total: 2min 12s Wall time: 2min 35s

*In[74]:*

*Out[74]:*

*In[51]:*

*Out[51]:*

Look at the blue lines and green arrows to see projections made by the model. The light blue line represents the predicted path…points matched to our map. For reference, here is what the points look like just mapped out:

leuvenmapmatching.readthedocs.io/en/latest/usage/latitudelongitude.html
The toolbox can deal with latitude-longitude coordinates directly. Map matching, however, requires a lot of repeated computations between points and latitude-longitude computations will be more expensive than Euclidean distances.

There are three different options how you can handle latitude-longitude coordinates: 1. Use Lat/Long Directly 2. Project Lat/Long to X/Y 3. Use Lat/Long as if they were X/Y

leuvenmapmatching.readthedocs.io/en/latest/usage/latitudelongitude.html#option-1-use-latitude-longitude-directly
Set the use_latlon flag in the Map to true.

For example to read in an OpenStreetMap file directly to a InMemMap object:

This is the same as above, so we will skip it

leuvenmapmatching.readthedocs.io/en/latest/usage/latitudelongitude.html#option-2-project-latitude-longitude-to-x-y

Latitude-Longitude coordinates can be transformed two a frame with two orthogonal axis.

*In[52]:*

route_latlon = [] # Add GPS locations route_xy = [map_con_xy.latlon2yx(latlon) for latlon in route_latlon]

*Out[52]:*

This can also be done directly using the pyproj toolbox. For example, using the Lambert Conformal projection to project the route GPS coordinates:

*In[66]:*

#Create Set of GPS Points Converted to GRS80 to_map = [] for x, y in zip(xs, ys): to_map.appendx, y

to_map

*Out[66]:*

leuvenmapmatching.readthedocs.io/en/latest/usage/latitudelongitude.html#option-3-use-latitude-longitude-as-if-they-are-x-y-points
A naive solution would be to use latitude-longitude coordinate pairs as if they are X-Y coordinates. For small distances, far away from the poles and not crossing the dateline, this option might work. But it is not advised.

For example, for long distances the error is quite large. In the image beneath, the blue line is the computation of the intersection using latitude-longitude while the red line is the intersection using Eucludean distances.