Category: Geek Hiking

In a previous article we showed how to extract a rudimentary databook from a route on a OSMnx graph. Now let us see how we can improve the output.

Reversing Simplify

One problem noticed since that article is that the OSMnx simplify_graph() function that removes 2nodes can cause trail names to be concatenated together, making it hard to know what trail you are supposed to be on. Another issue is that we want to list in our databook each time a trail crosses a road or another trail, and that information was lost when we simplified and removed 1nodes.

After simplify_graph(), here is an example of a ‘name’ attribute on a trail section (edge):

Cooney Canyon Trail (201),Cooney Canyon/Mineral Creek Trail,McKean                             

Or another example:

Bursum Road,Middle Fork Trail,Snow Canyon Trail #142,Forestry Road 652,Forest Road 1421,Loco Mountain Road   

These concatenated names are also displayed on our map when you hover over trails, and also on my tracks that I exported for use by a GPS app.

The simplify_graph() function has an argument that preserves all the node pairs in a simplified edge, saving them in the edge attribute ‘merged_edges‘. From the node pairs we should be able to go back to the original graph and “desimplify” to see all individual trail segments before they are concatenated together.

for u,v,k,d in J.edges(keys=True,data=True):
    if 'merged_edges' in d:
        merged_edges = d['merged_edges']
    else:
        merged_edges = [ [u,v] ]

    if len(merged_edges)==0:
        continue

    for s,t in merged_edges:
        if GTR.has_edge(s,t):
            for kg in GTR[s][t]:  #fix improve
                dat = GTR[s][t][kg]
                K.add_edge(s,t,key=None,**dat)
                break
        elif GTR.has_edge(t,s): # in case 1-way divided street around Silver City
            for kg in GTR[t][s]:  #fix improve
                dat = GTR[t][s][kg]
                K.add_edge(t,s,key=None,**dat)
                break
        else:
            print('FATAL ERROR WITH MERGED EDGES',s,t) 
            continue

The resulting graph is rather large, with edges typically 0.1 miles in length or smaller, showing each bend in the trail, going from 217 nodes in the simplified graph to over 40000 in the unsimplified.

[link to map full screen]

We will be able to use the desimplified graph to create a more detailed databook, which we explore in the next post.

Download source code here.

Related Posts:

• Trails and Graph Theory: The Series

In a previous post we imported elevation information into our graph and databook, using a JSON query to an open data website.

#https://stackoverflow.com/questions/68534454/python-obtaining-elevation-from-latitude-and-longitude-values/68540685#68540685
def get_elevation(lat, long):
    return 0
    query = ('https://api.open-elevation.com/api/v1/lookup'
             f'?locations={lat},{long}')
    r = requests.get(query).json()  # json object, various ways you can extract value
    elevation = pd.json_normalize(r, 'results')['elevation'].values[0]
    return elevation # returns in meters, not freedom units

Now, in preparation to revisiting our databook code, we will take advantage of the elevation module in OSMNX, reading in elevations for all nodes in the graph with one function call. (Even though the API call is add_node_elevations_google with ‘google‘ in the name, we are not required to use the Google service, for which I do not have a key).

ox.settings.elevation_url_template = 'https://api.open-elevation.com/api/v1/lookup?locations={locations}'

ox.elevation.add_node_elevations_google(J, api_key=None,
        batch_size=350,
        pause=1.0,
            )

(If you use the service at open-elevation.com, please throw them a donation.)

We can colorize our nodes and edges by elevation, making the crude approximation that an edge elevation is the mean elevation of its two nodes.

def colorize_elevation(Q):
    QR = Q.copy()
    for node, dat in Q.nodes(data=True):
        dat['color'] = dat['elevation']
        QR.add_node(node,**dat)                #replace node value

    for u, v, k, dat in Q.edges(keys=True,data=True):
        e1 = Q.nodes[u]['elevation']
        e2 = Q.nodes[v]['elevation']
        dat['color'] = (e1 + e2)/2.0
        QR.add_edge(u,v,key=k,**dat)
    QR.graph['color']=True
    return QR

A small change to our draw() function uses magma as our pre-defined cmap.

Because my graph is simplified, meaning many 2nodes are merged together, the elevation appears to have large steps. We will address this in the next post.

Download source code here.

Related Posts:

• Trails and Graph Theory: The Series

Now that we have been able to create GPX files and a databook for our trip, a logical improvement would be to include water sources. OpenStreetMaps has several tags for identifying water, which seem to have been added over time in an ad-hoc manner. Let us include as many relevant options as we can find.

water_tags = {'natural':['spring','water'],
              'water':['pond','lake','reservoir'],
              'landuse':['reservoir','basin'],
              'amenity':['watering_place','drinking_water','water'],
              'waterway':'stream',
              'man_made':'water_tank'
        }

### ###########  FIND WATER SOURCES    #############
bbox_gila = (34.35704, 32.77935, -107.65983, -108.94574) # rough box of Gila NF
wgdf = ox.features.features_from_bbox(bbox=bbox_gila, tags=water_tags)

The water features are imported as polygons. That is probably more detail than we need, so let’s reduce the polygon to a point, with the view of being able to import a list of waypoints.

# I really did not want to learn about GeoPandas, but there you go...
water_points = wgdf.copy()
water_points['geometry'] = water_points['geometry'].centroid

We understand that reducing from polygons to points is an approximation, and might not work so well for rivers and streams and large lakes. Still, we approximate for convenience, and see if the results are close enough to be useful.

We should be able to convert the GeoDataFrame back to a “network” with only nodes and no edges, but our OSMnx library functions do not work well with “networks” with no edges. We could add one dummy edge, or just stay with GeoDataFrames. Eventually we are able to plot results on our map.

Zoom in…

With more than 9k water sources in the Gila imported, this is too many to be practical. Some waypoints seem to be multiple locations along a creek. Other points are so numerous in areas, that we begin to doubt their veracity.

Let us filter out all water sources that are more than X distance from our route.

WJ = WG.copy()

E,D = ox.distance.nearest_edges(JD,X=water_points['x'],Y=water_points['y'],return_dist=True)

max_distance = 2000 #meters
EARTH_RADIUS_M = 6_371_009 # new convenience notation for big numbers 
max_distance = max_distance / EARTH_RADIUS_M

# Not clear from documentation, but nearest_edges returns distance in units of
# earth radius if projection is lat/long.
for node,d in zip(WG.nodes(),D):
    if d > max_distance:
        WJ.remove_node(node)

With the distance filter, our water sources are down to 77. Here they are on the map, and when you hover over a source, it shows some description of the water.

We can also plot the water sources with our route. Water sources appear as darker blue dots.

We could try to add these water sources to our databook, but perhaps the easiest approach is just to import water source waypoints to our favorite GPX app.

Download source code here.

Related Posts:

• Trails and Graph Theory: The Series