CODEWITHSUNDEEP
osmnx
rtphokie
rtphokie

Reputation: 686

osmnx returning .... unexpected... bearings

I'm trying to find all streets in a given city that run along a specific bearing, I'm building up a list of those edges like this:

G = ox.graph_from_place('Milpitas, CA')
G = ox.add_edge_bearings(ox.get_undirected(G))
sixtydegrees = []
uninteresting = []
for u, v, k, data in G.edges(keys=True, data=True):
    if np.isnan(data['bearing']):
        uninteresting.append((u, v, k))
    else:
         b = round(data['bearing'])
         if (b == 60 or b == 240):
             sixtydegrees.append((u, v, k))
         else:
             uninteresting.append((u, v, k))
G.remove_edges_from(G.edges - uninteresting)

However that produces edges like the one below that include multiple segments along multiple bearings. While the beginning and end nodes are 60º from each other, the edge is made up of 4 segments none of which run along 60º. It's imaginary line running through the park, connecting origin and destination nodes that runs along 60º.

enter image description here

osmnx.bearing.add_edge_bearings is (by design) a simple calculation of the angle between origin node to destination node.

To really find streets, which probably means sections of streets as described above given how OSM data is organized, that are along a given bearing, am I going to have to iterate through the xy pairs and calculate each bearing? Or is there another more elegant solution I'm missing?

Upvotes: 0

Views: 308

Answers (1)

rtphokie
rtphokie

Reputation: 686

Here's what I've got working, wondering if there is a better way.

for u, v, k, data in G.edges(keys=True, data=True):
    # iterate over edges
    if v==u:
        # loops don't have meaningful bearings
        uninteresting.append((u, v, k))
    else:
        lons, lats = data['geometry'].xy
        lat1 = None
        lon1 = None
        hit = False
        for lon, lat in zip(lons, lats):
            if lat1:
                b = round(ox.bearing.get_bearing((lat1, lon1),  (lat, lon)))
                d = round(ox.distance.great_circle_vec(lat1, lon1, lat, lon))
                if d > 250: #ignore segments shorter than 250m
                    hit = hit or (59 <= b <= 61)
            lat1 = lat
            lon1 = lon
        if hit:
            sixtydegrees.append((u, v, k))
        else:
            uninteresting.append((u, v, k))

Upvotes: 0

Related Questions