I found the way to do this in case anyone might want to do something similar. I used shapely package to work with geojson data.
Find the bounding box:
Since I wanted to work with the data inside my region, I needed the bounding box of the polygon. So I used bounds
function from shapely like so: Docs
with open('region.geojson', 'r') as r:
r_content = json.load(r)
coordinates = r_content['coordinates'][0]
region = Polygon(coordinates)
print(region.bounds) # output: (7.638, 49.574, 9.469, 50.699)
Use overpass turbo to filter the osm data and get what you need. since I wanted the data for all highways starting with A
in their ref in my boundary, I used the below query with regex: (highways like: ref=A 3
, ref=A 5
, ...)
remember to use the bounding box in the correct order for overpass Docs
[out:xml];
way["ref"~"^A"]
(49.574, 7.638, 50.699, 9.469);
out body;
>;
out skel qt;
This will give all the needed data in xml
format. I wanted to work with geojson
, so I exported it as geojson
from overpass turbo.


Find the intersection points:
Now we have the region and the highway's data, we can find their intersection with the help of shapely intersection
method.Docs
- Remember that the intersection between a Polygon and a LineString is a LineString, and the intersection between two LineString is a Point. So we need to convert the Polygon to a LineString first.
ring = LineString(list(region.exterior.coords)) # convert Polygon to Linestring
with open('A_data.geojson', 'rb') as a:
a_content = json.load(a)
for feature in a_content['features']:
# get all LineStrings which have the highway tag
if feature['geometry']['type'] == 'LineString' and feature['properties'].get('highway'):
line = LineString(Point(list(cor)[0], list(cor)[1]) for cor in feature['geometry']['coordinates'])
if ring.intersection(line):
point = (ring.intersection(line))
print(f"intersection point for {feature['properties']['ref']}: ", point)