Trails and Graph Theory 13: Longer Roadwalks

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 12: Mind the Gap

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')
        for node in Q.nodes:
            d =
            if d == 2: #Q is multidigraph, so this is an end-point of trail
                if distance_road(Q,node,QR) < road_distance_limit :
                    print('added node to trailhead list: ', node)
        pickle.dump(trailhead_nodes,open(trailheads_filename, 'wb'))
        print('finished calculating trailhead nodes')
    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')
        for node in trailhead_nodes:
            for node2 in trailhead_nodes:
                if node==node2:
                if (min(node,node2),max(node,node2)) in trailhead_pairs:
                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')
    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:
            Q = QT
    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 11: Add Connections

We have built tools to find a fairly long circuit in the Gila National Forest, without repeating any trail segments. We can make a few improvements.

  • Add forest roads, or even a cross-country route if necessary, to connect the Aldo Leopold Wilderness trail system to the rest of the Gila.
  • Add forest roads, paved roads, and BLM trails to get near resupply towns: Silver City, Pie Town, Alma/Glenwood, Gila Hot Springs (Doc Campbell’s), and perhaps others.
  • Include trails in the Gila Cliff Dwellings. Since this is a national monument, its trails were not added when we imported Gila National Forest trails.
  • When trails enter short sections of non-USFS land, such as private or state trust, add connections for those areas, since they would not be included in the original Gila import command.
  • Perhaps add part of Bursum Road, NM 159. It is one of my favorite road walks of the Grand Enchantment Trail, and connects several trails together.
  • Add the Continental Divide Wilderness Study Area, south of Pie Town, that contains a good span of the CDT and connects two northern protrusions of the Gila National Forest.
  • Add several brief road-walks that are part of the Continental Divide Trail and Grand Enchantment Trail, that would not be imported if we only add foot-paths.

The OSMnx library does not allow us a direct way to import GPX files of specific routes, since we need to use the node values already defined in OpenStreetMaps. As far as I know, OSMnx does not support importing specific trails or roads by name. We can import a rectangular region into a graph, and then merge it with our existing Gila trail graph.

As an example, Alma and Glenwood are near Mineral Creek Trail (quite an interesting route with a lot of history and geology) but this trail would normally be excluded from our loop because it dead-ends at a trailhead. By connecting Mineral Creek Trail with short road-walks to Alma (good cafe and small convenience store) and Glenwood (post office and maybe convenience store) and then to Little Whitewater Creek Trail (or Whitewater Creek Trail with the Catwalk), we can include these trails in our loop.

To add these roads, which might include other roads we do not care about, we can get away with a single rectangle.

# custom filter imports roads AND foot-paths
cf_add_routes = '["highway"~"path | footway | 
    motorway | motorway_link | trunk | primary | secondary | 
    tertiary | unclassified | residential | living_street | 
    residential | service | track"]'

n,s,e,w = 33.42443, 33.26350, -108.78571, -108.92392
GT = ox.graph_from_bbox(n,s,e,w,simplify=False,truncate_by_edge=True,
# merge with existing graph G
G2 = nx.compose(GT,G2)

Whenever we add a connection, we add code for a test to make sure any roads and trails that we add properly connect to each other.

description = 'glenwood mineral creek to little whitewater'
Point1 = (33.42070, -108.82503)
Point2 = (33.32292, -108.81121)

node1 = ox.nearest_nodes(T,Point1[1], Point1[0])
node2 = ox.nearest_nodes(T,Point2[1], Point2[0])
    path_length = nx.shortest_path_length(T,node1,node2,weight='length')
    print('no path found: ', description)
    print(description,' distance(miles): ', str(round(path_length*meters_to_mile,1)))
glenwood mineral creek to little whitewater  distance(miles):  14.7

For Silver City, the CDT connects to the city going northwest on a long road walk, and there is a network of trails about 5 miles to the east at Arenas Valley. Corre Caminos bus line goes from Silver City to Arenas Valley a few times each weekday, so we will include these connecting roads to pass through Silver City.

Several short roads connect between trails in the Aldo Leopold Wilderness and Gila Wilderness. When importing only trails into NetworkX, these paths would be interpreted as disconnected. Fixing this might include more of the Aldo. (However, the Aldo is blocked on the south and east, and to the north has few connected trails going out to the rest of the Gila, aside from the CDT, so the effect may be limited.)

We will not detail all the other route additions here, but they are documented in the source code, available for download below. Pie Town did not seem to have a loop opportunity, so that trail town will be added as an optional supply route.

An option was added to our draw() command to display all the connection rectangles we added. This was essential for debugging.

Finding and adding all the brief road connections for the CDT was very tedious and took many iterations, and I am still not finished. When one is making this many additions by hand, it is usually a sign that the problem should be solved another way.

Might it be possible to write a script that finds these gaps? And speaking of added connections, might it be possible to write a program to find forest roads that connect the ends of trails, and span less than 0.5 miles in length or so, and automagically add those road segments to our network?

That will be the subject of the next blog post.

Download the source code, available here. (For illustrative purposes only, because I am sure there is a better solution.)

Related Posts: