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')
        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 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,
    retain_all=True,custom_filter=cf_add_routes)
# 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])
try:
    path_length = nx.shortest_path_length(T,node1,node2,weight='length')
except:
    print('no path found: ', description)
else:
    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:

Trails and Graph Theory: The Series

Landing page for the entire series of posts on graph theory and trails.

It all started with a question:

What is the longest continuous route on existing trails in the Gila National Forest, that does not repeat any segment (but intersections are OK)?

And what is the longest continuous loop?

Can the result be mapped and documented as a Long Trail to share with hikers?Sagebrush

Graph Theory is a branch of mathematics that can help us analyze a trail network. A graph has nodes (or points or vertices) that are connected by edges, and can be represented visually: