EVR 628- Intro to Environmental Data Science

Working with spatial data in R

Juan Carlos Villaseñor-Derbez (JC)

Quick announcements

VIzUM 2025: Critical Cartography. Monday, November 10, 2025 | 4:00-7:00 PM

Introduction

GIS vs Geocomputation

GIS:

  • Powerful
  • Emphasize graphical user interface (GUI)
  • You can only do what the developers allow you to do
  • Difficult to reproduce

Geocomputation:

  • Just as powerful (same underlying code)
  • Command line interface (CLI)
  • Flexible
  • Reproducible

Academic research, software development and practical applications that use geographic data to solve problems, with a focus on reproducibility, flexibility and tool development. - Lovlace, Nowosad, Muenchow (GeocompR, 2024)

How do we get the world into a computer?

  • Imagine an attribute that occurs somewhere in the World:
    • SST
    • site identity
    • hurricane path
    • city boundaries
  • How do we represent it?
  • How do we work on it?
  • How do we store the data?

We need a model

Spatial data can be represented in a vector model, or a raster model

In today’s class:

  • Introduce the why behind each model
  • Discuss where and when they are used
  • Focus on vector:
    • It’s main attributes
    • It’s implementation in R

Learning Objectives

By the end of the lecture, you should be able to:

  • list the two main data models, and articulate their main characteristics
  • understand the principles of vector data
  • list the three main components of a spatial vector object
  • list the tree main types of a simple feature
  • given an sf object, plot it and inspect it

Two models

vector

  • World is represented using geometries: points, lines, and polygons
  • Discrete and well-defined borders = high precision
  • Computationally intensive

raster

  • Space is divided into equal-sized grid cells (or pixels)
  • Aggregate spatial features to a given resolution
  • You can lose some precision, but much faster to work with

Example: Rosenstiel

Code 
# Build the "spatial" part
rsmaes_pt_coords <- st_point(x = c(-80.1632879, 25.7328129),
                             dim = "XY") |> 
  st_sfc(crs = 4326)
# Combine with the "attributes" part
rsmaes_pt <- st_sf(id = "Rosenstiel (point)",
                   geometry = rsmaes_pt_coords)
# Visualize
m <- mapview(rsmaes_pt) 
m@map |> leaflet::setView(lng = -80.16328, lat = 25.73281, zoom = 17)
Code 
rsmaes_line_coords <- st_linestring(
  x = matrix(data = c(-80.163017, 25.733950,
                      -80.164236, 25.732816,
                      -80.163772, 25.732353,
                      -80.163924, 25.732148,
                      -80.163455, 25.731597,
                      -80.162187, 25.731605,
                      -80.160968, 25.732172,
                      -80.163017, 25.733950),
             ncol = 2,
             byrow = T),
  dim = "XY") |> 
  list() |> 
  st_sfc(crs = 4326)

rsmaes_line <- st_sf(id = "Rosenstiel (line)",
                     geometry = rsmaes_line_coords)

m <- mapview(rsmaes_line)
m@map |> leaflet::setView(lng = -80.16328, lat = 25.73281, zoom = 17)
Code 
rsmaes_poly_coords <- st_polygon(
  x = list(
    matrix(
      data = c(
        -80.163017, 25.733950,
        -80.164236, 25.732816,
        -80.163772, 25.732353,
        -80.163924, 25.732148,
        -80.163455, 25.731597,
        -80.162187, 25.731605,
        -80.160968, 25.732172,
        -80.163017, 25.733950),
      ncol = 2,
      byrow = T)),
  dim = "XY") |> 
  st_sfc(crs = 4326)

rsmaes_poly <- st_sf(id = "Rosenstiel (polygon)",
                     geometry = rsmaes_poly_coords)

m <- mapview(rsmaes_poly)
m@map |> leaflet::setView(lng = -80.16328, lat = 25.73281, zoom = 17)
Code 
rsmaes_rast <- rasterize(
  x = vect(rsmaes_poly),
  y = rast(resolution = 0.0001,
           crs = "EPSG:4326",
           val = 0,
           xmin = -80.168, xmax = -80.155,
           ymin = 25.730, ymax = 25.735),
  field = 1)

names(rsmaes_rast) <- "Rosenstiel"
rsmaes_rast[is.na(rsmaes_rast)] <- 0

m <- mapview(rsmaes_rast, legend = F)
m@map |> leaflet::setView(lng = -80.16328, lat = 25.73281, zoom = 17)

Common uses

vector

  • Common in social sciences:
    • Census tract boundaries
    • Political boundaries
    • Roads

raster

  • Common in natural sciences:
    • Sea surface temperature
    • Satellite imagery
    • Species distribution models

  • As environmental scientists and practitioners we will typically use both

  • The “correct” type of data varies by use case

“raster is faster, but vector is ‘correcter’”

Vector data model

World is represented using geometries: points, lines, and polygons

Q: What is a line? What is a polygon?

A: A sequence of points

Q: How do we represent a point?

A: Coordinates

Q: What are coordinates?

A: Two numbers representing distance from an origin: a coordinate reference system

What is where? And what does “where” mean?

Coordinates, Coordinate Reference Systems, and Attributes

What does “where” mean?

Consider the following two coordinates:

  • (-80.1632879, 25.7328129)

  • (1186764, 2863941)

Where are they?

The importance of attributes

  • Coordinates and a CRS are enough to place the point on a map
  • Also enough to perform spatial analysis
  • But not enough to convey information about the location
  • Often not enough to derive insights
  • We also need attributes

Note

The three main components of a spatial vector object are attributes (what), geometry, and CRS (where)

Attributes

  • Are not inherently spatial
  • Can be a discrete or categorical variable representing a characteristic of interest:
    • Identity
    • Sea Surface Temperature
    • Time
  • Note that you can have more than one attribute tied to a geometry
    • 2 attributes: vessel-tracking data or hurricane path (id and time)
    • 2+ attributes: Buoy weather data (id, time, and all variables measured)

Handling spatial data in a computer

How are these compotents handled in a computer?

We rely on four key low-level libraries for geocomputation:

  • GDAL: read, write, and manipulating geographic data
  • PROJ: coordinate system transformation
  • GEOS: planar geometry engine only to be used with data with a projected CRS)
  • S2: a spherical geometry engine, written in C++ (fast) and developed by Google (reliable)
  • None of these are R packages

All are incredibly powerful and open source, but not exactly user-friendly

Ways of interacting with these libraries

. . .

From R, our main point of access will be the sf package

The sf package

{sf}: a package that provides simple features access for R

{sf} provides:

  • A consistent command-line interface to communicate with the libraries
  • Tidyverse-like access to attributes data (dplyr and tidyr)
  • Visualization of attribute and spatial data (geom_sf())
  • Spatial and geometry operations

When we load sf, it automatically links to these low level libraries

library(sf)
Linking to GEOS 3.13.0, GDAL 3.10.0, PROJ 9.5.1; sf_use_s2() is TRUE

Introduction to Simple Features

Simple Features

  • Simple features is an official, open standard ISO 19125
  • Hierarchical data model to represent types of geometries (i.e. points, lines, polygons, multi-*, and collections)

The standard says:

“A simple feature is defined by the OpenGIS Abstract specification to have both spatial and non-spatial attributes. Spatial attributes are geometry valued, and simple features are based on 2D geometry with linear interpolation between vertices.”

—Pebesma, 2024

Quick example 1: inspecting sf objects

library(EVR628tools)
data("data_milton_sf") # Load
class(data_milton_sf) # Inspect it
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
dim(data_milton_sf) # Dimensions
[1] 62  6
head(data_milton_sf)
Simple feature collection with 6 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -95.5 ymin: 20.8 xmax: -94.6 ymax: 21.7
Geodetic CRS:  WGS 84
# A tibble: 6 × 6
  name   iso_time            wind_speed pressure sshs      geometry
  <chr>  <dttm>                   <dbl>    <dbl> <chr>  <POINT [°]>
1 MILTON 2024-10-04 18:00:00         30     1009 -3      (-94.6 21)
2 MILTON 2024-10-04 21:00:00         30     1009 -3    (-94.9 20.8)
3 MILTON 2024-10-05 00:00:00         30     1008 -3    (-95.2 20.8)
4 MILTON 2024-10-05 03:00:00         30     1008 -3      (-95.3 21)
5 MILTON 2024-10-05 06:00:00         30     1008 -3    (-95.4 21.4)
6 MILTON 2024-10-05 09:00:00         30     1008 -3    (-95.5 21.7)

Quick example 2: Visually inspecting sf objects

plot(data_milton_sf) # Plot it

Note you get one plot per column (variable), except for geometry

Quick example 3: An sf object is just a data.frame / tibble

You can quickly calculate summary statistics for any attribute variable (i.e. column)

summary(data_milton_sf$wind_speed)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  30.00   53.50   74.00   87.18  134.50  155.00

Keep only some columns

milton_light <- data_milton_sf |> 
  select(iso_time, wind_speed)
milton_light
Simple feature collection with 62 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -95.5 ymin: 20.8 xmax: -65.9 ymax: 29.5
Geodetic CRS:  WGS 84
# A tibble: 62 × 3
   iso_time            wind_speed     geometry
   <dttm>                   <dbl>  <POINT [°]>
 1 2024-10-04 18:00:00         30   (-94.6 21)
 2 2024-10-04 21:00:00         30 (-94.9 20.8)
 3 2024-10-05 00:00:00         30 (-95.2 20.8)
 4 2024-10-05 03:00:00         30   (-95.3 21)
 5 2024-10-05 06:00:00         30 (-95.4 21.4)
 6 2024-10-05 09:00:00         30 (-95.5 21.7)
 7 2024-10-05 12:00:00         30   (-95.5 22)
 8 2024-10-05 15:00:00         33 (-95.5 22.3)
 9 2024-10-05 18:00:00         35 (-95.5 22.5)
10 2024-10-05 21:00:00         35 (-95.5 22.6)
# ℹ 52 more rows

geometry columns are “sticky”

Quick example 4: Yes, sf objects play along with ggplot

Either for the spatial data

ggplot(data = milton_light) +
  geom_sf() +
  labs(x = "Lon",
       y = "Lat")

Or the attribute variables

ggplot(data = milton_light,
       mapping = aes(x = iso_time,
                     y = wind_speed)) +
  geom_line() +
  geom_point() +
  labs(x = "Date",
       y = "Wind speed (knots)")

Why simple features?

Why use simple features

  • STANDARD = transferability and reproducibility
  • Widely supported data model (e.g. QGIS and PostGIS)
  • Relatively quick I/O
  • Full R support: sf objects can be treated as data.frames (wrangling and visualization, pipe-friendly)

Visualizing simple features

(Just an overview for now)

Basic maps

You can specify you only want one panel

USA <- rnaturalearth::ne_states(geounit = "United States of America") |> 
  filter(!name %in% c("Hawaii", "Alaska"))

plot(USA, max.plot = 1)

Other parameters to plot()

Or that you specifically want a column (either by name or by position)

plot(USA["name"])

plot(USA[,18])

  • These are mainly helpful as check for ourselves
  • You won’t always get a legend

Adding layers

Use reset = FALSE and add = TRUE

plot(USA[,1], reset = F)
plot(FL, add = T, col = "orange")
plot(data_milton_sf, add = T, col = "red")
plot(MEX, add = T, col = "gray")

  • reset = FALSE needed when plotting the first layer produces a key OR the extent needs to be updated as new layers come in
  • See ?sf::plot_sf() for more details

Basic = limitations

Base plot() is great for:

  • Rapidly visualizing your data (just a few lines of code)
  • Checking correspondence between layers

Base plot() is terrible for:

  • Building complex visualizations (including formal maps)

We will cover visualization in more detail in a few weeks. In the meantime, two quick ways

Visualizing sf objects with ggplot

p <- ggplot(data = USA) +
  geom_sf(fill = "gray90")

p

p +
  geom_sf(data = FL, fill = "orange") +
  geom_sf(data = data_milton_sf,
          mapping = aes(color = wind_speed), size = 1) +
  geom_sf(data = MEX)

This is what I use to build 99% of my maps, including those for publications and proposals

Visualizing sf objects with mapview

  • Just like plot(), mapview() produces quick visualizations
  • Customization is painful, but worth it for communication
mapview(list(USA, FL, data_milton_sf, MEX))

This is what I use to interactively explore my data

Note

For now, recall that you can quickly visualize sf objects using plot(), ggplot(), or mapview()

Geometry types

The main types of simple features

  • POINT zero-dimensional geometry containing a single point
  • LINESTRING sequence of points connected by straight, non-self intersecting line pieces; one-dimensional geometry
  • POLYGON geometry with a positive area (two-dimensional); sequence of points form a closed, non-self intersecting ring; the first ring denotes the exterior ring, zero or more subsequent rings denote holes in this exterior ring
  • MULTIPOINT set of points
  • MULTILINESTRING set of linestrings
  • MULTIPOLYGON set of polygons
  • GEOMETRYCOLLECTION set of geometries of any type except GEOMETRYCOLLECTION

GeoCompR, 2024

sf objects

Remember the objects I built on the slide with tabs?

Point

class(rsmaes_pt)
[1] "sf"         "data.frame"
rsmaes_pt
Simple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -80.16329 ymin: 25.73281 xmax: -80.16329 ymax: 25.73281
Geodetic CRS:  WGS 84
                  id                   geometry
1 Rosenstiel (point) POINT (-80.16329 25.73281)

Polygon

class(rsmaes_poly)
[1] "sf"         "data.frame"
rsmaes_poly
Simple feature collection with 1 feature and 1 field
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -80.16424 ymin: 25.7316 xmax: -80.16097 ymax: 25.73395
Geodetic CRS:  WGS 84
                    id                       geometry
1 Rosenstiel (polygon) POLYGON ((-80.16302 25.7339...

Note

Both are sf objects, which are data.frames /tibbles with a special column called geometry

Let’s dive deeper into how an object of class sf is made

The sf class

An sf object

Two main parts: geometries (spatial) and attributes (non-geographic)

GeoCompR, 2024 - Rhomboids are functions, rectangles are objects
rsmaes_pt <- c(-80.1632879, 25.7328129) |>      # Coordinates
  st_point() |>                                 # Build sfg
  st_sfc(crs = "ESRI:4326") |>                  # Build sfc (an specify CRS)
  st_sf(                                        # Build sf
    data.frame(id = "Rosenstiel (point)",       # Add the data.frame portion
               foundation = 1943))        
rsmaes_pt
Simple feature collection with 1 feature and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -80.16329 ymin: 25.73281 xmax: -80.16329 ymax: 25.73281
CRS:           NA
                  id foundation
1 Rosenstiel (point)       1943
  st_sfc.st_point.c..80.1632879..25.7328129....crs....ESRI.4326..
1                                      POINT (-80.16329 25.73281)

Simple feature geometries (sfg)

  • These are the essential building blocks
  • They represent the type of simple features that can be handled in R
  • There is one function to create each type
  • You will rarely use them
  • They take numbers, matrices, or lists as inputs
  • Do not hold information on CRS

Simple feature columns (sfc)

  • A collection of 1 + sfg objects
  • Its contents may be more than one type of geometry (e.g. points and polygons)
  • Can contain a CRS, but must be the same across all geometries

Spherical geometries

What are spherical geometries?

  • Most software used to perform spatial operations assumes planar (e.g. projected) geometries
  • But our world is not flat
  • Therefore, when working with planar geometries, your choice of projection can influence your calculations
  • Using spherical geometries can bypass this
  • Spherical geometries are turned on by default (sf_use_s2(TRUE))

What do spherical geometries do?

S2 geometry assumes that straight lines between points on the globe are not formed by straight lines in the equirectangular projection, but by great circles: the shortest path over the sphere.

mex <- st_point(c(-99.165852, 19.376395)) |>
  st_sfc(crs = 4326)
ams <- st_point(c(4.646219, 52.38739)) |>
  st_sfc(crs = 4326)

# What is the distance between these points?
st_distance(mex, ams)
Units: [m]
        [,1]
[1,] 9206771
sf_use_s2(F) # Turn off spherical geometries
st_distance(mex, ams)
Units: [m]
        [,1]
[1,] 9220537

Why are they different?

One tries to calculate the distance along a straight (red) line. The other does so along the “curved” (black) line. On a sphere, the black line is shorter.

Which one is right?

In reality

Wrap up

Main points

  1. There are two models with which we represent spatial data (vector and raster)
  2. The vector model relies on simple features: points, lines, and polygons (coordinates)
  3. sf is the main R package for working with vector data
  4. an sf object should have three things:
    • attributes (what: data frame)
    • features (where: geometry column)
    • CRS (in relation to what?: in the header)
  5. sf objects can be treated like tibbles / data.frames
  6. Quickly visualize sf objects using plot(), ggplot(), or mapview()

BEFORE Next class

install.packages("sf")
install.packages("rnaturalearth")