Quickstart

In this page, we will quickly work through the Julia code necessary to load a network dataset, assign weights to it, and identify and score missing links. To run the code, you'll need to install a recent version of Julia.

I recommend running this code in an interactive, notebook style environment; I use Quarto and Visual Studio Code with the Julia and Quarto extensions. More details about all of the functions from the MissingLinks package are available in the API Documentation.

Load packages

First, we need to load some packages. The Missing Links code is distributed as a Julia package, so we need to load MissingLinks. We also need to load the GeoDataFrames package for spatial data manipulation, the Plots library for visualization, the Graphs library for working with graph data structures, the StatsBase library for statistics, and the Mmap package for working with memory-mapped files. Most of these packages are available from the Julia general registry, and can be installed with the Pkg-mode add command. Mmap is built into Julia and does not need to be installed.

Press ] to enter Pkg-mode.

add GeoDataFrames Plots Graphs StatsBase GeoFormatTypes MissingLinks DataFrames

You'll run this code in the REPL (in VSCode, choose View -> Command Palette and search for "Julia: Start REPL"). All of the other code you'll run should be in your notebook interface so you can save it. While it is optional, I do recommend creating a Julia environment to keep MissingLinks code separate from any other Julia projects you may work on.

Press backspace to exit Pkg-mode.

Once packages are installed, we are ready to load them.

using MissingLinks, GeoDataFrames, Plots, Mmap, Graphs, StatsBase, GeoFormatTypes, DataFrames

If you're using Quarto and VSCode, place this code between triple-backticks to create a Julia cell, then hover over it and click "run cell":

```{julia}
using MissingLinks, GeoDataFrames, Plots, Mmap, Graphs, StatsBase, GeoFormatTypes, DataFrames
```

Likewise, all the other code in this quickstart will also go in a cell.

Loading data

Next, we will read our data. To keep things fast for this example, we are working with a tiny segment of the Charlotte pedestrian network, in the Northlake area. We need data on the pedestrian network, and origins and destinations. In this case we are using residential parcels as origins, and a selection of commercial destinations in the Northlake area digitized from OpenStreetMap as destinations. The data are included with the MissingLinks package; running

MissingLinks.get_example_data()

"example_data"

will create a folder example_data in the current working directory, containing the data.

All of our layers are already projected to EPSG:32119 (State Plane North Carolina, meters). I recommend using a meter-based coordinate system. The system will not work correctly with a geographic coordinate system, and has not been tested with feet.

sidewalks = GeoDataFrames.read("example_data/Northlake.gpkg")
parcels = GeoDataFrames.read("example_data/Northlake_Parcels.gpkg")
destinations = GeoDataFrames.read("example_data/Northlake_Destinations.gpkg")

We can plot what we read in:

# sidewalks have an elevation component that we need to remove to plot
sidewalks.geom = MissingLinks.remove_elevation!.(sidewalks.geom)


plot(parcels.geom, color="#0f0", legend=false, aspect_ratio=:equal)
plot!(sidewalks.geom, color="#f00")
plot!(destinations.geometry, color="#f00")

Building the network

Next, we need to convert our network GIS layer into a mathematical "graph"—a data structure of nodes and edges representing intersections and streets. This network dataset is already "fully noded"—i.e. all intersections occur at the ends of line segments. The tool can also work with "semi-noded" data—where all intersections are at the end of one line segment, but may be in the middle of another. To use that kind of data, run the preprocessing function semi_to_fully_noded before building the graph.

To build the graph, we will use the graph_from_gdal function, which will convert GeoDataFrames into a graph.

graph = graph_from_gdal(sidewalks)

Meta graph based on a Graphs.SimpleGraphs.SimpleGraph{Int64} with vertex labels of type VertexID, vertex metadata of type Tuple{Float64, Float64}, edge metadata of type @NamedTuple{length_m::Float64, link_type::Union{Missing, String}, geom::ArchGDAL.IGeometry{ArchGDAL.wkbLineString}}, graph metadata given by nothing, and default weight Inf

It is possible to specify more than one layer (for example, if you have sidewalks and crosswalks in separate layers) by separating the layers with commas. It is also possible to set a tolerance for when nodes are considered coincident. However, I recommend leaving this at the default small value and using add_short_edges! to close such gaps, to avoid issues that are described in the linked documentation. Similarly, remove_tiny_islands can be used to remove small "islands" in the graph that may be data errors. There are also tools for troubleshooting the graph, namely find_dead_ends to find locations that are dead ends, and find_disconnected_crossings to find places where sidewalks cross without intersecting. While both of these things will be present in a correct graph, they may also represent data errors.

Alternate graph construction using OpenStreetMap data

It is also possible to construct a graph from OpenStreetMap data. First, obtain an OpenStreetMap data file in PBF format (e.g. from slice.openstreetmap.us). Then you can create the graph using the following syntax:

graph_from_osm(pbf, TraversalPermissionSettings(), projection)

pbf is a path to the PBF file. TraversalPermissionSettings define what ways should be included in the graph; see TraversalPermissionSettings for details on customization options.

MissingLinks always works with Euclidean coordinates, so coordinates need to be projected. projection is a GeoFormatTools coordinate system to project the OSM data to. e.g. for North Carolina I use GeoFormatTools.EPSG(32119) (you'll need to install and import GeoFormatTools).

Assigning weights to nodes

Next, we need to assign weights to the nodes. We will create two sets of weights, for origins and destinations. For origins, the weights are just the number of units on the parcel. The weight of each parcel will get divided among the graph edges within 20 meters, and then the weight for each edge will be evenly divided between the nodes it connects. [:,1] retrieves the first column of the weights (if multiple weight columns are specified instead of just units, there will be one column in the output for each column specified. This avoids re-doing computationally-intensive geographic calculations if there are more weight columns—for instance, if both our origins and our destinations were coded in the parcel layer, or if we wanted to measure accessibility from multiple types of origin).

origin_weights = create_graph_weights(graph, parcels, [:units], 20)[:,1]

1696-element Vector{Float64}:
 1.966666666666667
 2.4333333333333336
 0.6833333333333332
 1.4916666666666667
 3.225
 3.5833333333333335
 1.7833333333333337
 1.75
 2.341666666666667
 4.058333333333334
 ⋮
 0.0
 0.0
 2.999625468164794
 3.0329588014981272
 3.0329588014981272
 0.0
 1.8928571428571426
 0.0
 0.0

We can look at what proportion of the units are served by sidewalks and therefore were assigned to the network. Some units will not be assigned to any edge and therefore will not count towards origin_weights.

sum(origin_weights) / sum(parcels.units)

0.8947839046199701

We can similarly create our destination weights. There is no weight column here—all of the destinations are equally weighted. However, we need a weight column so we just create a column of all ones. Since destinations are further from the sidewalk network due to parking lots, we use a 100 meter distance threshold. In the paper we use parcels and a 20 meter threshold, on the assumption that the edges of parcels will be close to the road, but here we have point data on destinations.

destinations.one .= 1
destination_weights = create_graph_weights(graph, destinations, [:one], 100)[:, 1]
sum(destination_weights) / sum(destinations.one)

0.9333333333333333

Creating the distance matrix

A point-to-point distance matrix is used extensively throughout the algorithm for identify and scoring links. We store distances in this matrix as UInt16 integer meters, to save memory (each UInt16 can represent values from 0-65535, and takes only two bytes of memory, as opposed to four or eight for traditional integers and longs). For a network as small as the one we're working with here, it would be possible to store this distance matrix in memory. However, in real-world applications, often the matrix will be larger than available memory, so storing it on disk and using memory-mapping of the disk file is recommended. Memory-mapping lets the computer treat a portion of the disk as if it were memory, and the CPU and the operating system work to make sure that the data are available quickly when needed.

matrix_file = open("quickstart_distances.bin", "w+")
matrix = Mmap.mmap(matrix_file, Matrix{UInt16}, (nv(graph), nv(graph)))
# this just makes sure the matrix is filled with 0s before we begin. This should not make any real difference as all values should be overwritten but is just a good practice.
fill!(matrix, 0)

fill_distance_matrix!(graph, matrix)

1696-element Vector{Vector{UInt16}}:
 [0x0000, 0x0037, 0x0022, 0x0059, 0x00a7, 0x00db, 0x0072, 0x00ea, 0x006c, 0x0099  …  0x080c, 0x0714, 0x0742, 0x0b0c, 0x0afc, 0x0ae6, 0x0271, 0x01a9, 0x0336, 0x050a]
 [0x0037, 0x0000, 0x0059, 0x0022, 0x0070, 0x00a5, 0x003b, 0x00b4, 0x0036, 0x0063  …  0x07d5, 0x06dd, 0x070b, 0x0ad5, 0x0ac6, 0x0aaf, 0x023a, 0x0172, 0x02ff, 0x04d3]
 [0x0022, 0x0059, 0x0000, 0x007b, 0x00c9, 0x00fd, 0x0094, 0x010d, 0x008e, 0x00bc  …  0x082e, 0x0736, 0x0764, 0x0b2e, 0x0b1e, 0x0b08, 0x0293, 0x01cb, 0x0358, 0x052c]
 [0x0059, 0x0022, 0x007b, 0x0000, 0x004e, 0x0083, 0x0019, 0x0092, 0x0014, 0x0041  …  0x07b3, 0x06bb, 0x06e9, 0x0ab4, 0x0aa4, 0x0a8d, 0x0218, 0x0151, 0x02de, 0x04b1]
 [0x00a7, 0x0070, 0x00c9, 0x004e, 0x0000, 0x0035, 0x0035, 0x0044, 0x003a, 0x0067  …  0x07da, 0x06e2, 0x0710, 0x0ada, 0x0aca, 0x0ab4, 0x023e, 0x0177, 0x0304, 0x04d8]
 [0x00db, 0x00a5, 0x00fd, 0x0083, 0x0035, 0x0000, 0x0069, 0x000f, 0x006f, 0x009c  …  0x080f, 0x0716, 0x0744, 0x0b0f, 0x0aff, 0x0ae8, 0x0273, 0x01ac, 0x0339, 0x050c]
 [0x0072, 0x003b, 0x0094, 0x0019, 0x0035, 0x0069, 0x0000, 0x0079, 0x0006, 0x0033  …  0x07a5, 0x06ad, 0x06db, 0x0aa5, 0x0a96, 0x0a7f, 0x020a, 0x0142, 0x02cf, 0x04a3]
 [0x00ea, 0x00b4, 0x010d, 0x0092, 0x0044, 0x000f, 0x0079, 0x0000, 0x007e, 0x00ab  …  0x081e, 0x0725, 0x0754, 0x0b1e, 0x0b0e, 0x0af7, 0x0282, 0x01bb, 0x0348, 0x051b]
 [0x006c, 0x0036, 0x008e, 0x0014, 0x003a, 0x006f, 0x0006, 0x007e, 0x0000, 0x002d  …  0x07a0, 0x06a7, 0x06d5, 0x0aa0, 0x0a90, 0x0a79, 0x0204, 0x013d, 0x02ca, 0x049d]
 [0x0099, 0x0063, 0x00bc, 0x0041, 0x0067, 0x009c, 0x0033, 0x00ab, 0x002d, 0x0000  …  0x0772, 0x067a, 0x06a8, 0x0a73, 0x0a63, 0x0a4c, 0x01d7, 0x0110, 0x029d, 0x0470]
 ⋮
 [0x0714, 0x06dd, 0x0736, 0x06bb, 0x06e2, 0x0716, 0x06ad, 0x0725, 0x06a7, 0x067a  …  0x02f4, 0x0000, 0x00a7, 0x05f5, 0x05e5, 0x05ce, 0x04a3, 0x056b, 0x03de, 0x08cb]
 [0x0742, 0x070b, 0x0764, 0x06e9, 0x0710, 0x0744, 0x06db, 0x0754, 0x06d5, 0x06a8  …  0x0322, 0x00a7, 0x0000, 0x0623, 0x0613, 0x05fc, 0x04d1, 0x0599, 0x040c, 0x08f9]
 [0x0b0c, 0x0ad5, 0x0b2e, 0x0ab4, 0x0ada, 0x0b0f, 0x0aa5, 0x0b1e, 0x0aa0, 0x0a73  …  0x0300, 0x05f5, 0x0623, 0x0000, 0x0010, 0x0027, 0x089c, 0x0963, 0x07d6, 0x0cc4]
 [0x0afc, 0x0ac6, 0x0b1e, 0x0aa4, 0x0aca, 0x0aff, 0x0a96, 0x0b0e, 0x0a90, 0x0a63  …  0x02f0, 0x05e5, 0x0613, 0x0010, 0x0000, 0x0017, 0x088c, 0x0953, 0x07c6, 0x0cb4]
 [0x0ae6, 0x0aaf, 0x0b08, 0x0a8d, 0x0ab4, 0x0ae8, 0x0a7f, 0x0af7, 0x0a79, 0x0a4c  …  0x02da, 0x05ce, 0x05fc, 0x0027, 0x0017, 0x0000, 0x0875, 0x093d, 0x07af, 0x0c9d]
 [0x0271, 0x023a, 0x0293, 0x0218, 0x023e, 0x0273, 0x020a, 0x0282, 0x0204, 0x01d7  …  0x059b, 0x04a3, 0x04d1, 0x089c, 0x088c, 0x0875, 0x0000, 0x00c7, 0x00c6, 0x0428]
 [0x01a9, 0x0172, 0x01cb, 0x0151, 0x0177, 0x01ac, 0x0142, 0x01bb, 0x013d, 0x0110  …  0x0663, 0x056b, 0x0599, 0x0963, 0x0953, 0x093d, 0x00c7, 0x0000, 0x018d, 0x0361]
 [0x0336, 0x02ff, 0x0358, 0x02de, 0x0304, 0x0339, 0x02cf, 0x0348, 0x02ca, 0x029d  …  0x04d6, 0x03de, 0x040c, 0x07d6, 0x07c6, 0x07af, 0x00c6, 0x018d, 0x0000, 0x04ee]
 [0x050a, 0x04d3, 0x052c, 0x04b1, 0x04d8, 0x050c, 0x04a3, 0x051b, 0x049d, 0x0470  …  0x09c3, 0x08cb, 0x08f9, 0x0cc4, 0x0cb4, 0x0c9d, 0x0428, 0x0361, 0x04ee, 0x0000]

Link identification

We are now ready to identify missing links. The identify_potential_missing_links will do this and return a vector of places in the network that meet the specified criteria (in this case, no longer than 100 meters, and with a network distance of 1000 meters or more; you can change these values in the code below).

all_links = identify_potential_missing_links(graph, matrix, 100, 1000)

896-element Vector{MissingLinks.CandidateLink{UInt16}}:
 77m CandidateLink connecting edge VertexID(30, :node)-VertexID(32, :node)@48m to VertexID(704, :node)-VertexID(705, :node)@0m; network distance 65535m

 72m CandidateLink connecting edge VertexID(32, :node)-VertexID(685, :node)@28m to VertexID(704, :node)-VertexID(705, :node)@0m; network distance 65535m

 68m CandidateLink connecting edge VertexID(66, :node)-VertexID(67, :node)@0m to VertexID(868, :node)-VertexID(1538, :node)@34m; network distance 1074m

 27m CandidateLink connecting edge VertexID(66, :node)-VertexID(67, :node)@0m to VertexID(865, :node)-VertexID(1539, :node)@0m; network distance 1129m

 27m CandidateLink connecting edge VertexID(66, :node)-VertexID(67, :node)@0m to VertexID(864, :node)-VertexID(865, :node)@0m; network distance 1129m

 27m CandidateLink connecting edge VertexID(66, :node)-VertexID(67, :node)@0m to VertexID(863, :node)-VertexID(864, :node)@30m; network distance 1135m

 47m CandidateLink connecting edge VertexID(66, :node)-VertexID(67, :node)@0m to VertexID(1538, :node)-VertexID(1539, :node)@24m; network distance 1098m

 67m CandidateLink connecting edge VertexID(66, :node)-VertexID(68, :node)@36m to VertexID(866, :node)-VertexID(867, :node)@0m; network distance 1003m

 66m CandidateLink connecting edge VertexID(66, :node)-VertexID(68, :node)@36m to VertexID(866, :node)-VertexID(868, :node)@0m; network distance 1003m

 37m CandidateLink connecting edge VertexID(66, :node)-VertexID(68, :node)@36m to VertexID(868, :node)-VertexID(1538, :node)@34m; network distance 1037m

 ⋮
 96m CandidateLink connecting edge VertexID(1645, :node)-VertexID(1647, :node)@149m to VertexID(702, :node)-VertexID(703, :node)@2m; network distance 1224m

 87m CandidateLink connecting edge VertexID(1645, :node)-VertexID(1647, :node)@149m to VertexID(703, :node)-VertexID(716, :node)@37m; network distance 1187m

 91m CandidateLink connecting edge VertexID(1645, :node)-VertexID(1647, :node)@63m to VertexID(715, :node)-VertexID(716, :node)@0m; network distance 1042m

 88m CandidateLink connecting edge VertexID(1645, :node)-VertexID(1647, :node)@50m to VertexID(714, :node)-VertexID(715, :node)@23m; network distance 1008m

 87m CandidateLink connecting edge VertexID(1654, :node)-VertexID(1655, :node)@106m to VertexID(418, :node)-VertexID(419, :node)@30m; network distance 1076m

 92m CandidateLink connecting edge VertexID(1654, :node)-VertexID(1655, :node)@108m to VertexID(420, :node)-VertexID(421, :node)@31m; network distance 1075m

 87m CandidateLink connecting edge VertexID(1654, :node)-VertexID(1655, :node)@106m to VertexID(419, :node)-VertexID(421, :node)@0m; network distance 1076m

 98m CandidateLink connecting edge VertexID(1655, :node)-VertexID(1657, :node)@0m to VertexID(418, :node)-VertexID(419, :node)@30m; network distance 1120m

 98m CandidateLink connecting edge VertexID(1655, :node)-VertexID(1657, :node)@0m to VertexID(419, :node)-VertexID(421, :node)@0m; network distance 1120m

We can plot the links identified (shown here in blue):

all_links_geo = links_to_gis(graph, all_links)
plot(sidewalks.geom, color="#000", legend=false, aspect_ratio=:equal)
plot!(all_links_geo.geometry, color="#4B9CD3")

Link deduplication

Clearly, most of these links are essentially duplicates. The next step is to deduplicate them, by merging links where both ends are within 100m of one another (configurable by changing the number below).

links = deduplicate_links(graph, all_links, matrix, 100)

22-element Vector{MissingLinks.CandidateLink{UInt16}}:
 72m CandidateLink connecting edge VertexID(32, :node)-VertexID(685, :node)@28m to VertexID(704, :node)-VertexID(705, :node)@0m; network distance 65535m

 26m CandidateLink connecting edge VertexID(66, :node)-VertexID(68, :node)@12m to VertexID(865, :node)-VertexID(1539, :node)@3m; network distance 1113m

 45m CandidateLink connecting edge VertexID(114, :node)-VertexID(115, :node)@34m to VertexID(1586, :node)-VertexID(1587, :node)@173m; network distance 1116m

 43m CandidateLink connecting edge VertexID(114, :node)-VertexID(115, :node)@34m to VertexID(1590, :node)-VertexID(1591, :node)@138m; network distance 1081m

 98m CandidateLink connecting edge VertexID(196, :node)-VertexID(198, :node)@43m to VertexID(1124, :node)-VertexID(1138, :node)@12m; network distance 1348m

 86m CandidateLink connecting edge VertexID(266, :node)-VertexID(268, :node)@19m to VertexID(1124, :node)-VertexID(1138, :node)@12m; network distance 1266m

 86m CandidateLink connecting edge VertexID(267, :node)-VertexID(269, :node)@81m to VertexID(1142, :node)-VertexID(1143, :node)@12m; network distance 1479m

 98m CandidateLink connecting edge VertexID(295, :node)-VertexID(297, :node)@49m to VertexID(1139, :node)-VertexID(1141, :node)@77m; network distance 1532m

 23m CandidateLink connecting edge VertexID(354, :node)-VertexID(359, :node)@6m to VertexID(881, :node)-VertexID(1039, :node)@119m; network distance 1862m

 20m CandidateLink connecting edge VertexID(876, :node)-VertexID(882, :node)@12m to VertexID(1155, :node)-VertexID(1163, :node)@11m; network distance 2008m

 ⋮
 20m CandidateLink connecting edge VertexID(447, :node)-VertexID(1166, :node)@92m to VertexID(1158, :node)-VertexID(1159, :node)@29m; network distance 2475m

 20m CandidateLink connecting edge VertexID(487, :node)-VertexID(488, :node)@69m to VertexID(1366, :node)-VertexID(1367, :node)@25m; network distance 65535m

 7m CandidateLink connecting edge VertexID(487, :node)-VertexID(488, :node)@60m to VertexID(489, :node)-VertexID(490, :node)@67m; network distance 65535m

 28m CandidateLink connecting edge VertexID(676, :node)-VertexID(677, :node)@75m to VertexID(704, :node)-VertexID(705, :node)@0m; network distance 65535m

 91m CandidateLink connecting edge VertexID(629, :node)-VertexID(686, :node)@58m to VertexID(704, :node)-VertexID(705, :node)@0m; network distance 65535m

 26m CandidateLink connecting edge VertexID(667, :node)-VertexID(669, :node)@7m to VertexID(700, :node)-VertexID(701, :node)@35m; network distance 65535m

 80m CandidateLink connecting edge VertexID(674, :node)-VertexID(675, :node)@56m to VertexID(704, :node)-VertexID(705, :node)@0m; network distance 65535m

 85m CandidateLink connecting edge VertexID(703, :node)-VertexID(716, :node)@16m to VertexID(1641, :node)-VertexID(1647, :node)@126m; network distance 1231m

 88m CandidateLink connecting edge VertexID(714, :node)-VertexID(715, :node)@23m to VertexID(1645, :node)-VertexID(1647, :node)@50m; network distance 1008m

We can now plot these deduplicated links:

links_geo = links_to_gis(graph, links)
plot(sidewalks.geom, color="#000", legend=false, aspect_ratio=:equal)
# set line width wider to make the links easier to see
plot!(links_geo.geometry, color="#4B9CD3", lw=2)

Link scoring

The final step is to score the identified links, using the score_links function. For this we need to specify what our origin and destination weights are, as well as our distance decay function. Here, I use a 1-mile (1609 meter) cutoff. This will return a vector with the score for each link. We have to specify both the distance decay function (first argument) and the point at which that function goes to zero (last argument); for a smooth decay function, I recommend instead specifying a piecewise function that goes to zero when the result is small enough that additional destinations do not materially affect accessibility.

link_scores = score_links(d -> d < 1609, graph, links, matrix, origin_weights, destination_weights, 1609)

22-element Vector{Float64}:
    0.0
    0.0
    1.727099085103875
    1.727099085103875
  767.2176519158897
  741.2873119298721
  537.5906372737354
  573.9455530960228
 4248.055325817284
 3545.6235159582443
    ⋮
 1982.501313235136
  953.7111111111101
  951.6982204893737
    0.0
    0.0
    0.0
    0.0
    0.0
    0.2789844733029542

Saving links to GIS

Generally, you will want to export your results to GIS for further analysis, including the scores. You can create a GeoDataFrame including the scores like this:

links_geo = links_to_gis(graph, links, :score=>link_scores)

You can then write this to GIS file like this:

GeoDataFrames.write("links.gpkg", links_geo)

We can also use it in further plotting. For example, this plot shows the top five links (the competerank function just ranks values largest to smallest).

plot(sidewalks.geom, color="#000", legend=false, aspect_ratio=:equal)
# set line width wider to make the links easier to see
plot!(links_geo[competerank(links_geo.score, rev=true) .<= 5, :geometry], color="#4B9CD3", lw=2)

The new links that cross the road on the east are most valuable from an accessibility standpoint. This is likely in part due to the accessibility metric we chose; the links in the west are simply too far from the destinations to be meaningful for increasing access within 1 mile.

Link-level analyses

To understand the impact of an individual link, MissingLinks provides two functions: routing and regional access calculations (isochrone support is planned). Each one is a higher level of aggregation than the last.

First, we'll find the highest scoring link, though these same procedures would work with any link. We also plot it so we can see what we're working with:

best_link = links[argmax(link_scores)]

plot(sidewalks.geom, color="#000", legend=false, aspect_ratio=:equal)
# set line width wider to make the link easier to see
plot!(links_geo.geometry[argmax(link_scores)], color="#4B9CD3", lw=4)

The link is a crossing in the southeast part of the graph.

One-to-one routes

The lowest level of aggregation is the one-to-one route. This allows us to visualize the impact of the link on a particular trip from one location to another.

Next, we'll create a route from one point to another (specified as latitude/longitude):

# all routing functions require an edge index
# we can build it once and cache it; this is optional
# but vastly improves performance with large graphs
index = MissingLinks.index_graph_edges(graph)

# now, we can do the route without the link
route_without_link = route_one_to_one(
    graph,
    (35.3393, -80.8620), # origin
    (35.3408, -80.8555), # destination
    crs=GeoFormatTypes.EPSG(32119), # the projection of your graph
    index=index # the prebuilt index
)

MissingLinks.OneToOneRouteResult(2228.7137101155386, Geometry: LINESTRING (440357.998913449 177889.83504524,44036 ... 7515))

Including the link

route_one_to_one always routes only on the edges in the graph. To use a potential link, we need to create a new graph that has the link in it. We can use realize_graph to create a new graph with the link added (we can also specify multiple links, by separating them with commas inside the []):

graph_with_link = realize_graph(graph, [best_link])

Meta graph based on a Graphs.SimpleGraphs.SimpleGraph{Int64} with vertex labels of type VertexID, vertex metadata of type Tuple{Float64, Float64}, edge metadata of type @NamedTuple{length_m::Float64, link_type::Union{Missing, String}, geom::ArchGDAL.IGeometry{ArchGDAL.wkbLineString}}, graph metadata given by nothing, and default weight Inf

And then we can use that to create the route. If you are only doing one route, you could leave off the line that creates the index, and remove the index parameter to route_one_to_one; the index will then be automatically created.

index_with_link = MissingLinks.index_graph_edges(graph_with_link)

route_with_link = route_one_to_one(
    graph_with_link,
    (35.3393, -80.8620), # origin
    (35.3408, -80.8555), # destination
    crs=GeoFormatTypes.EPSG(32119), # the projection of your graph
    index=index_with_link # the prebuilt index
)

MissingLinks.OneToOneRouteResult(801.6475169585985, Geometry: LINESTRING (440357.998913449 177889.83504524,44036 ... 7515))

And we can plot the result, with the original route in blue and the new one in red:

plot(sidewalks.geom, color="#000", legend=false, aspect_ratio=:equal)
plot!(route_without_link.geom, color="#4B9CD3", lw=4)
plot!(route_with_link.geom, color="#a24040", lw=4)

We can also write the routes out for use in GIS (in this case, to routes.gpkg), by converting to a DataFrame and saving to a GIS file:

routes = DataFrame([route_without_link, route_with_link])
metadata!(routes, "geometrycolumns", (:geom,))
# Change the CRS code to the CRS of your graph so the routes display properly in GIS
GeoDataFrames.write("routes.gpkg", routes, crs=GeoFormatTypes.EPSG(32119))

"routes.gpkg"

Regional access

The regional access metric shows for every point around a link, how much the access from that point increases due to the link. Let's see an example:

access = regional_access(
    d -> d < 1609, # decay function, we consider locations accessible that are within 1 mile (1609 meters)
    graph,
    matrix,
    best_link,
    destination_weights, # using origin_weights here would map how many people can get to each location
    1609 # point at which our decay function is (essentially) zero - in this case, 1609 meters
)

┌ 66×66 Raster{Float64, 2} ┐
├──────────────────────────┴───────────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 439059.59136290743:50.0:442309.59136290743 ForwardOrdered Regular Points,
  → Y Sampled{Float64} 179763.1496767101:-50.0:176513.1496767101 ReverseOrdered Regular Points
├────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(X = (439059.59136290743, 442309.59136290743), Y = (176513.1496767101, 179763.1496767101))
└──────────────────────────────────────────────────────────────────────────────┘
 ↓ →       1.79763e5  1.79713e5  …  1.76613e5  1.76563e5  1.76513e5
 4.3906e5  0.0        0.0           0.0        0.0        0.0
 4.3911e5  0.0        0.0           0.0        0.0        0.0
 4.3916e5  0.0        0.0           0.0        0.0        0.0
 4.3921e5  0.0        0.0           0.0        0.0        0.0
 ⋮                               ⋱             ⋮          
 4.4216e5  0.0        0.0           0.0        0.0        0.0
 4.4221e5  0.0        0.0           0.0        0.0        0.0
 4.4226e5  0.0        0.0        …  0.0        0.0        0.0
 4.4231e5  0.0        0.0           0.0        0.0        0.0

# xlim and ylim just set the map extents, adjust or remove to map different areas
plot(access, xlim=(439200, 441500), ylim=(177500, 179500))
# add sidewalks on top of the map
plot!(sidewalks.geom, color="white")
# and the link on top of that
plot!(links_geo.geometry[argmax(link_scores)], color="#e84646ff", lw=4)

The largest increases in access happen on the west side of the link, presumably because most destinations are on the east side.

We can also save this to GIS (GeoTIFF format):

write("access.tif", access)

"access.tif"