+ - 0:00:00
Notes for current slide
Notes for next slide

Spatial data

Maximilian H.K. Hesselbarth

University of Michigan (EEB)

2022/10/24

1 / 19

Spatial data

2 / 19

Spatial data

  • All data that contains spatial information about location
2 / 19

Spatial data

  • All data that contains spatial information about location

  • More complex than just two extra columns (x,y coordinates)

2 / 19

Spatial data

  • All data that contains spatial information about location

  • More complex than just two extra columns (x,y coordinates)

  • Two major data types, namely raster and vector data

2 / 19

Spatial data

  • All data that contains spatial information about location

  • More complex than just two extra columns (x,y coordinates)

  • Two major data types, namely raster and vector data

2 / 19

Spatial data

  • All data that contains spatial information about location

  • More complex than just two extra columns (x,y coordinates)

  • Two major data types, namely raster and vector data

2 / 19

raster data

3 / 19

raster data

  • Cells of equal size represent data
3 / 19

raster data

  • Cells of equal size represent data

  • Closely related to satellite imagery

3 / 19

raster data

  • Cells of equal size represent data

  • Closely related to satellite imagery

  • terra packages provides SpatRaster, SpatVector and SpatExtent classes

    (successor of formerly raster packages)

3 / 19

raster data

  • Cells of equal size represent data

  • Closely related to satellite imagery

  • terra packages provides SpatRaster, SpatVector and SpatExtent classes

    (successor of formerly raster packages)

3 / 19

SpatRaster data

4 / 19

SpatRaster data

  • Multi-layer raster data
4 / 19

SpatRaster data

  • Multi-layer raster data

  • Stores number of columns/rows, spatial extent, and coordinate reference system

4 / 19

SpatRaster data

  • Multi-layer raster data

  • Stores number of columns/rows, spatial extent, and coordinate reference system

  • Values are held in RAM or on disk (if too large for memory)

4 / 19

SpatRaster data

  • Multi-layer raster data

  • Stores number of columns/rows, spatial extent, and coordinate reference system

  • Values are held in RAM or on disk (if too large for memory)

(x <- terra::rast(nrow = 50, ncol = 50, xmin = -50, xmax = 50, ymin = -50, ymax = 50, 
                  crs = NA))
## class       : SpatRaster 
## dimensions  : 50, 50, 1  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. :
4 / 19

SpatRaster data

  • Multi-layer raster data

  • Stores number of columns/rows, spatial extent, and coordinate reference system

  • Values are held in RAM or on disk (if too large for memory)

(x <- terra::rast(nrow = 50, ncol = 50, xmin = -50, xmax = 50, ymin = -50, ymax = 50, 
                  crs = NA))
## class       : SpatRaster 
## dimensions  : 50, 50, 1  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. :
terra::ext(x)
## SpatExtent : -50, 50, -50, 50 (xmin, xmax, ymin, ymax)
terra::crs(x)
## [1] ""
terra::hasValues(x)
## [1] FALSE
4 / 19

SpatRaster data

5 / 19

SpatRaster data

  • Adding values and information
vals <- runif(n = terra::nrow(x) * 
                terra::ncol(x))
terra::values(x) <- vals
x$lyr.2 <- vals * -1
terra::crs(x) <- "epsg:4326"
## class       : SpatRaster 
## dimensions  : 50, 50, 2  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : memory  
##               memory  
## names       :        lyr.1,         lyr.2 
## min values  : 0.0006893259, -0.9999481156 
## max values  : 0.9999481156, -0.0006893259
5 / 19

SpatRaster data

  • Adding values and information
vals <- runif(n = terra::nrow(x) * 
                terra::ncol(x))
terra::values(x) <- vals
x$lyr.2 <- vals * -1
terra::crs(x) <- "epsg:4326"
## class       : SpatRaster 
## dimensions  : 50, 50, 2  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : memory  
##               memory  
## names       :        lyr.1,         lyr.2 
## min values  : 0.0006893259, -0.9999481156 
## max values  : 0.9999481156, -0.0006893259
plot(x)

5 / 19

SpatRaster data

6 / 19

SpatRaster data

  • Access data as matrix or data.frame
x[]
terra::as.data.frame(x, xy = TRUE)
##          lyr.1      lyr.2
## [1,] 0.3201632 -0.3201632
## [2,] 0.5526989 -0.5526989
## [3,] 0.1201288 -0.1201288
## [4,] 0.7382276 -0.7382276
## [5,] 0.8010045 -0.8010045
## [6,] 0.1125737 -0.1125737
##     x  y     lyr.1      lyr.2
## 1 -49 49 0.3201632 -0.3201632
## 2 -47 49 0.5526989 -0.5526989
## 3 -45 49 0.1201288 -0.1201288
## 4 -43 49 0.7382276 -0.7382276
## 5 -41 49 0.8010045 -0.8010045
## 6 -39 49 0.1125737 -0.1125737
6 / 19

SpatRaster data

  • Access data as matrix or data.frame
x[]
terra::as.data.frame(x, xy = TRUE)
##          lyr.1      lyr.2
## [1,] 0.3201632 -0.3201632
## [2,] 0.5526989 -0.5526989
## [3,] 0.1201288 -0.1201288
## [4,] 0.7382276 -0.7382276
## [5,] 0.8010045 -0.8010045
## [6,] 0.1125737 -0.1125737
##     x  y     lyr.1      lyr.2
## 1 -49 49 0.3201632 -0.3201632
## 2 -47 49 0.5526989 -0.5526989
## 3 -45 49 0.1201288 -0.1201288
## 4 -43 49 0.7382276 -0.7382276
## 5 -41 49 0.8010045 -0.8010045
## 6 -39 49 0.1125737 -0.1125737
  • Use mathematical calculations with raster cells
x2 <- x["lyr.1"] + x["lyr.2"]
c(x, x2)
## class       : SpatRaster 
## dimensions  : 50, 50, 3  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : memory  
##               memory  
##               memory  
## names       :        lyr.1,         lyr.2, lyr.1 
## min values  : 0.0006893259, -0.9999481156,     0 
## max values  : 0.9999481156, -0.0006893259,     0
6 / 19

National Land Cover Dataset

7 / 19

National Land Cover Dataset

(nlcd <- terra::rast("data/nlcd_2019.img"))
## class       : SpatRaster 
## dimensions  : 104424, 161190, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : -2493045, 2342655, 177285, 3310005  (xmin, xmax, ymin, ymax)
## coord. ref. : Albers Conical Equal Area 
## source      : nlcd_2019.img 
## color table : 1 
## categories  : NLCD Land Cover Class, Histogram, Red, Green, Blue, Opacity 
## name        :        NLCD Land Cover Class 
## min value   :                 Unclassified 
## max value   : Emergent Herbaceous Wetlands
7 / 19

National Land Cover Dataset

(nlcd <- terra::rast("data/nlcd_2019.img"))
## class       : SpatRaster 
## dimensions  : 104424, 161190, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : -2493045, 2342655, 177285, 3310005  (xmin, xmax, ymin, ymax)
## coord. ref. : Albers Conical Equal Area 
## source      : nlcd_2019.img 
## color table : 1 
## categories  : NLCD Land Cover Class, Histogram, Red, Green, Blue, Opacity 
## name        :        NLCD Land Cover Class 
## min value   :                 Unclassified 
## max value   : Emergent Herbaceous Wetlands

7 / 19

Cropping raster

  • Crop to area of interest
bbox <- terra::ext(-84.023438, -83.320313, 42.032974, 42.553080)
raster_crop <- terra::project(terra::rast(bbox), terra::crs(nlcd))
nlcd_crop <- terra::crop(x = nlcd, y = raster_crop) |>
  terra::project(y = "epsg:4326")

8 / 19

Reclassification

  • Re-classify values using classification matrix
class_matrix <- cbind(c(21, 22, 23, 41, 42, 43), 
                      rep(x = c(1, 2), each = 3))
nlcd_class <- terra::classify(nlcd_crop, rcl = class_matrix, others = NA)

9 / 19

Vector data

10 / 19

Vector data

  • Based on points, lines, and polygons
10 / 19

Vector data

  • Based on points, lines, and polygons

  • Often allows realistic representation of spatial objects

10 / 19

Vector data

  • Based on points, lines, and polygons

  • Often allows realistic representation of spatial objects

  • sf packages provides main classes: POINT, LINESTRING, POLYGON (and MULTI*)

    (successor of formerly sp packages)

10 / 19

Vector data

  • Based on points, lines, and polygons

  • Often allows realistic representation of spatial objects

  • sf packages provides main classes: POINT, LINESTRING, POLYGON (and MULTI*)

    (successor of formerly sp packages)

10 / 19

Read simple features

11 / 19

Read simple features

cities <- sf::st_read("data/cities-mi.shp") |> 
  sf::st_transform(crs = terra::crs(nlcd))
## Reading layer `cities-mi' from data source 
##   `/Users/mhessel/Dropbox (University of Michigan)/06-Talks/Advanced-R-Workshop/slides/spatial/data/cities-mi.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1773 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -90.15955 ymin: 41.71626 xmax: -82.43374 ymax: 48.09634
## Geodetic CRS:  WGS84(DD)
11 / 19

Read simple features

cities <- sf::st_read("data/cities-mi.shp") |> 
  sf::st_transform(crs = terra::crs(nlcd))
## Reading layer `cities-mi' from data source 
##   `/Users/mhessel/Dropbox (University of Michigan)/06-Talks/Advanced-R-Workshop/slides/spatial/data/cities-mi.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1773 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -90.15955 ymin: 41.71626 xmax: -82.43374 ymax: 48.09634
## Geodetic CRS:  WGS84(DD)
states <- sf::st_read("data/provinces.shp") |>
  sf::st_transform(crs = terra::crs(nlcd))
## Reading layer `provinces' from data source 
##   `/Users/mhessel/Dropbox (University of Michigan)/06-Talks/Advanced-R-Workshop/slides/spatial/data/provinces.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 83 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -171.7911 ymin: 18.91619 xmax: -66.96466 ymax: 71.35776
## Geodetic CRS:  WGS 84
11 / 19

simple features

## Simple feature collection with 1773 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 448831.5 ymin: 2121120 xmax: 1096453 ymax: 2810620
## Projected CRS: Albers Conical Equal Area
## First 10 features:
##          name     type     county          gazetteer_ squaremile latitude
## 1     Merritt Township        Bay    Merritt township     31.653 43.51838
## 2      Sciota Township Shiawassee     Sciota township     26.770 42.90760
## 3    Kingston Township    Tuscola   Kingston township     35.843 43.44955
## 4       Grant Township      Mason      Grant township     48.830 44.13276
## 5  Bay De Noc Township      Delta Bay de Noc township     68.145 45.77043
## 6    Sheridan Township      Mason   Sheridan township     35.853 44.05481
## 7     Gilford Township    Tuscola    Gilford township     34.795 43.51942
## 8      Vernon  Village Shiawassee      Vernon village      1.022 42.93992
## 9      Ensign Township      Delta     Ensign township     58.441 45.87256
## 10    Cornell Township      Delta    Cornell township     60.023 45.91859
##    longitude                 geometry
## 1  -83.74868 POINT (983267.6 2344054)
## 2  -84.32242 POINT (945817.9 2270564)
## 3  -83.17502  POINT (1030069 2342515)
## 4  -86.36132 POINT (767460.5 2388076)
## 5  -86.91345 POINT (706236.7 2565161)
## 6  -86.09039 POINT (789859.7 2381650)
## 7  -83.63167 POINT (992590.2 2345386)
## 8  -84.03009 POINT (968917.9 2277091)
## 9  -86.86951 POINT (708553.8 2576777)
## 10 -87.29777 POINT (674947.3 2578747)
12 / 19

Work with simple features

  • Allow to use tidyverse syntax
(regions <- dplyr::filter(states, !name %in% c("Alaska", "Hawaii")) |>
  dplyr::select(name, region) |>
  dplyr::mutate(area_m = as.numeric(sf::st_area(geometry))) |>
  dplyr::group_by(region) |>
  dplyr::summarise(area_m = sum(area_m), area_km = area_m / 1000000) |>
  dplyr::arrange(-area_km))
## Simple feature collection with 4 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2357181 ymin: 311395.7 xmax: 2257301 ymax: 3167163
## Projected CRS: Albers Conical Equal Area
## # A tibble: 4 × 4
##   region     area_m  area_km                                            geometry
##   <chr>       <dbl>    <dbl>                                       <POLYGON [m]>
## 1 West      3.08e12 3075966. ((-1617641 1143814, -1746332 1223770, -1733511 124…
## 2 South     2.33e12 2326982. ((726613.2 841097.8, 654204.3 827934.5, 615792.6 8…
## 3 Midwest   2.13e12 2130036. ((-511400 1906969, -504303.7 2015504, -670478 2027…
## 4 Northeast 4.52e11  451826. ((1714981 2054464, 1708595 2047918, 1704096 203790…
13 / 19

Filter simple features

  • Easy to spatially filter objects
row_id <- sf::st_intersects(x = states, y = dplyr::filter(cities, type == "City"),
                            sparse = FALSE)
state_mi <- states[row_id, ]
14 / 19

Plotting sf objects

  • Can be plotted with ggplot2 geom
ggplot() + 
  geom_sf(data = regions, aes(fill = area_km), alpha = 0.5) +
  geom_sf(data = state_mi, fill = NA, color = "red") +
  scale_fill_viridis() +
  theme_classic()
15 / 19

raster-vector interactions

  • raster and vector data can be used together
states_mainland <- dplyr::filter(states, !name %in% c("Alaska", "Hawaii"))
plot(nlcd)
plot(states_mainland$geometry, add = TRUE)

16 / 19

raster-vector interactions (I)

  • crop cut outs parts using extent

  • mask changes all values outside the mask

nlcd_mi <- terra::crop(x = nlcd, y = state_mi) |> 
  terra::mask(mask = state_mi)

17 / 19

raster-vector interactions (II)

  • extract returns raster values at locations of vector object
nlcd_cities <- terra::extract(x = nlcd_mi, 
                              y = dplyr::filter(cities, type == "Village"), 
                              cells = TRUE)
head(nlcd_cities)
##   ID        NLCD Land Cover Class      cell
## 1  1     Developed, Low Intensity 421920227
## 2  2  Developed, Medium Intensity 417355557
## 3  3     Developed, Low Intensity 370187955
## 4  4 Emergent Herbaceous Wetlands 361982561
## 5  5     Developed, Low Intensity 368848346
## 6  6     Developed, Low Intensity 443220896
18 / 19

Thank you for your attention

Questions?

Further resources: https://mhesselbarth.github.io/advanced-r-workshop/resources

Exercise: https://mhesselbarth.github.io/advanced-r-workshop/exercise-spatial

mhessel@umich.edu

https://www.maxhesselbarth.com

@MHKHesselbarth

19 / 19

Spatial data

2 / 19
Paused

Help

Keyboard shortcuts

↑, ←, Pg Up, k Go to previous slide
↓, →, Pg Dn, Space, j Go to next slide
Home Go to first slide
End Go to last slide
Number + Return Go to specific slide
b / m / f Toggle blackout / mirrored / fullscreen mode
c Clone slideshow
p Toggle presenter mode
t Restart the presentation timer
?, h Toggle this help
Esc Back to slideshow