I am excited to introduce a new R package called flownet providing high-performance tools for transport modeling—network processing, route enumeration, and traffic assignment in R. As of now, the package implements the Path-Sized Logit (PSL) model (Ben-Akiva and Bierlaire, 1999)—a powerful and well-established method for stochastic traffic assignment accounting for route overlap—alongside a novel route enumeration algorithm and a highly efficient All-or-Nothing assignment solution. Furthermore, it provides powerful utility functions for (multimodal) network processing, including recursive graph consolidation/contraction, and/or simplification. These features, combined with an accessible representation of graphs and (sparse) origin-destination (OD) matrices using data frames and a parsimonious API, make it a compelling toolbox for both transport analytics and graph manipulation. As has come to be expected of my packages, flownet strives to be computationally efficient, building on a whole array of fastverse libraries including collapse and kit for fast data manipulation, igraph for efficient routing, geodist and leaderCluster for spatial distances and clustering, and mirai for efficient R-level parallelism. These dependencies are supplemented by custom C code, e.g., the PSL is implemented fully in C. There is still significant room for improvements, particularly if igraph can be replaced with internal C/C++ routing code, which could port the entire route enumeration algorithm to C/C++. Such innovation may be the subject of future updates, alongside new traffic assignment methods. The package supports directed and undirected graphs and multimodal networks.1 In the following I demonstrate flownet on two networks: The first is the African continental road network I built in my Job Market Paper. It is included in the package. This network is already consolidated and has 1379 nodes and 2344 edges representing 315,000 km of road network. It is thus a rather small network as far as the capabilities of flownet are concerned. The second network is the US Freight Analysis Framework (FAF5) network comprising 487,384 links (edges) and 348,498 nodes. This network is very large and requires careful tuning of the PSL method. It also provides an excellent test case to demonstrate flownet’s network simplification capabilities. Example 1: Africa’s Continental Road Network The africa_network, taken from Krantz (2024) Optimal Investments in Africa’s Road Network, is included in the package. It was generated by intersecting and simplifying fastest OSRM car routes between 453 African cities with population > 100,000 and international ports—provided in africa_cities_ports. The procedure is described in the paper, and further demonstrated in the introductory vignette using flownet’s network processing tools. For this blog post we will take the network as given and defer the network processing workflow to the second (FAF5) example below. The network includes 481 new links proposed by a network finding algorithm described in Section 4.2 of the paper. Below I plot the network with proposed links in grey and the 453 (port-)cities in black. library(fastverse) ## -- Attaching packages ------------------------------------------------------------------------------------------------------------- fastverse 0.3.4 -- ## v data.table 1.17.0 v kit 0.0.21 ## v magrittr 2.0.4 v collapse 2.1.6 fastverse_extend(flownet, sf, tmap) ## -- Attaching extension packages --------------------------------------------------------------------------------------------------- fastverse 0.3.4 -- ## v flownet 0.2.0 v tmap 4.2 ## v sf 1.0.20 # Plot network colored by travel speed tm_basemap("CartoDB.Positron", zoom = 4) + tm_shape(africa_network) + tm_lines(col = "speed_kmh", col.scale = tm_scale_continuous(values = "turbo", values.range = c(0.1, 0.9), value.na = "grey"), col.legend = tm_legend("Speed (km/h)", position = c("left", "bottom"), frame = FALSE, text.size = 0.8, title.size = 1, item.height = 2, na.show = FALSE), lwd = 1.5) + tm_shape(africa_cities_ports) + tm_dots(size = 0.1) + tm_layout(frame = FALSE, outer.margins = 0) Evidently, OSRM link speeds differ widely, with speeds along many central African roads around 40km/h, and higher speeds closer to urban centers. Below, I further describe the network after removing the proposed links. As detailed in Section 3 of the paper, I extracted border and documentary times and costs of exporting and importing from the last World Bank Doing Business Survey (data collected in 2019), and converted them to travel-time-equivalent frictions (in minutes) along border-crossing links. There is thus a frictionless duration (travel time) variable and border_time (which includes documentary time—median 76h), which added together give the total_time needed to trade across a link. # Removing potential new links africa_net cat(sep = ", ") ## from, to, from_ctry, to_ctry, FX, FY, TX, TY, sp_distance, distance, duration, speed_kmh, passes, gravity, gravity_rd, border_dist, total_dist, border_time, total_time, duration_100kmh, total_time_100kmh, rugg, pop_wpop, pop_wpop_km2, cost_km, upgrade_cat, ug_cost_km, geometry qsu(africa_net, cols = .c(from, to, speed_kmh, distance, duration, border_time, total_time)) ## N Mean SD Min Max ## from 2344 669.4347 391.7331 1 1377 ## to 2344 688.3921 395.8184 2 1379 ## speed_kmh 2344 65.2302 16.0805 5.1429 127.749 ## distance 2344 134411.13 115553.485 16968 1'801830.5 ## duration 2344 135.585 136.9445 12.05 2151.9 ## border_time 2344 687.2841 1910.1195 0 11871.9998 ## total_time 2344 822.869 1947.0973 12.05 12134.6998 The package also provides, in africa_trade, a dataset of bilateral trade flows between the 47 African countries that accommodate any of the 453 selected (port-)cities. It is derived from the CEPII BACI database, HS96 version, averaged over 2012-2022 by 21 HS sections. # Convert to tibble for compact printing of section names qTBL(africa_trade) ## # A tibble: 27,721 × 8 ## iso3_o iso3_d section_code section_name hs2_codes value value_kd quantity ## ## 1 AGO BDI 16 Machinery and mechanical appliances, electrical equipment, parts thereof, sound rec… 84 5.10e+0 5.10 9.87e-1 ## 2 AGO BEN 1 Live animals and animal products 3 6.81e+3 6527. 8.78e+3 ## 3 AGO BEN 4 Prepared foodstuffs, beverages, spirits and vinegar, tobacco and manufactured tobac… 19, 20, … 9.51e+0 9.10 2.07e+1 ## 4 AGO BEN 5 Mineral products 27 3.58e+3 3536. 8.37e+3 ## 5 AGO BEN 6 Product of the chemicals or allied industries 30, 32, … 9.44e-1 0.932 2.75e-2 ## 6 AGO BEN 7 Plastics and articles thereof, rubber and articles thereof 39, 40 7.13e+1 72.0 3.13e+1 ## 7 AGO BEN 8 Raw hides and skins, leather, furskins and articles thereof, saddlery and harness, … 42 1.18e+0 1.17 2.5 e-2 ## 8 AGO BEN 9 Wood and articles of wood, wood charcoal, cork and articles of cork, manufacturers … 44, 46 1.5 e-1 0.140 3.2 e-2 ## 9 AGO BEN 10 Pulp of wood or of other fibrous cellulosic material, recovered (waste and scrap) p… 48, 49 8.9 e-1 0.863 4.36e-1 ## 10 AGO BEN 11 Textile and textile articles 52, 57, … 3.00e+1 29.6 1.08e+1 ## # ℹ 27,711 more rows Our exercise will be to assign these trade flows, in particular the quantity in metric tons, to the network. While certainly not all intra-African trade travels by road, UNECA 2022 estimate, in their 2019 baseline scenario, that 76.7% of intra-African freight traveled by road, followed by 22.1% maritime, 0.9% air, and only 0.3% by rail. Thus, focusing on the road network alone is not a horrible simplification. For simplicity, we will aggregate the trade data across sections and consider assignments with or without trade frictions and the 481 proposed new links. # Aggregating the trade data across HS sections africa_trade_agg collap(quantity ~ iso3_o + iso3_d, fsum) The Africa network is already in a graph format with from and to node columns. Thus, all we need to do here to get the representation flownet expects is drop the geometry. Next, we need to map the (port-)city locations to the nearest network nodes. # Convert to graph graph join(city_pop |> add_stub("_d", FALSE), multiple = TRUE) |> fmutate(flow = quantity * pop_share_o * pop_share_d) |> frename(from = node_o, to = node_d) |> fsubset(flow > 0 & from != to) ## left join: africa_trade_agg[iso3_o] 1971/1971 (100%) y[iso3_o] 452/453 (99.8%) ## left join: x[iso3_d] 19685/19685 (100%) y[iso3_d] 452/453 (99.8%) nrow(od_matrix_trade) ## [1] 189020 od_matrix_trade |> fselect(from:flow) |> head(3) ## from city_o pop_share_o to city_d pop_share_d flow ## ## 1: 577 Luanda - Angola 0.53616303 962 Bujumbura - Burundi 1 0.52919291 ## 2: 550 Lubango - Angola 0.03794550 962 Bujumbura - Burundi 1 0.03745221 ## 3: 540 Cabinda - Angola 0.03202289 962 Bujumbura - Burundi 1 0.03160660 We are now ready to assign traffic using the PSL model, which is straightforward using run_assignment(). To increase performance, I set nthreads = 8 to use 8 cores on this M4 Pro chip. R-level parallelism using mirai is highly efficient, but requires additional memory. In this example we also save all intermediate results, including individual path alternatives and associated probabilities, by setting return.extra = "all". This can get very memory-intensive, especially if there are many plausible paths per OD pair. Thus, we limit the number of paths per OD-pair to 150 by setting npaths.max = 150. Another crucial step for PSL modelling is cost-scaling. Multinomial logits, of which the PSL is a penalized variant, do not deal well with very large path-costs. Thus, mean or median scaling of the generalized cost (dividing all data points by the mean or median) is good practice. We will later also need to calibrate the PSL parameter beta (default 1) to achieve the desired level of penalization of overlapping paths. But first, we conduct a basic dry run of the model and explore the outputs. # Median scaling cost column: recommended for PSL settfm(graph, duration_ms = fmedian(duration, TRA = "/")) # Run Traffic Assignment using the PSL model system.time({ result st_transform(4326) |> linestrings_to_graph() nrow(graph) ## [1] 487394 graph[1:5, 1:10] ## edge from FX FY to TX TY .length ID LENGTH ## 1 1 1 -132.9750 56.80636 2 -132.9750 56.80701 86.86074 [m] 1 0.05408278 ## 2 2 2 -132.9750 56.80701 4 -132.9762 56.80849 184.49595 [m] 2 0.11485356 ## 3 3 2 -132.9750 56.80701 3 -132.9746 56.80664 52.08027 [m] 3 0.03243200 ## 4 4 3 -132.9746 56.80664 6 -132.9601 56.80941 1004.86210 [m] 4 0.62628710 ## 5 5 4 -132.9762 56.80849 10 -132.6638 56.71238 24432.32444 [m] 5 15.22300434 # 'edge' is always unique any_duplicated(graph$edge) ## [1] FALSE # But there may be duplicate edges between the same from-to nodes sum(fduplicated(graph[, c("from", "to")])) ## [1] 679 Evidently, casting to LINESTRING did not increase the row count in this case. linestrings_to_graph() then removed the geometry (SHAPE) column and added 8 columns: a unique edge identifier, from and to node IDs with corresponding geodetic coordinates (WGS84), and .length computed directly from the shape in meters (which here matches the LENGTH column expressed in miles). We can again extract nodes with nodes_from_graph() to determine the nearest node to each of the 3796 FAF5 zones provided as origin and destination locations for commodity flows: zones fsubset(!is.na(CentroidID)) |> st_transform(4326) ## Reading layer `FAF5_Nodes' from data source `/Users/sebastiankrantz/Documents/CPCS/FAF5/Networks/Geodatabase Format/FAF5Network.gdb' using driver `OpenFileGDB' ## Simple feature collection with 974788 features and 15 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -167.0152 ymin: 18.91474 xmax: -66.98005 ymax: 71.3198 ## Geodetic CRS: NAD83 nrow(zones) ## [1] 3796 head(zones, 1) ## Simple feature collection with 1 feature and 15 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -132.9393 ymin: 56.83877 xmax: -132.9393 ymax: 56.83877 ## Geodetic CRS: WGS 84 ## ID DATA Entry_or_Exit Exit_Number Interchange Centroid CentroidID Facility_Type Facility_Name County State StateID StateName FAFID StateNameBak ## 1 14 47179393 1 2195 2 Alaska 20 Alaska ## SHAPE ## 1 POINT (-132.9393 56.83877) nodes ftransform(id_o = as.integer(substr(dms_orig, 1, 3)), id_d = as.integer(substr(dms_dest, 1, 3))) nrow(OD) ## [1] 16798 head(OD, 3) ## dms_orig dms_dest dms_mode thousand_tons_in_2024 id_o id_d ## ## 1: 011-Birmingham AL 011-Birmingham AL 1-Truck 35382.9943 11 11 ## 2: 011-Birmingham AL 012-Mobile AL 1-Truck 970.8396 11 12 ## 3: 011-Birmingham AL 019-Rest of AL 1-Truck 10178.2525 11 19 This data can be matched to zone nodes, but only at the coarser 3-digit FAF zone level: zones |> join(fslice(OD, id_o), on = c("FAFID" = "id_o")) |> invisible() ## left join: zones[FAFID] 3796/3796 (100%) y[id_o] 132/132 (100%) zones |> join(fslice(OD, id_d), on = c("FAFID" = "id_d")) |> invisible() ## left join: zones[FAFID] 3796/3796 (100%) y[id_d] 132/132 (100%) We therefore aggregate the centroids to the 3-digit level and map them to the nearest network nodes: # Aggregating zones and computing aggregate centroids zones_agg fgroup_by(FAFID) |> fsummarise(SHAPE = st_centroid(st_combine(SHAPE))) # Matching each aggregate zone centroid to the nearest network node zones_agg$NN st_drop_geometry() |> fselect(FAFID, NN), on = c("id_o" = "FAFID")) |> frename(from = NN) |> join(zones_agg |> st_drop_geometry() |> fselect(FAFID, NN), on = c("id_d" = "FAFID")) |> frename(to = NN) |> fmutate(flow = thousand_tons_in_2024) |> fsubset(is.finite(flow) / flow > 0, from, dms_orig, to, dms_dest, flow) ## left join: OD[id_o] 16798/16798 (100%) y[FAFID] 132/132 (100%) ## left join: x[id_d] 16798/16798 (100%) y[FAFID] 132/132 (100%) head(OD_mapped) ## from dms_orig to dms_dest flow ## ## 1: 6551 011-Birmingham AL 6551 011-Birmingham AL 35382.9943 ## 2: 6551 011-Birmingham AL 124499 012-Mobile AL 970.8396 ## 3: 6551 011-Birmingham AL 125253 019-Rest of AL 10178.2525 ## 4: 6551 011-Birmingham AL 1691 020-Alaska 0.0000 ## 5: 6551 011-Birmingham AL 15659 041-Phoenix AZ 50.0572 ## 6: 6551 011-Birmingham AL 5981 042-Tucson AZ 1.3969 Let us now proceed to a basic flows assignment using the efficient All-or-Nothing (AoN) method. Before running it, we need to further process the graph in two ways: Create a column representing the generalized cost of traversing a link. Ensure directed links are properly represented. For (1), I combine the free flow time (in minutes) with a truck toll estimate (average per-mile toll for a 5-axle vehicle $\times$ segment length), converted to minutes assuming a value of time (VOT) of 50 USD/hour. fnobs(graph$TRUCKTOLL) # Number of links with tolls ## [1] 3245 graph %% fmutate(TollTime = replace_na(TRUCKTOLL, 0) / 50 * 60, cost_AB = replace_na(AB_FreeFlowTime + TollTime, 1e4), cost_BA = replace_na(BA_FreeFlowTime + TollTime, 1e4)) graph %$% fquantile(setv(TollTime / cost_AB, 0, NA), Qprobs) ## 0% 1% 5% 25% 50% 75% 95% 99% 100% ## 0.09420283 0.11402505 0.13810940 0.22898990 0.28057551 0.84988865 0.98935454 0.99612581 0.99964115 Among the 3245 links with tolls, tolls account for ~28% of the generalized cost. A few links (~3800 in each direction) are missing free flow time; these are imputed with a very high cost of 10,000 minutes. The network provides free flow time estimates in both directions, and, to complicate things further, a DIR column indicates whether the link is bidirectional (0) or unidirectional (1 or -1). We thus need to convert the graph to a fully directed one by duplicating bidirectional links with swapped endpoints. table(graph$DIR) # See direction or undirectedness of links ## ## -1 0 1 ## 391 168215 318788 # Create fully directed graph representation graph_dir fselect(from, to, FX, FY, TX, TY, cost = cost_AB), graph |> fsubset(DIR == 0, from = to, to = from, FX = TX, FY = TY, TX = FX, TY = FY, cost = cost_BA) ) nrow(graph_dir) ## [1] 655609 We are now ready for directed All-or-Nothing assignment: system.time({ result_AoN

R-bloggers

flownet is a new R package for high-performance transport modeling, implementing the Path-Sized Logit (PSL) model for stochastic traffic assignment with route overlap correction. It provides network processing tools including graph consolidation and simplification, route enumeration algorithms, and efficient All-or-Nothing assignment. Built on fastverse libraries (collapse, kit, igraph, geodist, mirai) with custom C code, it handles both small networks (Africa's 1,379-node road network) and large ones (US FAF5's 487,384 links). The package uses data frames for graphs and OD matrices, supports multimodal networks, and enables probabilistic traffic flow assignments with customizable parameters.

Introducing flownet: Efficient Transport Modeling in R

Example 1: Africa’s Continental Road Network

Example 2: The US Freight Analysis Framework (FAF) Network