This tutorial will walk you through the fundamental principles of working with image raster data in R.
R Skill Level: Intermediate
After completing this activity, you will be able to:
- Import multiple image rasters and create a stack of rasters.
- Plot three band RGB images in R.
- Export single band and multiple band image rasters in R.
Things You’ll Need To Complete This Tutorial
You will need the most current version of
R and, preferably,
on your computer to complete this tutorial.
Install R Packages
These remote sensing data files provide information on the vegetation at the National 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.
This data download contains several files. You will only need the RGB .tif files for this tutorial. The path to this file is: NEON-DS-Field-Site-Spatial-Data/SJER/RGB/* . The other data files in the downloaded data directory are used for related tutorials. You should set your working directory to the SJER directory to follow the code exactly.
You may benefit from reviewing these related tutorials prior to this tutorial:
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. An color image raster is a bit different from other rasters in that it has multiple bands. Each band represents reflectance values for a particular color or spectra of light. If the image is RGB, then the bands are in the red, green and blue portions of the electromagnetic spectrum. These colors together create what we know as a full color image.
Work with Multiple Rasters
In a previous tutorial, we loaded a
single raster into R. We made sure we knew the
(coordinate reference system) and extent of the dataset among other key metadata
attributes. This raster was a Digital Elevation Model so there was only a single
raster that represented the ground elevation in each pixel. When we work with
color images, there are multiple rasters to represent each band. Here we’ll learn
to work with multiple rasters together.
A raster stack is a collection of raster layers. Each raster layer in the raster stack needs to have the same
- projection (CRS),
- spatial extent and
You might use raster stacks for different reasons. For instance, you might want to group a time series of rasters representing precipitation or temperature into one R object. Or, you might want to create a color images from red, green and blue band derived rasters.
In this tutorial, we will stack three bands from a multi-band image together to create a composite RGB image.
First let’s load the R packages that we need:
# load the raster, sp, and rgdal packages library(raster) library(sp) library(rgdal) # set the working directory to the data #setwd("pathToDirHere")
Next, let’s create a raster stack with bands representing
- blue: band 19, 473.8nm
- green: band 34, 548.9nm
- red; band 58, 669.1nm
This can be done by individually assigning each file path as an object.
# import tiffs band19 <- "RGB/band19.tif" band34 <- "RGB/band34.tif" band58 <- "RGB/band58.tif"
Note that if we wanted to create a stack from all the files in a directory (folder)
you can easily do this with the
list.files() function. We would use
full.names=TRUE to ensure that R will store the directory path in our list of
# create list of files to make raster stack rasterlist1 <- list.files('RGB', full.names=TRUE)
Or, if your directory is consistes of some .tif files and other file types you don’t want in your stack, you can ask R to only list those files with a .tif extension.
rasterlist2 <- list.files('RGB', full.names=TRUE, pattern="tif")
Back to creating our raster stack with three bands. We only want three of the
bands in the RGB directory and not the fourth
band90, so will create the stack
from the bands we loaded individually. We do this with the
# create raster stack rgbRaster <- stack(band19,band34,band58) # example syntax for stack from a list #rstack1 <- stack(rasterlist1)
This has now created a stack of rasters three thick. Let’s view them.
# check attributes rgbRaster ## class : RasterStack ## dimensions : 502, 477, 239454, 3 (nrow, ncol, ncell, nlayers) ## resolution : 1, 1 (x, y) ## extent : 256521, 256998, 4112069, 4112571 (xmin, xmax, ymin, ymax) ## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 ## names : band19, band34, band58 ## min values : 84, 116, 123 ## max values : 13805, 15677, 14343 # plot stack plot(rgbRaster)
From the attributes we see the CRS, resolution, and extent of all three rasters. The we can see each raster plotted. Notice the different shading between the different bands. This is because the landscape relects in the red, green, and blue spectra differently.
Check out the scale bars. What do they represent?
This reflectance data is radiance corrected for atmospheric effects. The data are typically unitless and ranges from 0-1. NEON Airborne Observation Platform data, where these rasters come from, has a scale factor of 10,000.
Plot an RGB Image
You can plot a composite RGB image from a raster stack. You need to specify the order of the bands when you do this. In our raster stack, band 19, which is the blue band, is first in the stack, whereas band 58, which is the red band, is last. Thus the order for a RGB image is 3,2,1 to ensure the red band is rendered first as red.
Thinking ahead to next time: If you know you want to create composite RGB images, think about the order of your rasters when you make the stack so the RGB=1,2,3.
We will plot the raster with the
rgbRaster() function and the need these
- R object to plot
- which layer of the stack is which color
- stretch: allows for increased contrast. Options are “lin” & “hist”.
Let’s try it.
# plot an RGB version of the stack plotRGB(rgbRaster,r=3,g=2,b=1, stretch = "lin")
Note: read the
raster package documentation for other arguments that can be
scale) to improve or modify the image.
Explore Raster Values - Histograms
You can also explore the data. Histograms allow us to view the distrubiton of values in the pixels.
# view histogram of reflectance values for all rasters hist(rgbRaster) ## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = ## main[y[i]], : 42% of the raster cells were used. 100000 values used. ## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = ## main[y[i]], : 42% of the raster cells were used. 100000 values used. ## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = ## main[y[i]], : 42% of the raster cells were used. 100000 values used.
Note about the warning messages: R defaults to only showing the first 100,000 values in the histogram so if you have a large raster you may not be seeing all the values. This saves your from long waits, or crashing R, if you have large datasets.
You can crop all rasters within a raster stack the same way you’d do it with a single raster.
# determine the desired extent rgbCrop <- c(256770.7,256959,4112140,4112284) # crop to desired extent rgbRaster_crop <- crop(rgbRaster, rgbCrop) # view cropped stack plot(rgbRaster_crop)
Challenge: Plot Cropped RGB
Plot this new cropped stack as an RGB image.
Raster Bricks in R
rgbRaster object we have a list of rasters in a stack. These rasters
are all the same extent, CRS and resolution. By creating a raster brick we
will create one raster object that contains all of the rasters so that we can
use this object to quickly create RGB images. Raster bricks are more efficient
objects to use when processing larger datasets. This is because the computer
doesn’t have to spend energy finding the data - it is contained within the object.
# create raster brick rgbBrick <- brick(rgbRaster) # check attributes rgbBrick ## class : RasterBrick ## dimensions : 502, 477, 239454, 3 (nrow, ncol, ncell, nlayers) ## resolution : 1, 1 (x, y) ## extent : 256521, 256998, 4112069, 4112571 (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 : band19, band34, band58 ## min values : 84, 116, 123 ## max values : 13805, 15677, 14343
While the brick might seem similar to the stack (see attributes above), we can see that it’s very different when we look at the size of the object.
- the brick contains all of the data stored in one object
- the stack contains links or references to the files stored on your computer
object.size() to see the size of an R object.
# view object size object.size(rgbBrick) ## 5759448 bytes object.size(rgbRaster) ## 40408 bytes # view raster brick plotRGB(rgbBrick,r=3,g=2,b=1, stretch = "Lin")
Notice the faster plotting? For a smaller raster like this the difference is slight, but for larger raster it can be considerable.
Write to GeoTIFF
We can write out the raster in GeoTIFF format as well. When we do this it will copy the CRS, extent and resolution information so the data will read properly into a GIS program as well. Note that this writes the raster in the order they are in. In our case, the blue (band 19) is first but most programs expect the red band first (RGB).
One way around this is to generate a new raster stack with the rasters in the proper order - red, green and blue format. Or, just always create your stacks R->G->B to start!!!
# Make a new stack in the order we want the data in orderRGBstack <- stack(rgbRaster$band58,rgbRaster$band34,rgbRaster$band19) # write the geotiff # change overwrite=TRUE to FALSE if you want to make sure you don't overwrite your files! writeRaster(orderRGBstack,"rgbRaster.tif","GTiff", overwrite=TRUE)
Import A Multi-Band Image into R
You can import a multi-band image into R too. To do this, you import the file as a stack rather than a raster (which brings in just one band). Let’s import the raster than we just created above.
# import multi-band raster as stack multiRasterS <- stack("rgbRaster.tif") # import multi-band raster direct to brick multiRasterB <- brick("rgbRaster.tif") # view raster plot(multiRasterB)