Earth analytics
earth-analytics Home

Lesson 6. Clip Raster in R

Leah Wasser

Learning Objectives

After completing this tutorial, you will be able to:

  • Crop a raster dataset in R using a vector extent object derived from a shapefile.
  • Open a shapefile in R.

What You Need

You need R and RStudio to complete this tutorial. Also you should have an earth-analytics directory set up on your computer with a /data directory with it.

R Libraries to Install:

  • raster: install.packages("raster")
  • rgdal: install.packages("rgdal")
    • sf: install.packages("sf")

If you have not already downloaded the week 3 data, please do so now. Download Week 3 Data (~250 MB)

In this lesson, you will learn how to crop a raster dataset in R. Previously, you reclassified a raster in R, however the edges of your raster dataset were uneven. In this lesson, you will learn how to crop a raster - to create a new raster object / file that you can share with colleagues and / or open in other tools such as QGIS.

Load Libraries

# load the raster and rgdal libraries
library(raster)
library(rgdal)
# if you want to use sf. you will use sf for future lessons!

Open Raster and Vector Layers

First, you will use the raster() function to open a raster layer. Let’s open the canopy height model that you created in the previous lesson

# open raster layer
lidar_chm <- raster("data/week-03/BLDR_LeeHill/outputs/lidar_chm.tif")

# plot CHM
plot(lidar_chm,
     col = rev(terrain.colors(50)))

lidar chm plot

Open Vector Layer

Next, let’s open up a vector layer that contains the crop extent that you want to use to crop your data. To open a shapefile you use the readOGR() function.

# import the vector boundary
crop_extent <- readOGR("data/week-03/BLDR_LeeHill/clip-extent.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/root/earth-analytics/data/week-03/BLDR_LeeHill/clip-extent.shp", layer: "clip-extent"
## with 1 features
## It has 1 fields
## Integer64 fields read as strings:  id

# plot imported shapefile
# notice that you use add = T to add a layer on top of an existing plot in R.
plot(crop_extent,
     main = "Shapefile imported into R - crop extent",
     axes = TRUE,
     border = "blue")

shapefile crop extent plot

The spatial extent of a shapefile or R spatial object represents the geographic 'edge' or location that is the furthest north, south east and west.
The spatial extent of a shapefile or `R` spatial object represents the geographic "edge" or location that is the furthest north, south east and west. Thus is represents the overall geographic coverage of the spatial object. Image Source: Colin Williams, NEON.

Now that you have imported the shapefile. You can use the crop() function in R to crop the raster data using the vector shapefile.

# crop the lidar raster using the vector extent
lidar_chm_crop <- crop(lidar_chm, crop_extent)
plot(lidar_chm_crop, main = "Cropped lidar chm")

# add shapefile on top of the existing raster
plot(crop_extent, add = TRUE)

lidar chm cropped with vector extent on top

Optional Challenge: Crop Change Over Time Layers

In the previous lesson, you created 2 plots:

  1. A classified raster map that shows positive and negative change in the canopy height model before and after the flood. To do this you will need to calculate the difference between two canopy height models.
  2. A classified raster map that shows positive and negative change in terrain extracted from the pre and post flood Digital Terrain Models before and after the flood.

Create the same two plots except this time CROP each of the rasters that you plotted using the shapefile: data/week-03/BLDR_LeeHill/crop_extent.shp

For each plot, be sure to:

  • Add a legend that clearly shows what each color in your classified raster represents.
  • Use proper colors.
  • Add a title to your plot.

You will include these plots in your final report due next week.

Check out this cheatsheet for more on colors in R.

Or type: grDevices::colors() into the r console for a nice list of colors!

Classify a Raster

Leave a Comment

You May Also Enjoy

Plot Data With Matplotlib

1 minute read

This Week

Calculate Seasonal Summary Values from Climate Data Variables Stored in NetCDF 4 Format: Work With MACA v2 Climate Data in Python

26 minute read

Learn how to calculate seasonal summary values for MACA 2 climate data using xarray and region mask in open source Python.

Calculate Summary Values Using Spatial Areas of Interest (AOIs) including Shapefiles for Climate Data Variables Stored in NetCDF 4 Format: Work With MACA v2 Climate Data in Python

32 minute read

Climate datasets stored in netcdf 4 format often cover the entire globe or an entire country. Learn how to subset climate data spatially and by time slices u...

How to Open and Process NetCDF 4 Data Format in Open Source Python

26 minute read

Historic and projected climate data are most often stored in netcdf 4 format. Learn how to open and process MACA version 2 climate data for the Continental U...