R {raster}: Introduction and basic plotting

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)

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 8 (NIR) Plot R

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)
#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

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!


About This Author

Martin was born in Czech Republic and studied at the University of Natural Resources and Life Sciences, Vienna. He is currently working at GeoVille - an Earth Observation Company based in Austria, specialised in Land Monitoring. His main interests are: Open-source applications like R, (geospatial) statistics and data-management, web-mapping and visualization. He loves travelling, geocaching, photography and sports.


You can post comments in this post.

  • Thanks for sharing! This is very helpful

    Liza le Roux 8 years ago Reply

  • Your code and contents are very clear and it is very helpful.

    If the data is hyperspectral(HS), then plotRGB() doesnt work I think so ,

    How the false color composite will work in HS dataset? Please assist me..

    Nandhini K 7 years ago Reply

    • Hey Nandhini! The function still works. You just need to specify which band to plot in red,blue and green. If you have a 4 band raster stack ‘ras’ and you want to plot the bands 4,3,2, this is how the function call looks like: plotRGB(ras, r=4, g=3, b=2)
      Cheers Martin

      Martin 7 years ago Reply

      • Thank you Martin…. it works well …. 🙂

        Nandhini K 7 years ago Reply

  • Thank you for this information.

    I have a question:

    I have a raster stack e.g x.
    In R I accidentally use the function plot(x,1,1) where the first 1 means first layer in stack x. But I cant figure out what the second 1 means. I played around with many numbers after the first 1 from 20 to 10000. It changes the appearance of the map but I am confused how is it working and what it means. Can someone please help out.

    Mun 6 years ago Reply

    • Dear Mun,
      thats a good question. If you load the raster package and type ‘?plot’ you will get to the plot help pages. Click on “Plot a Raster* object” and you will get to the raster plotting method. There you will find the answer to your question: The third argument is the “maxpixel” argument: “integer > 0. Maximum number of cells to use for the plot.”
      If you have troubles finding the help page form R, you can also use this link: https://artax.karlin.mff.cuni.cz/r-help/library/raster/html/plot.html

      Hope it helps!


      Martin 6 years ago Reply

  • Hi Martin
    I fond your blog very useful. I have raster stack with 84 individual raster layers. When I plot the raster stack only a few are plotted. Do you know why it might have happened?

    Dev 6 years ago Reply

    • Hey! There is an argument you can use inside the plot function:
      maxnl = integer. Maximum number of layers to plot (for a multi-layer object)
      when rasterStack is your stack, then:
      plot(rasterStack, maxnl=nlayers(rasterStack))
      should work.

      Good luck!

      Martin 6 years ago Reply

      • That worked. thanks a lot Martin

        Dev 6 years ago Reply

  • Hi Martin,

    Thanks for your blog which is very helpful. I’m trying to carry out a layer stack using Landsat 8, band 10 included since I want to perform land cover analysis and compare this band with the rest. However in R band 10 dimensions are twice the size of the other bands. I resampled it but pixel values are different in the end, though it’s dimensions match with the other bands using:
    b10resampled <- projectRaster(stack10,stack1_7,method = 'ngb')

    Is there a way I can stack them or I have to analyze it separately?

    Patricia 6 years ago Reply

  • How do you add coordinates and titles to the plot?

    Gorden Jiang 5 years ago Reply

  • How do you add coordinates and titles to the true colour images?

    Gorden Jiang 5 years ago Reply

Post A Reply