Introduction

For this in-class exercise, we’ll practice working with raster data. Specifically, we’ll learn:

  • How to crop rasters using a polygon
  • How to make a raster mask from a polygon
  • How to rasterize a polygon using a template raster
  • How to use levelplot() to make a [relatively] nice map of a raster of factor data

Data

You will need the following datasets, which you can download from the course GitHub site:

  • temperature.zip A raster of average annual temperatures for the Contiguous United States (CONUS)
  • wy_countes.zip A vector data set of Wyoming county borders.

Download and unzip them to subdirectories of your data folder.

Load the data into an R session

I’ll provide the code to read the temperature data, you can figure out the WY counties code:

temperature = raster(here("data", "temperature", "temperature.tif"))

Make sure your temperature raster and county borders are in the same projection.

Given what we know about transforming rasters, you’ll need to select which of the following is best:

  1. Reproject the temperature raster into the WY counties projection.
  2. Reproject the WY data into the coordinate system of the temperature raster

Once you’ve got your datasets loaded, you can create a quick plot as a reality check to verify that the projections match:

Cropping a raster

Rasters are often huge. and therefore slow to work with and plot. Sometimes only need a small part of a larger raster, which can speed up our processing time.

It’s usually a good idea to crop a larger raster as one of the first steps in your analysis.

The crop() function in the raster package is your friend!

Even though it’s in the raster package, it also knows how to work with vector spatial types, like SpatialPolygonsDataFrame objects.

Try using crop() on the temperature data. You should end up with something that looks like this:

If you zoom in closely, you might notice that there are some small areas just inside the WY state border that are white, i.e. they don’t have any raster cells.

It’s a good practice to do a small buffer around your polygon before you crop.

Remember our friends raster::buffer() and rgeos::gBuffer()?

Try adding a 10km buffer around the edges of WY, then cropping:

## Loading required namespace: rgeos

Raster Masks

It might be hard to see in the cropping example above (since Wyoming is pretty rectangular), but the crop() function doesn’t crop to the borders of the, but rather to the extent or bounding box.

It’s easier to see if we crop to one of the more irregular counties:

carbon_cnty = subset(wy_cnty, COUNTYNAME == "Carbon")
carbon_cnty_rst = crop(tmp_wy_buf, buffer(carbon_cnty, width = 10000))
plot(carbon_cnty_rst)
plot(carbon_cnty, add = T, lwd = 2)

res(carbon_cnty_rst)
## [1] 800 800
carbon_yst = projectRaster(
  carbon_cnty_rst, 
  res =800,
  crs = projection(ystone))

carbon_cnty_rst
## class      : RasterLayer 
## dimensions : 293, 284, 83212  (nrow, ncol, ncell)
## resolution : 800, 800  (x, y)
## extent     : -12024457, -11797257, 5001861, 5236261  (xmin, xmax, ymin, ymax)
## crs        : +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs 
## source     : memory
## names      : temperature 
## values     : -1.53, 7.94  (min, max)
carbon_yst
## class      : RasterLayer 
## dimensions : 236, 230, 54280  (nrow, ncol, ncell)
## resolution : 800, 800  (x, y)
## extent     : 718141.8, 902141.8, -372903.7, -184103.7  (xmin, xmax, ymin, ymax)
## crs        : +proj=lcc +lat_0=44.25 +lon_0=-109.5 +lat_1=49 +lat_2=45 +x_0=600000 +y_0=0 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : temperature 
## values     : -1.371816, 7.944703  (min, max)
plot(carbon_cnty_rst)

plot(carbon_yst)

proj4string(carbon_yst)
## [1] "+proj=lcc +lat_0=44.25 +lon_0=-109.5 +lat_1=49 +lat_2=45 +x_0=600000 +y_0=0 +datum=NAD83 +units=m +no_defs"
proj4string(wy_cnty)
## Warning in proj4string(wy_cnty): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
projection(wy_cnty)
## [1] "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"

Some things to note:

  • I first buffered the county by 10km.
  • I used the border of Carbon county, the 6th county in the dataset.
  • My variable tmp_wy_buf is the temperature raster cropped to a 10km buffer of the Wyoming border.

If we want to keep only the cells within the border, we can create a raster mask.

It’s easy with the mask function in the raster package:

plot(mask(carbon_cnty_rst, carbon_cnty))
plot(carbon_cnty, add = T)

Practice cropping the temperature raster to a different Wyoming county, and creating a mask.

Rasterize a polygon

Sometimes we need to turn a polygons into a raster layer.

It’s fairly easy to do with the rasterize() function.

Try using this funciton to rasterize the Wyoming counties. Hint: you’ll need a template raster.

When you’re successful, you should see something like this:

Note that this operation may be a little slow.