Today

  • terra (raster) intro
    • classes
    • dummy rasters
    • reading and writing
    • agg/disagg, masking
    • assigning values
    • zonal
    • focal
    • global
    • terrain

terra objects

  • SpatRaster object - 1 layer to many layers, from multiple or single file sources
  • SpatVector object - Vector data

accessing bands

  • use list [[ ... ]] notation
chirps <- rast(system.file("extdata/chirps.tif", package = "geospaar"))
chirps[[3]] ## access by index
#chirps[["Y16301"]] ## access by band name
plot(chirps[[3]])

Creating raster from scratch

  • how does extent work?
library(terra)
r <- rast(ext(30, 31, -14, -13), resolution = 0.1, crs = "EPSG:4326")
values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
# values(r) <- 1:15
# plot(r)

par(mar = c(0, 0, 0, 0))
plot(districts %>% st_geometry)
plot(r, add = TRUE)

# plot(vect(districts))
# plot(districts)
  • Creating multi-layer stack
s <- c(r, log10(r))
names(s) <- c("dummy", "log10dummy")
plot(s)
print(class(s))
  • Creating multi-layer brick
  • what are the dimensions of this raster? resolution?
  • how are the values filled?
b2 <- lapply(1:10, function(x) {
  r <- rast(ext(30, 31, -14, -13), resolution = 0.1, 
            crs = "EPSG:4326")
  set.seed(x)
  values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
  r
}) %>% do.call(c, .)
names(b2) <- paste0("b", 1:10)

plot(chirps[[1:10]])
# plot(b2)

A large brick

  • how does this differ from above?
b3 <- lapply(1:10, function(x) {
  r <- rast(ext(30, 31, -14, -13), resolution = 0.001, 
              crs = "EPSG:4326")
  set.seed(x)
  values(r) <- sample(1:10, size = ncell(r), replace = TRUE)
  r
}) %>% do.call(c, .)
names(b3) <- paste0("b", 1:10)
plot(b3)
plot(b3[[1]])
plot(b3[[10]])
plot(b3[[c(1, 10)]])

read and write

## write raster
writeRaster(r, filename = file.path(tempdir(), "mydummy.tif"), overwrite = TRUE)
dir(tempdir())

## read in raster
r <- rast(file.path(tempdir(), "mydummy.tif"))
plot(r)

read and write

  • brick function can be used both to read and write
## write file
writeRaster(s, filename = file.path(tempdir(), "mybrick.tif"), overwrite = TRUE)

## read file
newbrick <- rast(file.path(tempdir(), "mybrick.tif"))

Projecting

  • project
  • similar to st_transform(), use crs of existing layer to project
  • need to define resolution with res argument
roads <- read_sf(system.file("extdata/roads.geojson", package = "geospaar"))
chirpsz_alb <- project(x = chirps, y = crs(roads), res = 5000) 
chirpsz_alb %>% print()
plot(chirpsz_alb[[1:10]])

Vector <–> Raster

  • rasterize converts to raster. Need to pass argument that gives raster dimensions.
  • In below example, we use chirpsz for dimensions
distsr <- rasterize(districts, chirpsz, field = "distName")
distsr %>% plot_noaxes
str(distsr)
# stack(distsr, chirpsz[[1:2]]) %>% plotRGB(stretch = "lin")

Vector <–> Raster

  • rasterToPolygons, use dissolve to merge common raster values
distsr_pol <- as.polygons(distsr, dissolve = TRUE)
distsr_pol %>% 
  st_as_sf %>% 
  st_geometry %>% 
  plot
distsr_pol %>% 
  st_as_sf %>% 
  slice(49) %>% 
  plot(add = TRUE)

Plotting

  • run examples below individually
plot(b)
plot(b2)
plot_noaxes(b)
plot_noaxes(b2, main = paste("Random", 1:10))
plot(b2, main = paste("Random", 1:10))
plot_noaxes(b, legend = FALSE)
plot_noaxes(b[[1]], legend = FALSE)
legend("right", legend = 1:10, fill = terrain.colors(10), border = FALSE, 
       bty = "n")
# plotting with larger datasets
plot_noaxes(b3, main = paste("Random", 1:10))
plot_noaxes(b3, nc = 3, nr = 4, main = paste("Random", 1:10))
  • rasterVis package includes additional plot options
rasterVis::levelplot(b3)#, #names.attr = paste("Random", 1:10))
rasterVis::levelplot(b2, names.attr = paste("Random", 1:10))
rasterVis::levelplot(b[[1]])
rasterVis::gplot(b[[1]]) + geom_tile(aes(fill = value)) + 
  scale_fill_viridis_c()
  • ggplot() works with geom_raster()
ggplot2::ggplot(as.data.frame(b[[1]], xy = TRUE)) + 
  # geom_tile(aes(x = x, y = y, fill = dummy)) + scale_fill_viridis_c()  
  geom_raster(aes(x = x, y = y, fill = dummy)) + scale_fill_viridis_c()
  • stars package
stars::st_as_stars(b) %>% plot(col = viridis::viridis(10))
stars::st_as_stars(b) %>% plot()

Raster to other types

plot(b)
b[1:10]
b[1:ncell(b)]
as.matrix(b)
as.data.frame(b) %>% as_tibble()

Pre-processing

  • Aggregating/disaggregating
  • what is default method for aggregating?
plot(b, main = paste0("Brick ", 1:2) )
aggregate(b, fact = 2) %>% 
  plot(., main = paste0("Brick agg", 1:2) )
# aggregate(b, fun = min, fact = 2) %>% plot
# aggregate(b, fun = sd, fact = 2) %>% plot

Pre-processing

  • what is default method for disaggregating?
plot(b, main = paste0("Brick ", 1:2) )
disaggregate(b, fact = 2) %>% 
  plot(., main = paste0("Brick disagg", 1:2) )
disaggregate(b, fact = 2, method = 'bilinear') %>% 
  plot(., main = paste0("Brick disagg bilinear", 1:2) )
  • Masking
  • mask makes values NA outside of mask area.
library(geospaar)
data(chirps)
plot_noaxes(chirps[[1]])
chirpsz <- mask(chirps, districts)
plot_noaxes(chirpsz[[1]])
plot(st_geometry(districts), add = TRUE)
# rasterVis::levelplot(chirpsz[[1:5]])