Compute edge headings for start and end

[nreinicke]

In our current generate_dataset function, we use the edge bearings generated from osnmx which computes the bearing from the origin node to the destination node. This is a good approximation but RouteE Compass allows us to specify a bearing for both the beginning and end of a road link. This is useful for computing the difference in angle between two links that might have significant curvature.

Imagine a (somewhat contrived) example where we have two links with curvature connected to each other like this:

If we just use the bearing between the start and end nodes, it will appear as if the vehicle is not turning at all between the two links but in reality the vehicle is making an ~90 degree turn.

In order to support this we should use the geometry from the link and then take the first two points in the linestring to compute the arrival_heading and the last two points in the linestring to compute the departure_heading.

[nreinicke]

Adding a bit more context to this:

We'll need to first extract the edge geometries from the incoming networkx graph. We can get the edge geometries in a geopandas dataframe like this:

edges_gdf = ox.graph_to_gdfs(g, edges=True, nodes=False)
edges_gdf["geometry"]

Then, we need to compute the arrival and departure heading for each edge in that dataframe. Here are a couple of functions that we could use:

import math 

def compass_heading(point1, point2):
    lon1, lat1 = point1
    lon2, lat2 = point2

    lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])

    dlon = lon2 - lon1

    x = math.sin(dlon) * math.cos(lat2)
    y = math.cos(lat1) * math.sin(lat2) - (
        math.sin(lat1) * math.cos(lat2) * math.cos(dlon)
    )

    initial_bearing = math.atan2(x, y)

    initial_bearing = math.degrees(initial_bearing)
    compass_bearing = (initial_bearing + 360) % 360

    return compass_bearing


def calculate_bearings(geom):
    if len(geom.coords) < 2:
        raise ValueError("Geometry must have at least two points")
    if len(geom.coords) == 2:
        # start and end heading is equal
        heading = int(compass_heading(geom.coords[0], geom.coords[1]))
        return (heading, heading)
    else:
        start_heading = int(compass_heading(geom.coords[0], geom.coords[1]))
        end_heading = int(compass_heading(geom.coords[-2], geom.coords[-1]))
        return (start_heading, end_heading)