09-2: R as GIS: Raster Data Basics

Before you start


Learning objectives

Learn how to handle raster datasets using R.


Table of contents

Tips to make the most of the lecture notes


Interactive navigation tools

  • Click on the three horizontally stacked lines at the bottom left corner of the slide, then you will see table of contents, and you can jump to the section you want

  • Hit letter “o” on your keyboard and you will have a panel view of all the slides


Running and writing codes

  • The box area with a hint of blue as the background color is where you can write code (hereafter referred to as the “code area”).
  • Hit the “Run Code” button to execute all the code inside the code area.
  • You can evaluate (run) code selectively by highlighting the parts you want to run and hitting Command + Enter for Mac (Ctrl + Enter for Windows).
  • If you want to run the codes on your computer, you can first click on the icon with two sheets of paper stacked on top of each other (top right corner of the code chunk), which copies the code in the code area. You can then paste it onto your computer.
  • You can click on the reload button (top right corner of the code chunk, left to the copy button) to revert back to the original code.

The terra and raster packages: Basics

  • The most popular R packages to handle raster data is the raster and terra packages. They are authored by the same person and terra is the successor of raster.

  • The terra package is now mature and does everything faster than raster

  • There are lots of packages that still depends on raster and do not work well with R object class defined by terra.

  • We primarily learns how to use the terra package

The terra and raster packages use different R object class to represent raster data:

  • terra
    • SpatRaster
  • raster (collectively referred to as Raster*)
    • RasterLayer
    • RasterStack
    • RasterBrick

Note

  • We will learn how to covert SpatRaster to Raster* and vice versa.
  • This is important as there are (and will be) packages that accepts only Raster*, especially the ones that are old.

Let’s take a look at a SpatRaster object, reflec_blue.

Check the class first:


What is inside?


Here are the explanations of the information provided:

  • class
  • dimensions
    • nrow: number of rows
    • ncol: number of columns
    • nlyr: number of layers
  • resolution:
    • x: how long the top and bottom lines of the grid are
    • y: how long the left and right lines of the grid are
  • extent: bounding box (just like what you get with st_bbox(sf))
  • coord. ref.: CRS
  • name: name of the attribute

You can use plot() to make a map quickly:

It automatically color the grids by the value of the attribute (blue).

You can simply use c() function to create a multi-layer SpatRaster just like you create a vector as long as all the layers have exactly the same dimensions, extent, and resolution.

Notice that nlyr is 3 now and you see three attribute names in names.


plot() create maps for all the attributes.

The raster package differentiates single-layer and multi-layer raster data.

  • RasterLayer: single-layer
  • RasterStack: multi-layer
  • RasterBrick: multi-layer

You can convert a SpatRaster to a Raster* object using as(SpatRaster, "Raster").

Since, reflec_blue is a single-layer SpatRaster, it was converted into a RasterLayer.


Since, reflec_all is a multi-layer SpatRaster, it was converted into a RasterBrick.


Note

You can convert a RasterBrick to a RasterStack by applying stack() to the RasterBrick if you want. But, you do not need to.

You can convert an Raster* object to SpatRaster using terra::rast() function.

Basic operations

This is very useful. As we will see later, when interacting two spatial objects (e.g., extracting values from a raster data to sf) some functions require that the two spatial objects has the same CRS. You can use terra::crs() to get the CRS of the raster data and apply it to another spatial object.

This is like st_bbox() for sf.

You can access the cell values using [] just like a vector. Note that head() is there to avoid a very long vector presented on the console.


Yes, there are so many NAs in this raster data. Let’s look at the value of 10100th through 10200th cells:


Note

  • We use this function to check our arithmetic operations are successful.
  • I have never had to use this function before in practice.

You can use terra::xyFromCell() to get the geographic coordinates of the cell centers.


Syntax

terra::xyFromCell(SpatRaster, cell)


  • cell: cell numbers


Example

You can access specific layers using subset().


Syntax

subset(SpatRaster, subset)
  • subset: layer names or corresponding integers


Examples


Raster data input and output

Raster data files can come in numerous different formats.

  • The most common format is GeoTiff (.tif as the file extension)
  • PRPISM weather data comes in the Band Interleaved by Line (BIL) format
  • Some of the Daymet data comes in netCDF format.
  • Other popular formats include SAGA, ENVI, and many others.

You can read data of almost all the existing file formats with the terra package.

You can use terra::rast() to read raster data files.


single raster data file

drone_blue_sr <- terra::rast(path to the file)


multiple raster data files

drone_blue_sr <- terra::rast(c(path to the file 1, path to the file 1, ...))

Instruction

  • download reflec_blue.tif, reflec_red.tif, and reflec_green.tif from this link.
  • find the path to the three files


single file

It looks like this for me:

reflec_blue <- terra::rast("data-files-for-participants/reflec_blue.tif")


multiple files

It looks like this for me:

reflec_all <- 
  terra::rast(
    c(
      "data-files-for-participants/reflec_blue.tif",
      "data-files-for-participants/reflec_red.tif",
      "data-files-for-participants/reflec_green.tif"
    )
  )

Instruction

  • download the folder named PRISM_ppt_stable_4kmD2_20120801_bil from this link.
  • find the path to the file named PRISM_ppt_stable_4kmD2_20120801_bil.bil, which is PRISM precipitation data observed on 08-01-2012.


It has a different file extension of .bil. Well it does no matter. Just use terra::rast() with path to the file inside it just like you did with the GeoTiff files.


This is what the code looks like for me:

prism_precip <- terra::rast("data-files-for-participants/PRISM_ppt_stable_4kmD2_20120801_bil/PRISM_ppt_stable_4kmD2_20120801_bil.bil")

Many weather data comes in the netCDF format with the file extension of .nc. (e.g., gridMET)

Instruction

  • download the file named gm_precip_2018.nc from this link.
  • find the path to the file


Let’s read the data now. This is what the code looks like for me:

prism_precip <- terra::rast("data-files-for-participants/gm_precip_2018.nc")

You can use terra::writeRaster() to write raster data to a data file.


Syntax

terra::writeRaster(SpatRaster, path to the file)


Example (does not run)

writeRaster(reflec_blue, "./data/reflec_blue.tif", overwrite = TRUE)

  • This code saves reflec_blue (a SpatRaster object) as a GeoTiff file.

  • writeRaster() infers the correct format from the extension of the file name, which if .tif here.

  • The overwrite = TRUE option is necessary if a file with the same name already exists and you are overwriting it.


Note

No distinction is necessary for single-layer and multi-layer SpatRaster objects.

single-layer

Write reflec_blue on your computer. Mine looks like this:

terra::writeRaster(reflec_blue, "data/reflec_blue.tif")


Then read it back:

reflec_blue_re_read <- terra::rast("data/reflec_blue.tif")


Confirm that reflec_blue_re_read is the same as reflec_blue.


multi-layer

Write reflec_all on your computer. Mine looks like this:

terra::writeRaster(reflec_all, "data/reflec_all.tif")


Then read it back:

reflec_all_re_read <- terra::rast("data/reflec_all.tif")


Confirm that reflec_all_re_read is the same as reflec_all.

While you can use terra::writeRaster() to write to a netCDF file, you will get a note saying that you should consider using terra::writeCDF().

Write reflec_all to a netCDF file. Mine looks like this.

terra::writeCDF(reflec_all, "data/reflec_all.nc")

Raster Data Operations

  • You can do basic arithmetic operations (addition, subtraction, division, etc) using raster layers as long as they share the same spatial extent and resolution

  • You can also apply a function like log() to transform the value of the cells

  • Raster arithmetic operations are done element-by-element (cell-by-cell) just like vector arithmetic operations.

For example, when two RasterLayers are added, then the two values associated with the same cell are added and the resulting value becomes the new value for the cell in the newly created RasterLayer.


Did it work?


Yes, looks like it did. Look at different cells yourself. Note that the name of the attribute in reflec_b_plus_g inherited the name of the attribute in reflec_blue (first SpatRaster in the addition above).

Multiply reflec_blue with reflec_red and add square root of reflec_green:


Look at several cells to confirm that multiplication was successful.


Answer:

Code 
reflect_temp <- reflec_blue * reflec_red + sqrt(reflec_green)

#--- confirm if the calculations were done right ---#
reflect_temp[10100] == reflec_blue[10100] * reflec_red[10100] + sqrt(reflec_green[10100])

Using NIR and RED, calculate NDVI and name it NDVI.

\(NDVI = (NIR-RED)/(NIR+RED)\)


Then, plot it to see what it looks like.


Answer:

Code 
NDVI <- (NIR - RED) / (NIR + RED)

#--- confirm if the calculations were done right ---#
plot(NDVI)

Sometimes, you want to make your raster data have a lower resolution. For example, satellite image is often very fine with spatial resolution of say 30cm. When trying to create a map using the data, it takes long time for R to render a plot and also the size of the figure can be very large.


Syntax

aggregate(SpatRaster, fact)
  • fact:

Example


Let’s compare before and after. After aggregating by factor of 5, the map is visibly more coarse. Maybe this was too much.

Aggregate NIR by factor of 2.


Then, check if its plot is acceptable or not.


Answer:

Code 
NIR_ag_2 <- aggregate(NIR, fact = 2)

plot(NIR_ag_2)

  • You have two raster layers that differ in any of the dimension and resolution.

  • You want to assign a value from one layer to each of the cells in the other layer so that you have consistent observation units for the variables from the two layers.

Fake precipitation data:

Fake soil data:

  • Map below shows grids from precip (red border) and from soil (blue border).

  • As you can see, the grids from the two layers are not nicely aligned.

  • See the the top left grid with green fill color from the precip layer. Resampling will find a single value to the grid based on the value of the nearby grids from the soil layer.

  • Even though the name “resample” sounds like it is a random process, there is no randomness.

Syntax

resample(sr_1, sr_2, method)


  • st_1: SpatRaster to be resampled
  • st_2: SpatRaster with the geometry that st_1 should be resampled to
  • method: method of assigning values
    • “near”: nearest neighbor (Default)
    • “cubicspline”: cubic-spline interpolation
    • “bilinear”: bi-linear interpolation
    • others (run ?terra::resample to see all the method options)

Since we are resampling values from the soil layer to the grids of precip,


Try yourself

  • Change the method and see how the resampling results change.
  • Many of the spatial datasets tend to exhibit a positive spatial correlation and the resampling outcomes are not as sensitive to method as you see here where cell values are completely independent.

Sometimes, you have two or more raster layers that have different spatial coverages. In such a case, you might want to merge them into a single raster layer.

For demonstration purpose, we will use two SpatRaster objects: prism_saunders and prism_douglas. They are PRISM maximum temperature observed on 08/01/2012 in the Saunders and Douglas counties in Nebraska, which are adjacent to each other.

Note

Note that this is different from combining multiple single-layer raster data of the same spatial extent and resolution into multi-layer raster data.

You can use the terra::merge() function to merge two raster datasets into one.


You can check the result of the merging below:


Note

Remember that raster object has to be perfectly rectangular. The result of merging will construct a rectangle that encompasses both prism_saunders and prism_douglas. All the cells that are not covered by prism_saunders and prism_douglas will be assigned NA.

When merging two SpatRaster objects and when they have spatial overlaps, the value of the first SpatRaster object will be respected.

Let’s run a little experiment. We will assign a high value to all the cells in prism_saunders. This will make this phenomenon easy to detect.


prism_saunders first


prism_douglas first