Find intersection between highways and GeoJSON polygon

Question

I have a GeoJSON which shows an area like so:

There are some highways, crossing this GeoJSON. I need to output the coordinates of each highway entering and exiting this polygon with Python.

First, I need to find the data for all the highways in this area, and then find a way to output the intersection between them. I could find some data using overpass turbo like so:

[out:json];
way["highway"~"motorway|trunk"]
(49.8194,8.3249,50.3218,9.7202);
out meta;
>;
out skel qt;

but I have no idea if this data is useful in anyway, or if it is, how I can use it. I'm new to all of this.

UPDATE:

I also tried this query which turns nodes to coordinates in the resulting json. But I still have no idea how I may be able to use this.

[out:json];
way["highway"~"motorway|trunk"]
(49.8194,8.3249,50.3218,9.7202);
convert item ::=::,::geom=geom(),_osm_type=type();
out geom;

Answer 1

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)

Find the data from Overpass turbo:

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)