Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Map matched route does not perfectly align with road network #5960

Closed
ghost opened this issue Feb 9, 2021 · 20 comments
Closed

Map matched route does not perfectly align with road network #5960

ghost opened this issue Feb 9, 2021 · 20 comments

Comments

@ghost
Copy link

ghost commented Feb 9, 2021

I've map matched a series of GPS points to the OSM road network (downloaded from geofabrik) using the following options:

overview=full
geometries=polyline6
tidy=true

Screen Shot 2021-02-09 at 13 55 04

According to the documentation, these options should yield the highest precision in terms of OSRM route geometry. However, when zooming in you can see that the map matched route (red) does not exactly align with OSM's road network (green). I measured gaps ranging from 1 millimeter to 50 millimeters.

Screen Shot 2021-02-09 at 13 55 18

My goal is to provide insight in where but more importantly how frequently the vehicle has been per road segment in the road network. If the map matched route would have perfectly aligned with the road network, I could have counted the number of overlapping (route) lines per street segment. Is it possible to increase the precision of the map matched route so that it perfectly aligns with the OSRM road network? For example by using other options than mentioned above? Or does someone have a better workflow for what I'm trying to achieve?

@ghost ghost changed the title Map matched route does not perfectly algin with the road network Map matched route does not perfectly align with road network Feb 9, 2021
@danpat
Copy link
Member

danpat commented Feb 9, 2021

OSRM only stores 6 decimal places of fixed precision for coordinates (about 10cm precision at the equator). OSM (and the PBF you downloaded) uses a 7-decimal places fixed precision format (about 1cm at the equator). So there's already up to 9cm of error introduced by that truncation.

OSRM saves space by not bothering with that extra decimal place - it was deemed unnecessary for navigation purposes. OSM itself is hardly accurate to 1cm, well, anywhere.

What you should be doing for this comparison is using the annotations=nodes request parameter. This will return the OSM node ID for each point in the matched response, and you can then just count the pairs that are returned in the API response. This will be unambiguous, and not dependent on flakey decimal comparison.

@ghost
Copy link
Author

ghost commented Feb 10, 2021

Thanks, that sounds like a much better solution. Where exactly do I find the node ID's in the OSM road network? I expected them to find in the Points layer of the PBF file I downloaded from geofabrik.

@danpat
Copy link
Member

danpat commented Feb 10, 2021

What tool are you using to inspect the PBF? Every coordinate in the file has an associated node ID - in XML format, they look like this:

  <node id="298884269" lat="54.0901746" lon="12.2482632" user="SvenHRO" uid="46882" visible="true" version="1" changeset="676636" timestamp="2008-09-21T21:37:45Z"/>

It'll be up to your tool to decide whether to display the node ID or not.

@ghost
Copy link
Author

ghost commented Feb 10, 2021

I'm using QGIS. What tool do you recommend?

@danpat
Copy link
Member

danpat commented Feb 10, 2021

Not QGIS - it throws away the node IDs.

What exactly is your desired output of this process? You will be able to get all the pairs and geometry from the OSRM API response alone. This would allow you to tabulate and rank the most common segments, and display them visually, no joining to a PBF required.

@ghost
Copy link
Author

ghost commented Feb 10, 2021

My initial idea was to (1) extract all node id's from the OSRM API (2) add them to a pandas dataframe and count each occurence (3) joining the number of occurences/hits of the node id's with the corresponding line segments in QGIS. But apparently there is another way to tabulate and rank the most common segments and visualize them. What would you recommend?

@ghost
Copy link
Author

ghost commented Feb 10, 2021

My desired output is (1) an attribute table of the road network with an extra column that shows the number of times the vehicle drove through the street, (2) a map visualisation of the flow of traffic through the network, with the width of each line segment representing the amount of traffic (i.e. based on the frequency values in the attribute table). Something like this:
Screen Shot 2021-02-10 at 16 17 08

@danpat
Copy link
Member

danpat commented Feb 10, 2021

How are you defining the rows in your road network table? A vehicle can turn off any given road name at any intersection - I assume you'd want to know if some sections of a named road have more visits than others?

@danpat
Copy link
Member

danpat commented Feb 10, 2021

Here's a script that will match an array of traces, collate all the number of times each piece of geometry is visited, and output GeoJSON features with a count property that you can style on a map.

import requests
import json

traces = [[[ -122.29175806045532, 37.809512710254495 ],
           [ -122.29154348373412, 37.810233196634634 ],
           [ -122.29129672050475, 37.810996056903875 ],
           [ -122.29098558425902, 37.811801289748544 ],
           [ -122.29082465171813, 37.81219966485621 ],
           [ -122.2902774810791, 37.812106428321556 ],
           [ -122.28985905647279, 37.81201319166914 ] ]]


paircounts = {}

for trace in traces:
  s = ";".join(map(lambda y: ",".join(str(x) for x in y), trace))
  url = "https://router.project-osrm.org/match/v1/car/%s?overview=full&geometries=geojson" % s
  response = requests.get(url)

  distinct_pairs = set()
  for matching in response.json()["matchings"]:
    coordinates = list(map(lambda c: (c[0],c[1]), matching["geometry"]["coordinates"]))
    for pair in zip(coordinates[:-1],coordinates[1:]):
      distinct_pairs.add(pair)

  for pair in distinct_pairs:
    if pair in paircounts:
      paircounts[pair] += 1
    else:
      paircounts[pair] = 1

geojson = { "type": "FeatureCollection", "features": []}
for pair,count in paircounts.items():
  feature = { "type": "Feature", "geometry": { "coordinates": [ [pair[0][0],pair[0][1]],[pair[1][0],pair[1][1]]], "type": "LineString"}, "properties": { "count": count }}
  geojson["features"].append(feature)

print(json.dumps(geojson))

You can probably do a pass over this geojson to join directly-connected lines together into contiguous features to get nicer visualization.

@Rougree
Copy link

Rougree commented Feb 11, 2021

@danpat I'm a little bit a noob with OSRM and plotting a route in an OSRM map. I'm testing with the route service How can you visualize a route (response of route request) into a map of OSRM. I want to see what kind of route OSRM returns to some route requests. I'm doing requests in python btw.

@jcoupey
Copy link

jcoupey commented Feb 11, 2021

@Rougree you may want to start with osrm-frontend https://github.com/Project-OSRM/osrm-frontend/ (see docs over there on how to point it to your own OSRM instance).

@ghost
Copy link
Author

ghost commented Feb 15, 2021

@danpat thanks for sharing your script. This is probably the method I should use, since I would like to use a custom road network for the analysis in the end. Although I have been able to convert a custom road network from .shp to .osm using ogr2osm (and successfully map match the GPS points to this network), the original node id's from the shapefile appear to have lost after conversion. I was hoping ogr2osm would be able to transfer the node id's from the .shp to the .osm file, as the original shapefile contains a column "source" and a column "target" to indicate the node id's. Instead, new node id's have been chronologically numbered in the .osm file (see below), which makes it useless to count the node id pairs.

Screen Shot 2021-02-12 at 15 06 55

Btw, I tried to rerun your script with a different trace:

trace = [[-122.292423248291, 37.807779275984046],
 [-122.29196190834045, 37.80918636999258],
 [-122.29170441627501, 37.810038242190245],
 [-122.29131281375885, 37.81067396130711],
 [-122.29036331176758, 37.81044934117778],
 [-122.29056179523468, 37.809410993963944],
 [-122.29086220264433, 37.80893207787961],
 [-122.29172050952911, 37.80897022175241],
 [-122.29196727275848, 37.80940251760008],
 [-122.29151666164398, 37.81054257980499],
 [-122.29074418544768, 37.81069515185008],
 [-122.2902774810791, 37.81024591103703],
 [-122.29073882102965, 37.8088981721984],
 [-122.2934317588806, 37.8094703284838],
 [-122.29287922382353, 37.8110045330848],
 [-122.29044377803802, 37.8104959605061],
 [-122.29091584682465, 37.8089363160887],
 [-122.28980541229247, 37.80862692618813],
 [-122.29021847248076, 37.807380877026844]]

This trace visits some streets a couple of times. However, this isn't translate to the count, which is 1 for all geometries. Do you get the same result?

@danpat
Copy link
Member

danpat commented Feb 15, 2021

I cannot reproduce your results without your map, I only have access to OpenStreetmap.

I would suggest you check the results of your ogr2ogr conversion and check that everything is connected as you expect. A few things are required:

  1. Generated OSM ways are connected at intersections
  2. Proper OSM-like tags are placed on ways in the .osm file so that OSRM knows what road classes/routing rules can be applied.

If (1) and (2) haven't been done, then the network you're trying to match is basically junk. I'm not sure how much magic ogr2ogr will do out of the box, but I suspect at least (2) will need some attention. Nothing will match against the OSM file if OSRM doesn't understand how everything is connnected together.

Try opening up your .osm file in JOSM and see if things are connected as you'd expect.

@ghost
Copy link
Author

ghost commented Feb 16, 2021

Sorry I maybe was a bit unclear.

The custom shapefile that I converted to an OSM file (using ogr2osm, not ogr2ogr btw) works totally fine; I have been able to match a GPS trace to this .osm network without any problems. The only drawback was that the original node id's are lost in the OSM file (see screenshot in my previous reply). Selected nodes in JOSM also don't show an associated id (see below). What I tried to explain was that your initial solution (counting the node id's pairs that are returned from the annotations=nodes request parameter) won't work, knowing that the original node id's are lost. This only leaves, as I understood it, the other option: the script that matches an array of traces, collates the number of times each piece of geometry is visited, and outputs a GeoJSON features with a count property. To test it out, I reran your script with a different trace (one that visits some streets a couple of times) using OpenStreetMap, not my own map. However, I got unexpected results. Could you try to reproduce it?
Screen Shot 2021-02-16 at 11 45 15

@datendelphin
Copy link

There is no way that a node has no id. Ways reference nodes only by id and have no coordinates of their own. But the node looks like a "new" node in josm. Maybe it's just a negative id which is usually used for new elements in editors.

@ghost
Copy link
Author

ghost commented Feb 16, 2021

That's correct, I meant that the original node id's are lost. In the screenshot of the .osm file you can see that ogr2osm created new negative id's

@danpat
Copy link
Member

danpat commented Feb 16, 2021

The script I posted, if you read it carefully, is only working with the coordinates, not the node IDs.

The critical line for detecting repeated visits to the same pair is this one:

    if pair in paircounts:

My Python is a bit rusty - I'm assuming that two separately constructed tuples containing the same numeric values will cause the above to evaluate to True. You didn't paste the output - if you look closely at the output list, are there any duplicates?

@ghost
Copy link
Author

ghost commented Feb 17, 2021

I understand, but given the trace I posted above (see also map below), you would expect a couple more counts for some coordinate pairs (see output at bottom). My python is also not very developed, so I might miss the reason why this is happening. It seems it has to do with the fact that only unique coordinates pairs are counted

Screen Shot 2021-02-17 at 13 28 28

{((37.811073, -122.292686), (37.811004, -122.292384)): 1,
 ((37.810166, -122.293224), (37.81038, -122.293146)): 1,
 ((37.810259, -122.290337), (37.808912, -122.290803)): 1,
 ((37.808961, -122.291398), (37.809081, -122.291923)): 1,
 ((37.810639, -122.291392), (37.810688, -122.291374)): 1,
 ((37.811004, -122.292384), (37.810955, -122.292166)): 1,
 ((37.808692, -122.290221), (37.808632, -122.289963)): 1,
 ((37.808179, -122.29225), (37.808449, -122.292152)): 1,
 ((37.808449, -122.292152), (37.808739, -122.292047)): 1,
 ((37.808833, -122.29083), (37.808692, -122.290221)): 1,
 ((37.810518, -122.290245), (37.810402, -122.290287)): 1,
 ((37.810668, -122.290908), (37.810636, -122.290766)): 1,
 ((37.809372, -122.291822), (37.810525, -122.291431)): 1,
 ((37.811016, -122.29293), (37.81112, -122.292894)): 1,
 ((37.810738, -122.291223), (37.810683, -122.290979)): 1,
 ((37.808965, -122.291965), (37.809081, -122.291923)): 1,
 ((37.810765, -122.291345), (37.810738, -122.291223)): 1,
 ((37.81054, -122.293092), (37.811016, -122.29293)): 1,
 ((37.810518, -122.290245), (37.81043, -122.290277)): 1,
 ((37.810935, -122.292078), (37.810907, -122.291953)): 1,
 ((37.808739, -122.292047), (37.808965, -122.291965)): 1,
 ((37.808632, -122.289963), (37.808598, -122.28991)): 1,
 ((37.810603, -122.290623), (37.810558, -122.290422)): 1,
 ((37.810402, -122.290287), (37.810259, -122.290337)): 1,
 ((37.808833, -122.29083), (37.80886, -122.290934)): 1,
 ((37.809105, -122.292033), (37.809153, -122.292256)): 1,
 ((37.808598, -122.28991), (37.808598, -122.28991)): 1,
 ((37.80967, -122.293396), (37.810166, -122.293224)): 1,
 ((37.809182, -122.291887), (37.810017, -122.291603)): 1,
 ((37.809081, -122.291923), (37.809182, -122.291887)): 1,
 ((37.80886, -122.290934), (37.808961, -122.291398)): 1,
 ((37.80917, -122.291891), (37.809182, -122.291887)): 1,
 ((37.810955, -122.292166), (37.810935, -122.292078)): 1,
 ((37.810603, -122.290623), (37.810518, -122.290245)): 1,
 ((37.81112, -122.292894), (37.811073, -122.292686)): 1,
 ((37.810558, -122.290422), (37.810518, -122.290245)): 1,
 ((37.810907, -122.291953), (37.81088, -122.291837)): 1,
 ((37.808545, -122.289896), (37.807402, -122.29031)): 1,
 ((37.81088, -122.291837), (37.810858, -122.291739)): 1,
 ((37.810402, -122.290287), (37.808912, -122.290803)): 1,
 ((37.81038, -122.293146), (37.81054, -122.293092)): 1,
 ((37.810402, -122.290287), (37.809425, -122.290625)): 1,
 ((37.807773, -122.292398), (37.80781, -122.292385)): 1,
 ((37.810668, -122.290908), (37.810603, -122.290623)): 1,
 ((37.810636, -122.290766), (37.810603, -122.290623)): 1,
 ((37.809153, -122.292256), (37.809433, -122.293492)): 1,
 ((37.810525, -122.291431), (37.810639, -122.291392)): 1,
 ((37.808961, -122.291398), (37.80903, -122.291699)): 1,
 ((37.810858, -122.291739), (37.810838, -122.291651)): 1,
 ((37.810683, -122.290979), (37.810668, -122.290908)): 1,
 ((37.810838, -122.291651), (37.810765, -122.291345)): 1,
 ((37.810688, -122.291374), (37.810765, -122.291345)): 1,
 ((37.810639, -122.291392), (37.810765, -122.291345)): 1,
 ((37.808598, -122.28991), (37.808545, -122.289896)): 1,
 ((37.809081, -122.291923), (37.80917, -122.291891)): 1,
 ((37.809081, -122.291923), (37.809105, -122.292033)): 1,
 ((37.808912, -122.290803), (37.808833, -122.29083)): 1,
 ((37.80781, -122.292385), (37.808179, -122.29225)): 1,
 ((37.808919, -122.2908), (37.808833, -122.29083)): 1,
 ((37.810017, -122.291603), (37.810639, -122.291392)): 1,
 ((37.809433, -122.293492), (37.80948, -122.293473)): 1,
 ((37.809425, -122.290625), (37.808919, -122.2908)): 1,
 ((37.809182, -122.291887), (37.809372, -122.291822)): 1,
 ((37.81043, -122.290277), (37.810402, -122.290287)): 1,
 ((37.80948, -122.293473), (37.80967, -122.293396)): 1,
 ((37.80903, -122.291699), (37.809081, -122.291923)): 1}

@danpat
Copy link
Member

danpat commented Feb 17, 2021

If I plot all those coordinate pairs on a map, this is what it looks like:

Screen Shot 2021-02-17 at 6 25 34 AM

which is basically nonsense. Have you actually looked at the match geometry from your API call? Does the matched path look like the one you expect?

@ghost
Copy link
Author

ghost commented Feb 17, 2021

Interesting. The returned matched path looks as expected:
Screen Shot 2021-02-17 at 15 58 21

This issue was closed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants