Hey there,

today we are going to have a look at the {raster} package for R. I will show you some of the functions of this package applied on Landsat 8 satellite imagery downloaded from the USGS Earthexplorer .

1. Install package and read in data

First we have to install the package and load the raster image we would like to process:

install.packages('raster') #install
library(raster) #load
setwd('./Landsat8') #set your working directory
imagelist <- list.files(pattern='TIF') #list with all .TIF-files in current folder

2. Stack layers (bands) and plot a single band

Now the time has arrived to use the first function from the {raster} packge: stack()

This basically just creates a collection of layers with the same spatial extent and resolution.

imagestack <- stack(imagelist)
imagestack

If you call your created image stack you get information about the number of layers, extent, projection, resolution, min/max values, ect.:

class       : RasterStack 
dimensions  : 7931, 7811, 61949041, 7  (nrow, ncol, ncell, nlayers)
resolution  : 30, 30  (x, y)
extent      : 434385, 668715, 5135685, 5373615  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names       : LC81900272013217LGN00_B1, LC81900272013217LGN00_B2, LC81900272013217LGN00_B3, LC81900272013217LGN00_B4, LC81900272013217LGN00_B5, LC81900272013217LGN00_B6, LC81900272013217LGN00_B7 
min values  : 0,0,0,0,0,0,0 
max values  : 65535,65535, 65535,65535,65535,65535,65535 

Lets do some simple plotting! Calling the base plot() function will allow us to generate a plot of a single band. Lets have a look at the Band 5 (Near Infrared):

plot(imagestack[[5]]) #double brackets needed!

As a result we get:

Band 5 (NIR) Basic Plot

The plot above already looks quite nice, for example we can clearly see low NIR values in the upper right corner of the image, indicating a lake (Neusiedlersee, Austria). But lets also have a look at multiple bands in one plot.

3. Plot multispectral information

For this purpose we will use the plotRGB() function, which can be tricky sometimes and cause some Errors. I.e. this:

Error in rgb(RGB[, 1], RGB[, 2], RGB[, 3], alpha = alpha, max = scale) 

A solution to this Error sometimes is to simply specify the scale parameter of the layers. The scale parameter stands for the maximum value in all three bands. In our case thats 65535 (look above at the imagestack information output).

I created four plots, two True Color  Composites of the bands 4,3,2 (one with a linear stretch) and two Near-Infrared False Color Composites 5,4,3 (one with a histogram stretch). With the parameter colNA we are able to specify the color of the NAvalues, in our case thats also the background color.

#Layout (you can skip this line)
par(mfrow=c(2,2))
#True Color RGB 432
plotRGB(imagestack, 4,3,2, scale=65535, colNA='white')
plotRGB(imagestack, 4,3,2, scale=65535, stretch='lin', colNA='white')
#False Color RGB 543
plotRGB(imagestack, 5,4,3, scale=65535, colNA='white')
plotRGB(imagestack, 5,4,3, scale=65535, stretch='hist', colNA='white')

The results look like this:

rgbPlot: False and True Color Compsites

Lets interpret the last NIR plot above:

We can clearly see larger, turquoise regions with low NIR values (in the upper right corner) surrounding the lake, which in this case represent soil on fields and some smaller settlements. As we move further to the west we can see more intensive redish colors, this is where the Wiener Wald (Vienna forest) starts.
For 3-4 lines of code, in my opninion, this is a quite impressive example of how powerful the {raster} package is for plotting raster images. In one of the following posts we will be looking on methods to filter raster images, handle NAvalues as well as write new raster images with R.
Fell free to ask any questions!

Cheers