Jim, Sagebrush – Page 40 – Sagebrush Trails

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.