2

I have a GeoJSON which shows an area like so:

enter image description here

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;

enter image description here

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;

1 Answer 1

2

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.

enter image description here

enter image description here

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.