Disclaimer: The contents of this document come from Chapter 5. Geometry operations of Geocomputation with R(Lovelace, Nowosad, & Muenchow, 2019). This document is prepared for CP6521 Advanced GIS, a graduate-level city planning elective course at Georgia Tech in Spring 2019. For any question, contact the instructor, Yongsung Lee, Ph.D. via yongsung.lee(at)gatech.edu.
This document is also published on RPubs.
install.packages("tidyverse", dependencies = TRUE)
install.packages("stringr", dependencies = TRUE)
install.packages("sf", dependencies = TRUE)
install.packages("raster", dependencies = TRUE)
install.packages("spData", dependencies = TRUE)
install.packages("devtools") # for this, you need Rtools installed on your machine 
devtools::install_github("Nowosad/spDataLarge")
# install.packages("sp", dependencies = TRUE)
# install.packages("mapview", dependencies = TRUE)
# install.packages("tmap", dependencies = TRUE)
install.packages("geosphere", dependencies = TRUE)
library(tidyverse)
library(stringr)     # for working with strings (pattern matching)
library(sf)          # classes and functions for vector data
library(raster)      # classes and functions for raster data
library(spData)      # load geographic data
library(spDataLarge) # load larger geographic data
# library(sp)
# library(mapview)
# library(tmap)
library(geosphere)

Learning Objectives

What we do:
1. Learn geometry operations on raster data
- Spatial intersection between rasters
- Incraese/decrease raster resolutions

  1. Learn how to use raster and vector together

5.3 Geometric operations on raster data

5.3.1 Geometric intersections
data("elev", package = "spData")
clip = raster(xmn = 0.9, xmx = 1.8, ymn = -0.45, ymx = 0.45, 
              res =0.3, vals = rep(1, 9))
elev[clip, drop = FALSE]
5.3.2 Extent and origin
data(elev, package = "spData") 
elev_2 = extend(elev, c(1, 2), value = 1000) # add one row and two columns to each side of the raster
plot(elev_2)

What if the extents of two rasters differ? R will return an error.

elev_3 = elev + elev_2 
# Raster objects have different extents. Result for their intersection is returned
# extend the first raster to be fitted to the second raster 
# note that the default value for new raster cells are NA here. 
elev_4 = extend(elev, elev_2)

Change the origin of a raster (the point in a raster, which is the closest to (0.0) based on the size of its resolution. That is, changing the resoltion of a raster may change its origin.)

origin(elev_4)
#> [1] 0 0
# change the origin
origin(elev_4) = c(0.25, 0.25)
plot(elev_4)
# and add the original raster
plot(elev, add = TRUE)
5.3.3. Aggregation and disaggregation

To match resolutions of rasters,
- decrease by aggregate()
- increase by disaggregate()

By default, output values are the mean of input values

data("dem", package = "RQGIS")
dem_agg = aggregate(dem, fact = 5, fun = mean)

disaggregate() requires a method that computes the values of smaller cells. By default, method="" returns values taken from immediate neighbors of those cells. method = "bilinear" interpolates using the values of four nearest (larger) cells and their distances to the new (smaller) cell.

dem_disagg = disaggregate(dem_agg, fact = 5, method = "bilinear")
identical(dem, dem_disagg)

resample() computes values for new pixel locations, and by default it uses bilinear interpolation.

# add 2 rows and columns, i.e. change the extent
dem_agg = extend(dem_agg, 2)
plot(dem)
plot(dem_agg)

# resample() takes the first argeument as an input and 
# returns a raster with the same extent and resolution as the second argument 
# below, the added rows/columns are removed after resampling 
dem_disagg_2 = resample(dem_agg, dem)

5.4 Raster-vector interactions

5.4.1 Raster cropping

crop() returns a raster in a new (smaller) extent defined by a cropping object.

srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))   # target object 
zion = st_read(system.file("vector/zion.gpkg", package = "spDataLarge")) # cropping object 
zion = st_transform(zion, projection(srtm))
srtm_cropped = crop(srtm, zion) # the extent becomes slightly smaller

mask() returns NA for those raster cells outside of the cropping object.

srtm_masked = mask(srtm, zion)
srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
5.4.2 Raster extraction

raster::extract() takes values from a ‘target’ raster for a given geographic ‘selector’ (typically vector: e.g., points, lines, or polygons).

data("zion_points", package = "spDataLarge") # 30 random points in the Zion National Park 
zion_points$elevation = raster::extract(srtm, zion_points) # create a new column 
summary(zion_points$elevation)

buffer argument generates a list of vectors. Each vector in the list contains raster values (elevation) inside the buffer of a given point.

raster::extract(srtm, zion_points, buffer = 1000)

extract() works with lines.

# create a line 
zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
  st_linestring() %>%
  st_sfc(crs = projection(srtm)) %>%
  st_sf()

# extract elevation along the line 
# returns a list containing a sinlge matrix: first column - cell id, second column - elevation 
# if multiple lines, then it will return a list of multiple matrices, one list element for each line 
transect = raster::extract(srtm, zion_transect, 
                           along = TRUE, cellnumbers = TRUE)
as.data.frame(transect[[1]])

purrr::map_dfr(): https://purrr.tidyverse.org/reference/map.html

transect_df = purrr::map_dfr(transect, as_data_frame, .id="ID") 
transect_coords = xyFromCell(srtm, transect_df$cell) 
transect_df$dist = c(0, cumsum(geosphere::distGeo(transect_coords)))

An example of raster::extract() work with a polygon. This helps generate summary statistics for raster cells, grouped by polygons (boundaries).

zion_srtm_values = raster::extract(x = srtm, y = zion, df = TRUE)
group_by(zion_srtm_values, ID) %>% 
  summarize_at(vars(srtm), list(~min, ~mean, ~max))

Besides continuous values such as elevation, raster::extract() also handles categorical values by counting the number of cells in each category given a geographic selector.

zion_nlcd <- raster::extract(nlcd, zion, df = TRUE, factors = TRUE)
dplyr::select(zion_nlcd, ID, levels) %>%
  tidyr::gather(key, value, -ID) %>%
  group_by(ID, key, value) %>%
  tally() %>%
  tidyr::spread(value, n, fill = 0)
5.4.3 Rasterization

Conversion of vector objects into their represeation in raster objects.
- the first argument of rasterize() is a vector object
- the seond argument is a raster object that works as a template for extent, resolution, crs, etc.

cycle_hire_osm_projected = st_transform(cycle_hire_osm, 27700) 
raster_template = raster(extent(cycle_hire_osm_projected), resolution = 1000, 
                         crs = st_crs(cycle_hire_osm_projected)$proj4string)

First, create a present/absence rater. filed defines values for all non-empty cells.

ch_raster1 = rasterize(cycle_hire_osm_projected, raster_template, field = 1)

fun specifies summary statistics, which convert multiple observations in close proximity into associated cells in the raster object: e.g., fun = "count" counts the number of cycle hire points in each grid cell here.

ch_raster2 <- rasterize(cycle_hire_osm_projected, raster_template, 
                        field = 1, fun = "count")

By specifying field and fun, we can create a raster presenting the capacity of cycle hire points across London.

ch_raster3 <- rasterize(cycle_hire_osm_projected, raster_template, 
                        field = "capacity", fun = sum)
5.4.4 Spatial vectorization

Conversion of spatialy continuous raster data into spatially discrete vector data such as points, lines, or polygons.

First, raster to points

elev_point <- rasterToPoints(elev, spatial = TRUE) %>%
  st_as_sf()

plot(elev)
plot(elev_point, add = TRUE)

Second, raster to lines

data(dem, package = "RQGIS")
cl = rasterToContour(dem)
plot(dem, axes = FALSE)
plot(cl, add = TRUE)

Better visualization options

# create hillshade
hs = hillShade(slope = terrain(dem, "slope"), aspect = terrain(dem, "aspect"))
plot(hs, col = gray(0:100 / 100), legend = FALSE)
# overlay with DEM
plot(dem, col = terrain.colors(25), alpha = 0.5, legend = FALSE, add = TRUE)
# add contour lines
contour(dem, col = "white", add = TRUE)

Third, raster to polygons

grain_poly <- rasterToPolygons(grain) %>% # each cell converted to separate polygons 
  st_as_sf()
plot(grain)
plot(grain_poly)

grain_poly2 <- grain_poly %>% # cells are "dissolved"; in total there are three multipolygon features 
  group_by(layer) %>%
  summarize()
plot(grain_poly2)

5.5 Exercises

Some of the exercises use a raster dataset (ndvi) from the RQGIS package. They also use a polygonal ‘convex hull’ derived from the vector dataset (ch) to represent the area of interest:

install.packages("RQGIS", dependencies = TRUE)
library(RQGIS)
data(random_points) #for questions before 7. 
data(ndvi)
ch = st_combine(random_points) %>% 
  st_convex_hull()
  1. Crop the ndvi raster using (1) the random_points dataset and (2) the ch dataset. Are there any differences in the output maps? Next, mask ndvi using these two datasets. Can you see any difference now? How can you explain that?
plot(ndvi)
plot(st_geometry(random_points), add = TRUE)
plot(ch, add = TRUE)

  1. Firstly, extract values from ndvi at the points represented in random_points. Next, extract average values of ndvi using a 90 buffer around each point from random_points and compare these two sets of values. When would extracting values by buffers be more suitable than by points alone?

  2. Subset points higher than 3100 meters in New Zealand (the nz_height object) and create a template raster with a resolution of 3 km. Using these objects:
    • Count numbers of the highest points in each grid cell.
    • Find the maximum elevation in each grid cell.
nz_height3100 <- nz_height %>% filter(elevation>3100)
new_graticule <- st_graticule(nz_height3100, datum = 2193)
plot(st_geometry(nz_height3100), graticule = new_graticule, axes = TRUE)

  1. Aggregate the raster counting high points in New Zealand (created in the previous exercise), reduce its geographic resolution by half (so cells are 6 by 6 km) and plot the result.
    Resample the lower resolution raster back to a resolution of 3 km. How have the results changed?
    Name two advantages and disadvantages of reducing raster resolution.

  2. Polygonize the grain dataset and filter all squares representing clay.
    • Name two advantages and disadvantages of vector data over raster data.
    • At which points would it be useful to convert rasters to vectors in your work?