Category: Geek Hiking

After adding a few short roads to fill in gaps to a trail network, my algorithm (for finding the longest continuous loop without repeating any trail section) hit a brick wall, and stopped giving improving results. It seems my code has some edge cases that were not being handled correctly with the new expanded network, and sleuthing is required.

I will go through fixing edge cases in the next post, but here we show some ways to visualize the trail graph, which will aid in debugging.

Color Nodes and Edges

Important in our algorithm is whether nodes have three edges or four, and we can color nodes and edges to help see what is going on. Here we use the cmap argument in the map drawing function explore.

geopandas.GeoDataFrame.explore(column,cmap,folium.Map, tiles, ... , **kwargs)

The cmap takes some attribute and translates it into a color. The particular color is calculated using the maximum/minimum range for the attribute, and calculates a color within the cmap spectrum. Several pre-built cmap spectrums are already defined in the matplotlib library. I do not particularly care what colors are chosen, as long as they are visually distinct. Here are a couple of cmaps I have used.

Let us name the attribute we add to nodes and edges ‘color‘, and calculate colors for nodes, like so:

    for node, dat in Q.nodes(data=True):
        d = Q.degree(node)
        dat['color'] = d

The way the cmap is used in the explore function is like this:

    m = edges.explore(column='color', cmap = 'tab10', legend=False, tiles="openstreetmap")

And when we draw a trail map, nodes and edges appear to help us diagnose trail connections.

For this map, nodes with odd number of edges are in blue, and nodes with 4 edges are in turquoise. Trails are purple, and road sections are turquoise.

We can add extra colorize functions to explore other properties of the graph. A bridge on a graph is an edge that, if cut, divides the graph into separate subgraphs. Bridges are a measure of the connectivity of an edge. When we start eliminating edges so that all nodes are even, we probably do not want to break up the graph very much, so in the future we may weight these edges to minimize the chance of deleting them. Find bridges in a graph with the Networkx function:

nx.bridges(T)

And we write a function colorize_edge_list to highlight these edges:

bridge_list = list(nx.bridges(T))
print('length of bridges:', len(bridge_list))
T = colorize_edge_list(T, bridge_list)
draw(T,'bridges on simplified graph', show_nodes=True)

During debugging, one can imagine writing several special-purpose colorize_xx() functions, for highlighting matched nodes, displaying unmatched nodes, showing all 4nodes that are connected directly to another 4node, and so forth.

While working with our draw() function, we might as well add or fix a few other features.

Tool Tip (Hover Text)

Previously when hovering over edges on our map, all attributes were displayed, which is rather ugly.

But the explore function has argument that selects which attribute, or column, do display as a tool tip, and ‘name’ is a good option for OSMnx plots:

m = edges.explore(column='color', cmap='tab10', legend=False, tiles="openstreetmap", tooltip='name')

Town Markers

We had previously used Markers to display features we were debugging, such as small gaps in our trail network. Now we choose to use markers to display potential trail resupply towns. The display code is already written, located in the draw() function several blog posts ago, so the only thing left is to add the array of town GPS coordinates.

town_markers = [(34.29963, -108.13268, 'Pie Town'),
         (33.19593, -108.20834, 'Gila Hot Springs'),
         (32.77981, -108.27533, 'Silver City'),
         (33.37939, -108.90369, 'Alma/Glenwood'),
         (33.03022, -108.16864, 'Lake Roberts'),
         (32.91685, -107.70593, 'Kingston')]

Three Dimensions

Looking towards the far future, assuming my tools are debugged to find a good long distance hiking route, it would be really great if we could view our trail in 3D, for planning a hike. OpenStreetMap.org does not currently have a 3-D display option, but height is included in OSM data. A few open-source projects do display OSM in 3D. OSM Buildings is one effort, but seems oriented towards city buildings.

Another project is Streets.gl, which looks better for what we want.

As far as I can discover, the website only displays roads, and not trails, but it is open source, so this might be a great project for any of you trail geeks.

Source code for coloring will be included in the next post, which should come out soon.

Related Posts:

• Trails and Graph Theory: The Series

In the previous post we learned how to automagically close short gaps in our trail network, and now we will develop a more refined way of adding longer road-walks without using a long and tedious list of rectangles to define import boundaries as in this post.

Our method involves importing a graph of the Gila that includes both trails and roads, with variable name GTR. With nodes node1 and node2, we find the shortest distance between these nodes in GTR

        path = nx.shortest_path(GTR,node1,node2,weight='length')

and then add that path to our original trails graph G. We choose node1 and node2 so that the shortest path between them is very likely to be the road we want to import, and not some alternate path. , and then import this path into our trail graph. The trail graph, G, and the road-and-trail graph GTR share the same node numbers, which is the unique OSMid for that location, which makes this method of selectively importing roads from GTR into G possible.

It is fairly easy to look at our map and come up with a list of roads we want to include. (I use the CalTopo website map because its user interface is rather good for finding longitudes and latitudes, and measuring distances.)

road_list = [
    ( (32.82663, -108.34204), (32.79465, -108.18348), 'silver city cottage san, silver heights, arenas' ),
    ( (33.01687, -107.94504), (32.93080, -108.01376), 'kelly mesa connector' ),
    ( (33.02712, -108.13412), (33.04191, -108.22340), 'lake roberts' ),
    ( (32.91072, -107.70901), (32.92754, -107.75264), 'kingston' ),
    ( (33.41948, -108.82821), (33.31934, -108.85576), 'alma/glenwood' ),
    ( (33.39806, -108.58210), (33.37838, -108.76716), 'bursum road near willow creek campground' ),
    ( (33.19050, -108.18217), (33.23077, -108.26747), 'gila hot springs' ),

    #cdt middle fork alt
    ( (33.42041, -108.49880), (33.45350, -108.49304), 'snow lake' ),
    ( (33.49220, -108.48115), (33.74251, -108.47609), 'fr 652, fr94' ),
    ( (33.74230, -108.47677), (33.74605, -108.48022), 'fr3070' ),

    ###### cdt wilderness study area
    ( (33.69325, -108.24452), (33.68589, -108.18571), 'cdt study area near pelona mtn' ),
    ( (33.67700, -108.05966), (33.68507, -108.02283), 'cdt study area' ),

    # aldo
    ( (33.55441, -107.81135), (33.57859, -107.80749), 'fr 4052P aldo' ),
    ( (33.51417, -107.83057), (33.52641, -107.80248), 'fr 4052T aldo' ),
    ( (33.34820, -107.81522), (33.33652, -107.85583), 'near lookout mtn' ),
    ( (33.33497, -107.82658), (33.34089, -107.84467), 'fr226 aldo' ),
    ( (33.11088, -107.73054), (33.11153, -107.76046), 'north seco creek' ),
    ( (33.29079, -108.05481), (33.27013, -108.15834), 'fr225 aldo' ),

    ]

We also do the mind-the-gap routine from the previous post after adding the longer roadwalks, because these additions fix some of the short gaps for free. We also have refactored our mind-the-gap code to utilize what we learned in this post, so it runs much faster and makes fewer errors in closing gaps.

The resulting trail graph is now pleasingly connected.

I might make one or two tweaks, but the graph finally looks in good shape. What happens when we run our gradient fuzz search routine now? Actually, the results are not very improved. My current theory is that my trail search algorithm does not correctly handle the case where two 4-nodes are neighbors with each other. With this better-connected graph, there are more 4-nodes, and a higher probability of 4-nodes as neighbors. I will attempt to fix that edge-case next time.

Download source code here.

Related Posts:

• Trails and Graph Theory: The Series

In our last post, we tried to add short connecting roads between trails in our quest to improve the longest non-repeating trail loop in the Gila National Forest. Several trails, the CDT in particular, will have some miles of footpath, then a mile of forest road, and then transitions back to footpath.

Trying to find all these gaps by hand, and debugging mistakes, seems like the wrong approach. We will try to modify our procedure, writing code to automagically find these gaps and import short spans of road into our trail network.

The first change is to replace the command that imports a graph from within the Gila National Forest boundary, and change it to a box roughly approximating the forest. This allows us to include the Continental Divide Study area, the Gila Cliff Dwellings National Monument, and several small private land inholdings and wildlife study areas.

n,s,e,w = 33.92485, 32.73704, -107.65983, -108.94574 # rough box of Gila NF
G =  ox.graph_from_bbox(n,s,e,w,simplify=False,retain_all=True,custom_filter=cf)

(Writing and debugging the program took many iterations, and several steps took a loooonnnnnggg time to execute, so for each step we save a Pickle file and are able to load it next time to skip ahead to the next step.)

Our first task is to make a list of all trailheads, that is, nodes with only one neighbor in our trail graph, that are really close to a node in our road graph.

# Because our program development needs a great deal of iteration,
# store intermediate results in pickle files
def find_trailheads(Q,QR,road_distance_limit = 30.0):
    trailheads_filename = 'trailheads.pkl'
    trailhead_nodes = set()
    if os.path.exists(trailheads_filename):
        trailhead_nodes = pickle.load(open(trailheads_filename, 'rb'))
        print('finished loading trailhead nodes')
        elapsed_time(start_time)
    else:
        for node in Q.nodes:
            d = Q.degree(node)
            if d == 2: #Q is multidigraph, so this is an end-point of trail
                if distance_road(Q,node,QR) < road_distance_limit :
                    trailhead_nodes.add(node)
                    print('added node to trailhead list: ', node)
                    continue
        pickle.dump(trailhead_nodes,open(trailheads_filename, 'wb'))
        print('finished calculating trailhead nodes')
        elapsed_time(start_time)
    return trailhead_nodes

number of trailhead nodes:  238

Now that we have a list of trailheads, we can identify pairs of trailheads that are close to each other, say perhaps 2 kilometers, and make a list of these pairs.

def find_trailhead_pairs(Q, trailhead_nodes, distance_limit = 2000.0):
    trailhead_pairs = set()
    trailhead_pairs_filename = 'trailhead_pairs.pkl'
    if os.path.exists(trailhead_pairs_filename):
        trailhead_pairs = pickle.load(open(trailhead_pairs_filename, 'rb'))
        print('finished loading trailhead pairs')
    else:            
        for node in trailhead_nodes:
            for node2 in trailhead_nodes:
                if node==node2:
                    continue
                if (min(node,node2),max(node,node2)) in trailhead_pairs:
                    continue
                d= distance(Q, node,node2)
                #print('distance between ',node, node2, ' = ', d)
                if d < distance_limit:
                    trailhead_pairs.add((min(node, node2),max(node, node2)))
                    print('trailhead pairs: ', node, node2)
        pickle.dump(trailhead_pairs,open(trailhead_pairs_filename, 'wb'))
        print('finished calculating trailhead pairs')
    elapsed_time(start_time)
    print('length of trailhead_pairs: ',len(trailhead_pairs))
    return trailhead_pairs            

length of trailhead_pairs:  326

Now, using these pairs of trailheads, we can attempt to find a road that goes between the trailhead pairs, and add that road to our graph.

def mind_the_gap(Q,QR,distance_limit = 2000.0):
    print('start mind_the_gap function')
    
    trailhead_nodes = set()
    trailhead_nodes = find_trailheads(Q,QR)
    print('number of trailhead nodes: ', len(trailhead_nodes))
    
    trailhead_pairs = find_trailhead_pairs(Q,trailhead_nodes,distance_limit)

    joined_trailhead_boxes = []
    not_joined_trailhead_boxes = []
    for node, node2 in trailhead_pairs:
        result, QT, box = connect_trailheads_with_road(Q,node,node2)
        if result:
            joined_trailhead_boxes.append(box)
            Q = QT
        else:
            not_joined_trailhead_boxes.append(box)
    
    return Q, joined_trailhead_boxes, not_joined_trailhead_boxes

length of added boxes:  174
length of not-added boxes:  152

We keep a list of boxes that show a successful bridging of trailheads, and a list of boxes that do not successfully bridge trailheads with road sections, to display on our map for troubleshooting purposes. The resulting map looks pretty good, and the process was much easier than manually identifying trailheads and manually adding bounding boxes to add short roadwalks, as done in our previous post. Boxes showing road connections are in orange, and boxes that were not able to connect are in green. Sometimes an orange and green box overlap, which is still a successful connection.

Download source code here.

In our next post, we will show a better way to import longer roadwalks to our graph, for going into trail towns and a road-walks longer than 2km.

Related Posts:

• Trails and Graph Theory: The Series