This lesson is being piloted (Beta version)

Work With Multi-Band Rasters in R

Overview

Teaching: 40 min
Exercises: 20 min
Questions
  • How can I visualize individual and multiple bands in a raster object?

Objectives
  • Identify a single vs. a multi-band raster file.

  • Import multi-band rasters into R using the raster package.

  • Plot multi-band color image rasters in R using the ggplot package.

Things You’ll Need To Complete This Episode

See the lesson homepage for detailed information about the software, data, and other prerequisites you will need to work through the examples in this episode.

We introduced multi-band raster data in an earlier lesson. This episode explores how to import and plot a multi-band raster in R.

Getting Started with Multi-Band Data in R

In this episode, the multi-band data that we are working with is imagery collected using the NEON Airborne Observation Platform high resolution camera over the NEON Harvard Forest field site. Each RGB image is a 3-band raster. The same steps would apply to working with a multi-spectral image with 4 or more bands - like Landsat imagery.

If we read a RasterStack object into R using the raster() function, it only reads in the first band.

RGB_band1_HARV <- raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")

We need to convert this data to a data frame in order to plot it with ggplot.

RGB_band1_HARV_df  <- as.data.frame(RGB_band1_HARV, xy = TRUE)
ggplot() +
  geom_raster(data = RGB_band1_HARV_df,
              aes(x = x, y = y, alpha = HARV_RGB_Ortho)) + 
  coord_quickmap()

plot of chunk harv-rgb-band1

Challenge

View the attributes of this band. What are its dimensions, CRS, resolution, min and max values, and band number?

Solution

RGB_band1_HARV
class      : RasterLayer 
band       : 1  (of  3  bands)
dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /home/travis/build/datacarpentry/r-raster-vector-geospatial/_episodes_rmd/data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif 
names      : HARV_RGB_Ortho 
values     : 0, 255  (min, max)

Notice that when we look at the attributes of this band, we see: band: 1 (of 3 bands)

This is R telling us that this particular raster object has more bands (3) associated with it.

Data Tip

The number of bands associated with a raster object can also be determined using the nbands() function: syntax is nbands(RGB_band1_HARV).

Image Raster Data Values

As we saw in the previous exercise, this raster contains values between 0 and 255. These values represent degrees of brightness associated with the image band. In the case of a RGB image (red, green and blue), band 1 is the red band. When we plot the red band, larger numbers (towards 255) represent pixels with more red in them (a strong red reflection). Smaller numbers (towards 0) represent pixels with less red in them (less red was reflected). To plot an RGB image, we mix red + green + blue values into one single color to create a full color image - similar to the color image a digital camera creates.

Import A Specific Band

We can use the raster() function to import specific bands in our raster object by specifying which band we want with band = N (N represents the band number we want to work with). To import the green band, we would use band = 2.

RGB_band2_HARV <-  raster("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif", band = 2)

We can convert this data to a data frame and plot the same way we plotted the red band:

RGB_band2_HARV_df <- as.data.frame(RGB_band2_HARV, xy = TRUE)
ggplot() +
  geom_raster(data = RGB_band2_HARV_df,
              aes(x = x, y = y, alpha = HARV_RGB_Ortho)) + 
  coord_equal()

plot of chunk rgb-harv-band2

Challenge: Making Sense of Single Band Images

Compare the plots of band 1 (red) and band 2 (green). Is the forested area darker or lighter in band 2 (the green band) compared to band 1 (the red band)?

Solution

We’d expect a brighter value for the forest in band 2 (green) than in band 1 (red) because the leaves on trees of most often appear “green” - healthy leaves reflect MORE green light than red light.

Raster Stacks in R

Next, we will work with all three image bands (red, green and blue) as an R RasterStack object. We will then plot a 3-band composite, or full color, image.

To bring in all bands of a multi-band raster, we use thestack() function.

RGB_stack_HARV <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")

Let’s preview the attributes of our stack object:

RGB_stack_HARV
class      : RasterStack 
dimensions : 2317, 3073, 7120141, 3  (nrow, ncol, ncell, nlayers)
resolution : 0.25, 0.25  (x, y)
extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : HARV_RGB_Ortho.1, HARV_RGB_Ortho.2, HARV_RGB_Ortho.3 
min values :                0,                0,                0 
max values :              255,              255,              255 

We can view the attributes of each band in the stack in a single output:

RGB_stack_HARV@layers
[[1]]
class      : RasterLayer 
band       : 1  (of  3  bands)
dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /home/travis/build/datacarpentry/r-raster-vector-geospatial/_episodes_rmd/data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif 
names      : HARV_RGB_Ortho.1 
values     : 0, 255  (min, max)


[[2]]
class      : RasterLayer 
band       : 2  (of  3  bands)
dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /home/travis/build/datacarpentry/r-raster-vector-geospatial/_episodes_rmd/data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif 
names      : HARV_RGB_Ortho.2 
values     : 0, 255  (min, max)


[[3]]
class      : RasterLayer 
band       : 3  (of  3  bands)
dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /home/travis/build/datacarpentry/r-raster-vector-geospatial/_episodes_rmd/data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif 
names      : HARV_RGB_Ortho.3 
values     : 0, 255  (min, max)

If we have hundreds of bands, we can specify which band we’d like to view attributes for using an index value:

RGB_stack_HARV[[2]]
class      : RasterLayer 
band       : 2  (of  3  bands)
dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
resolution : 0.25, 0.25  (x, y)
extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /home/travis/build/datacarpentry/r-raster-vector-geospatial/_episodes_rmd/data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif 
names      : HARV_RGB_Ortho.2 
values     : 0, 255  (min, max)

We can also use the ggplot functions to plot the data in any layer of our RasterStack object. Remember, we need to convert to a data frame first.

RGB_stack_HARV_df  <- as.data.frame(RGB_stack_HARV, xy = TRUE)

Each band in our RasterStack gets its own column in the data frame. Thus we have:

str(RGB_stack_HARV_df)
'data.frame':	7120141 obs. of  5 variables:
 $ x               : num  731999 731999 731999 731999 732000 ...
 $ y               : num  4713535 4713535 4713535 4713535 4713535 ...
 $ HARV_RGB_Ortho.1: num  0 2 6 0 16 0 0 6 1 5 ...
 $ HARV_RGB_Ortho.2: num  1 0 9 0 5 0 4 2 1 0 ...
 $ HARV_RGB_Ortho.3: num  0 10 1 0 17 0 4 0 0 7 ...

Let’s create a histogram of the first band:

ggplot() +
  geom_histogram(data = RGB_stack_HARV_df, aes(HARV_RGB_Ortho.1))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

plot of chunk rgb-harv-hist-band1

And a raster plot of the second band:

ggplot() +
  geom_raster(data = RGB_stack_HARV_df,
              aes(x = x, y = y, alpha = HARV_RGB_Ortho.2)) + 
  coord_quickmap()

plot of chunk rgb-harv-plot-band2

We can access any individual band in the same way.

Create A Three Band Image

To render a final three band, colored image in R, we use the plotRGB() function.

This function allows us to:

  1. Identify what bands we want to render in the red, green and blue regions. The plotRGB() function defaults to a 1=red, 2=green, and 3=blue band order. However, you can define what bands you’d like to plot manually. Manual definition of bands is useful if you have, for example a near-infrared band and want to create a color infrared image.
  2. Adjust the stretch of the image to increase or decrease contrast.

Let’s plot our 3-band image. Note that we can use the plotRGB() function directly with our RasterStack object (we don’t need a dataframe as this function isn’t part of the ggplot2 package).

plotRGB(RGB_stack_HARV,
        r = 1, g = 2, b = 3)

plot of chunk plot-rgb-image

The image above looks pretty good. We can explore whether applying a stretch to the image might improve clarity and contrast using stretch="lin" or stretch="hist".

Image Stretch

When the range of pixel brightness values is closer to 0, a darker image is rendered by default. We can stretch the values to extend to the full 0-255 range of potential values to increase the visual contrast of the image.

Image Stretch light

When the range of pixel brightness values is closer to 255, a lighter image is rendered by default. We can stretch the values to extend to the full 0-255 range of potential values to increase the visual contrast of the image.

plotRGB(RGB_stack_HARV,
        r = 1, g = 2, b = 3,
        scale = 800,
        stretch = "lin")

plot of chunk plot-rbg-image-linear

plotRGB(RGB_stack_HARV,
        r = 1, g = 2, b = 3,
        scale = 800,
        stretch = "hist")

plot of chunk plot-rgb-image-hist

In this case, the stretch doesn’t enhance the contrast our image significantly given the distribution of reflectance (or brightness) values is distributed well between 0 and 255.

Challenge - NoData Values

Let’s explore what happens with NoData values when working with RasterStack objects and using the plotRGB() function. We will use the HARV_Ortho_wNA.tif GeoTIFF file in the NEON-DS-Airborne-Remote-Sensing/HARVRGB_Imagery/ directory.

  1. View the files attributes. Are there NoData values assigned for this file?
  2. If so, what is the NoData Value?
  3. How many bands does it have?
  4. Load the multi-band raster file into R.
  5. Plot the object as a true color image.
  6. What happened to the black edges in the data?
  7. What does this tell us about the difference in the data structure between HARV_Ortho_wNA.tif and HARV_RGB_Ortho.tif (R object RGB_stack). How can you check?

Answers

1) First we use the GDALinfo() function to view the data attributes.

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")
rows        2317 
columns     3073 
bands       3 
lower left origin.x        731998.5 
lower left origin.y        4712956 
res.x       0.25 
res.y       0.25 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
file        data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif 
apparent band summary:
   GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1 Float64           TRUE       -9999          1       3073
2 Float64           TRUE       -9999          1       3073
3 Float64           TRUE       -9999          1       3073
apparent band statistics:
  Bmin Bmax     Bmean      Bsd
1    0  255 107.83651 30.01918
2    0  255 130.09595 32.00168
3    0  255  95.75979 16.57704
Metadata:
AREA_OR_POINT=Area 

2) From the output above, we see that there are NoData values and they are assigned the value of -9999.

3) The data has three bands.

4) To read in the file, we will use the stack() function:

HARV_NA <- stack("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_Ortho_wNA.tif")

5) We can plot the data with the plotRGB() function:

plotRGB(HARV_NA,
        r = 1, g = 2, b = 3)

plot of chunk harv-na-rgb

6) The black edges are not plotted. 7) Both data sets have NoData values, however, in the RGB_stack the NoData value is not defined in the tiff tags, thus R renders them as black as the reflectance values are 0. The black edges in the other file are defined as -9999 and R renders them as NA.

GDALinfo("data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif")
rows        2317 
columns     3073 
bands       3 
lower left origin.x        731998.5 
lower left origin.y        4712956 
res.x       0.25 
res.y       0.25 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
file        data/NEON-DS-Airborne-Remote-Sensing/HARV/RGB_Imagery/HARV_RGB_Ortho.tif 
apparent band summary:
   GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1 Float64           TRUE   -1.7e+308          1       3073
2 Float64           TRUE   -1.7e+308          1       3073
3 Float64           TRUE   -1.7e+308          1       3073
apparent band statistics:
  Bmin Bmax Bmean Bsd
1    0  255   NaN NaN
2    0  255   NaN NaN
3    0  255   NaN NaN
Metadata:
AREA_OR_POINT=Area 

Data Tip

We can create a RasterStack from several, individual single-band GeoTIFFs too. We will do this in a later episode, Raster Time Series Data in R.

RasterStack vs RasterBrick in R

The R RasterStack and RasterBrick object types can both store multiple bands. However, how they store each band is different. The bands in a RasterStack are stored as links to raster data that is located somewhere on our computer. A RasterBrick contains all of the objects stored within the actual R object. In most cases, we can work with a RasterBrick in the same way we might work with a RasterStack. However a RasterBrick is often more efficient and faster to process - which is important when working with larger files.

More Resources

You can read the help for the brick() function by typing ?brick.

We can turn a RasterStack into a RasterBrick in R by using brick(StackName). Let’s use the object.size() function to compare RasterStack and RasterBrick objects. First we will check the size of our RasterStack object:

object.size(RGB_stack_HARV)
44320 bytes

Now we will create a RasterBrick object from our RasterStack data and view its size:

RGB_brick_HARV <- brick(RGB_stack_HARV)

object.size(RGB_brick_HARV)
170897168 bytes

Notice that in the RasterBrick, all of the bands are stored within the actual object. Thus, the RasterBrick object size is much larger than the RasterStack object.

You use the plotRGB() function to plot a RasterBrick too:

plotRGB(RGB_brick_HARV)

plot of chunk plot-brick

Challenge: What Functions Can Be Used on an R Object of a particular class?

We can view various functions (or methods) available to use on an R object with methods(class=class(objectNameHere)). Use this to figure out:

  1. What methods can be used on the RGB_stack_HARV object?
  2. What methods can be used on a single band within RGB_stack_HARV?
  3. Why do you think there is a difference?

Answers

1) We can see a list of all of the methods available for our RasterStack object:

methods(class=class(RGB_stack_HARV))
  [1] !              !=             [              [[            
  [5] [[<-           [<-            %in%           ==            
  [9] $              $<-            addLayer       adjacent      
 [13] aggregate      all.equal      animate        approxNA      
 [17] area           Arith          as.array       as.character  
 [21] as.data.frame  as.integer     as.list        as.logical    
 [25] as.matrix      as.vector      atan2          bbox          
 [29] boxplot        brick          calc           cellFromRowCol
 [33] cellFromXY     cellStats      clamp          click         
 [37] coerce         colFromCell    colFromX       colSums       
 [41] Compare        coordinates    corLocal       cover         
 [45] crop           crosstab       crs<-          cut           
 [49] cv             density        dim            dim<-         
 [53] disaggregate   dropLayer      extend         extent        
 [57] extract        flip           freq           getValues     
 [61] getValuesBlock getValuesFocal head           hist          
 [65] image          interpolate    intersect      is.factor     
 [69] is.finite      is.infinite    is.na          is.nan        
 [73] isLonLat       KML            labels         length        
 [77] levels         levels<-       log            Logic         
 [81] mask           match          Math           Math2         
 [85] maxValue       mean           merge          minValue      
 [89] modal          mosaic         names          names<-       
 [93] ncell          ncol           ncol<-         nlayers       
 [97] nrow           nrow<-         origin         origin<-      
[101] overlay        pairs          persp          plot          
[105] plotRGB        predict        print          proj4string   
[109] proj4string<-  quantile       raster         rasterize     
[113] readAll        readStart      readStop       reclassify    
[117] res            resample       rotate         rowColFromCell
[121] rowFromCell    rowFromY       rowSums        sampleRandom  
[125] sampleRegular  scale          select         setMinMax     
[129] setValues      shift          show           spplot        
[133] stack          stackSelect    subs           subset        
[137] Summary        summary        t              tail          
[141] text           trim           unique         unstack       
[145] values         values<-       weighted.mean  which.max     
[149] which.min      whiches.max    whiches.min    writeRaster   
[153] xFromCell      xFromCol       xmax           xmin          
[157] xres           xyFromCell     yFromCell      yFromRow      
[161] ymax           ymin           yres           zonal         
[165] zoom          
see '?methods' for accessing help and source code

2) And compare that with the methods available for a single band:

methods(class=class(RGB_stack_HARV[1]))
 [1] [             [<-           anyDuplicated as_tibble     as.data.frame
 [6] as.raster     as.tbl_cube   bind          boxplot       brick        
[11] coerce        coordinates   determinant   duplicated    edit         
[16] extent        extract       head          initialize    isSymmetric  
[21] Math          Math2         Ops           raster        rasterize    
[26] relist        subset        summary       surfaceArea   tail         
[31] trim          unique        weighted.mean writeValues  
see '?methods' for accessing help and source code

3) There are far more things one could or want to ask of a full stack than of a single band.

Key Points

  • A single raster file can contain multiple bands or layers.

  • Use the stack() function to load all bands in a multi-layer raster file into R.

  • Individual bands within a stack can be accessed, analyzed, and visualized using the same functions as single bands.