Key packages

Vector data

  • sf — simple features
    • points, lines, polygons
    • all functions prefixed st_
    • reads (GIS) shape files

Raster data

  • terra — modern raster (and vector)
    • successor to raster
    • faster, actively developed
  • raster — older, still widely used, still has features not available in terra

Visualization

  • mapview — quick, interactive leaflet maps (one-liner)
  • ggplot2::geom_sf() — static maps in the grammar of graphics
  • maps — base R country/state outlines

Older vector packages (sp, rgdal, rgeos) are retired as of 2023, but will come up often … especially as dependencies in other packages.

Simple Features (sf)

What is a “simple feature”?

A simple feature is a “standard” for 2D vector geometry. Every sf object is a data frame with a geometry column:

  • POINT — e.g., a GPS fix, a nest, a city
  • LINESTRING — a track, a river, a road
  • POLYGON — a home range, a lake, a park
  • … plus MULTI* versions

Each row is one feature; the geometry column holds its coordinates. Every sf function starts with st_ (spatial / temporal).

source: r-spatial.github.io/sf

Data frame to sf

Four Syracuse landmarks:

syr <- data.frame(
  name = c("Great Law of Peace Center", "SUNY-ESF", "Barry Park", "Onondaga Lake"),
  lon  = c(-76.1967, -76.1356, -76.1149, -76.2175),
  lat  = c( 43.0936,  43.0354, 43.0267,  43.0967)
)
syr
##                        name      lon     lat
## 1 Great Law of Peace Center -76.1967 43.0936
## 2                  SUNY-ESF -76.1356 43.0354
## 3                Barry Park -76.1149 43.0267
## 4             Onondaga Lake -76.2175 43.0967

The key function is st_as_sf:

require(sf)
syr_sf <- st_as_sf(syr, coords = c("lon", "lat"), crs = 4326)

crs = 4326 → WGS84 lon/lat (default for GPS data). coords = specifies the column names with lon-lat (NOT lat-lon!).

syr_sf
## Simple feature collection with 4 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.2175 ymin: 43.0267 xmax: -76.1149 ymax: 43.0967
## Geodetic CRS:  WGS 84
##                        name                 geometry
## 1 Great Law of Peace Center POINT (-76.1967 43.0936)
## 2                  SUNY-ESF POINT (-76.1356 43.0354)
## 3                Barry Park POINT (-76.1149 43.0267)
## 4             Onondaga Lake POINT (-76.2175 43.0967)

Note the geometry list-column.

st_geometry and st_coordinates

st_geometry() extracts just the geometry column (a list of POINTs, a.k.a. an sfc):

st_geometry(syr_sf)
## Geometry set for 4 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -76.2175 ymin: 43.0267 xmax: -76.1149 ymax: 43.0967
## Geodetic CRS:  WGS 84

st_coordinates() returns a plain numeric matrix of coordinates, handy for base R or for passing to non-spatial functions:

st_coordinates(syr_sf)
##             X       Y
## [1,] -76.1967 43.0936
## [2,] -76.1356 43.0354
## [3,] -76.1149 43.0267
## [4,] -76.2175 43.0967

Compare: st_drop_geometry() goes the other way — returns the attribute table only.

Base plot

plot(st_geometry(syr_sf), pch = 19, col = "darkred", cex = 2)
text(st_coordinates(syr_sf), labels = syr_sf$name, pos = 4, cex = 0.8)

ggplot (geom_sf)

require(ggplot2)
ggplot(syr_sf) + 
  geom_sf(color = "darkred", size = 3) +
  geom_sf_text(aes(label = name), nudge_y = 0.001, size = 3)

geom_sf() auto-handles the CRS and aspect ratio.

mapview

One line, fully interactive (pan, zoom, basemap, popups):

require(mapview)
mapview(syr_sf, label = "name", col.regions = "darkred", cex = 6)

Works inline in any RMarkdown.

Rasters

What is a raster?

A raster is a regular grid of cells, each holding a value.

  • Continuous — elevation, temperature, NDVI
  • Categorical — land cover, soil type, habitat class

Defined by:

  • extent (xmin, xmax, ymin, ymax)
  • resolution (cell size)
  • CRS (projection)
  • values (one or more layers / bands)

terra vs. raster

  • terra is the modern replacement — faster, C++ backend, unified vector + raster
  • Core class is SpatRaster (vs. RasterLayer in raster)
  • Most function names carry over: rast(), crop(), extract(), project(), plot()

Loading and plotting

require(terra)
syr_nlcd <- rast("syr_nlcd.tif")
syr_nlcd
## class       : SpatRaster 
## size        : 802, 548, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 1588575, 1605015, 2383635, 2407695  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83 / Conus Albers (EPSG:5070) 
## source      : syr_nlcd.tif 
## color table : 1 
## categories  : NLCD Land Cover Class, Histogram, Red, Green, Blue, Opacity 
## name        :         NLCD Land Cover Class 
## min value   :                    Open Water 
## max value   : Emergent Herbaceuous Wetlands

Note the CRS: Albers Equal Area Conic / NAD83 — the native NLCD projection (units: meters), not lon/lat.

plot(syr_nlcd, plg = list(cex = 0.6))

National Land Cover Database — 30 m categorical land cover for the U.S., updated every 2–3 years.

Interactive view with mapview

require(mapview)
mapview(syr_nlcd)

mapview reprojects on the fly to web-mercator (EPSG:3857) and overlays on an interactive basemap — good for sanity-checking extent and alignment.

Extracting values: align CRS first

terra::extract() pulls raster values at vector locations — but the CRS must match. Our points are in WGS84 lon/lat (EPSG:4326); the raster is in Albers NAD83.

st_crs(syr_sf)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         MEMBER["World Geodetic System 1984 (G2296)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
crs(syr_nlcd, describe = TRUE)$name
## [1] "NAD83 / Conus Albers"

Transform the points into the raster’s CRS:

syr_sf <- st_transform(syr_sf, crs(syr_nlcd))

Coordinates are now in meters (Albers), matching the raster.

Extract land cover at each point:

syr_sf$nlcd <- extract(syr_nlcd, syr_sf)[,2]
syr_sf
## Simple feature collection with 4 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 1590369 ymin: 2393706 xmax: 1600121 ymax: 2399643
## Projected CRS: NAD83 / Conus Albers
##                        name                geometry                      nlcd
## 1 Great Law of Peace Center POINT (1592088 2399643)  Developed, Low Intensity
## 2                  SUNY-ESF POINT (1598278 2394308) Developed, High Intensity
## 3                Barry Park POINT (1600121 2393706)          Deciduous Forest
## 4             Onondaga Lake POINT (1590369 2399633)                Open Water

Overlay on the raster

plot(syr_nlcd, plg = list(cex = 0.6))
plot(st_geometry(syr_sf), add = TRUE, pch = 21, bg = "yellow")

Some spatial take-aways

  • Easy and powerful and fast.
  • When you combine vector data with raster data … ALWAYS transform the vector, never the raster.