Mapping: Creating lines from points

A common task for mapping humanities data is to create lines from sets of points to demonstrate the movement of people, objects, or ideas. To do this we need to get the data in the right format, create point data using the sf package as described in Mapping with R, and then connect the points into lines using a group_by() and summarize() pipeline.

Two different ways to organize the data are demonstrated here, resulting in different types of line data. This page also demonstrates the usefulness of joins between multiple data frames and working with geocoded data as demonstrated in Wrangling data in the tidyverse.

Download

Download a version of this Quarto document to follow along.

We start, as always, by loading the packages we will be using. The main packages are tidyverse and sf, but we will also be creating interactive maps with leaflet and using paletteer for a color palette.

library(tidyverse)
library(sf)
library(leaflet)
library(paletteer)

The data

The data is available through the datasets folder on the class syllabus.

movements <- read_csv("https://raw.githubusercontent.com/jessesadler/vt5444s26/refs/heads/main/datasets/movements.csv")
locations <- read_csv("https://raw.githubusercontent.com/jessesadler/vt5444s26/refs/heads/main/datasets/locations.csv")

The movements data shows the date and location of siblings from two different families at the end of the 16th century.

movements

# A tibble: 74 × 3
   person                 place     date      
   <chr>                  <chr>     <date>    
 1 Andries van der Meulen Antwerp   1584-07-01
 2 Andries van der Meulen Bremen    1585-09-01
 3 Anna van der Meulen    Cologne   1584-07-01
 4 Anna van der Meulen    Bremen    1588-06-01
 5 Anna van der Meulen    Stade     1588-08-01
 6 Anna van der Meulen    Bremen    1592-03-01
 7 Sara van der Meulen    Antwerp   1584-07-01
 8 Sara van der Meulen    Bremen    1585-09-01
 9 Sara van der Meulen    Frankfurt 1586-03-01
10 Sara van der Meulen    Cologne   1586-05-01
# ℹ 64 more rows

The locations data has the latitude and longitude values for the locations that these individuals visited in their travels. This data could have been created by geocoding the locations with tidygeocoder as shown in Mapping with R. It is best practice to save this output to a data folder so you are not re-geocoding the locations every time you run the code.

locations

# A tibble: 21 × 3
   location          lat   lng
   <chr>           <dbl> <dbl>
 1 Antwerp          51.2  4.40
 2 Bremen           53.1  8.81
 3 Cologne          50.9  6.96
 4 Stade            53.6  9.48
 5 Frankfurt        50.1  8.68
 6 Utrecht          52.1  5.13
 7 Leiden           52.2  4.48
 8 Geertruidenberg  51.7  4.85
 9 Bergen op Zoom   51.5  4.29
10 The Hague        52.1  4.31
# ℹ 11 more rows

Let’s convert the locations data to an sf object to make it truly spatial data. This data can be used to make the points on the map.

locations_sf <- st_as_sf(locations,
                         coords = c("lng", "lat"),
                         crs = 4326)

The workflow

We can now work on creating lines from the data. First, we will take the data as it is, creating one line for each individuals in the movements data. Secondly, we will create separate lines for each journey between cities, which will entail a bit more wrangling of the data.

The basic workflow is:

Join the movements data to the latitude and longitude values in locations using a left_join().
Convert the resulting data frame into an sf object with st_as_sf().
Group and summarize the data and the cast from points to lines using st_cast().

Note

It is possible to combine steps one and two by doing a join on locations_sf. However, the sf object must be the first data frame in a left_join() in order to create an sf object, and this puts the geometry column on the left side of the data frame, which I find aesthetically displeasing.

One line for each individual

Let’s start by doing steps one and two. This is a common workflow when working with geographic data, especially vector data in the form of points, joining geocoded data to attribute data and then turning it into an sf object.

movements_geo <- left_join(movements, locations,
                           by = join_by(place == location))
movements_sf <- st_as_sf(movements_geo,
                         coords = c("lng", "lat"),
                         crs = 4326)

Now we can move onto the third step, which is the most crucial. This uses the group by and summarize workflow that is familiar from data wrangling in the tidyverse, but there are some aspects when working with sf objects. summarise() works on the geometry column by default so it can be left blank. However, we set do_union = FALSE to have the points be arranged in order from top to bottom of the data frame. This is why the pipeline begins with arrange(). It is st_cast() that does the actual transformation from points to lines. In this case, it transforms the geometry column from "MULTIPOINT" to "LINESTRING".

movements_lines <- movements_sf |> 
1  arrange(date) |>
2  group_by(person) |>
3  summarise(do_union = FALSE) |>
4  st_cast("LINESTRING")

1: Arrange the rows from earliest to latest date to ensure lines are made in the right order.
2: Group by person so that each person will have one line.
3: Summarize geometry column, turning POINT data to MULTIPOINT.
4: Cast MULTIPOINT to LINESTRING.

movements_lines

Simple feature collection with 11 features and 1 field
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -0.1276474 ymin: 50.11064 xmax: 10.00991 ymax: 53.59979
Geodetic CRS:  WGS 84
# A tibble: 11 × 2
   person                                                               geometry
   <chr>                                                        <LINESTRING [°]>
 1 Andries van der Meulen                 (4.399708 51.22111, 8.807165 53.07582)
 2 Anna della Faille                       (4.64356 52.38371, 4.399708 51.22111)
 3 Anna van der Meulen    (6.959974 50.93836, 8.807165 53.07582, 9.475438 53.59…
 4 Carlo della Faille     (4.686789 51.76895, 4.604206 51.64844, 4.399708 51.22…
 5 Daniel van der Meulen  (4.399708 51.22111, 3.613696 51.49973, 4.362724 51.99…
 6 Hester della Faille    (4.64356 52.38371, 8.807165 53.07582, 4.481109 52.151…
 7 Jacques della Faille   (4.64356 52.38371, 10.00991 53.54376, 7.2058 53.36704…
 8 Jan della Faille       (4.399708 51.22111, 4.481109 52.15182, 4.846743 51.69…
 9 Marten della Faille    (4.399708 51.22111, 4.289741 51.49659, 4.311346 52.07…
10 Sara van der Meulen    (4.399708 51.22111, 8.807165 53.07582, 8.682092 50.11…
11 Steven della Faille    (4.686789 51.76895, 4.399708 51.22111, 4.351697 50.84…

You can see that each individual has a single row. Their movements are represented by a single line, but those with many movements have more stopping points in the line. We can see the results by creating a leaflet map, using paletteer to help create a palette to represent the individuals.

# Create a palette
pal <- colorFactor(
  palette = as.character(paletteer_d("colorBlindness::paletteMartin")),
  domain = movements_lines$person
)

leaflet() |> 
  addTiles() |>
  addPolylines(data = movements_lines,
               color = ~pal(person),
               opacity = 0.8,
               label = ~person) |> 
  addCircleMarkers(data = locations_sf,
                   fillOpacity = 0.5,
                   radius = 5,
                   label = ~location) |> 
  addLegend("bottomleft", pal = pal, values = movements_lines$person)

Figure 1: A single line representing the movements of individuals.

Note

There is currently a bug in how Quarto renders leaflet legends. You can fix it with this solution. For an example, see how it is done on this page.

One line for each journey

This get’s us the visualization we wanted, but it loses the data about the date of the journey and the direction. This can be rectified by making one line for each journey. However, to do this we need to do a bit of data wrangling. I am sure that there are many ways to do this process, but this is how I would proceed.

First, we need have start and destination locations on each row, so that each row represents a single journey. The key to this is the lag() function that fills in the previous value combined with grouping by person. To see how lag() works run lag(1:5).

movements_wide <- movements |> 
  arrange(date) |> 
1  group_by(person) |>
2  mutate(start = lag(place), .after = person) |>
3  ungroup() |>
4  rename(destination = place) |>
5  filter(!is.na(start)) |>
6  rowid_to_column("journey_id")

1: It is important to group the data so that lag() looks to the previous value for each person, not just from the previous row.
2: Create a new column named start and place it after the person column for convenience.
3: Ungroup the data frame since we are done acting on the grouped data.
4: Rename the place column to destination to better represent the nature of the data in the column.
5: Remove the NAs introduced into the start column by lag() where there was no previous location.
6: Add an id column to keep track of each journey. This will be used for grouping in the next step.

movements_wide

# A tibble: 63 × 5
   journey_id person                 start      destination date      
        <int> <chr>                  <chr>      <chr>       <date>    
 1          1 Daniel van der Meulen  Antwerp    Middelburg  1584-08-01
 2          2 Daniel van der Meulen  Middelburg Delft       1584-09-01
 3          3 Jan della Faille       Antwerp    Leiden      1584-11-01
 4          4 Daniel van der Meulen  Delft      Haarlem     1584-12-01
 5          5 Daniel van der Meulen  Haarlem    Delft       1585-02-01
 6          6 Andries van der Meulen Antwerp    Bremen      1585-09-01
 7          7 Sara van der Meulen    Antwerp    Bremen      1585-09-01
 8          8 Anna della Faille      Haarlem    Antwerp     1585-09-01
 9          9 Steven della Faille    Dordrecht  Antwerp     1585-09-01
10         10 Jacques della Faille   Haarlem    Hamburg     1585-10-01
# ℹ 53 more rows

movements_wide accomplishes one goal in turning each row into a single journey, but now the task of adding our spatial data is more difficult. We have two locations per row, but we only want one geometry column. The solution is to transform our data from its current wide format to a longer format in which each journey is split into two rows, one for the start and one for the destination. An example of this is shown in Wrangling data in the tidyverse.

movements_long <- movements_wide |> 
  pivot_longer(
1    cols = start:destination,
2    names_to = "type",
3    values_to = "location"
  )

1: The columns to pivot.
2: The name of the column for the column names that will be pivoted.
3: The name of the column for the values in the pivoted columns.

The meaning of the arguments are more clear when you see the outcome.

movements_long

# A tibble: 126 × 5
   journey_id person                date       type        location  
        <int> <chr>                 <date>     <chr>       <chr>     
 1          1 Daniel van der Meulen 1584-08-01 start       Antwerp   
 2          1 Daniel van der Meulen 1584-08-01 destination Middelburg
 3          2 Daniel van der Meulen 1584-09-01 start       Middelburg
 4          2 Daniel van der Meulen 1584-09-01 destination Delft     
 5          3 Jan della Faille      1584-11-01 start       Antwerp   
 6          3 Jan della Faille      1584-11-01 destination Leiden    
 7          4 Daniel van der Meulen 1584-12-01 start       Delft     
 8          4 Daniel van der Meulen 1584-12-01 destination Haarlem   
 9          5 Daniel van der Meulen 1585-02-01 start       Haarlem   
10          5 Daniel van der Meulen 1585-02-01 destination Delft     
# ℹ 116 more rows

Now we can do the same as before, joining the locations data by the single location column, converting it to an sf object and then creating lines.

movements_long_sf <- movements_long |> 
  left_join(locations, by = "location") |> 
  st_as_sf(coords = c("lng", "lat"), crs = 4326) |> 
  group_by(journey_id) |> 
  summarise(do_union = FALSE) |> 
  st_cast("LINESTRING")

The result, as before, is a sf data frame with two columns: journey_id and geometry. We can get back the attribute columns for the journeys by joining movements_long_sf and movements_wide, being careful to place movements_long_sf as the first data frame so that the result is an sf object.

journeys_sf <- left_join(movements_long_sf, movements_wide, by = "journey_id")
journeys_sf

Simple feature collection with 63 features and 5 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: -0.1276474 ymin: 50.11064 xmax: 10.00991 ymax: 53.59979
Geodetic CRS:  WGS 84
# A tibble: 63 × 6
   journey_id                       geometry person start destination date      
        <int>               <LINESTRING [°]> <chr>  <chr> <chr>       <date>    
 1          1 (4.399708 51.22111, 3.613696 … Danie… Antw… Middelburg  1584-08-01
 2          2 (3.613696 51.49973, 4.362724 … Danie… Midd… Delft       1584-09-01
 3          3 (4.399708 51.22111, 4.481109 … Jan d… Antw… Leiden      1584-11-01
 4          4 (4.362724 51.99946, 4.64356 5… Danie… Delft Haarlem     1584-12-01
 5          5 (4.64356 52.38371, 4.362724 5… Danie… Haar… Delft       1585-02-01
 6          6 (4.399708 51.22111, 8.807165 … Andri… Antw… Bremen      1585-09-01
 7          7 (4.399708 51.22111, 8.807165 … Sara … Antw… Bremen      1585-09-01
 8          8 (4.64356 52.38371, 4.399708 5… Anna … Haar… Antwerp     1585-09-01
 9          9 (4.686789 51.76895, 4.399708 … Steve… Dord… Antwerp     1585-09-01
10         10 (4.64356 52.38371, 10.00991 5… Jacqu… Haar… Hamburg     1585-10-01
# ℹ 53 more rows

Now we can plot the data, but first let’s create a label column that we can use in the leaflet map using paste() to create a character vector and <br/> to implement line breaks in the HTML used for the popups.

journeys_sf <- journeys_sf |> 
  mutate(label = paste(
    person, "<br/>",
    start, "to", destination, "<br>",
    day(date), month(date, label = TRUE), year(date)))

leaflet() |> 
  addTiles() |>
  addPolylines(data = journeys_sf,
               color = ~pal(person),
               opacity = 0.8,
               popup = ~label) |> 
  addCircleMarkers(data = locations_sf,
                   fillOpacity = 0.5,
                   radius = 5,
                   label = ~location) |> 
  addLegend("bottomleft", pal = pal, values = journeys_sf$person)

Figure 2: Multiple line representing the journeys of individuals.

Bonus: Creating great circles

The lines that are created this way are rhumb lines, straight lines at a constant bearing. But another way to represent the lines is with great circles, which represent the shortest distance between points on a spherical representation of the earth.

Happily, there is a function in sf designed specifically for this, st_segmentize(). This function adds segments to the line, making it appear curved. The function depends on the lwgeom packge, which needs to be installed separately to use the function. st_segmentize() also takes advantage of the units package, which allows you to specify the type of units a numeric value represents. This is used for calculating the longest distance between segments of the line. Let’s see how this works with journeys_sf.

journeys_gc <- journeys_sf |> 
  st_segmentize(units::set_units(20, km))

If you compare the geometry column of journeys_gc to journeys_sf, you will notice how many more segments there are in each line in journeys_gc. Substitute journeys_gc for journeys_sf in the above map to see the difference. This difference is more noticeable the longer the distance of the line.

blacksburg_london <- tibble(
  location = c("Blacksburg", "London"),
  lat = c(37.22966, 51.50732),
  lng = c(-80.41368, -0.1276474)
) |> 
  st_as_sf(coords = c("lng", "lat"), crs = 4326)

rhumb <- blacksburg_london |> 
  # do not need group by because want to summarize all rows (2) into one
  summarise(do_union = FALSE) |> 
  st_cast("LINESTRING")

great_circle <- rhumb |> 
  st_segmentize(units::set_units(20, km))

comparison <- bind_rows(rhumb, great_circle) |> 
  add_column(type = c("Rhumb line", "Great circle"))

Let’s see what the difference looks like.

# Create a palette
pal <- colorFactor(
  palette = as.character(paletteer_d("colorblindr::OkabeIto")),
  domain = comparison$type
)

leaflet() |> 
  addTiles() |>
  addPolylines(data = comparison,
               opacity = 1,
               label = ~type,
               color = ~pal(type)) |> 
  addCircleMarkers(data = blacksburg_london,
                   fillOpacity = 0.5,
                   radius = 5,
                   label = ~location)

Figure 3: Rhumb line vs great circle.