Here’s how you can find route intersections using GeoDesk and display them on a map.
First, you need to do a little setup:
-
Install Python and a Java (if your system doesn’t have them already)
-
Install the GeoDesk Python module: pip install geodesk
-
Download and install the GOL Tool
-
Download the OSM data for your region (list of extract services)
-
Turn the .osm.pbf file into a GOL): gol build drammen extract (Where extract is the name of the OSM data file)
Now we can write the actual script:
from geodesk import *
drammen = Features("drammen")
This gives us access to all the features in drammen.gol.
Next, we’ll define two subsets: cycling routes, and streets (any highway, except footpaths and minor ways like alleys and driveways).
routes = drammen("r[route=bicycle]")
streets = drammen("w[highway][highway != path,footway,pedestrian,service]")
We’re going to mark each point where a cycling route intersects with another street. First, we’ll create a Map object:
map = Map("cycle-routes")
We’ll process each route with a simple for-loop, then traverse each of its member ways and mark them in blue.
for route in routes:
print(f"{route} {route.name} ({route.ref}):")
count = 0
for member in route.members.ways:
map.add(member, color="blue", opacity=0.3, weight=20)
Inside this loop, we’ll traverse the nodes of each route segment and check their parent ways (excluding ways that are part of the route itself). If the node belongs to any street other than one on the route, we’ll mark it in red.
for node in member.nodes:
has_intersection = False
for way in streets.parents_of(node):
if not route in way.parents:
print(f" - Intersects with {way.highway} {way.name}")
has_intersection = True
if has_intersection:
map.add(node, color="red")
count += 1
For each route, we’ll print the total number of intersections:
print(f" {count} intersections.")
And finally, we’ll display our map in a browser:
map.show()
The output will look like this:
relation/6723229 Løkkeberget (S15):
- Intersects with tertiary Engene
- Intersects with tertiary Engene
- Intersects with unclassified Rådhusgata
- Intersects with unclassified Thornegata
- Intersects with residential Cappelens gate
...
28 intersections.
Here’s the full script:
cycleways.py
from geodesk import *
drammen = Features("c:\\geodesk\\tests\\drammen")
routes = drammen("r[route=bicycle]")
streets = drammen("w[highway][highway != path,footway,pedestrian,service]")
map = Map("cycle-routes")
for route in routes:
name = route.name if route.name else "Unnamed route"
print(f"{route} {name} ({route.ref}):")
count = 0
for member in route.members.ways:
map.add(member, color="blue" opacity=0.3, weight=20)
for node in member.nodes:
has_intersection = False
for way in streets.parents_of(node):
if not route in way.parents:
name = way.name if way.name else "without name"
print(f" - Intersects with {way.highway} {name}")
has_intersection = True
if has_intersection:
map.add(node, color="red")
count += 1
print(f" {count} intersections.")
map.show()
You can run it with
python cycleways.py
(If we wanted to get fancy, we could analyze the intersections further to determine if they are “sign-worthy” – if the route follows a main road, and a dirt track veers off to the right, there’s no need to post signage. Likewise, if the route goes straight through an intersection of two residential roads, it probably does not need a sign, either. When I get a chance, I’ll flesh this out into a full-fledged example.)
