9  OD pairs and desire lines

Desire lines are a way to represent the flow of people or goods between two points. They are often used in transport planning to represent the flow of trips between zones.

There are many ways to create desire lines, connecting origins and destinations. For instance:

Examples of representation od desire lines: a) multiple origins and 1 destination; b) transport zones; c) od jittering

The stplanr package is a collection of functions for sustainable transport planning with R, and it is built on top of the sf package (Lovelace and Ellison 2018).

In this workshop, we will use the od package, a lightweight package with a few functions from stplanr, namely the ones to create desire lines from origin-destination (OD) pairs (Lovelace and Morgan 2024).

9.1 Desire lines with od_2_sf

To create desire lines, we need a dataset with OD pairs and other dataset with the corresponding transport zones (spatial data).

The TRIPSmode.Rds dataset includes origins, destinations and number of trips between municipalities.

Code
TRIPSmode = readRDS("data/TRIPSmode.Rds")

The Municipalities_geo.gpkg dataset includes the geometry of the transport zones.

Code
library(sf)
Municipalities_geo = st_read("data/Municipalities_geo.gpkg", quiet = TRUE) # supress mesage

Then, we need to load the od package. We will use the od_to_sf() function to create desire lines from OD pairs.

Code
# install.packages("od")
library(od)

TRIPSdlines = od_to_sf(TRIPSmode, z = Municipalities_geo) # z for zones

For this magic to work smoothly, the first two columns of the TRIPSmode dataset must be the origin and destination zones, and these zones need to correspond to the first column of the Municipalities_geo dataset (with an associated geometry).

See more options with the ?stplanr::od2line function.

Now we can visualize the desire lines using the mapview function.

Code
library(mapview)
mapview(TRIPSdlines, zcol = "Total")

As you can see, this is too much information to be able to understand the flows.

9.1.1 Filtering desire lines

Filter intrazonal trips.

Code
library(dplyr)

TRIPSdlines_inter = TRIPSdlines |> 
  filter(Origin != Destination) |> # remove intrazonal trips
  filter(Total > 5000) # remove noise

mapview(TRIPSdlines_inter, zcol = "Total", lwd = 5)

Filter trips with origin or destination not in Lisbon.

Code
TRIPSdl_noLX = TRIPSdlines_inter |> 
  filter(Origin != "Lisboa", Destination != "Lisboa")

mapview(TRIPSdl_noLX, zcol = "Total", lwd = 8) # larger line width

Try to replace the Total with other variables, such as Car, PTransit, and see the differences.

9.2 Oneway desire lines

Note that the od_to_sf() function creates bidirectional desire lines. This can be not the ideal for visualization / representation purposes, as you will have 2 lines overlaping.

The function od_oneway() aggregates oneway lines to produce bidirectional flows.

By default, it returns the sum of each numeric column for each bidirectional origin-destination pair.

Code
nrow(TRIPSdlines)
[1] 315
Code
TRIPSdlines_oneway = od_oneway(TRIPSdlines)
nrow(TRIPSdlines_oneway)
[1] 168

Note that for the last municipalities you have less combinations now. Nevertheless, all the possible combinations are represented.

Code
head(TRIPSdlines_oneway[,c(1,2)]) # just the first 2 columns
Simple feature collection with 6 features and 2 fields
Attribute-geometry relationships: identity (2) (with 3 geometries empty)
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -9.229502 ymin: 38.63562 xmax: -8.91588 ymax: 38.75981
Geodetic CRS:  WGS 84
          o         d                       geometry
1 Alcochete Alcochete               LINESTRING EMPTY
2 Alcochete    Almada LINESTRING (-8.91588 38.735...
3    Almada    Almada               LINESTRING EMPTY
4 Alcochete   Amadora LINESTRING (-8.91588 38.735...
5    Almada   Amadora LINESTRING (-9.193069 38.63...
6   Amadora   Amadora               LINESTRING EMPTY
Code
tail(TRIPSdlines_oneway[,c(1,2)])
Simple feature collection with 6 features and 2 fields
Attribute-geometry relationships: identity (2) (with 1 geometry empty)
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -9.357651 ymin: 38.4949 xmax: -8.806644 ymax: 38.92211
Geodetic CRS:  WGS 84
                      o                   d                       geometry
163             Palmela Vila Franca de Xira LINESTRING (-8.806644 38.61...
164              Seixal Vila Franca de Xira LINESTRING (-9.108801 38.60...
165            Sesimbra Vila Franca de Xira LINESTRING (-9.120124 38.49...
166             Setúbal Vila Franca de Xira LINESTRING (-8.887481 38.51...
167              Sintra Vila Franca de Xira LINESTRING (-9.357651 38.82...
168 Vila Franca de Xira Vila Franca de Xira               LINESTRING EMPTY

Example of visualization with Public Transit trips in both ways.

Code
TRIPSdlines_oneway_noLX = TRIPSdlines_oneway |> 
  filter(o != d) |> # remove intrazonal trips
  filter(PTransit > 5000) # reduce noise

mapview(TRIPSdlines_oneway_noLX, zcol = "PTransit", lwd = 8)

9.3 Using population centroids

The od_to_sf() function uses the geometric center of the zones to create the desire lines. But if we replace those zones by the weighted centroids, we can have a more realistic representation of the flows.

Code
# Centroid_pop = st_read("data/Centroid_pop.gpkg")

TRIPSdlines_pop = od_to_sf(TRIPSmode, z = Centroid_pop) |>  # works the same way
  od_oneway() # oneway

Check differences of lines with trips from/to Lisbon:

Code
TRIPSdlines_geo_LX = TRIPSdlines_oneway |> 
  filter(o == "Lisboa" | d == "Lisboa") # or condition
TRIPSdlines_pop_LX = TRIPSdlines_pop |> 
  filter(o == "Lisboa" | d == "Lisboa")

mapview(TRIPSdlines_geo_LX, color = "blue") + mapview(TRIPSdlines_pop_LX, color = "red")

The next step will be estimating the euclidean distances between these centroids, and compare them with the routing distances.


  1. See (Lovelace, Félix, and Carlino 2022).↩︎