3.3 Manipulating raster objects

In contrast to the vector data model underlying simple features (which represents points, lines and polygons as discrete entities in space), raster data represent continuous surfaces.This section shows how raster objects work by creating them from scratch, building on Section 2.3.1.Because of their unique structure, subsetting and other operations on raster datasets work in a different way, as demonstrated in Section 3.3.1.

The following code recreates the raster dataset used in Section 2.3.3, the result of which is illustrated in Figure 3.2.This demonstrates how the raster() function works to create an example raster named elev (representing elevations).

  1. elev = raster(nrows = 6, ncols = 6, res = 0.5,
  2.               xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
  3.               vals = 1:36)

The result is a raster object with 6 rows and 6 columns (specified by the nrow and ncol arguments), and a minimum and maximum spatial extent in x and y direction (xmn, xmx, ymn, ymax).The vals argument sets the values that each cell contains: numeric data ranging from 1 to 36 in this case.Raster objects can also contain categorical values of class logical or factor variables in R.The following code creates a raster representing grain sizes (Figure 3.2):

  1. grain_order = c("clay", "silt", "sand")
  2. grain_char = sample(grain_order, 36, replace = TRUE)
  3. grain_fact = factor(grain_char, levels = grain_order)
  4. grain = raster(nrows = 6, ncols = 6, res = 0.5, 
  5.                xmn = -1.5, xmx = 1.5, ymn = -1.5, ymx = 1.5,
  6.                vals = grain_fact)

raster objects can contain values of class numeric, integer, logical or factor, but not character.To use character values, they must first be converted into an appropriate class, for example using the function factor().The levels argument was used in the preceding code chunk to create an ordered factor:clay < silt < sand in terms of grain size.See the Data structures chapter of Wickham (2014a) for further details on classes.

raster objects represent categorical variables as integers, so grain[1, 1] returns a number that represents a unique identifier, rather than “clay”, “silt” or “sand”.The raster object stores the corresponding look-up table or “Raster Attribute Table” (RAT) as a data frame in a new slot named attributes, which can be viewed with ratify(grain) (see ?ratify() for more information).Use the function levels() for retrieving and adding new factor levels to the attribute table:

  1. levels(grain)[[1]] = cbind(levels(grain)[[1]], wetness = c("wet", "moist", "dry"))
  2. levels(grain)
  3. #> [[1]]
  4. #>   ID VALUE wetness
  5. #> 1  1  clay     wet
  6. #> 2  2  silt   moist
  7. #> 3  3  sand     dry

This behavior demonstrates that raster cells can only possess one value, an identifier which can be used to look up the attributes in the corresponding attribute table (stored in a slot named attributes).This is illustrated in command below, which returns the grain size and wetness of cell IDs 1, 11 and 35:

  1. factorValues(grain, grain[c(1, 11, 35)])
  2. #>   VALUE wetness
  3. #> 1  sand     dry
  4. #> 2  silt   moist
  5. #> 3  clay     wet

Raster datasets with numeric (left) and categorical values (right).
Figure 3.2: Raster datasets with numeric (left) and categorical values (right).

3.3.1 Raster subsetting

Raster subsetting is done with the base R operator [, which accepts a variety of inputs:

  • Row-column indexing
  • Cell IDs
  • Coordinates
  • Another raster object
    Here, we only show the first two options since these can be considered non-spatial operations.If we need a spatial object to subset another or the output is a spatial object, we refer to this as spatial subsetting.Therefore, the latter two options will be shown in the next chapter (see Section 4.3.1 in the next chapter).

The first two subsetting options are demonstrated in the commands below —both return the value of the top left pixel in the raster object elev (results not shown):

  1. # row 1, column 1
  2. elev[1, 1]
  3. # cell ID 1
  4. elev[1]

To extract all values or complete rows, you can use values() and getValues().For multi-layered raster objects stack or brick, this will return the cell value(s) for each layer.For example, stack(elev, grain)[1] returns a matrix with one row and two columns — one for each layer.For multi-layer raster objects another way to subset is with raster::subset(), which extracts layers from a raster stack or brick. The [[ and $ operators can also be used:

  1. r_stack = stack(elev, grain)
  2. names(r_stack) = c("elev", "grain")
  3. # three ways to extract a layer of a stack
  4. raster::subset(r_stack, "elev")
  5. r_stack[["elev"]]
  6. r_stack$elev

Cell values can be modified by overwriting existing values in conjunction with a subsetting operation.The following code chunk, for example, sets the upper left cell of elev to 0:

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

Leaving the square brackets empty is a shortcut version of values() for retrieving all values of a raster.Multiple cells can also be modified in this way:

  1. elev[1, 1:2] = 0

3.3.2 Summarizing raster objects

raster contains functions for extracting descriptive statistics for entire rasters.Printing a raster object to the console by typing its name returns minimum and maximum values of a raster.summary() provides common descriptive statistics (minimum, maximum, interquartile range and number of NAs).Further summary operations such as the standard deviation (see below) or custom summary statistics can be calculated with cellStats().

  1. cellStats(elev, sd)

If you provide the summary() and cellStats() functions with a raster stack or brick object, they will summarize each layer separately, as can be illustrated by running: summary(brick(elev, grain)).

Raster value statistics can be visualized in a variety of ways.Specific functions such as boxplot(), density(), hist() and pairs() work also with raster objects, as demonstrated in the histogram created with the command below (not shown):

  1. hist(elev)

In case a visualization function does not work with raster objects, one can extract the raster data to be plotted with the help of values() or getValues().

Descriptive raster statistics belong to the so-called global raster operations.These and other typical raster processing operations are part of the map algebra scheme, which are covered in the next chapter (Section 4.3.2).

Some function names clash between packages (e.g., select(), as discussed in a previous note). In addition to not loading packages by referring to functions verbosely (e.g., dplyr::select()), another way to prevent function names clashes is by unloading the offending package with detach(). The following command, for example, unloads the raster package (this can also be done in the package tab which resides by default in the right-bottom pane in RStudio): detach(“package:raster”, unload = TRUE, force = TRUE). The force argument makes sure that the package will be detached even if other packages depend on it. This, however, may lead to a restricted usability of packages depending on the detached package, and is therefore not recommended.