introduction

In this lab, we’ll explore the basics of working with raster data, including attribute, spatial, and geometry operations. This lab follows chapters 3, 4, and 5 of Geocomputation with R by Robin Lovelace.

prerequisites

library(terra)
## terra 1.7.39
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(spData)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
library(spDataLarge)
## Warning: package 'spDataLarge' was built under R version 4.3.1
library(tmap)
## 
## Attaching package: 'tmap'
## The following object is masked from 'package:datasets':
## 
##     rivers
library(geodata)

manipulating raster objects

Raster data represents continuous surfaces, as opposed to the discrete features represented in the vector data model. Here we’ll learn how to create raster data objects from scratch and how to do basic data manipulations.

Let’s create a SpatRaster object using a digitial elevation model for Zion National Park.

raster_filepath <- system.file("raster/srtm.tif", package = "spDataLarge") #establish file path
my_rast <- rast(raster_filepath) # create raster object

class(my_rast) # test class of raster object
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
my_rast # gives summary information
## class       : SpatRaster 
## dimensions  : 457, 465, 1  (nrow, ncol, nlyr)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : srtm.tif 
## name        : srtm 
## min value   : 1024 
## max value   : 2892
plot(my_rast)

We can also create rasters from scratch using the rast() function. Here we create 36 cells centerd around (0,0). By default the CRS is set to WGS84, but we could change this with the “crs” argument. Because we are working in WGS84, the resolution is in units of degrees. rast() fills the values of the cells row-wise starting in the upper left corner.

new_raster = rast(nrows = 6, ncols = 6, resolution = 0.5,
                  xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
                  vals = 1:36)

tm_shape(new_raster) +
  tm_raster()

The SpatRaster class can also handle multiple layers.

multi_raster_file <- system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast <- rast(multi_raster_file)
multi_rast
## class       : SpatRaster 
## dimensions  : 1428, 1128, 4  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 12N (EPSG:32612) 
## source      : landsat.tif 
## names       : landsat_1, landsat_2, landsat_3, landsat_4 
## min values  :      7550,      6404,      5678,      5252 
## max values  :     19071,     22051,     25780,     31961
nlyr(multi_rast) # test number of layers in raster object
## [1] 4

We can subset layers using either the layer number or name

multi_rast3 <- subset(multi_rast, 3)
multi_rast4 <- subset(multi_rast, "landsat_4")

We can combine SpatRaster objects into one, using the c function

multi_rast34 <- c(multi_rast3, multi_rast4)

Let’s create an example raster for elevation

elev <- rast(nrows = 6, ncols = 6, resolution = 0.5,
            xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
            vals = 1:36)

plot(elev)

Rasters can also hold categorical data. Let’s create an example raster for soil types.

grain_order <- c("clay", "silt", "sand") # set soil types
grain_char <- sample(grain_order, 36, replace = TRUE) # randomly create character string of soil types
grain_fact <- factor(grain_char, levels = grain_order) # convert to factors

grain <- rast(nrows = 6, ncols = 6, resolution = 0.5,
             xmin = -1.5, xmax = 1.5, ymin = -1.5, ymax = 1.5,
             vals = grain_fact)

plot(grain)

raster subsetting

We can index rasters using, row-column indexing, cell IDs, coordinates, other spatial objects.

# row 1, column 1
elev[1, 1]
##   lyr.1
## 1     1
# cell ID 1
elev[1]
##   lyr.1
## 1     1

If we had a two layered raster, subsetting would return the values in both layers.

two_layers <- c(grain, elev)
two_layers[1]
##   lyr.1 lyr.1
## 1  silt     1

We can also modify/overwrite cell values.

elev[1, 1] = 0
elev[]
##       lyr.1
##  [1,]     0
##  [2,]     2
##  [3,]     3
##  [4,]     4
##  [5,]     5
##  [6,]     6
##  [7,]     7
##  [8,]     8
##  [9,]     9
## [10,]    10
## [11,]    11
## [12,]    12
## [13,]    13
## [14,]    14
## [15,]    15
## [16,]    16
## [17,]    17
## [18,]    18
## [19,]    19
## [20,]    20
## [21,]    21
## [22,]    22
## [23,]    23
## [24,]    24
## [25,]    25
## [26,]    26
## [27,]    27
## [28,]    28
## [29,]    29
## [30,]    30
## [31,]    31
## [32,]    32
## [33,]    33
## [34,]    34
## [35,]    35
## [36,]    36

Replacing values in multi-layer rasters requires a matrix with as many columns as layers and rows as replaceable cells.

two_layers[1] <- cbind(c(1), c(4))
two_layers[]
##       lyr.1 lyr.1
##  [1,]     1     4
##  [2,]     2     2
##  [3,]     2     3
##  [4,]     2     4
##  [5,]     3     5
##  [6,]     1     6
##  [7,]     2     7
##  [8,]     2     8
##  [9,]     1     9
## [10,]     1    10
## [11,]     1    11
## [12,]     1    12
## [13,]     3    13
## [14,]     1    14
## [15,]     2    15
## [16,]     1    16
## [17,]     1    17
## [18,]     1    18
## [19,]     3    19
## [20,]     2    20
## [21,]     2    21
## [22,]     2    22
## [23,]     1    23
## [24,]     1    24
## [25,]     2    25
## [26,]     3    26
## [27,]     3    27
## [28,]     1    28
## [29,]     2    29
## [30,]     3    30
## [31,]     1    31
## [32,]     3    32
## [33,]     3    33
## [34,]     1    34
## [35,]     2    35
## [36,]     1    36

summarizing raster objects

We can get info on raster values just by typing the name or using the summary function.

elev
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     0 
## max value   :    36
summary(elev)
##      lyr.1      
##  Min.   : 0.00  
##  1st Qu.: 9.75  
##  Median :18.50  
##  Mean   :18.47  
##  3rd Qu.:27.25  
##  Max.   :36.00

We can get global summaries, such as standard deviation.

global(elev, sd)
##             sd
## lyr.1 10.58432

Or we can use freq() to get the counts with categories.

freq(grain)
##   layer value count
## 1     1  clay    15
## 2     1  silt    13
## 3     1  sand     8
hist(elev)

spatial subsetting

We can move from subsetting based on specific cell IDs to extract info based on spatial objects.

To use coordinates for subsetting, we can “translate” coordinates into a cell ID with the terra function cellFromXY() or terra::extract().

id <- cellFromXY(elev, xy = matrix(c(0.1, 0.1), ncol = 2))
elev[id]
##   lyr.1
## 1    16
# the same as
terra::extract(elev, matrix(c(0.1, 0.1), ncol = 2))
##   lyr.1
## 1    16

Raster objects can also subset with another raster object. Here we extract the values of our elevation raster that fall within the extent of a masking raster.

clip <- rast(xmin = 0.9, xmax = 1.8, ymin = -0.45, ymax = 0.45,
            resolution = 0.3, vals = rep(1, 9))

elev[clip]
##   lyr.1
## 1    18
## 2    24
# we can also use extract
terra::extract(elev, ext(clip))

In the previous example, we just got the values back. In some cases, we might want the output to be the raster cells themselves. We can do this use the “[” operator and setting “drop = FALSE”

This example returns the first 2 cells of the first row of the “elev” raster.

elev[1:2, drop = FALSE]
## class       : SpatRaster 
## dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, -0.5, 1, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     0 
## max value   :     2

Another common use of spatial subsetting is when we use one raster with the same extent and resolution to mask the another. In this case, the masking raster needs to be composed of logicals or NAs.

# create raster mask of the same resolution and extent
rmask <- elev 
# randomly replace values with NA and TRUE to use as a mask
values(rmask) <- sample(c(NA, TRUE), 36, replace = TRUE) 

# spatial subsetting
elev[rmask, drop = FALSE]   # with [ operator
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     0 
## max value   :    36
mask(elev, rmask)           # with mask()
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     3 
## max value   :    36

We can also use a similar approach to replace values that we suspect are incorrect.

elev[elev < 20] = NA

map algebra

Here we define map algebra as the set of operations that modify or summarize raster cell values with reference to surrounding cells, zones, or statistical functions that apply to every cell.

local operations

Local operations are computed on each cell individually. We can use oridinary arithemetic or logical statements.

elev + elev
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :    40 
## max value   :    72
elev^2
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :   400 
## max value   :  1296
log(elev)
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        :    lyr.1 
## min value   : 2.995732 
## max value   : 3.583519
elev > 5
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :  TRUE 
## max value   :  TRUE

We can also classify intervals of values into groups. For example, we could classify a DEM into low, middle, and high elevation cells

rcl <- matrix(c(0, 12, 1, 12, 24, 2, 24, 36, 3), ncol = 3, byrow = TRUE)
rcl
##      [,1] [,2] [,3]
## [1,]    0   12    1
## [2,]   12   24    2
## [3,]   24   36    3
# we then use this matrix to reclassify our elevation matrix
recl <- classify(elev, rcl = rcl)
recl
## class       : SpatRaster 
## dimensions  : 6, 6, 1  (nrow, ncol, nlyr)
## resolution  : 0.5, 0.5  (x, y)
## extent      : -1.5, 1.5, -1.5, 1.5  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 
## source(s)   : memory
## name        : lyr.1 
## min value   :     2 
## max value   :     3

For more efficient processing, we can use a set of map algebra functions. - app() applies a function to each cell of a raster to summarize the values of multiple layers into one layer - tapp() is an extension of “app()” that allows us to apply on operation on a subset of layers - lapp() allows us to apply a function to each cell using layers as arguments

We can use the lapp() function to compute the Normalized Difference Vegetation Index (NDVI). Let’s calculate NDVI for Zion National Park using multispectral satellite data.

multi_raster_file <- system.file("raster/landsat.tif", package = "spDataLarge")
multi_rast <- rast(multi_raster_file)

We need to define a function to calculate NDVI.

ndvi_fun = function(nir, red){
  (nir - red) / (nir + red)
}

So now we can use lapp() to calculate NDVI in each raster cell. To do so, we just need the NIR and red bands.

ndvi_rast <- lapp(multi_rast[[c(4, 3)]], fun = ndvi_fun)

tm_shape(ndvi_rast) +
  tm_raster()
## SpatRaster object downsampled to 1126 by 889 cells.
## Variable(s) "col" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

focal operations

Local operations operate on one cell, though from multiple layers. Focal operations take into account a central (focal) cell and its neighbors. The neighborhood (or kernel, moving window, filter) can take any size or shape. A focal operation applies an aggregation function to all cells in the neighborhood and updates the value of the central cell before moving on to the next central cell

We can use the focal() function to perform spatial filtering. We define the size, shape, and weights of the moving window using a matrix. Here we find the minimum.

elev <- rast(system.file("raster/elev.tif", package = "spData"))

r_focal <- focal(elev, w = matrix(1, nrow = 3, ncol = 3), fun = min)

plot(elev)

plot(r_focal)

zonal operations

Similar to focal operations, zonal operations apply an aggregation function to multiple cells. However, instead of applying operations to neighbors, zonal operations aggregate based on “zones”. Zones can are defined using a categorical raster and do not necessarily have to be neighbors

For example, we could find the average elevation for different soil grain sizes.

zonal(elev, grain, fun = "mean")
##   lyr.1     elev
## 1  clay 19.26667
## 2  silt 14.76923
## 3  sand 23.12500

merging rasters

In some cases, data for a region will be stored in multiple, contiguous files. To use them as a single raster, we need to merge them.

In this example, we download elevation data for Austria and Switzerland and merge the two rasters into one.

aut = geodata::elevation_30s(country = "AUT", path = tempdir())
ch = geodata::elevation_30s(country = "CHE", path = tempdir())
aut_ch = merge(aut, ch)

geometric operations

When merging or performing map algebra, rasters need to match in their resolution, projection, origin, and/or extent

In the simplest case, two images differ only in their extent. Let’s start by increasing the extent of a elevation raster.

elev = rast(system.file("raster/elev.tif", package = "spData"))
elev_2 = extend(elev, c(1, 2)) # add one row and two columns

plot(elev)

plot(elev_2)

Performing algebraic operations on objects with different extents doesn’t work.

elev + elev_2

We can align the extent of the 2 rasters using the extend() function. Here we extend the elev object to the extent of elev_2 by adding NAs.

elev_4 <- extend(elev, elev_2)

the origin function returns the coordinates of the cell corner closes to the coordinates (0,0). We can also manually change the origin.

origin(elev_4)
## [1] 0 0
origin(elev_4) <- c(0.25, 0.25)
origin(elev_4)
## [1] 0.25 0.25

aggregation and disaggregation

Faster datasets can also differ in their resolution to match resolutions we can decrease the resolution by aggregating or increase the resolution by disaggregating.

Let’s start by changing the resolution of a DEM by a factor of 5, by taking the mean.

dem <- rast(system.file("raster/dem.tif", package = "spDataLarge"))
dem_agg <-  aggregate(dem, fact = 5, fun = mean)

plot(dem)

plot(dem_agg)

We have some choices when increasing the resolution. Here we try the bilinear method.

dem_disagg <- disagg(dem_agg, fact = 5, method = "bilinear")
identical(dem, dem_disagg)
## [1] FALSE
plot(dem_disagg)

resampling

Aggregation/disaggregation work when both rasters have the same origins. What do we do in the case where we have two or more rasters with different origins and resolutions? Resampling computes values for new pixel locations based on custom resolutions and origins.

In most cases, the target raster would be an object you are already working with, but here we define a target raster.

target_rast <- rast(xmin = 794600, xmax = 798200,
                   ymin = 8931800, ymax = 8935400,
                   resolution = 150, crs = "EPSG:32717")

dem_resampl <- resample(dem, y = target_rast, method = "bilinear")

plot(dem)

plot(dem_resampl)