Trails and Graph Theory 21: More Databook

In a previous article we sketched out how to create a databook for our long hike, a text summary of the route, which can include waypoints, elevation, trail names, and intersections. Our first result was somewhat crude, so we are circling back to make a better effort.

In a recent post we described how to break apart trail sections that had been concatenated together as part of our algorithm to find a maximum length route. We also found a different way to add elevation to our network, using built-in OSMNx functions. We have so many nodes in our desimplified network graph that it is better to take steps to only add elevation to important nodes.

def add_elevation(Q):
    H = nx.MultiGraph()
    H.graph['crs'] = ox.settings.default_crs

    for node, dat in Q.nodes(data=True):
        if 'street_count' in dat:
            street_count = dat['street_count']
            if street_count > 2:
                H.add_node(node,**dat)
    ox.settings.elevation_url_template = \
        'https://api.open-elevation.com/api/v1/lookup?locations={locations}'
    ox.elevation.add_node_elevations_google(H, api_key=None,
        batch_size=350,
        pause=1.0,
            )
    print('Added elevation data for ',H.number_of_nodes(), ' nodes')
    for node,dat in H.nodes(data=True):
        Q.add_node(node,**dat)
    return Q

Another improvement we should make is to document any road crossing or trail crossing on the route. Along with our trail graph G, we use the graph of all trails and roads GTR. Whenever a G node has a ‘street_count’ attribute greater than 2, we expect a crossing, and go back to GTR to attempt to find the name of the road or trail that our route crosses.

def report_path(Q,Q_orig, txt='',freedom_units=False,lat_long=False):
    meters_to_mile = 0.0006213712
    meters_to_feet = 3.28084
    meters_to_km = 0.001
    
    if nx.is_eulerian(Q):
        total_length = 0
        data = []
        
        for u,v,k in nx.eulerian_circuit(Q, keys=True): #fix: use predefined circuit, not arbitrary
             data_line = []
             if freedom_units:
                 data_line.append(str(round(total_length*meters_to_mile,1)))
             else:
                 data_line.append(str(round(total_length*meters_to_km,1)))
                 
             d = Q.get_edge_data(u,v)[k]
            # print(d)
             length=0                 
             if 'length' in d:
                 length = d['length'] #meters
                 total_length += length
             node_attr = Q.nodes(data=True)[v]
            # print( node_attr )
            
             street_count = node_attr['street_count']
             if street_count == 2:
                 continue
             
             long = node_attr['x']
             lat = node_attr['y']
             if lat_long:
                 data_line.append(lat)
                 data_line_append(long)
             #elevation = get_elevation(lat,long) #in meters
             if 'elevation' in node_attr:
                 elevation = node_attr['elevation']
                 if freedom_units:
                     data_line.append(str(round(elevation*meters_to_feet)))
                 else:
                     data_line.append(str(round(elevation)))
             else:
                 data_line.append(' ')

             trail_name=''
             if 'name' in d:
                 trail_name = d['name']
                 if not isinstance(trail_name, str):
                     trail_name = ','.join(trail_name) #in case is a list
             
             data_line.append(trail_name)                     
             
             crossing = set()
             crossing_txt = ''
             for x,neighbor,k,d in Q_orig.edges(v,data=True, keys=True):
                 if 'name' in d:
                     name = d['name']
                     if not isinstance(name, str): #FIX with a foreach
                         name = ','.join(name) #in case is a list
                     if name != trail_name:
                         crossing.add(name)
                 elif 'highway' in d:
                     name = d['highway']
                     name = '<' + name + '>'
                     crossing.add(name)
             if len(crossing)>0:
                 crossing_txt = ','.join(crossing)
             data_line.append(crossing_txt)
                              
             data.append(data_line)                 
        header = []
        if freedom_units:
            header.append('MILE')
        else:
            header.append('KM')
        if lat_long:
            header.append('LAT')
            header.append('LONG')
        header.append('ELEVATION')
        header.append('TRAIL')
        header.append('INTERSECTION')
        table = columnar(data, header, no_borders=True,terminal_width=132)
        print(table)     
    else:
        print('The graph is not eulerian')

Ideally, it might be useful to automatically add nearby water sources. We did document a way to find water sources, but inserting them into our databook looks challenging, so perhaps we will add that feature at a later date. Algorithmically adding other features, such as post offices, stores, or established campsites, are left as an exercise for the reader. 😁

A short example excerpt from the databook shows the improvements:

  MILE   ELEVATION  TRAIL                                    INTERSECTION                                                      

234.6 8373 West Fork Trail #151 Iron Creek Trail #172
234.8 8376 West Fork Trail #151 Iron Creek Mesa Trail #171
237.4 7831 West Fork Trail #151 Cooper Trail #141
240.3 7841 Cooper Trail #141 Clayton Mesa Trail #175
242.6 7106 Clayton Mesa Trail #175 Middle Fork Trail #157
242.7 7090 Middle Fork Trail #157 Iron Creek Mesa Trail #171
246.0 7326 Middle Fork Trail #157 Snow Canyon Trail #142
246.6 7467 Snow Canyon Trail #142 <path>
246.8 7356 Snow Canyon Trail #142 <path>
246.9 7359 Snow Canyon Trail #142 <service>

Download source code and improved databook here.

Related Posts:

Trails and Graph Theory 20: Desimplify

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 19: Elevation

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: