class: center, middle, inverse, title-slide .title[ # Spatial Data in R ] .subtitle[ ## vectors, rasters, and maps ] .author[ ### Elie Gurarie ] .date[ ### 2026-04-20 ] --- ## Key packages .pull-left[ ### Vector data - **`sf`** — *simple features* - points, lines, polygons - all functions prefixed `st_` - reads (GIS) shape files ] .pull-right[ ### 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 > .footnotesize[ Older vector packages (`sp`, `rgdal`, `rgeos`) are **retired** as of 2023, but will come up often ... esp. as dependencies in other packages ] --- class: inverse, center, middle # Simple Features (`sf`) --- ## What is a "simple feature"? .pull-left[ 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*). ] .pull-right[ <br> <img src='https://r-spatial.github.io/sf/articles/sf_xfig.png' width='100%'/> .footnotesize[*source: r-spatial.github.io/sf*] ] --- .pull-left[ ## Data frame to `sf` Four Syracuse landmarks: .small[ ``` r 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 ``` ] ] .pull-right[ Key function: `st_as_sf`: ``` r 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!) ``` r 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) ``` .darkred[Note: the **`geometry`** list-column.] ] --- ## `st_geometry` and `st_coordinates` .pull-left[ **`st_geometry()`** — extract just the geometry column (a list of POINTs, a.k.a. an `sfc`): ``` r 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 ``` ] .pull-right[ **`st_coordinates()`** — get a plain numeric matrix of coordinates, handy for base R or for passing to non-spatial functions: ``` r 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. ] --- .pull-left-50[ ## Base plot .small[ ``` r 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) ``` <img src="spatial-data_files/figure-html/unnamed-chunk-8-1.png" alt="" style="display: block; margin: auto;" /> ]] .pull-right[ ### `ggplot` (`geom_sf`) .small[ ``` r require(ggplot2) ggplot(syr_sf) + geom_sf(color = "darkred", size = 3) + geom_sf_text(aes(label = name), nudge_y = 0.001, size = 3) ``` <img src="spatial-data_files/figure-html/unnamed-chunk-9-1.png" alt="" width="80%" style="display: block; margin: auto;" /> ] `geom_sf()` auto-handles the CRS and aspect ratio. ] --- .pull-left-40[ ## `mapview` One line, fully interactive (pan, zoom, basemap, popups): ``` r require(mapview) mapview(syr_sf, label = "name", col.regions = "darkred", cex = 6) ``` .red[Works inline in any RMarkdown] ] .pull-right-50[
] --- class: inverse, center, middle # Rasters --- ## What is a raster? .pull-left[ A **raster** is a regular grid of cells, each holding a value. - **Continuous** — elevation, temperature, NDVI - **Categorical** — land cover, soil type, habitat class ] .pull-right[ 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()` --- .pull-left-60[ ## Loading and plotting ``` r 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. ] .pull-right-40[ ``` r plot(syr_nlcd, plg = list(cex = 0.6)) ``` <img src="spatial-data_files/figure-html/unnamed-chunk-13-1.png" alt="" width="95%" style="display: block; margin: auto;" /> ] .footnotesize[**National Land Cover Database** — 30 m categorical land cover for the U.S., updated every 2–3 years.] --- .pull-left[ ## Interactive view with `mapview` ``` r 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. ] .pull-right[
] --- ## Extracting values: align CRS first .pull-left-50[ `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**. ``` r 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]] ``` ``` r crs(syr_nlcd, describe = TRUE)$name ``` ``` ## [1] "NAD83 / Conus Albers" ``` Transform the points into the raster's CRS: ``` r syr_sf <- st_transform(syr_sf, crs(syr_nlcd)) ``` Coordinates are now in meters (Albers), matching the raster. ] .pull-right-50[ Extract land cover at each point: ``` r 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 ``` ] --- .pull-left-40[ ### Overlay on the raster: ``` r plot(syr_nlcd, plg = list(cex = 0.6)) plot(st_geometry(syr_sf), add = TRUE, pch = 21, col = "yellow") ``` ] .pull-right-60[ <img src="spatial-data_files/figure-html/unnamed-chunk-20-1.png" alt="" style="display: block; margin: auto;" /> ] --- class: inverse ## Some spatial take-aways .box-green[ - Easy and powerful and fast - When you combine vector data with raster data ... ALWAYS transform the *vector*, never the *raster* ]