Data Science II

Session 2: sf Objects and Basic Mapping

Bogdan G. Popescu

bogdan.popescu@johncabot.edu

John Cabot University

What You’ll Learn Today

What an sf object actually is, built from scratch
The three classes: sfg, sfc, and sf
Taking a layer apart — geometry, attributes, coordinates
Turning a table of coordinates into a map — the skill you will use most
Reading and writing the common vector file formats

By the end, you will map 6,460 real health establishments in Yucatán from a CSV.

Where We Left Off

Last week: coordinates, CRS, and reading a shapefile with st_read.

We said an sf object is a dataframe with a geometry column, and left it there.

Today we open it up. Knowing what is inside is what lets you fix it when it breaks.

Before We Start

Same folder shape as last week — code reads ../data/...:

ciencia_datos_2/
├── data/
│   ├── gadm41_MEX_shp/
│   ├── denue/denue_salud_yucatan.csv
│   └── mex_formats/
└── session2/
    └── your_notes.qmd

Everything today runs from these two lines:

library(sf); library(ggplot2); library(dplyr)
sf_use_s2(FALSE)

The complete, ordered code for this session is in session2_code.R.

Part A: What an `sf` Object Is

Vector File Formats

The most common vector files are:

Type	Format	File extension
Binary	ESRI Shapefile	`.shp`, `.shx`, `.dbf`, `.prj`
	GeoPackage (GPKG)	`.gpkg`
Plain Text	GeoJSON	`.json` or `.geojson`
	GPS Exchange Format (GPX)	`.gpx`
	Keyhole Markup Language (KML)	`.kml`
Spatial Databases	PostgreSQL / PostGIS

A Word on `sp`

The first R package for vector data was sp (2005–2023), with rgdal and rgeos.

They dominated spatial analysis in R for fifteen years. rgdal was retired in 2023.

Warning

Do not use them. But recognise them — most Stack Overflow answers you find were written for sp, and code using SpatialPolygonsDataFrame will not run today.

The `sf` Package

sf is the standard package for working with vector data in R.

It relies on three external components: GDAL (reading formats), GEOS (geometry operations), and PROJ (projections).

It contains three classes, and they nest:

sfg — a single geometry
sfc — a geometry column: many sfg plus a CRS
sf — a layer: an sfc inside a dataframe of attributes

Everything today is one of these three.

The `sf` Class, Illustrated

A polygon layer with three features and six non-spatial attributes:

The `sf` Class, Mapped

The same layer, drawn:

Attributes and geometry are two views of one object.

Geometry Types

sfg is a single geometry, and it can be one of several types:

st_point
st_multipoint
st_linestring
st_multilinestring
st_polygon
st_multipolygon
st_geometrycollection

The Seven Simple Feature Types

Building a Point

Let us create one point — the centre of Mérida:

library(sf)
library(ggplot2)
merida <- st_point(c(-89.6237, 20.9674))

Printing it gives the WKT (Well-Known Text) representation:

print(merida)

Warning

Longitude first. c(-89.6237, 20.9674) is x then y. Swap them and Mérida lands in the Indian Ocean.

What Class Is It?

merida <- st_point(c(-89.6237, 20.9674))
class(merida)

[1] "XY"    "POINT" "sfg"

XY — two-dimensional geometry

POINT — the geometry type

sfg — the general class: a simple feature geometry

Note what is missing: there is no CRS. A bare sfg is just numbers.

Building a Polygon

a <- st_polygon(list(cbind(c(0, 0, 7.5, 7.5, 0),
                           c(0, -1, -1, 0, 0))))
print(a)

The first and last coordinate are the same — that is what closes the ring. Omit it and sf errors.

Plotting the Polygon

a <- st_polygon(list(
  cbind(c(0, 0, 7.5, 7.5, 0),
        c(0, -1, -1, 0, 0))))

plot(a)

The Same Polygon with `ggplot`

ggplot() +
  geom_sf(data = a) +
  theme_bw()

geom_sf accepts a bare sfg. You rarely do this in practice — but it shows the geometry is the only thing ggplot needs.

A More Complex Polygon

b <- st_polygon(list(
  cbind(c(0,1,2,3,4,5,6,7,7,0),
        c(1,0,0.5,0,0,0.5,-0.5,-0.5,1,1))))

plot(b)

Same function, more vertices. A municipio boundary is this with a few thousand more.

Combining Geometries

We combine geometries with c():

ab <- c(a, b)
class(ab)

[1] "XY"           "MULTIPOLYGON" "sfg"

The result is a single MULTIPOLYGON — one geometry composed of several shapes.

This is exactly what a state like Baja California Sur is: one feature, many islands.

Plotting the Combination

ab <- c(a, b)
plot(ab)

Intersecting Geometries

New geometries can be calculated from existing ones:

i <- st_intersection(a, b)
class(i)

[1] "XY"                 "GEOMETRYCOLLECTION" "sfg"

Note the difference:

c(a, b) gives a MULTIPOLYGON — everything in either
st_intersection(a, b) gives what is in both

Intersection is the workhorse of Session 5. Here we just meet it.

Combination and Intersection Compared

#Step1: Two panels side by side
par(mfrow = c(1, 2))

#Step2: Everything in either shape, then only what is in both
plot(ab, main = "c(a, b) --- either", col = "#dbe6f0", border = "#33475b")
plot(i,  main = "st_intersection(a, b) --- both", col = "#c2d6e8", border = "#33475b")

Part B: Geometry Columns and Layers

The Geometry Column: `sfc`

A single geometry is rarely useful. We collect them into a geometry column with st_sfc.

Three points in Yucatán:

#Step1: Three bare geometries -- longitude first, then latitude
merida     <- st_point(c(-89.6237, 20.9674))
valladolid <- st_point(c(-88.2020, 20.6896))
progreso   <- st_point(c(-89.6636, 21.2820))

An sfc — unlike an sfg — carries a CRS.

Four Ways to Specify a CRS

an EPSG code — 4326
a crs object, as returned by st_crs()
a PROJ4 string — '+proj=longlat +datum=WGS84 +no_defs'
a WKT string

Use the EPSG code. It is shortest and least error-prone.

Warning

If you omit the crs argument, the CRS is set to NA — undefined. The object still works, still plots, and then fails the moment you combine it with anything else.

Building an `sfc`

combined <- st_sfc(merida, valladolid, progreso, crs = 4326)
combined

Geometry set for 3 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -89.6636 ymin: 20.6896 xmax: -88.202 ymax: 21.282
Geodetic CRS:  WGS 84

Notice what printing gives you now: geometry type, dimension, bounding box, and CRS. None of that existed on the bare sfg.

Plotting the `sfc`

ggplot() +
  geom_sf(data = combined, size = 3) +
  theme_bw()

Three points, no attributes. Useful, but we cannot yet say which is which.

From `sfc` to `sf`: Adding Attributes

#Step1: An ordinary attribute table, one row per geometry
names_df <- data.frame(
  city = c("Mérida", "Valladolid", "Progreso"),
  population = c(995129, 85460, 63400))

#Step2: Weld the attributes to the geometry column
layer <- st_sf(names_df, combined)
layer

Simple feature collection with 3 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -89.6636 ymin: 20.6896 xmax: -88.202 ymax: 21.282
Geodetic CRS:  WGS 84
        city population                 combined
1     Mérida     995129 POINT (-89.6237 20.9674)
2 Valladolid      85460  POINT (-88.202 20.6896)
3   Progreso      63400  POINT (-89.6636 21.282)

st_sf() welds an attribute table to a geometry column. That is all an sf object is.

Mapping with `plot()`

plot(layer)

Base plot() draws one panel per attribute — here, city and population.

Mapping with `ggplot()`

ggplot() +
  geom_sf(data = layer,
          aes(size = population),
          colour = "#33475b") +
  theme_bw()

`plot()` or `ggplot()`?

plot()

Instant, zero typing
Shows every attribute at once
Good for checking data

ggplot()

Full control of colour, size, legend
Layers stack
Good for producing maps

Tip

Use plot() while you work, ggplot() for anything anyone else will see. Getting into the habit of plot()-ing a layer the moment you read it catches most data problems in two seconds.

Part C: Taking a Layer Apart

Extracting Layer Components

Sometimes you need to go the other way and pull out:

the geometry
the attribute table
the coordinates

The Geometry

st_geometry(layer)

Geometry set for 3 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -89.6636 ymin: 20.6896 xmax: -88.202 ymax: 21.282
Geodetic CRS:  WGS 84

This gives you back the sfc — geometry and CRS, no attributes.

The Attribute Table

st_drop_geometry(layer)

        city population
1     Mérida     995129
2 Valladolid      85460
3   Progreso      63400

A plain dataframe. Use this whenever you want to inspect, summarise, or export attributes without dragging thousands of coordinates along.

The Coordinates

st_coordinates(layer)

            X       Y
[1,] -89.6237 20.9674
[2,] -88.2020 20.6896
[3,] -89.6636 21.2820

A matrix. Turn it into a dataframe when you need it as data:

head(data.frame(st_coordinates(layer)), 2)

         X       Y
1 -89.6237 20.9674
2 -88.2020 20.6896

Part D: From a Table to a Map

The Skill You Will Use Most

Most data does not arrive as a shapefile. It arrives as a spreadsheet with a latitude column and a longitude column.

Turning that into a map is one function: st_as_sf().

The Data: DENUE

DENUE — Directorio Estadístico Nacional de Unidades Económicas — is INEGI’s register of every economic establishment in Mexico, with coordinates.

Today’s extract: every health establishment in Yucatán.

6,460 establishments
106 municipios
clinics, dental practices, laboratories, care homes

https://www.inegi.org.mx/app/mapa/denue/

Reading the Table

#Step1: Read the DENUE extract -- an ordinary CSV, no geometry yet
df <- read.csv("../data/denue/denue_salud_yucatan.csv")
dim(df)

[1] 6460   15

names(df)

 [1] "id"         "nom_estab"  "codigo_act" "nombre_act" "per_ocu"   
 [6] "cve_ent"    "entidad"    "cve_mun"    "municipio"  "localidad" 
[11] "telefono"   "www"        "tipoUniEco" "latitud"    "longitud"

The two columns that matter: longitud and latitud.

`st_as_sf()`: Table to Layer

#Step1: Turn the two coordinate columns into a geometry column
#       x first (longitud), then y (latitud); always give a crs
pts <- st_as_sf(df, coords = c("longitud", "latitud"), crs = 4326)
class(pts)

[1] "sf"         "data.frame"

dim(df)

[1] 6460   15

dim(pts)

[1] 6460   14

Fifteen columns became fourteen: longitud and latitud were consumed and replaced by one geometry column.

Warning

coords = c("longitud", "latitud") — x first, then y. And always pass crs = 4326; without it the layer is undefined.

Mapping the Points

ggplot() +
  geom_sf(data = pts, size = 0.2,
          alpha = 0.4) +
  theme_bw()

Six thousand dots. You can see Mérida, and you can guess at the coast — but only if you already know Yucatán.

Points Need Context

#Step1: Read municipio polygons (GADM level 2)
mex <- st_read(
  "../data/gadm41_MEX_shp/gadm41_MEX_2.shp",
  quiet = TRUE)

#Step2: Keep only Yucatán
yuc <- mex[mex$NAME_1 == "Yucatán", ]

#Step3: Polygons first, points on top
ggplot() +
  geom_sf(data = yuc,
          fill = "#eef3f8",
          colour = "#9db2c6",
          linewidth = 0.15) +
  geom_sf(data = pts, size = 0.2,
          alpha = 0.4,
          colour = "#33475b") +
  theme_bw()

Same points. Now the concentration in Mérida, and the emptiness of the east, are visible facts rather than guesses.

`sf` Properties: Rows and Columns

nrow(pts)

[1] 6460

ncol(pts)

[1] 14

dim(pts)

[1] 6460   14

All the base R you know still applies. The geometry column counts as one column.

`sf` Properties: Bounding Box

st_bbox(pts)

     xmin      ymin      xmax      ymax 
-90.40105  20.06068 -87.56142  21.59810

The extreme coordinates: westernmost, southernmost, easternmost, northernmost.

Tip

st_bbox is your first sanity check. If the numbers are not roughly where Yucatán should be (about -90 to -87, 20 to 22), something is wrong — usually swapped coordinates.

`sf` Properties: CRS

st_crs(pts)$input

[1] "EPSG:4326"

st_crs(pts)$units_gdal

[1] "degree"

Degrees — because we set 4326. To measure anything in metres, st_transform(pts, 6372) first.

Subsetting by Attributes

Only the establishments that list a telephone number:

#Step1: Keep only rows that list a telephone number
with_phone <- subset(pts, telefono != "" & !is.na(telefono))

#Step2: How many did we lose?
nrow(pts)

[1] 6460

nrow(with_phone)

[1] 3584

Ordinary subset(). The geometry follows the rows automatically — you never subset the geometry yourself.

Comparing the Subsets

#Step1: patchwork lets us place two ggplots side by side with |
library(patchwork)

Comparing the Subsets

f1 <- ggplot() +
  geom_sf(data = yuc, fill = "#eef3f8", colour = "#9db2c6", linewidth = 0.12) +
  geom_sf(data = pts, size = 0.18, alpha = 0.4, colour = "#33475b") +
  ggtitle(paste0("All establishments (", nrow(pts), ")")) +
  theme_bw(base_size = 9)
f2 <- ggplot() +
  geom_sf(data = yuc, fill = "#eef3f8", colour = "#9db2c6", linewidth = 0.12) +
  geom_sf(data = with_phone, size = 0.18, alpha = 0.4, colour = "#33475b") +
  ggtitle(paste0("With a telephone (", nrow(with_phone), ")")) +
  theme_bw(base_size = 9)
f1 | f2

Just over half. And the ones that drop out are not randomly distributed — a reminder that missing data is itself spatial.

Subsetting Columns

Selecting columns works as it does on any dataframe — with one difference:

small <- pts[, c("nom_estab", "municipio")]
names(small)

[1] "nom_estab" "municipio" "geometry"

You asked for two columns and got three. The geometry column is sticky — it comes along unless you explicitly drop it.

Dropping the Geometry

plain <- st_drop_geometry(small)
names(plain)

[1] "nom_estab" "municipio"

class(plain)

[1] "data.frame"

Now it is an ordinary dataframe. Do this before write.csv, or before any operation that does not need the geometry — it is dramatically faster.

Part E: Reading and Writing Files

Reading a Shapefile

states_shp <- st_read("../data/gadm41_MEX_shp/gadm41_MEX_1.shp", quiet = TRUE)
nrow(states_shp)

[1] 32

Remember: a shapefile is four files minimum. st_read points at the .shp, but silently needs the others beside it.

Reading a GeoJSON

states_json <- st_read("../data/mex_formats/mex_states.geojson", quiet = TRUE)
nrow(states_json)

[1] 32

One file, plain text. Open it in a text editor and you can read the coordinates.

Reading a GeoPackage

states_gpkg <- st_read("../data/mex_formats/mex_states.gpkg", quiet = TRUE)
nrow(states_gpkg)

[1] 32

One file, binary, and it can hold many layers at once. st_layers() tells you what is inside.

Which Format Should You Use?

	Use it when
Shapefile	someone insists — it is the old standard
GeoJSON	small data, web maps, version control
GeoPackage	your default — one file, no size limit, multiple layers

Warning

Shapefile limitations that will bite you: column names truncated to 10 characters, a 2 GB size cap, and no proper support for accented characters. Nombre_del_Municipio becomes Nombre_del.

Writing Files

st_write(with_phone, "salud_yucatan.gpkg", delete_dsn = TRUE)
st_write(with_phone, "salud_yucatan.geojson", delete_dsn = TRUE)

st_write picks the format from the file extension. delete_dsn = TRUE overwrites; without it, a second run errors.

What Have We Learned?

An sf object is three nested things: sfg → sfc → sf
A bare geometry has no CRS; an sfc does
st_sf() welds attributes to geometry — that is the whole trick
st_geometry, st_drop_geometry, st_coordinates take a layer apart
st_as_sf(df, coords = c("lon", "lat"), crs = 4326) turns any table into a map
The geometry column is sticky — it survives subsetting until you drop it
GeoPackage should be your default format, not shapefile

Next session: choropleths — joining data to boundaries, classification, colour, and interactive maps with leaflet.

Data Science II

What You’ll Learn Today

Where We Left Off

Before We Start

Part A: What an sf Object Is

Vector File Formats

A Word on sp

The sf Package

The sf Class, Illustrated

The sf Class, Mapped

Geometry Types

The Seven Simple Feature Types

The Seven Simple Feature Types

The Seven Simple Feature Types

The Seven Simple Feature Types

The Seven Simple Feature Types

The Seven Simple Feature Types

The Seven Simple Feature Types

Building a Point

What Class Is It?

Building a Polygon

Plotting the Polygon

The Same Polygon with ggplot

A More Complex Polygon

Combining Geometries

Plotting the Combination

Intersecting Geometries

Combination and Intersection Compared

Part B: Geometry Columns and Layers

The Geometry Column: sfc

Four Ways to Specify a CRS

Building an sfc

Plotting the sfc

From sfc to sf: Adding Attributes

Mapping with plot()

Mapping with ggplot()

plot() or ggplot()?

Part C: Taking a Layer Apart

Extracting Layer Components

The Geometry

The Attribute Table

The Coordinates

Part D: From a Table to a Map

The Skill You Will Use Most

The Data: DENUE

Reading the Table

st_as_sf(): Table to Layer

Mapping the Points

Points Need Context

sf Properties: Rows and Columns

sf Properties: Bounding Box

sf Properties: CRS

Subsetting by Attributes

Comparing the Subsets

Comparing the Subsets

Subsetting Columns

Dropping the Geometry

Part E: Reading and Writing Files

Reading a Shapefile

Reading a GeoJSON

Reading a GeoPackage

Which Format Should You Use?

Writing Files

What Have We Learned?

Part A: What an `sf` Object Is

A Word on `sp`

The `sf` Package

The `sf` Class, Illustrated

The `sf` Class, Mapped

The Same Polygon with `ggplot`

The Geometry Column: `sfc`

Building an `sfc`

Plotting the `sfc`

From `sfc` to `sf`: Adding Attributes

Mapping with `plot()`

Mapping with `ggplot()`

`plot()` or `ggplot()`?

`st_as_sf()`: Table to Layer

`sf` Properties: Rows and Columns

`sf` Properties: Bounding Box

`sf` Properties: CRS