Workshop: Using R for Geodata Processing

A Guide for the Digital Humanist

Author
Affiliation
Notice also…

While this page attempts to make the map in {leaflet}, I have now added another page in which I make it with the {ggplot2} package. Check it out.

1 Cheatsheet for working in R

1.1 Where to work

  1. Go to www.posit.cloud
  2. Open a free account using a Google email (e.g. @utexas.edu)
  3. Go into your Posit workspace
  4. Create a new RStudio project, not from a template
  5. Rearrange the panes in your project
    1. Go to Tools > Global Options > Pane Layout
    2. Place Source at top left, Console at top right, Environment at bottom left, and Files at bottom right.

Posit vs. RStudio
Posit is the new name for RStudio. It’s still mostly the same thing: a graphical and highly useful wrapper around R, except now you can also use it for Python.

1.2 Essential keyboard shortcuts

Action Shortcut
Insert pipe %>%. Shift + Cmd + M
Insert assignment operator <- Alt + -
Run chunk in which your cursor is Cmd + Enter
Run all (the whole script) Option + Cmd + R

2 Setup

2.1 Setup, loading packages, getting data ready

To get your working session in R started, you need to think about three things:

  1. Setting up a new RStudio project
  2. Loading some “packages”, i.e. software modules that expand the capabilities of R
  3. Creating a new, empty R script file by clicking File > New File > R Script. Hit “save” and give it a name; store it in the project folder (the default option).

Start your script by copy/pasting the following chunk at the top of the script.

These steps are a great way to start a new project in R/RStudio. Re-use this beginning code chunk if you wish.

# If the {pacman} package is not yet installed in your RStudio 
# environment, uncomment the following line (and re-comment once 
# installed)
# install.packages("pacman")
library(pacman)

# Load some general-use packages
# (i.e., begin every script with these)
p_load(rio, tidyverse, janitor, conflicted)

# Also load packages that are more specific to the demands
# of the type of data and analytical task you're working
# with, e.g. the geomapping package {leaflet} for our 
# mapping exercise
p_load(leaflet)

# The following statement resolves a potential conflict 
# between two different functions that exist in different parts of 
# R which both have the same name
conflict_prefer("filter", "dplyr")

Our sample data for this exercise will be the Archestratos data. Download it (link below) and save the CSV file in your project folder. You should be able to see it in the “Files” pane, bottom right.

Now load the data using import() from the {rio} package. We’ll call our dataframe df.

The package name {rio} stands for “R in/out”. It contains a number of functions and capabilities for data reading/loading as well as writing to different file formats.

# Read the data file, then clean colnames and format
# as a tibble.
df <- import("archestratos_rmaps.csv") %>% 
  clean_names() %>% 
  tibble()

To preview your data, you have a number of options.

  1. Find the df data object in the “Environment” pane, where it should now exist. Click on it: a preview should open in RStudio in a new tab next to your script file.
  2. Enter just the name of the data object - df - as code in your script and “Run” the statement. In the “Console” you should see a tibble version of the data.
  3. Use the glimpse() function: as a practice, run the following code to see this extremely useful preview function: df %>% glimpse(). The output produced by this function is a list of all columns in your dataset.

A tibble is a special kind of table, a way that the {tidyverse} formats datasets. It has certain useful features like the information about the number of rows and columns in it, and about the type of data in each column (numeric vs. character-type data).

I’ll use the second option here, so we can look at a basic version of the data.

df
# A tibble: 85 × 8
   fragment_line food_type     food_specific location pleiades_uri     lon   lat
           <dbl> <chr>         <chr>         <chr>    <chr>          <dbl> <dbl>
 1          44   fish          gray mullet   Abdera   http://pleiad…  25.0  41.0
 2          56.1 invertebrates cuttle-fish   Abdera   http://pleiad…  25.0  41.0
 3           7.1 shellfish     oysters       Abydos   http://pleiad…  26.4  40.2
 4          43.1 fish          gray mullet   Aigina   http://pleiad…  23.4  37.8
 5           7.3 shellfish     scallops      Ambracia http://pleiad…  21.0  39.2
 6          16.1 fish          boar-fish     Ambracia http://pleiad…  21.0  39.2
 7          26.3 shellfish     large shrimp  Ambracia http://pleiad…  21.0  39.2
 8          31.2 fish          chromis       Ambracia http://pleiad…  21.0  39.2
 9          45.5 fish          sea-bass      Ambracia http://pleiad…  21.0  39.2
10          55.2 invertebrates squid         Ambracia http://pleiad…  21.0  39.2
# ℹ 75 more rows
# ℹ 1 more variable: latlon <chr>

3 Make a leaflet map from our data

3.1 Map with standard (modern) tiles

The following is the complete map code presented in class on 9/17/25. Notice that this code creates the map and assigns it to a variable name (map1). After that, to display the map, the variable name is called.

map1 <- df %>% 
  leaflet() %>%
  addTiles() %>%  # Add default OpenStreetMap tiles
  addMarkers(
    ~lon, 
    ~lat, 
    label = ~food_type,  # Display food_type on the flag
    popup = ~paste0(
      "<b>Food Specific:</b> ", food_specific, "<br>",
      "<b>Website:</b> <a href='", pleiades_uri, "' target='_blank'>Go to Website</a>"
    ),  # Create a clickable link in the popup
    labelOptions = labelOptions(noHide = TRUE)  # Keep the labels always visible
  )

map1

3.2 Map with ancient tiles

I tried to get this to display on a map of the ancient world. Following instructions I found online I downloaded an old map of the Mediterranean (of 238 BC) from Wikimedia. It was just an image this one, so it needed to be “warped”.

Warping maps seems to be an indexing process where a map gets squeezed into a digital idea of what the world is like. You can get this done by using a website called Map Warper.

Having tried various methods of getting this scanned-and-warped map to appear in the leaflet map, the one that worked (sort of) is the one that uses the GeoTiff version.

The image is initially too large in file size, so must be reduced.

p_load(png, raster)

# Load the GeoTIFF file
raster_file <- raster("./97142.tif")

# Reduce the resolution by aggregating (factor = 2 
# combines 2x2 pixels into 1)
raster_resampled <- aggregate(raster_file, fact = 2)  

After that, we can run kind of the same map code as before, but insert the GeoTiff as an overlay.

map2 <- df %>% 
  leaflet() %>%
  addTiles() %>%
  addRasterImage(
    raster_resampled,
    #colors = pal,  # Apply the custom color palette
    opacity = 0.9
  ) %>%
  addMarkers(
    ~lon, 
    ~lat, 
    label = ~food_type,  # Display food_type on the flag
    popup = ~paste0(
      "<b>Food Specific:</b> ", food_specific, "<br>",
      "<b>Website:</b> <a href='", pleiades_uri, 
      "' target='_blank'>Go to Website</a>"
    ),  # Create a clickable link in the popup
    labelOptions = labelOptions(noHide = TRUE)  
  ) %>% 
  setView(
    lng = 15.0, lat = 37.0,  # Center on Medit.
    zoom = 4             
  )
map2
Two steps forward, one step back?

Clearly we’re not all the way there yet, but something is working. Love how carefully the warper folded the square map sheet around the imaginary globe. - What is up with these colors?