Logjam

WORLD_LIMITS

A constant defining the geographical limits for a world map projection.

• The first tuple specifies the longitude limits in degrees: (-180, 180).
• The second tuple specifies the latitude limits in degrees: (-75, 75).

US_LIMITS

A constant defining the geographical limits for a map projection of the contiguous United States.

• The first tuple specifies the longitude limits in degrees: (-180, -65).
• The second tuple specifies the latitude limits in degrees: (15, 72).

CUS_LIMITS

A constant defining the geographical limits for a map projection of the contiguous United States (excluding Alaska and Hawaii).

• The first tuple specifies the longitude limits in degrees: (-125, -65).
• The second tuple specifies the latitude limits in degrees: (24, 50).

makemap(x::Union{Nothing, AbstractVector{<:Real}, NTuple{2, <:Real}} = nothing,
        y::Union{Nothing, AbstractVector{<:Real}, NTuple{2, <:Real}} = nothing;
        region::Symbol = :World, backend::Symbol = :CairoMakie,
        xexpand::Real = 0.3, yexpand::Real = 0.1, doRoadbkgd::Bool = true, maxroadlatspan::Real = 30.0) -> Figure, GeoAxis, Vector, Tuple

Creates map visualization for predefined or user-defined region of interest.

The map can focus on different predefined regions (the world, U.S., or continental U.S.) or a user-defined region of interest that contains a set of longitude-latitude points. Provides Mercator projection of geographical features such as country borders, U.S. state borders, and roads using GeoMakie.

Arguments

• x: Optional set of at least two longitudes to define the region of interest; if nothing, defaults to a pre-defined region based on the region parameter.
• y: Optional set of at least two latitudes to define the region of interest; if nothing, defaults to a pre-defined region based on the region parameter.
• region::Symbol: Specifies the region to focus on. Options include:
  • :World: Default. Focuses on the entire world.
  • :US: Focuses on the United States.
  • :CUS: Focuses on the continental U.S. without showing country borders.
• backend::Symbol: Specifies the rendering backend. Options are:
  • :CairoMakie: Default. Uses CairoMakie for rendering.
  • :GLMakie: Uses GLMakie for interactive rendering. Requires using GLMakie before calling.
• xexpand::Float64: Expansion factor for the x-axis limits. Default is 0.3.
• yexpand::Float64: Expansion factor for the y-axis limits. Default is 0.1.
• doRoadbkgd::Bool: Whether to include roads as background features if maximum latitude span is less than maxroadlatspan. Default is true.
• maxroadlatspan::Float64: Maximum latitude span for displaying roads. Default is 30.0° (allows continental US coverage with FAF5 interstate network).
• showgrid::Bool: Whether to display grid lines and coordinate labels. Default is false.

Returns

• fig::Figure: The figure object containing the map.
• ax::GeoAxis: The axis object where the map is drawn.
• hborders::Vector: A vector of handles for the lines plotted on the map in the following order:
  • hborders[1]: Interstate roads, if used (derived from FAF5: https://geodata.bts.gov/datasets/usdot::freight-analysis-framework-faf5-network-links/about).
  • hborders[2]: U.S. state borders, if used (derived from: https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json)).
  • hborders[3]: Country borders, if used (derived from https://github.com/PublicaMundi/MappingAPI/blob/master/data/geojson/countries.geojson?short_path=b27f2ec)).
• limits::Tuple: The geographic limits (bounding box) used for the map.

Behavior

• Automatically selects the appropriate region and borders based on the provided x, y, and region parameters.
• Chooses the rendering backend and activates it accordingly.
• Draws country borders, U.S. state borders, and FAF5 Interstate Highways depending on the specified options and region.
• If x and y are provided, calculates the bounding box with optional expansion and adjusts the map view accordingly. Expansion allows for better visualization around x and y points.

Examples

# Create a world map using CairoMakie
fig, ax, hborders, limits = makemap()

# Create a U.S. map
fig, ax, hborders, limits = makemap(region=:US)

# Create a map focused on a specific region with expanded limits
using GeoMakie   # Required for scatter! function
x = [-84.0, -83.0, -82.0]
y = [41.0, 42.0, 43.0]
fig, ax, hborders, limits = makemap(x, y)
scatter!(ax, x, y, markersize=12, color=:red)
fig

mapbbox(x::Union{AbstractVector{<:Real}, Tuple{Vararg{Real}}},
        y::Union{AbstractVector{<:Real}, Tuple{Vararg{Real}}};
        xexpand::Real=0.0, yexpand::Real=0.0) -> Tuple, Tuple

Calculates the bounding box for a set of geographic coordinates, with optional expansion along the x and y axes.

Arguments

• x: A vector or tuple of x-coordinates (longitude values).
• y: A vector or tuple of y-coordinates (latitude values).
• xexpand: A Float64 value (default = 0.0) specifying the fractional expansion of the bounding box along the x-axis. For example, xexpand=0.1 expands the bounding box by 10% on each side.
• yexpand: A Float64 value (default = 0.0) specifying the fractional expansion of the bounding box along the y-axis.

Returns

• A tuple of x-limits and y-limits after applying any expansions, in the form ((xmin, xmax), (ymin, ymax)).
• A tuple of the original x-limits and y-limits without any expansion.

Details

• The function first calculates the minimum and maximum values of x and y, ignoring any NaN values.
• It then applies the specified xexpand and yexpand to enlarge the bounding box.
• The x-limits are clamped to the range [-180, 180] to ensure valid longitude values.
• The y-limits are clamped to slightly above -90 and slightly below 90 to ensure valid latitude values and avoid issues with map projections.
• Throws ArgumentError if x or y contains no non-NaN values.

aligntext(x::Union{Real, AbstractVector{<:Real}, Tuple{Vararg{Real}}},
          y::Union{Real, AbstractVector{<:Real}, Tuple{Vararg{Real}}};
          offsetamt::Real=1, mindistratio::Real=1.5, idx=nothing) -> Pair, Pair

Determines text alignment and offset positions for given points.

This function attempts to calculate the best alignment and offset positions for text labels based on the spatial arrangement of the points provided. It is particularly useful for positioning labels or annotations on a plot, ensuring that they do not overlap and remain readable. The function can handle various input formats for the points, including scalars, vectors, and tuples, and adjusts the text position to try to avoid collisions with nearby labels or graphical elements.

Arguments

• x: Scalar, vector, or tuple representing the x-coordinates for the points.
• y: Scalar, vector, or tuple representing the y-coordinates for the points.
• offsetamt: Scalar value specifying the amount of offset to apply to the text labels. This controls the distance by which the text is shifted away from the point. Default is 1.
• mindistratio: Scalar value that sets the minimum distance ratio used to decide the best alignment for text labels relative to adjacent points. Default is 1.5.
• idx: Optional index or index vector. When provided, alignment is computed using all points but only results for the specified indices are returned. Useful for labeling a subset of points (e.g., a depot) while considering all points for placement.

Returns

• :align => alignout: A Pair where :align is associated with an array of 2-tuples representing horizontal and vertical alignment symbols (e.g., (:left, :bottom), (:center, :top)) corresponding to each point.
• :offset => offsetout: A Pair where :offset is associated with an array of 2-tuples representing the x and y offsets to be applied to the text labels for each point.

Behavior

• Single Point: If only one point is provided, the function returns a default alignment (:left, :bottom) with the specified offsetamt.
• Two Points: If two points are provided, the function calculates the angle between the points and determines the best alignment and offset in both directions.
• Three or More Points: For three or more points, the function uses Delaunay triangulation to determine the optimal alignment by analyzing the angles and distances between adjacent points. It ensures that labels do not overlap and are well-positioned relative to each other.

Example

# Example: Cities in North Carolina with populaions over 100,000
using GeoMakie, DataFrames
df = filter(r -> (r.STFIP == st2fips(:NC)) && (r.POP > 100_000), usplace())
x, y, name = df.LON, df.LAT, df.NAME
fig, ax = makemap(x, y)
scatter!(ax, x, y)
text!(ax, x, y, text=name; aligntext(x, y)...)  # Note ";" and "..." for splatting
display(fig)

isptinbbox(pt, bbox::Tuple{Union{Tuple{<:Real, <:Real}, AbstractVector{<:Real}},
                           Union{Tuple{<:Real, <:Real}, AbstractVector{<:Real}}}) -> Bool

Determines whether a given point lies within a specified bounding box.

Arguments

• pt: A tuple or vector of exactly two elements representing the coordinates of the point (x, y).
• bbox: A tuple of two tuples or arrays, each containing two elements representing the bounding box. The first tuple/array defines the x-limits (xmin, xmax) and the second tuple/array defines the y-limits (ymin, ymax).

Returns

• A Bool value:
  • true if the point pt lies within the bounding box bbox.
  • false otherwise.

Example

julia> bbox = ((0, 10), (0, 15));

julia> isptinbbox((5, 10), bbox)  # inside
true

julia> isptinbbox((15, 10), bbox)  # outside
false

alloclines(W, hub_xy, spoke_xy; tol=sqrt(eps())) -> (X, Y)

Convert allocation matrix to NaN-separated line segments for visualization.

Creates line segments connecting hubs to their allocated spokes. Returns one vector of coordinates per hub, enabling per-hub formatting (e.g., different colors).

Arguments

• W: n×m allocation matrix where W[i,j] indicates allocation weight from hub i to spoke j.
• hub_xy: n×2 matrix of hub coordinates [lon, lat] or [x, y].
• spoke_xy: m×2 matrix of spoke coordinates [lon, lat] or [x, y].
• tol: Threshold for nonzero allocation (default: √eps ≈ 1.5e-8).

Returns

• (X, Y): Tuple of Vector{Vector{Float64}}, each of length n (one per hub). X[i] and Y[i] contain NaN-separated coordinates for hub i's allocation lines.

Example

k = [100.0, 100.0, 150.0]
C = [0 3 7 10; 3 0 4 8; 7 4 0 5]
y, TC, W = ufl(k, C; verbose=false)

hubs = [-80.0 35.0; -78.0 36.0; -79.0 35.5]
spokes = [-80.5 35.2; -78.5 35.8; -79.2 36.1; -78.0 35.0]

X, Y = alloclines(W, hubs, spokes)

plotroads!(ax::GeoAxis, dfN::DataFrame, dfL::DataFrame;
           show_connectors::Bool = false) -> Vector{Lines}

Overlay road networks on GeoAxis with FCLASS-based styling inspired by OSM Carto.

Renders roads from DataFrames with differentiation by FCLASS (functional class). FAF5 and OSM roads at the same FCLASS receive identical visual treatment. Uses muted OSM Carto hues with adaptive zoom-based rendering: two-pass casing at close zoom, single-line with hue tints at medium zoom, minimal at wide zoom.

Arguments

• ax::GeoAxis: Geographic axis from makemap() or manual creation.
• dfN::DataFrame: Nodes with required IDX (col 1), LON (col 2), LAT (col 3).
• dfL::DataFrame: Links with required SRC (col 1), DST (col 2); optional SOURCE, FCLASS.
• show_connectors::Bool: Whether to render CONNECTOR links (default: false).

Returns

• Dict{Symbol, Any}: Named handles to plotted line objects. Keys: :fill_1 through :fill_5 (by FCLASS tier), :casing_1 through :casing_5 (close zoom only), :connector (if shown). Only non-empty tiers are included.

Styling

Roads are styled by FCLASS tier with zoom-adaptive rendering:

• latspan > 20° (CONUS-scale): Minimal single-line, faint (lw 0.25–0.05)
• 10° < latspan ≤ 20° (Regional): Single-line with subtle hue tints (lw 0.7–0.15)
• latspan ≤ 10° (City/metro): Casing + fill with muted OSM Carto hues (lw 1.0–0.3)

FCLASS tiers: 1=Interstate (blue), 2=Freeway (green), 3=Arterial (warm yellow), 4=Collector (pale yellow), 5+=Local (white/gray). Colors are muted to serve as background beneath overlaid data.

Data Requirements

• dfL must have at least 2 columns: SRC, DST (node IDs as integers)
• dfN must have at least 3 columns: IDX, LON, LAT (coordinates in WGS84)
• Optional dfL.SOURCE: "FAF5", "OSM", or "CONNECTOR" (missing treated as "FAF5")
• Optional dfL.FCLASS: Integer functional class (1-7, per FHWA/OSM mapping)

Examples

using GeoMakie

# Basic FAF5 plot
dfN, dfL = faf5nodes(), faf5links()
fig, ax = makemap(region=:CUS)
handles = plotroads!(ax, dfN, dfL)
display(fig)

# Customize interstate fill color
handles[:fill_1].color = :darkblue

# Check available handles
keys(handles)  # e.g., [:fill_1, :fill_2, :casing_1, :casing_2, ...]

Notes

• FCLASS-based styling applies uniformly to FAF5 and OSM roads
• Returns Dict{Symbol, Any} with keys :fill_N, :casing_N (N=1-5), :connector
• Casing keys only present at close zoom (latspan ≤ 10°)
• Connectors are hidden by default to avoid visual clutter from synthetic edges
• Node lookup uses Dict to handle non-sequential OSM node IDs efficiently
• Compatible with both CairoMakie and GLMakie backends

plotroute!(ax, lx, ly; kwargs...) → Vector

Base rendering method. Plots a route from coordinate vectors. Requires using CairoMakie, GeoMakie.

usplace() -> DataFrame

Returns DataFrame containing U.S. place data.

Geographic and population data for each place in the U.S., where each place is a city, town, or census-designated place (CDP). The latitude-longitude of each place represents a central location interior to the place and not its center of population. Does not include U.S. territories.

Columns

• STFIP: Integer representing state FIPS (Federal Information Processing Standards) code.
• PLFIP: Integer representing place FIPS code.
• NAME: String containing name of place (city, town, or CDP).
• ST: Symbol representing state abbreviation (e.g., :AL for Alabama).
• LAT: Float representing an interior latitude of place.
• LON: Float representing an interior longitude of place.
• POP: Integer representing population of place.
• ALAND: Float representing land area of place in square miles.
• AWATER: Float representing water area of place in square miles.
• LSAD: Integer representing legal/statistical area description code (e.g., 25 for a place).
• FUNCSTAT: String representing functional status of place (e.g., A for active).
• CBSA: Integer or None representing Core-Based Statistical Area code associated with place.
• ISCUS: Boolean indicating whether place is within continental U.S. (true or false).

Sources

Geographic data derived from [1], population data from [2], and CBSA from [3]. ICUS determined from LAT and LON.

1. U.S. Census Bureau, 2020 Gazetteer Files, 2020Gazplace_national.txt

2. U.S. Census Bureau, 2020 Census Demographic and Housing Characteristics File (DHC), DECENNIALDHC2020.P1

3. U.S. Census Bureau, Principal cities of metropolitan and micropolitan statistical areas, list2_2023.xls

uscounty() -> DataFrame

Returns DataFrame containing U.S. county-level data.

Geographic and population data for each U.S. county, including latitude-longitude coordinates representing the center of population of the county. Does not include U.S. territories.

Columns

• STFIP: Integer representing state FIPS (Federal Information Processing Standards) code.
• COFIP: Integer representing county FIPS code.
• NAME: String containing name of county.
• ST: Symbol representing state abbreviation (e.g., :AL for Alabama).
• LAT: Float representing latitude of county center of population.
• LON: Float representing longitude of county center of population.
• POP: Integer representing population of county.
• ALAND: Float representing land area of county in square miles.
• AWATER: Float representing water area of county in square miles.
• CBSA: Integer or None representing Core-Based Statistical Area code associated with county.

Sources

Area data from [1], population and center of population data from [2], and CBSA data from [3].

1. U.S. Census Bureau, 2020 Gazetteer Files, 2020Gazcounty_national.txt

2. U.S. Census Bureau, Centers of Population, CenPop2020MeanCO.txt

3. U.S. Census Bureau, Core based statistical areas (CBSAs), metropolitan divisions, and combined statistical areas (CSAs), list1_2023.xls

uscentract() -> DataFrame

Returns DataFrame containing U.S. census tract-level data.

Geographic and population data for each U.S. census tract, including latitude-longitude coordinates representing the center of population of the tract. Does not include U.S. territories.

Columns

• STFIP: Integer representing state FIPS (Federal Information Processing Standards) code.
• COFIP: Integer representing county FIPS code.
• TRFIP: Integer representing census tract FIPS code.
• ST: Symbol representing state abbreviation (e.g., :AL for Alabama).
• LAT: Float representing latitude of census tract center of population.
• LON: Float representing longitude of census tract center of population.
• POP: Integer representing population of census tract.
• ALAND: Float representing land area of census tract in square miles.
• AWATER: Float representing water area of census tract in square miles.
• ISCUS: Boolean indicating whether census tract is within continental U.S. (true or false).

Sources

Area data from [1]. Population and center of population data from [2].

1. U.S. Census Bureau, 2020 Gazetteer Files, 2020Gaztract_national.txt

2. U.S. Census Bureau, Centers of Population, CenPop2020MeanTR.txt

uscenblkgrp() -> DataFrame

Returns DataFrame containing U.S. census block group-level data.

Geographic and population data for each U.S. census block group, including latitude-longitude coordinates representing the center of population of the block group. Does not include U.S. territories.

Columns

• STFIP: Integer representing state FIPS (Federal Information Processing Standards) code.
• COFIP: Integer representing county FIPS code.
• TRFIP: Integer representing census tract FIPS code.
• BGFIP: Integer representing census block group FIPS code.
• LAT: Float representing latitude of block group center of population.
• LON: Float representing longitude of block group center of population.
• POP: Integer representing population of block group.
• ALAND: Float representing land area of block group in square miles.
• AWATER: Float representing water area of block group in square miles.

Sources

Area data from [1]. Population and center of population data from [2].

1. U.S. Census Bureau, TIGER/Line Shapefiles for 2020 Census Block Groups, [https://www2.census.gov/geo/tiger/TIGER2020/BG/]

2. U.S. Census Bureau, Centers of Population, CenPop2020MeanBG.txt

uszcta5() -> DataFrame

Returns DataFrame containing U.S. 5-digit ZIP Code Tabulation Area (ZCTA5) data.

Geographic and population data for each U.S. 5-digit ZIP Code Tabulation Area (ZCTA5). The latitude-longitude of each ZCTA5 represents a central location within the ZCTA5 and not its center of population. Does not include U.S. territories.

Columns

• ZCTA5: Integer (<= 5 digits) representing the ZIP Code Tabulation Area (ZCTA5) code.
• LAT: Float representing an interior latitude of ZCTA5.
• LON: Float representing an interior longitude of ZCTA5.
• POP: Integer representing population of ZCTA5.
• ALAND: Float representing land area of ZCTA5 in square miles.
• AWATER: Float representing water area of ZCTA5 in square miles.
• ISCUS: Boolean indicating whether ZCTA5 is within continental U.S. (true or false).

Sources

Geographic data derived from [1]. Population data derived from [2]. ISCUS determined from LAT and LON.

1. U.S. Census Bureau, 2023 Gazetteer Files, 2023Gazzcta_national.txt

2. U.S. Census Bureau, 2020 Census Demographic and Housing Characteristics File (DHC), DECENNIALDHC2020.P1

uszcta3() -> DataFrame

Returns DataFrame containing U.S. 3-digit ZIP Code Tabulation Area (ZCTA3) data.

Geographic and population data for each U.S. 3-digit ZIP Code Tabulation Area (ZCTA3). The latitude-longitude of each ZCTA3 represents its approximate center of population. ZCTA3s are approximated by aggregating ZCTA5s: summing population and areas, and approximating the center of population by calculating the population-weighted centroid of the ZCTA5s' interior locations. Does not include U.S. territories.

Columns

• ZCTA3: Integer (<= 3 digits) representing the 3-digit ZIP Code Tabulation Area (ZCTA3) code.
• LAT: Float representing latitude of ZCTA3 approximate center of population.
• LON: Float representing longitude of ZCTA3 approximate center of population.
• POP: Integer representing population of ZCTA3.
• ALAND: Float representing land area of ZCTA3 in square miles.
• AWATER: Float representing water area of ZCTA3 in square miles.
• ISCUS: Boolean indicating whether ZCTA3 is within continental U.S. (true or false).

Sources

All data derived from uszcta5.

uscbsa() -> DataFrame

Returns DataFrame containing U.S. Core-Based Statistical Area (CBSA) data.

Geographic and population data for each U.S. CBSA. The latitude-longitude of each CBSA represents its center of population. CBSA geographic and populations are determiend by aggregating its constituent counties: summing population and areas, and determining the center of population by calculating the population-weighted centroid of the county centers of population. Does not include U.S. territories.

Columns

• CBSA: Integer representing Core-Based Statistical Area code.
• NAME: String containing name of CBSA.
• LAT: Float representing latitude of CBSA center of population.
• LON: Float representing longitude of CBSA center of population.
• POP: Integer representing population of CBSA.
• ALAND: Float representing land area of CBSA in square miles.
• AWATER: Float representing water area of CBSA in square miles.
• M_MSA: String indicating whether CBSA is a Metropolitan Statistical Area or a Micropolitan Statistical Area.
• CSA: Integer or None representing Combined Statistical Area code if CBSA is part of a CSA.
• ISCUS: Boolean indicating whether CBSA is within continental U.S. (true or false).

Sources

CBSA delineations and classifications from [1]. Geographic and population data from uscounty().

1. U.S. Census Bureau, Core based statistical areas (CBSAs), metropolitan divisions, and combined statistical areas (CSAs), list1_2023.xls

uscsa() -> DataFrame

Returns DataFrame containing U.S. Combined Statistical Area (CSA) data.

Geographic and population data for each U.S. CSA. The latitude-longitude of each CSA represents its center of population. CSA geographic and populations are determined by aggregating its constituent CBSAs: summing population and areas, and determining the center of population by calculating the population-weighted centroid of the CBSA centers of population. Does not include U.S. territories.

Columns

• CSA: Integer representing Combined Statistical Area code.
• NAME: String containing name of CSA.
• LAT: Float representing latitude of CSA center of population.
• LON: Float representing longitude of CSA center of population.
• POP: Integer representing population of CSA.
• ALAND: Float representing land area of CSA in square miles.
• AWATER: Float representing water area of CSA in square miles.

Sources

CSA delineations and classifications from [1]. Geographic and population data from uscbsa().

1. U.S. Census Bureau, 2023 Combined Statistical Area (CSA) Codes, list2_2023.xls

st2fips(state::Symbol) -> Integer

Convert a two-character symbol for US states and territories to its corresponding FIPS code.

Valid two-character symbols of the state or territory are AK, AL, AR, AS, AZ, CA, CO, CT, DC, DE, FL, GA, GU, HI, IA, ID, IL, IN, KS, KY, LA, MA, MD, ME, MI, MN, MO, MP, MS, MT, NC, ND, NE, NH, NJ, NM, NV, NY, OH, OK, OR, PA, PR, RI, SC, SD, TN, TX, UT, VA, VI, VT, WA, WI, WV, WY

Examples

julia> st2fips(:NC)
37

julia> st2fips.([:NC, :NY])
2-element Vector{Int64}:
 37
 36

fips2st(fips::Integer) -> Symbol

Convert a FIPS code to its corresponding two-character state or territory symbol.

Valid FIPS codes are: 1, 2, 4, 5, 6, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49, 50, 51, 53, 54, 55, 56, 60, 66, 69, 72, 78

Arguments

• fips: An integer representing the FIPS code.

Returns

• A Symbol representing the two-character state or territory abbreviation.

Throws

• ArgumentError if the FIPS code is not valid.

Examples

julia> fips2st(37)
:NC

julia> fips2st(36)
:NY

julia> fips2st.(37:39)
3-element Vector{Symbol}:
 :NC
 :ND
 :OH

loc2lonlat(s::AbstractString; state=nothing, country=:US,
           cache_dir=".", force_download=false)

Geocode a single location string. Accepts addresses, city names, postal codes.

Returns a named tuple (lon, lat, source, uncert, status).

Examples

loc2lonlat("123 Main St, Raleigh, NC 27601")
loc2lonlat("Raleigh", state=:NC)
loc2lonlat("27601")

loc2lonlat(v::Vector{<:AbstractString}; state=nothing, country=:US,
           cache_dir=".", force_download=false)

Geocode a vector of location strings. Auto-detects postal codes vs place names.

Returns a DataFrame with LON, LAT, GCSOURCE, GCUNCERT, GC_STATUS columns.

Examples

loc2lonlat(["Raleigh", "Durham", "Chapel Hill"], state=:NC)
loc2lonlat(["27601", "27708", "27514"])

loc2lonlat(df::DataFrame; street=:STREET, city=:CITY, state=:STATE,
           postalcode=:POSTALCODE, county=:COUNTY, country=:US,
           cache_dir=".", force_download=false)

Geocode a DataFrame of addresses. Inspects which columns exist to determine resolution tier per row.

Returns the input DataFrame with LON, LAT, GCSOURCE, GCUNCERT, GC_STATUS appended.

Examples

df = DataFrame(CITY=["Raleigh", "Durham"], STATE=[:NC, :NC])
loc2lonlat(df)

df2 = DataFrame(addr=["123 Main St"], town=["Raleigh"], st=["NC"])
loc2lonlat(df2; street=:addr, city=:town, state=:st)

lonlat2loc(lon::Real, lat::Real, df::DataFrame; threshold=nothing)

Find nearest place to a (lon, lat) coordinate pair (reverse geocoding).

Returns a named tuple with fields: name, st, dist, bearing, dir, desc.

Example

cities = usplace()
lonlat2loc(-78.6382, 35.7796, cities)

lonlat2loc(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, df::DataFrame; threshold=nothing)

Find nearest places to coordinate vectors (reverse geocoding).

Returns DataFrame with columns: idx, name, st, dist, bearing, dir, desc.

Example

cities = usplace()
lonlat2loc([-78.6382, -80.8431], [35.7796, 35.2271], cities)

lonlat2loc(xy::AbstractVector, df::DataFrame; threshold=nothing)

Find nearest place to a [LON, LAT] coordinate (reverse geocoding).

Uses area-based "in" radius by default: if within √(ALAND/π) of centroid, reported as "in City"; otherwise "X mi DIR of City".

Returns a named tuple with fields: name, st, dist, bearing, dir, desc.

Example

cities = usplace()
lonlat2loc([-78.6382, 35.7796], cities)

lonlat2loc(XY::AbstractMatrix, df::DataFrame; threshold=nothing)

Find nearest places to n×2 matrix of [LON LAT] coordinates.

Returns DataFrame with columns: idx, name, st, dist, bearing, dir, desc.

lonlat2loc(df_in::DataFrame, df_ref::DataFrame;
           lon_col=:LON, lat_col=:LAT, threshold=nothing)

Reverse geocode coordinates from a DataFrame.

Returns the input DataFrame with location description columns appended.

faf5nodes() -> DataFrame

Returns DataFrame containing FAF5 road network nodes.

Physical road network nodes from the Freight Analysis Framework version 5 (FAF5). Centroid nodes (artificial FAF zone connectors) are excluded.

Columns

• IDX: Integer node identifier (non-sequential due to centroid removal).
• LON: Float representing longitude of node.
• LAT: Float representing latitude of node.
• ENTRY_EXIT: Entry/exit indicator for interchanges.
• EXIT_NUM: Exit number (if applicable).
• INTERCHANGE: Interchange name (if applicable).
• FACILITY_TYPE: Facility type code (port, airport, etc.).
• FACILITY_NAME: Facility name (if applicable).
• STATEID: Integer state FIPS code.
• FAFID: Integer FAF zone identifier.

Sources

U.S. Department of Transportation, Bureau of Transportation Statistics, Freight Analysis Framework (FAF5) Network, https://geodata.bts.gov/datasets/usdot::freight-analysis-framework-faf5-network-nodes/about

faf5links() -> DataFrame

Returns DataFrame containing FAF5 road network links.

Physical road network links from the Freight Analysis Framework version 5 (FAF5). Links connecting to centroid nodes are excluded.

Columns

• SRC: Integer source node ID.
• DST: Integer destination node ID.
• DIST: Float link distance in miles.
• STFIP: Integer state FIPS code.
• COFIP: Integer county FIPS code.
• SPEED: Integer posted speed limit (mph).
• NHS: Integer National Highway System code.
• SIGN: String signed route (e.g., "I 40", "US 1").
• NAME: String road name.
• FCLASS: Integer functional classification code.
• URBAN: String urban area code.
• AB_SPEED: Float final speed A→B direction (mph).
• BA_SPEED: Float final speed B→A direction (mph).
• AB_TIME: Float free-flow travel time A→B (minutes).
• BA_TIME: Float free-flow travel time B→A (minutes).
• AB_LANES: Integer lane count A→B direction.
• BA_LANES: Integer lane count B→A direction.
• DIR: Integer direction code.
• STRAHNET: Integer Strategic Highway Network code.
• NHFN: Integer National Highway Freight Network code.
• TOLL_TYPE: Integer toll classification.
• TOLL_LINK: Integer toll road indicator (0/1).
• BORDER_LINK: Integer border crossing indicator (0/1).
• FAFZONE: Integer FAF zone identifier.
• STATUS: Integer road status code.

Sources

U.S. Department of Transportation, Bureau of Transportation Statistics, Freight Analysis Framework (FAF5) Network, https://geodata.bts.gov/datasets/usdot::freight-analysis-framework-faf5-network-links/about

faf5interstate() -> Tuple{Vector{Float64}, Vector{Float64}}

Returns polyline vectors for FAF5 interstate road network.

Returns a tuple (x, y) of coordinate vectors with NaN separators between segments. This format is directly compatible with Makie's lines!() function for plotting road network backgrounds.

Usage

x, y = faf5interstate()
lines!(ax, x, y, color=(:steelblue, 0.3), linewidth=0.5)

Sources

Derived from FAF5 links where SIGN starts with "I " (interstate routes).

prt(X; rows, cols, title, fracdig, allfrac, row_title)

Pretty-print a matrix, vector, or DataFrame with smart formatting.

Formatted display with right-aligned columns. For matrices, applies smart per-column digit detection inspired by MATLOG's mdisp: integer columns show 0 decimal places, all-fractional columns (all values < 1) show allfrac decimal places, and mixed columns show fracdig decimal places. NaN values display as blank. Numbers ≥ 1,000 include comma separators. The row_title keyword places a header in the upper-left corner of the table (the row-label column header).

Matrix method

prt(X::AbstractMatrix; rows=1:size(X,1), cols=1:size(X,2), title="", fracdig=2, allfrac=4, row_title="", str=false)

Vector method

prt(v::AbstractVector; rows=1:length(v), title="", row_title="", str=false)

DataFrame method

prt(df::DataFrame; title="", str=false)

With str=true, returns the formatted table as a string (separator lines stripped) instead of printing. Default behavior prints to stdout and returns nothing.

Examples

julia> prt([1000 0.1234; 2000 0.5678])
       1        2
────────────────────
  1  1,000   0.1234
  2  2,000   0.5678

julia> prt([1 2; 3 4]; rows=["a","b"], cols=["X","Y"], row_title="ID")
  ID    X    Y
──────────────
   a    1    2
   b    3    4

dcf()

Display current figure. Requires using CairoMakie.

mat2df(X::AbstractMatrix, cols::AbstractVector; rows=nothing) -> DataFrame

Convert a matrix to a labeled DataFrame.

Arguments

• X: Matrix of values.
• cols: Column names (length must equal size(X, 2)).
• rows: Optional row labels. If provided (length must equal size(X, 1)), inserted as the first column with a blank header.

Example

julia> mat2df([1 2; 3 4], ["A", "B"])
2×2 DataFrame
 Row │ A      B
     │ Int64  Int64
─────┼──────────────
   1 │     1      2
   2 │     3      4

julia> mat2df([1 2; 3 4], ["A", "B"]; rows=["r1", "r2"])
2×3 DataFrame
 Row │         A      B
     │ String  Int64  Int64
─────┼──────────────────────
   1 │ r1          1      2
   2 │ r2          3      4

snapvals(x; atol=1e-8) -> Array
snapvals(x, n, m; atol=1e-8) -> Matrix

Snap near-integer and near-zero floating-point values in an array.

Designed for post-processing mathematical programming solutions where solvers introduce small numerical artifacts (e.g., 2.9999999997 instead of 3, or 1.2e-12 instead of 0). Each element is rounded to the nearest integer if within atol, and values that round to zero are set to exactly 0.0.

Caveat: Snapping may cause constraint violations in some models. When snapped values are used in further computation, verify feasibility — particularly for models with big-M constraints where rounding a near-zero binary variable to 0 can violate the associated constraint by a large amount.

Arguments

• x: Array of numeric values (typically from value.(x) after JuMP solve).
• n, m: Optional dimensions for reshaping the result.
• atol: Tolerance for snapping (default: 1e-8).

Example

julia> snapvals([2.9999999997, 1e-12, 0.5])
3-element Vector{Float64}:
 3.0
 0.0
 0.5

julia> snapvals([1.0, 2.0, 3.0, 4.0], 2, 2)
2×2 Matrix{Float64}:
 1.0  3.0
 2.0  4.0

dgc(xy₁, xy₂; unit=:mi) -> Float64

Calculate the great circle distance between two points.

Uses the haversine formula to compute the shortest distance over the Earth's surface between two points specified by longitude-latitude coordinates.

Arguments

• xy₁: Tuple or vector of (longitude, latitude) for the first point.
• xy₂: Tuple or vector of (longitude, latitude) for the second point.
• unit: Distance unit, either :mi (miles, default) or :km (kilometers).

Returns

• Great circle distance in the specified unit.

Example

# Distance from Raleigh to Charlotte
dgc((-78.6382, 35.7796), (-80.8431, 35.2271))  # ≈ 130 miles

d1(x₁, x₂) -> Float64

Calculate rectilinear (Manhattan, L₁) distance between two points.

Arguments

• x₁: First point (vector or tuple).
• x₂: Second point (same dimension as x₁).

Returns

• Sum of absolute differences: Σ|x₁ᵢ - x₂ᵢ|.

Example

d1([0, 0], [3, 4])  # Returns 7.0
d1([1, 2, 3], [4, 6, 2])  # Returns 8.0

d2(x₁, x₂) -> Float64

Calculate Euclidean (L₂) distance between two points.

Arguments

• x₁: First point (vector or tuple).
• x₂: Second point (same dimension as x₁).

Returns

• Euclidean distance: √(Σ(x₁ᵢ - x₂ᵢ)²).

Example

d2([0, 0], [3, 4])  # Returns 5.0
d2([1, 2, 3], [4, 6, 2])  # Returns 6.0

dists(X1, X2[, p]) -> Matrix{Float64}

Compute distance matrix between two point sets using specified metric.

Replaces: Dgc with unified interface supporting all metrics.

Arguments

• X1: m×n matrix of m points in n dimensions
• X2: k×n matrix of k points in n dimensions
• p: Distance metric (default: 2)
  • 1: Rectilinear (Manhattan) distance
  • 2: Euclidean distance (default)
  • :mi, :km, :rad: Great circle distance (requires n=2, lon-lat coordinates)

Returns

• D: m×k matrix where D[i,j] = distance from X1[i,:] to X2[j,:]

Examples

# Euclidean (default)
X1 = [0.0 0.0; 10.0 0.0]
X2 = [5.0 0.0; 5.0 5.0]
D = dists(X1, X2)  # 2×2 matrix

# Manhattan distance
D = dists(X1, X2, 1)

# Great circle distance (lon-lat coordinates)
cities = [-78.64 35.78; -122.42 37.77]  # Raleigh, SF
dc = [-77.04 38.91]                      # Washington DC
D = dists(cities, dc, :mi)               # Statute miles
D = dists(cities, dc, :km)               # Kilometers

See also: dgc, d1, d2

prune_reindex(dfN::DataFrame, dfL::DataFrame;
              src_col=1, dst_col=2, node_col=1) -> Tuple{DataFrame, DataFrame}

Prune network to common vertices and reindex node IDs sequentially.

Removes any nodes that don't appear in links and any links that reference non-existent nodes. Then reindexes all node IDs to be sequential starting from 1.

Arguments

• dfN: Nodes DataFrame with columns [IDX, ...].
• dfL: Links DataFrame with columns [SRC, DST, ...].
• src_col: Source-node column in dfL (index or symbol, default 1).
• dst_col: Destination-node column in dfL (index or symbol, default 2).
• node_col: Node-ID column in dfN (index or symbol, default 1).

Returns

• Tuple of (prunednodes, prunedlinks) DataFrames with sequential node IDs.

Example

dfN, dfL = faf5nodes(), faf5links()
# Filter to North Carolina
dfL_nc = filter(r -> r.STFIP == 37, dfL)
dfN_nc = filter(r -> r.STATEID == 37, dfN)
# Prune and reindex
dfN_nc, dfL_nc = prune_reindex(dfN_nc, dfL_nc)

thin(dfN::DataFrame, dfL::DataFrame;
     agg::Dict{String, Function} = Dict(),
     must_match::Vector{String} = String[],
     keep_index::Bool = false,
     verbose::Bool = false)

Remove degree-2 nodes from a road network, merging their incident links.

When a node has exactly two neighbors, it can often be removed by merging the two incident links into a single link. This simplifies the network while preserving connectivity and total distance.

Arguments

• dfL: Links DataFrame with columns [SRC, DST, DIST, ...].
• dfN: Nodes DataFrame with columns [IDX, LON, LAT, ...].
• agg: Dict mapping column names to aggregation functions (default: keep first value). Example: Dict("SPEED" => minimum, "NHS" => maximum). Note: DIST is always summed regardless of this parameter.
• must_match: Vector of column names that must match for merging to occur. If attributes differ, the degree-2 node is preserved. Example: ["FCLASS", "STFIP"].
• keep_index: If true, returns merge log tracking which original links were combined.
• verbose: If true, prints thinning statistics (nodes, links, distance reduction).

Returns

• If keep_index=false: (thinned_nodes, thinned_links) tuple.
• If keep_index=true: (thinned_nodes, thinned_links, merge_log) tuple where merge_log is a Dict mapping new link index to vector of original link indices.

Example

# Load and crop network for a small region
dfN, dfL = faf5nodes(), faf5links()
dfN, dfL = cropnetwork(dfN, dfL, [-79.5, -78.5], [35.5, 36.0])

# Simple thinning (distance only)
dfN_thin, dfL_thin = thin(dfN, dfL)

# Conservative: require functional class to match
dfN_thin, dfL_thin = thin(dfN, dfL; must_match=["FCLASS"])

# Custom aggregation with verbose output
dfN_thin, dfL_thin = thin(dfN, dfL;
    agg=Dict("SPEED" => minimum, "NHS" => maximum),
    must_match=["STFIP"], verbose=true)

addconnectors(dfN::DataFrame, dfL::DataFrame, x_prime::Vector, y_prime::Vector;
              circuity::Real=1.3, add_nf_nf::Bool=true,
              src_col=1, dst_col=2, dist_col=3,
              node_col=1, x_col=2, y_col=3) -> Tuple{DataFrame, DataFrame}

Add connector links from demand points to a road network.

Uses Delaunay triangulation to efficiently find nearby network nodes for each demand point, then creates connector links with estimated distances. When add_nf_nf is enabled, also creates direct connectors between neighboring demand points using a separate Delaunay triangulation of the demand points.

Demand points are assigned node IDs 1:n, and existing network nodes are shifted by n (i.e., old node 1 becomes n+1).

Arguments

• dfN: Nodes DataFrame with columns [IDX, LON, LAT, ...].
• dfL: Links DataFrame with columns [SRC, DST, DIST, ...].
• x_prime: Vector of demand point longitudes.
• y_prime: Vector of demand point latitudes.
• circuity: Circuity factor for connector distances (default 1.3).
• add_nf_nf: If true (default), add direct connectors between neighboring demand points.
• src_col, dst_col, dist_col: Link columns for source, destination, and distance.
• node_col, x_col, y_col: Node columns for node ID, longitude, and latitude.

Returns

• Tuple of (newnodes, newlinks) DataFrames with connectors added.
• All input columns are preserved. Connector rows populate SRC/DST/DIST (and set DIR/ONEWAY to 0 when present); additional columns are missing.

Example

dfN, dfL = faf5nodes(), faf5links()
# Define warehouse locations
warehouses_lon = [-78.9, -80.2, -79.5]
warehouses_lat = [35.8, 35.9, 36.1]
# Add connectors (includes NF-to-NF connections)
dfN_conn, dfL_conn = addconnectors(dfN, dfL, warehouses_lon, warehouses_lat)

links2graph(dfL::DataFrame; weight=3, ab_weight=nothing, ba_weight=nothing,
            dir_col=:DIR, oneway_val=1, reindex=:auto, return_map=false,
            sparse_ratio=10, max_id_threshold=500_000,
            src_col=1, dst_col=2) -> SimpleWeightedDiGraph or (SimpleWeightedDiGraph, Dict, Dict)

Convert a links DataFrame to a directed weighted graph, with optional auto-reindexing.

This method can reindex sparse or non-sequential node IDs to avoid creating very large graphs with many unused vertices.

Note: When reindex=:auto triggers reindexing, a warning is emitted: "links2graph reindexed node IDs. Pass returnmap=true to get the ID mapping." Use `returnmap=trueto retrieve the old↔new ID mappings, orreindex=false` to suppress reindexing (requires integer IDs).

Arguments

• dfL: Links DataFrame with columns [SRC, DST, weight_col, ...].
• weight, ab_weight, ba_weight, dir_col, oneway_val: Same as core graph conversion.
• reindex: :auto (default), true, or false.
  • :auto reindexes if node IDs are non-integer, or if maximum(node_id) > max_id_threshold, or maximum(node_id) > sparse_ratio * n_unique_ids.
  • true always reindexes.
  • false never reindexes (requires integer node IDs).
• return_map: If true, returns (graph, id_map, inv_map).
  • id_map: Dict{ID,Int} mapping old IDs to new IDs.
  • inv_map: Dict{Int,ID} mapping new IDs to old IDs.
  • When reindexing is not performed, these maps are identity maps.
• sparse_ratio: Threshold for :auto reindexing. Default 10.
• max_id_threshold: Absolute max node ID threshold for :auto reindexing. Default 500_000.
• src_col, dst_col: Source and destination columns (index or symbol). Defaults to 1 and 2.

Returns

• SimpleWeightedDiGraph, or a tuple with mapping dictionaries if return_map=true.

x2ln(g, x) -> Vector{Float64}

Convert graph edges to NaN-separated line coordinates for plotting.

Given a graph g and coordinate vector x, produces a vector suitable for plotting all edges as lines with NaN separators.

Arguments

• g: A Graphs.jl graph object.
• x: Vector of coordinates (e.g., longitudes or latitudes) indexed by vertex.

Returns

• Vector with pattern [x[src], x[dst], NaN, x[src], x[dst], NaN, ...].

Example

using Graphs, CairoMakie
g = SimpleGraph(3)
add_edge!(g, 1, 2)
add_edge!(g, 2, 3)
lon = [-78.0, -79.0, -80.0]
lat = [35.0, 37.0, 36.0]
x_line = x2ln(g, lon)
y_line = x2ln(g, lat)
fig, ax = makemap(x_line, y_line)
lines!(ax, x_line, y_line)
fig

cropnetwork(nodes::DataFrame, links::DataFrame, x, y;
            xexpand=0.1, yexpand=0.1) -> Tuple{DataFrame, DataFrame}

Crop a road network to a bounding box defined by coordinates.

Filters nodes to those within the expanded bounding box of the given coordinates, then filters links to those connecting remaining nodes, and reindexes.

Arguments

• nodes: Nodes DataFrame with columns [IDX, LON, LAT, ...].
• links: Links DataFrame with columns [SRC, DST, ...].
• x: Vector of longitudes defining the region of interest.
• y: Vector of latitudes defining the region of interest.
• xexpand: Expansion factor for longitude bounds (default 0.1).
• yexpand: Expansion factor for latitude bounds (default 0.1).

Returns

• Tuple of (croppednodes, croppedlinks) with sequential node IDs.

Example

nodes, links = faf5nodes(), faf5links()
cities = filter(r -> r.ST == :NC && r.POP > 100_000, usplace())
nodes_nc, links_nc = cropnetwork(nodes, links, cities.LON, cities.LAT)

shortestpaths(g, weights, n::Int) -> Tuple{Matrix{Float64}, Vector{Vector{Int}}}

Compute shortest path distances and parent pointers from the first n nodes.

Uses Dijkstra's algorithm to compute shortest paths from each of the first n nodes to all other nodes in the graph. Returns both the distance matrix and parent vectors for path reconstruction.

Arguments

• g: A Graphs.jl graph object.
• weights: Edge weight matrix (sparse or dense) where weights[i,j] is the cost from i to j.
• n: Number of source nodes to compute paths from (typically the number of demand points).

Returns

• D: n×n distance matrix where D[i,j] is the shortest distance from node i to node j.
• P: Vector of parent vectors for path reconstruction. P[i][j] gives the predecessor of node j on the shortest path from node i.

Example

using Graphs, SparseArrays

# Load network and add 3 demand points
nodes, links = faf5nodes(), faf5links()
nodes, links = cropnetwork(nodes, links, [-79.0, -78.0], [35.5, 36.0])
locs_lon, locs_lat = [-78.5, -78.7, -78.3], [35.7, 35.8, 35.6]
nodes, links = addconnectors(nodes, links, locs_lon, locs_lat)

# Build graph and explicit weight matrix
g = links2graph(links)
n = Graphs.nv(g)
weights = sparse(links.SRC, links.DST, links.DIST, n, n)
weights = weights + weights'  # Make symmetric

# Compute paths from the 3 demand points
D, P = shortestpaths(g, weights, 3)

Notes

• The first n nodes are typically demand points added via addconnectors.
• Parent vectors can be used with rte2lines to reconstruct physical paths.

shortestpaths(g::SimpleWeightedDiGraph, n::Int) -> Tuple{Matrix{Float64}, Vector{Vector{Int}}}

Compute shortest path distances and parent pointers from the first n nodes.

This method works with SimpleWeightedDiGraph from links2graph, using the graph's internal edge weights directly.

Arguments

• g: A SimpleWeightedDiGraph from links2graph.
• n: Number of source nodes to compute paths from (typically the number of demand points).

Returns

• D: n×n distance matrix where D[i,j] is the shortest distance from node i to node j.
• P: Vector of parent vectors for path reconstruction. P[i][j] gives the predecessor of node j on the shortest path from node i.

Example

# Load network and add demand points
nodes, links = faf5nodes(), faf5links()
nodes, links = cropnetwork(nodes, links, [-79.0, -78.0], [35.5, 36.0])
locs_lon, locs_lat = [-78.5, -78.7, -78.3], [35.7, 35.8, 35.6]
nodes, links = addconnectors(nodes, links, locs_lon, locs_lat)

# Build weighted graph and compute paths in two lines
g = links2graph(links)
D, P = shortestpaths(g, 3)  # Paths from the 3 demand points

Notes

• The first n nodes are typically demand points added via addconnectors.
• Parent vectors can be used with rte2lines to reconstruct physical paths.
• This method uses the weights embedded in the graph, no separate matrix needed.

tracepath(parents::Vector{Int}, origin::Int, dest::Int) → Vector{Int}

Reconstruct a shortest-path node sequence from a Dijkstra parent vector.

Walk parents backwards from dest to origin and return the forward-ordered node sequence. Throws an ArgumentError if dest is unreachable from origin.

Example

using Graphs

# Build a small network graph
dfN, dfL = faf5nodes(), faf5links()
dfN, dfL = cropnetwork(dfN, dfL, [-79.5, -78.5], [35.5, 36.0])
g = links2graph(dfL)
ds = Graphs.dijkstra_shortest_paths(g, 1)
path = tracepath(ds.parents, 1, 5)

osm_roads(bbox; cache_dir=".", force_download=false) → (dfN, dfL)

Download OpenStreetMap drivable roads for a bounding box and return Logjam-compatible DataFrames. Requires LightOSM and NearestNeighbors to be loaded (using LightOSM, NearestNeighbors).

Arguments

• bbox: Tuple (xmin, xmax, ymin, ymax) where x = longitude, y = latitude, as returned by mapbbox.

Keywords

• cache_dir::String=".": Directory for CSV cache files.
• force_download::Bool=false: If true, bypass cache and re-download from Overpass.

Returns

• dfN::DataFrame: Nodes with columns IDX (Int), LON (Float64), LAT (Float64), SOURCE (String).
• dfL::DataFrame: Links with columns SRC (Int), DST (Int), DIST (Float64, miles), SPEED (Int, mph), FCLASS (Int), DIR (Int), NAME (String), SOURCE (String).

Speed defaults by OSM highway class (mph)

Speeds are typical US defaults used when OSM maxspeed tags are missing (~70% of US roads lack maxspeed in OSM):

Notes

• Results are cached as CSV files in cache_dir. Subsequent calls with the same bounding box skip the Overpass download.
• Large bounding boxes (>~200 sq mi) emit a warning; the Overpass API may be slow or fail for very large regions.
• Coordinates are stored as (LON, LAT) per Logjam convention.

stitchnetworks(dfN_base, dfL_base, dfN_osm, dfL_osm; tolerance_m=2000, fclass_max=99) → (dfN, dfL)

Stitch an OSM regional network to an existing combined network (FAF5 or FAF5 + prior OSM regions) by creating connector edges between nearby nodes.

Arguments

• dfN_base::DataFrame: Nodes of the existing (base) network.
• dfL_base::DataFrame: Links of the existing (base) network.
• dfN_osm::DataFrame: Nodes of the OSM region to attach.
• dfL_osm::DataFrame: Links of the OSM region to attach.

Keywords

• tolerance_m::Real=2000: Maximum connection distance in meters (~1.2 miles).
• fclass_max::Int=99: Maximum FCLASS on base-network links for eligible connection nodes (default 99 = all road classes).

Returns

• dfN::DataFrame: Combined nodes (base + OSM + reindexed).
• dfL::DataFrame: Combined links (base + OSM + connector edges).

Details

The function:

1. Identifies eligible base-network nodes (on links with FCLASS ≤ fclass_max).
2. Identifies candidate OSM nodes near the region boundary.
3. Builds a KDTree over eligible base nodes and finds matches within tolerance_m.
4. Filters for grade separation (OSM layer tags) and respects directionality.
5. Creates connector edges (SOURCE="CONNECTOR") with distance via great-circle calculation and speed set to the lower of the two connected roads.
6. Combines all DataFrames and runs prune_reindex for sequential node IDs.
7. Validates topology (warns if no cross-boundary path exists).

Requires LightOSM and NearestNeighbors to be loaded.

segcost(loc, C) -> Vector

Calculate the cost of each segment in a location sequence.

Returns a vector of costs for traveling between consecutive locations.

Arguments

• loc: Vector of location indices representing a sequence of stops.
• C: Cost matrix where C[i,j] is the cost from location i to location j.

Returns

• Vector of segment costs of length length(loc) - 1.

Example

C = [0 10 20; 10 0 15; 20 15 0]  # 3x3 cost matrix
loc = [1, 2, 3, 1]               # Tour visiting all nodes
segcost(loc, C)                  # Returns [10, 15, 20]

rteTC(rte, sh, C, tr=(b=[], e=[])) -> Float64

Calculate the total cost of a route.

Converts the route to a location sequence and sums all segment costs.

Arguments

• rte: Route as a vector of shipment indices (each appears twice: pickup and delivery).
• sh: Shipments table with columns b (begin/pickup) and e (end/delivery) node indices.
• C: Cost matrix where C[i,j] is the cost from location i to location j.
• tr: Named tuple with depot locations (b=start_nodes, e=end_nodes). Default empty.

Returns

• Total route cost as the sum of all segment costs.

Example

using DataFrames

# Define 6 locations as (lon, lat) pairs
locs = [-78.6 35.8; -80.8 35.2; -79.0 36.1; -77.5 35.5; -81.0 35.0; -78.0 36.0]

# Create cost matrix using great circle distances
C = dists(locs, locs, :mi)

# Define 3 shipments: pickup at locs 1,2,3 and deliver to locs 4,5,6
sh = DataFrame(b=[1, 2, 3], e=[4, 5, 6])

# Route visiting each shipment sequentially (pickup then deliver)
rte = [1, 1, 2, 2, 3, 3]
cost = rteTC(rte, sh, C)

isorigin(rte) -> Vector{Bool}

Determine which positions in a route are origins (first occurrence).

For routes with pickup-delivery pairs, each shipment index appears twice. This function identifies which occurrence is the pickup (origin).

Arguments

• rte: Route as a vector of shipment indices.

Returns

• Boolean vector where true indicates first occurrence (origin/pickup).

Example

rte = [1, 2, 1, 3, 2, 3]  # Pickup 1, pickup 2, deliver 1, pickup 3, deliver 2, deliver 3
isorigin(rte)             # Returns [true, true, false, true, false, false]

rte2loc(rte, sh, tr=(b=[], e=[])) -> Vector

Convert a route to a location sequence.

Maps shipment indices in a route to their actual node locations, using pickup location for origins and delivery location for destinations.

Arguments

• rte: Route as a vector of shipment indices.
• sh: Shipments table with columns b (begin/pickup) and e (end/delivery) node indices.
• tr: Named tuple with depot locations (b=start_nodes, e=end_nodes). Default empty.

Returns

• Vector of node indices representing the physical locations visited.

Example

using DataFrames
sh = DataFrame(b=[10, 20], e=[15, 25])  # 2 shipments
rte = [1, 2, 1, 2]                       # Pickup both, then deliver both
rte2loc(rte, sh)                         # Returns [10, 20, 15, 25]

mincostinsert(idx, rte, rteTCh) -> Tuple{Vector, Float64}

Insert a shipment into a route at the minimum cost position.

For pickup-delivery problems, both the pickup and delivery must be inserted. This function tries all valid insertion positions and returns the best.

Arguments

• idx: Shipment index to insert.
• rte: Current route as a vector of shipment indices.
• rteTCh: Function that calculates route cost: rteTCh(route) -> cost.

Returns

• Tuple of (newroute, totalcost).

Example

using DataFrames

# Define locations and cost matrix
locs = [-78.6 35.8; -80.8 35.2; -79.0 36.1; -77.5 35.5]
C = dists(locs, locs, :mi)

# Define 2 shipments
sh = DataFrame(b=[1, 3], e=[2, 4])
rteTCh = rte -> rteTC(rte, sh, C)

# Insert shipment 2 into a route containing only shipment 1
rte = [1, 1]
rte, cost = mincostinsert(2, rte, rteTCh)

Notes

• Returns (route, Inf) if no feasible insertion exists.
• The pickup must occur before the delivery in the route.

pairwisesavings(rteTCh, sh) -> Tuple{Vector{Int}, Vector{Int}, Vector{Float64}}

Calculate savings from combining pairs of single-shipment routes.

The savings value represents cost reduction from serving two shipments on one route instead of two separate routes.

Arguments

• rteTCh: Function that calculates route cost: rteTCh(route) -> cost.
• sh: Shipments DataFrame (used to determine number of shipments via nrow).

Returns

• Tuple of (iindices, jindices, savings_values), sorted by savings (descending).

Example

using DataFrames

# Define locations and cost matrix
locs = [-78.6 35.8; -80.8 35.2; -79.0 36.1; -77.5 35.5]
C = dists(locs, locs, :mi)

# Define 2 shipments
sh = DataFrame(b=[1, 3], e=[2, 4])
rteTCh = rte -> rteTC(rte, sh, C)

# Calculate pairwise savings (pairs with highest savings are first)
iˢ, jˢ, sˢ = pairwisesavings(rteTCh, sh)

Notes

• Only returns pairs with positive savings.
• Savings = cost([i,i]) + cost([j,j]) - cost(combined route).

savings(rteTCh, sh) -> Vector{Vector{Int}}

Construct routes for a Pickup and Delivery Problem (PDP) using a savings-ordered insertion heuristic.

Unlike the standard Clarke-Wright algorithm (which concatenates routes at endpoints), this function prioritizes merges based on savings values but executes them by inserting all shipments from a source route into optimal positions within a target route.

Algorithm

1. Initialize a separate route for each shipment [i, i] (pickup, delivery).
2. Compute pairwise savings to prioritize merge attempts.
3. Iterate through pairs with highest savings:
   • Attempt to merge routes by inserting shipments from the shorter route into optimal positions in the longer route using mincostinsert.
   • If insertion satisfies constraints (finite cost), finalize the merge.

Arguments

• rteTCh: Function that calculates route cost: rteTCh(route) -> cost. Should return Inf for infeasible routes (capacity, time windows, etc.).
• sh: Shipments DataFrame where each row is a shipment requiring pickup and delivery.

Returns

• Vector of routes, where each route is a vector of shipment indices.

Example

using DataFrames

# Define locations as (lon, lat) pairs
locs = [-78.6 35.8; -80.8 35.2; -79.0 36.1; -77.5 35.5]  # 4 locations

# Create cost matrix using great circle distances
C = dists(locs, locs, :mi)

# Define 2 shipments: pickup at loc 1 deliver to loc 2, pickup at loc 3 deliver to loc 4
sh = DataFrame(b=[1, 3], e=[2, 4])

# Build routes using savings heuristic
rteTCh = rte -> rteTC(rte, sh, C)
routes = savings(rteTCh, sh)

Notes

• This is a hybrid of savings-based prioritization and greedy insertion.
• Constraints are enforced through rteTCh returning Inf.
• Each shipment appears exactly twice in its route (pickup and delivery).

twoopt(r, rTCh) -> Tuple{Vector, Float64}

Improve a route using the 2-opt neighborhood search.

The 2-opt algorithm iteratively reverses segments of the route to reduce total cost. It continues until no improving moves are found.

Arguments

• r: Initial route as a vector (can be TSP tour or shipment sequence).
• rTCh: Function that calculates route cost: rTCh(route) -> cost.

Returns

• Tuple of (improvedroute, totalcost).

Example

C = [0 10 25 20; 10 0 15 30; 25 15 0 10; 20 30 10 0]
rTCh = r -> sum(C[r[i], r[i+1]] for i in 1:length(r)-1)
initial = [1, 2, 3, 4, 1]
improved, cost = twoopt(initial, rTCh)

Notes

• For TSP: route should start and end at same node (depot).
• For pickup-delivery: ensure feasibility is maintained by the cost function.

ufladd(k, C; y=Int[], p=nothing) -> (Vector{Int}, Float64, SparseMatrixCSC)

Greedy ADD construction heuristic for uncapacitated facility location.

Iteratively adds facilities that provide the greatest cost reduction until no improvement is possible or p facilities are selected.

Arguments

• k: Fixed costs. Scalar (same cost for all) or vector (one per site).
• C: n×m cost matrix where C[i,j] is cost of serving customer j from facility i.
• y: Initial facility set (default: empty, start from scratch).
• p: Maximum facilities to select (default: nothing, no limit).

Returns

• (y, TC, W): Selected facility indices, total cost, and n×m sparse allocation matrix where W[i,j]=1 if facility i serves customer j.

Example

k = [10, 10, 10]
C = [2 5 4; 4 1 3; 5 4 2]
y, TC, W = ufladd(k, C)

References

• M.G. Kay, Facility Location (course notes), NC State University

ufldrop(k, C; y=nothing, p=nothing) -> (Vector{Int}, Float64, SparseMatrixCSC)

Greedy DROP construction heuristic for uncapacitated facility location.

Starts with all facilities open and iteratively drops the one providing the greatest cost reduction until no improvement or p facilities remain.

Arguments

• k: Fixed costs. Scalar (same cost for all) or vector (one per site).
• C: n×m cost matrix where C[i,j] is cost of serving customer j from facility i.
• y: Initial facility set (default: all facilities).
• p: Target number of facilities (default: nothing, drop until no improvement).

Returns

• (y, TC, W): Selected facility indices, total cost, and n×m sparse allocation matrix where W[i,j]=1 if facility i serves customer j.

Example

k = [10, 10, 10]
C = [2 5 4; 4 1 3; 5 4 2]
y, TC, W = ufldrop(k, C)

References

• M.G. Kay, Facility Location (course notes), NC State University

uflxchg(k, C, y) -> (Vector{Int}, Float64, SparseMatrixCSC)

Pairwise EXCHANGE improvement heuristic for uncapacitated facility location.

Performs steepest descent pairwise swaps (close one facility, open another) until local optimum is reached.

Arguments

• k: Fixed costs. Scalar (same cost for all) or vector (one per site).
• C: n×m cost matrix where C[i,j] is cost of serving customer j from facility i.
• y: Initial facility set.

Returns

• (y, TC, W): Improved facility indices, total cost, and n×m sparse allocation matrix where W[i,j]=1 if facility i serves customer j.

Example

k = [10, 10, 10]
C = [2 5 4; 4 1 3; 5 4 2]
y, TC, W = uflxchg(k, C, [1, 3])

References

• M.G. Kay, Facility Location (course notes), NC State University

ufl(k, C; verbose=true) -> (Vector{Int}, Float64, SparseMatrixCSC)

Hybrid UFL heuristic combining ADD, DROP, and EXCHANGE procedures.

Iterates through ADD → EXCHANGE → (ADD vs DROP, take best) until no improvement. Typically produces high-quality solutions for uncapacitated facility location problems.

Arguments

• k: Fixed costs. Scalar (same cost for all) or vector (one per site).
• C: n×m cost matrix where C[i,j] is cost of serving customer j from facility i.
• verbose: Print iteration costs (default: true).

Returns

• (y, TC, W): Best facility indices, total cost, and n×m sparse allocation matrix where W[i,j]=1 if facility i serves customer j.

Example

k = [10, 10, 15]
C = [0 3 7; 3 0 4; 7 4 0]
y, TC, W = ufl(k, C)  # Prints: Add: 13.0, Xchg: 13.0

References

• M.G. Kay, Facility Location (course notes), NC State University

pmedian(p, C; verbose=true) -> (Vector{Int}, Float64, SparseMatrixCSC)

p-median facility location (fixed number of facilities, no fixed costs).

Selects exactly p facilities to minimize total transportation cost. Uses ufladd with k=0 to select p facilities, then uflxchg to improve the solution.

Arguments

• p: Number of facilities to select.
• C: n×m cost matrix where C[i,j] is cost of serving customer j from facility i.
• verbose: Print iteration costs (default: true).

Returns

• (y, TC, W): Selected facility indices, total cost, and n×m sparse allocation matrix where W[i,j]=1 if facility i serves customer j.

Example

C = [0 3 7; 3 0 4; 7 4 0]
y, TC, W = pmedian(2, C; verbose=false)

References

• M.G. Kay, Facility Location (course notes), NC State University

randX(P::AbstractMatrix, n::Int=1) -> Matrix

Generate n random points within the bounding box of point set P.

Useful for generating random facility locations for optimization problems.

Arguments

• P: m×d matrix of m points in d dimensions.
• n: Number of random points to generate (default: 1).

Returns

• n×d matrix of random points, uniformly distributed within min/max of each dimension.

Example

P = [0 0; 2 0; 2 3]
X = randX(P, 5)  # 5 random points in [0,2] × [0,3]

rate_ltl(q, s, d; ppi=104.2) -> Float64 or Vector{Float64}

Estimate LTL (less-than-truckload) transportation rate.

Uses nonlinear regression model fitted to industry tariff tables (CzarLite) covering 55,800 O-D pairs across the continental US. Returns Inf for out-of-bounds inputs.

Arguments

• q: Shipment weight (tons). Clamped to minimum 0.075 tons (150 lb).
• s: Shipment density (lb/ft³).
• d: Shipment distance (miles). Clamp to minimum 37 miles.
• ppi: LTL Producer Price Index (default: 104.2, the 2004 baseline).

Returns

• Rate in /ton-mile. Returns Inf if q > 5, d > 3354, or 2000q/s > 650.

Formula

r_LTL = PPI_LTL × [(s²/8 + 14) / ((q^(1/7) × d^(15/29) - 7/2) × (s² + 2s + 14))]

Valid Ranges

• q ∈ [0.075, 5] tons (150 lb to 10,000 lb)
• d ∈ [37, 3354] miles
• 2000q/s ≤ 650 ft³ (cube-out limit)

Example

rate_ltl(0.5, 8.0, 250.0)  # 0.5 ton, 8 lb/ft³, 250 miles → ~0.12 $/ton-mi
rate_ltl([0.5, 1.0], [8.0, 10.0], [250.0, 500.0])  # Vectorized

References

• M.G. Kay (2023), Freight Transport (course notes), Section 1.5.2, Eq (1.5), NC State University
• Kay & Warsing (2009), "Estimating LTL rates using publicly available empirical data," Int. J. Logistics Research and Applications, 12(3):165–193 doi:10.1080/13675560802392415

charge_tl(q, d, s; r=2.00, Kwt=25.0, Kcu=2750.0, ppi=102.7) -> Float64

Calculate TL transport charge (combining distance cost and minimum charge).

Arguments

• q: Shipment weight (tons).
• d: Distance (miles).
• s: Shipment density (lb/ft³).
• r: TL rate per loaded truck-mile (/mi). Default: 2.00.
• Kwt: Truck weight capacity (tons). Default: 25.0.
• Kcu: Truck effective cube capacity (ft³). Default: 2750.0.
• ppi: TL Producer Price Index. Default: 102.7.

Returns

• Transport charge in dollars.

Formula

c_TL = ⌈q/q_max⌉ × max(r × d, MC_TL)
where q_max = min(Kwt, s × Kcu / 2000)

Example

charge_tl(10.0, 500.0, 8.0)  # 10 tons, 500 mi, 8 lb/ft³ → ~1000
charge_tl(25.0, 500.0, 8.0)  # Full truckload → ~1000
charge_tl(30.0, 500.0, 8.0)  # Needs 2 trucks → ~2000

References

• M.G. Kay (2023), Freight Transport (course notes), Section 1.5.1, Eq (1.1-1.2), NC State University

charge_ltl(q, d, s; ppi=104.2) -> Float64

Calculate LTL transport charge (rate × weight × distance, or minimum charge).

Arguments

• q: Shipment weight (tons).
• d: Distance (miles).
• s: Shipment density (lb/ft³).
• ppi: LTL Producer Price Index. Default: 104.2.

Returns

• Transport charge in dollars.

Formula

c_LTL = max(r_LTL × q × d, MC_LTL)

Example

charge_ltl(0.5, 250.0, 8.0)  # 0.5 ton, 250 mi → ~15-20
charge_ltl(2.0, 500.0, 10.0) # 2 tons, 500 mi → ~120-150

References

• M.G. Kay (2023), Freight Transport (course notes), Section 1.5.2 and 1.5.3, Eq (1.5) and (1.8), NC State University

mincharge_tl(r::Real=2.00; ppi=102.7) -> Float64

Calculate TL minimum charge (independent of distance).

Represents fixed costs of loading/unloading at origin and destination. The constant 45 is an empirically derived value reflecting terminal handling costs that do not vary with distance, fitted from industry tariff data using methods similar to Kay & Warsing (2009).

Arguments

• r: TL revenue per loaded truck-mile (/mi). Default: 2.00 (2004 baseline).
• ppi: TL Producer Price Index (default: 102.7, 2004 baseline).

Returns

• Minimum charge in dollars.

Formula

MC_TL = (r/2) × 45  or  MC_TL = (ppi/102.7) × 45

Example

mincharge_tl()  # Uses default r=2.00 → 45.0
mincharge_tl(2.11; ppi=108.6)  # 2005 rates → 47.6

References

• M.G. Kay (2023), Freight Transport (course notes), Section 1.5.3, Eq (1.7), NC State University. Empirical derivation methodology: Kay & Warsing (2009), Int. J. Logistics Research and Applications, 12(3):165–193

mincharge_ltl(d::Real; ppi=104.2) -> Float64

Calculate LTL minimum charge.

The LTL minimum charge includes both a fixed terminal handling cost (45) and a distance-dependent component that accounts for the shipment being loaded/unloaded at multiple terminals along its route. The parameters (28/19 exponent, 1625 divisor) were empirically fitted from industry tariff data using nonlinear regression methods similar to those described in Kay & Warsing (2009).

Arguments

• d: Distance (miles). Must be > 0 and ≤ 3354.
• ppi: LTL Producer Price Index (default: 104.2, 2004 baseline).

Returns

• Minimum charge in dollars.

Formula

MC_LTL = (ppi/104.2) × [45 + (d^(28/19) / 1625)]

where 28/19 ≈ 1.474 (economies of scale exponent) and 1625 is the scaling constant.

Example

mincharge_ltl(0.0)    # 0.0 (no shipment)
mincharge_ltl(250.0)  # 250 miles → ~67.3
mincharge_ltl(500.0)  # 500 miles → ~90.8

References

• M.G. Kay (2023), Freight Transport (course notes), Section 1.5.3, Eq (1.8), NC State University. Empirical derivation methodology: Kay & Warsing (2009), Int. J. Logistics Research and Applications, 12(3):165–193

maxpayld(s, Kwt, Kcu) -> Float64 or Vector{Float64}

Determine maximum payload limited by weight or cube capacity.

Arguments

• s: Shipment density (lb/ft³).
• Kwt: Truck weight capacity (tons). Typical: 25.0.
• Kcu: Truck cube capacity (ft³). Typical: 2750.0.

Returns

• Maximum payload in tons: min(Kwt, s·Kcu/2000).

Example

maxpayld(8.0, 25.0, 2750.0)  # 8 lb/ft³ → 11.0 tons (cube-limited)
maxpayld(15.0, 25.0, 2750.0) # 15 lb/ft³ → 20.6 tons (cube-limited)
maxpayld(25.0, 25.0, 2750.0) # 25 lb/ft³ → 25.0 tons (weight-limited)

References

• M.G. Kay (2023), Freight Transport (course notes), Section 1.5.1, Eq (1.2), NC State University

totlogcost(q, c, f, a, v, h) -> Float64
totlogcost(q, c, params) -> Float64

Calculate total logistics cost (transport + inventory).

TLC = TC + IC, where:

• TC (transport cost) = c·f/q (cost per shipment × annual shipments)
• IC (inventory cost) = q·a·v·h (avg inventory × value × holding rate)

Arguments

• q: Shipment size (tons).
• c: Transport charge per shipment ($).
• f: Annual demand (tons/year).
• a: Inventory fraction (0-D correction, typically 0.5).
• v: Product value (/ton).
• h: Annual holding cost rate (fraction, e.g., 0.25 = 25%/year).
• params: Any object with fields .f, .a, .v, .h (e.g., NamedTuple from aggshmt or a DataFrameRow).

Returns

• Total logistics cost (/year).

Example

# Scalar form
totlogcost(5.0, 450.0, 100.0, 0.5, 1000.0, 0.25)
# Returns: TC = 9000, IC = 625, TLC = 9625

# Using aggshmt output
using DataFrames
products = DataFrame(f=[100, 200], s=[8, 10], v=[1000, 1500], h=[0.25, 0.25], a=[0.5, 0.5])
agg = aggshmt(products)
q, d = 2.0, 300.0
c = charge_ltl(q, d, agg.s)
totlogcost(q, c, agg)

aggshmt(df::DataFrame) -> NamedTuple

Aggregate multiple shipments into a single equivalent shipment.

Weighted harmonic mean for density, weighted arithmetic mean for value/holding cost/inventory fraction. Weights are annual demands (f) or shipment sizes (q).

Arguments

• df: DataFrame with columns :f (demand), :s (density), :v (value), :h (holding rate), :a (inventory fraction). If :f missing, uses :q.

Returns

• NamedTuple with aggregated fields: (f=..., s=..., v=..., h=..., a=...).

Example

using DataFrames
df = DataFrame(
    f = [100, 200, 150],
    s = [8, 10, 6],
    v = [1000, 1500, 800],
    h = [0.25, 0.25, 0.25],
    a = [0.5, 0.5, 0.5]
)
agg = aggshmt(df)
# agg.f = 450, agg.s ≈ 8.18 (harmonic mean), agg.v ≈ 1200 (weighted avg)

transport_costs(shipments::DataFrame; mode=:auto, kwargs...) -> DataFrame

Calculate transport costs for a collection of shipments, automatically selecting TL vs LTL.

Arguments

• shipments: DataFrame with columns :weight (tons), :density (lb/ft³), :distance (miles).
• mode: :auto (select min cost), :tl, or :ltl.
• kwargs: Additional parameters passed to charge functions (e.g., ppi, Kwt, Kcu, r).

Returns

• DataFrame with original columns plus :cost ($), `:mode` (`:tl` or `:ltl`), `:rate` ($/ton-mi).

Example

using DataFrames
shipments = DataFrame(
    weight = [0.5, 2.0, 15.0, 30.0],
    density = [8.0, 10.0, 12.0, 15.0],
    distance = [250.0, 500.0, 800.0, 1200.0]
)
results = transport_costs(shipments)
# Automatically selects LTL for small shipments, TL for large