This activity will walk you through the fundamental principles of working with raster data in R.
R Skill Level: Intermediate - you’ve got the basics of
Goals / ObjectivesAfter completing this activity, you will know:
- What a raster dataset is and its fundamental attributes.
- How to import rasters into `R` using the raster library.
- How to perform raster calculations in `R`.
Things You'll Need To Complete This Lesson
R Libraries to Install:
Tools To InstallPlease be sure you have the most current version of `R` and preferably R studio to write your code.
Data to DownloadNational Ecological Observatory Network's San Joaquin Experimental Range and Soaproot Saddle field sites. This data is intended for educational purposes, for access to all the data for research purposes visit the NEON Data Portal.
Recommended Pre-Lesson Reading
About Raster Data
Raster or “gridded” data are data that are saved in pixels. In the spatial world, each pixel represents an area on the Earth’s surface. For example in the raster below, each pixel represents a particular land cover class that would be found in that location in the real world. More on rasters here.
To work with rasters in R, we need two key libraries,
To install the raster library you can use
When you install the raster library,
sp should also install. Also install the
install.packages('rgdal'). Among other things,
allow us to export rasters to geotiff format.
#load the raster, sp, and rgdal packages library(raster) library(sp) library(rgdal) #Set your working directory to the folder where your data for this workshop #are stored. NOTE: if you created a project file in R studio, then you don't #need to set the working directory as it's part of the project. #setwd("~/yourWorkingDirectoryHere")
Next, let’s load a raster containing elevation data into R.
#load raster in an R object called 'DEM' DEM <- raster("DigitalTerrainModel/SJER2013_DTM.tif") # look at the raster attributes. DEM ## class : RasterLayer ## dimensions : 5060, 4299, 21752940 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 254570, 258869, 4107302, 4112362 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## data source : C:\Users\lwasser\Documents\1_Workshops\ESA2015\DigitalTerrainModel\SJER2013_DTM.tif ## names : SJER2013_DTM
Notice a few things about this raster.
- Nrow, Ncol: the number of rows and columns in the data (imagine a spreadsheet or a matrix).
- Ncells: the total number of pixels or cells that make up the raster.
- Resolution: the size of each pixel (in meters in this case). 1 meter pixels means that each pixel represents a 1m x 1m area on the earth’s surface.
- Extent: the spatial extent of the raster. This value will be in the same coordinate units as the coordinate reference system of the raster.
- Coord ref: the coordinate reference system string for the raster. This raster is in UTM (Universal Trans mercator) zone 11 with a datum of WGS 84. More on UTM here.
Defining Min/Max Values
By default the raster doesn’t have the min / max values associated with it’s attributes Let’s change that.
#calculate and save the min and max values of the raster to the raster object DEM <- setMinMax(DEM) #view raster attributes DEM ## class : RasterLayer ## dimensions : 5060, 4299, 21752940 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 254570, 258869, 4107302, 4112362 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## data source : C:\Users\lwasser\Documents\1_Workshops\ESA2015\DigitalTerrainModel\SJER2013_DTM.tif ## names : SJER2013_DTM ## values : 228.1, 518.66 (min, max)
You can also view the rasters min and max values and the range of values containes within the pixels.
#Get min and max cell values from raster #NOTE: this code may fail if the raster is too large cellStats(DEM, min) ##  228.1 cellStats(DEM, max) ##  518.66 cellStats(DEM, range) ##  228.10 518.66
Working with Rasters in R
Now that we have the raster loaded into R, let’s grab some key raster attributes.
#view coordinate reference system DEM@crs ## CRS arguments: ## +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 ## +towgs84=0,0,0 #view raster extent DEM@extent ## class : Extent ## xmin : 254570 ## xmax : 258869 ## ymin : 4107302 ## ymax : 4112362 #plot the raster #note that this raster represents a small region of the NEON SJER field site plot(DEM, main="Digital Elevation Model, NEON Field Site")
We can also create a histogram to view the distribution of values in our raster.
Note that the max number of pixels that R will plot by default is 100,000. We
can tell it to plot more using the
#we can look at the distribution of values in the raster too hist(DEM, main="Distribution of elevation values", col= "purple", maxpixels=21752940)
Basic Raster Math
We can also perform calculations on our raster. For instance, we could multiply all values within the raster by 2.
#multiple each pixel in the raster by 2 DEM2 <- DEM * 2 DEM2 ## class : RasterLayer ## dimensions : 5060, 4299, 21752940 (nrow, ncol, ncell) ## resolution : 1, 1 (x, y) ## extent : 254570, 258869, 4107302, 4112362 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## data source : in memory ## names : SJER2013_DTM ## values : 456.2, 1037.32 (min, max) #plot the new DEM plot(DEM2, main="DEM with all values doubled")
Plotting Raster Data
R has an
image function that allows you to control the way a raster is
rendered on the screen. The
plot command in R has a base setting for the number
of pixels that it will plot (100,000 pixels). The image command thus might be
better for rendering larger rasters.
#create a plot of our raster image(DEM)
#specify the range of values that you want to plot in the DEM #just plot pixels between 250 and 300 m in elevation image(DEM, zlim=c(250,300))
#we can specify the colors too col <- terrain.colors(5) image(DEM, zlim=c(250,300), main="Digital Elevation Model (DEM)", col=col)
In the above example.
R to create a palette of colors within the
terrain.colors color ramp. There are other palettes that you can use as well include
What happens if you change the number of colors in the
What happens if you change the
zlim values in the image command?
What are the other attributes that you can specify when using the
Breaks and Colorbars in R
A digital elevation model (DEM) is an example of a continuous raster. It contains elevation values for a range. For example, elevations values in a DEM might include any set of values between 200 m and 500 m. Given this range, you can plot DEM pixels using a gradient of colors.
By default, R will assign the gradient of colors uniformly across the full range of values in the data. In our case, our DEM has values between 250 and 500. However, we can adjust the “breaks” which represent the numeric locations where the colors change if we want too.
#add a color map with 5 colors col=terrain.colors(5) #add breaks to the colormap (6 breaks = 5 segments) brk <- c(250, 300, 350, 400,450,500) plot(DEM, col=col, breaks=brk, main="DEM with more breaks")
# Expand right side of clipping rect to make room for the legend par(xpd = FALSE,mar=c(5.1, 4.1, 4.1, 4.5)) #DEM with a custom legend plot(DEM, col=col, breaks=brk, main="DEM with a Custom (buf flipped) Legend",legend = FALSE) #turn xpd back on to force the legend to fit next to the plot. par(xpd = TRUE) #add a legend - but make it appear outside of the plot legend( par()$usr, 4110600, legend = c("lowest", "a bit higher", "middle ground", "higher yet", "Highest"), fill = col)
Notice that the legend is in reverse order in the previous plot. Let’s fix that.
# Expand right side of clipping rect to make room for the legend par(xpd = FALSE,mar=c(5.1, 4.1, 4.1, 4.5)) #DEM with a custom legend plot(DEM, col=col, breaks=brk, main="DEM with a Custom Fixed Legend",legend = FALSE) #turn xpd back on to force the legend to fit next to the plot. par(xpd = TRUE) #add a legend - but make it appear outside of the plot legend( par()$usr, 4110600, legend = c("Highest", "Higher yet", "Middle","A bit higher", "Lowest"), fill = rev(col))
We can add a custom color map with fewer breaks as well.
#add a color map with 4 colors col=terrain.colors(4) #add breaks to the colormap (6 breaks = 5 segments) brk <- c(200, 300, 350, 400,500) plot(DEM, col=col, breaks=brk, main="DEM with fewer breaks")
A discrete dataset has a set of unique categories or classes. One example could be landuse classes. The example below shows elevation zones generated using the same DEM.
Cropping Rasters in R
You can crop rasters in R using different methods. You can crop the raster directly drawing a box in the plot area. To do this, first plot the raster. Then define the crop extent by clicking twice:
- Click in the UPPER LEFT hand corner where you want the crop box to begin.
- Click again in the LOWER RIGHT hand corner to define where the box ends.
You’ll see a red box on the plot. NOTE that this is a manual process that can be used to quickly define a crop extent.
#plot the DEM plot(DEM) #Define the extent of the crop by clicking on the plot cropbox1 <- drawExtent() #crop the raster, then plot the new cropped raster DEMcrop1 <- crop(DEM, cropbox1) #plot the cropped extent plot(DEMcrop1) You can also manually assign the extent coordinates to be used to crop a raster. We'll need the extent defined as (`xmin`, `xmax`, `ymin` , `ymax`) to do this. This is how we'd crop using a GIS shapefile (with a rectangular shape) #define the crop extent cropbox2 <-c(255077.3,257158.6,4109614,4110934) #crop the raster DEMcrop2 <- crop(DEM, cropbox2) #plot cropped DEM plot(DEMcrop2)
Open up file:
it to an R object called
DSM. Convert the raster data to feet.
DSM in feet that you created above. Use a custom
color map. Create numeric breaks that make sense given the distribution of the data.
Hint: your breaks might represent