In this chapter we describe how to access and explore attributes of remote sensing images with R. We also show how to plot (make maps).

We will primarily use a spatial subset of Landsat 8 scene collected on June 14, 2017. The subset covers the area between Concord and Stockton, in California, USA.

All Landsat image scenes have a unique product ID and metadata. You can find the information on Landsat sensor, satellite, location on Earth (WRS path, WRS row) and acquisition date from the product ID. For example, the product identifier of the data we will use is ‘LC08_044034_20170614’. Based on this guide, you can see that the Sensor-Satellite is OLI/TIRS combined Landsat 8, WRS Path 44, WRS Row 34 and collected on June 14, 2017. Landsat scenes are most commonly delivered as zipped file, which contains separate files for the reflectance values at each bandwidth (wavelength).

We will start by exploring and visualizing the data (See [chapter 1] for data downloading instructions if you have not already done so).

View image properties

Create RasterLayer objects for single Landsat layers (bands)

## Loading required package: sp
# Blue
b2 <- raster('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- raster('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- raster('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- raster('data/rs/LC08_044034_20170614_B5.tif')

Print the variables to check. E.g.

## class      : RasterLayer
## dimensions : 1245, 1497, 1863765  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source     : c:/github/rspatial/rspatial-web/source/rs/_R/data/rs/LC08_044034_20170614_B2.tif
## names      : LC08_044034_20170614_B2
## values     : 0.0748399, 0.7177562  (min, max)

You can see the spatial resolution, extent, number of layers, coordinate reference system and more.

Image information and statistics

The below shows how you can access various properties from a Raster* object (this is the same for any raster data set).

# coordinate reference system (CRS)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# Number of rows, columns, or cells
## [1] 1863765
## [1] 1245 1497    1
# spatial resolution
## [1] 30 30
# Number of bands
## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
## [1] TRUE

You can create a RasterStack (an object with multiple layers) from the existing RasterLayer (single band) objects.

landsatRGB <- stack(b4, b3, b2)
landsatFCC <- stack(b5, b4, b3)
# Check the properties of the RasterStack
## class      : RasterStack
## dimensions : 1245, 1497, 1863765, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names      : LC08_044034_20170614_B4, LC08_044034_20170614_B3, LC08_044034_20170614_B2
## min values :              0.02084067,              0.04259216,              0.07483990
## max values :               0.7861769,               0.6924697,               0.7177562

Notice the order of the layers. These are suitable for plotting different color composites. You can learn more about color composites in remote sensing here and also in the section below.

You can also create the RasterStack using the filenames.

# first create a list of raster layers to use
filenames <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
##  [1] "data/rs/LC08_044034_20170614_B1.tif"
##  [2] "data/rs/LC08_044034_20170614_B2.tif"
##  [3] "data/rs/LC08_044034_20170614_B3.tif"
##  [4] "data/rs/LC08_044034_20170614_B4.tif"
##  [5] "data/rs/LC08_044034_20170614_B5.tif"
##  [6] "data/rs/LC08_044034_20170614_B6.tif"
##  [7] "data/rs/LC08_044034_20170614_B7.tif"
##  [8] "data/rs/LC08_044034_20170614_B8.tif"
##  [9] "data/rs/LC08_044034_20170614_B9.tif"
## [10] "data/rs/LC08_044034_20170614_B10.tif"
## [11] "data/rs/LC08_044034_20170614_B11.tif"
landsat <- stack(filenames)
## class      : RasterStack
## dimensions : 1245, 1497, 1863765, 11  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names      : LC08_044034_20170614_B1, LC08_044034_20170614_B2, LC08_044034_20170614_B3, LC08_044034_20170614_B4, LC08_044034_20170614_B5, LC08_044034_20170614_B6, LC08_044034_20170614_B7, LC08_044034_20170614_B8, LC08_044034_20170614_B9, LC08_044034_20170614_B10, LC08_044034_20170614_B11
## min values :            9.641791e-02,            7.483990e-02,            4.259216e-02,            2.084067e-02,            8.457669e-04,           -7.872183e-03,           -5.052945e-03,            3.931751e-02,           -4.337332e-04,             2.897978e+02,             2.885000e+02
## max values :              0.73462820,              0.71775615,              0.69246972,              0.78617686,              1.01243150,              1.04320455,              1.11793602,              0.82673049,              0.03547901,             322.43139648,             317.99530029

Above we created a RasterStack with 11 layers. The layers represent reflection intensity in the following wavelengths: Ultra Blue, Blue, Green, Red, Near Infrared (NIR), Shortwave Infrared (SWIR) 1, Shortwave Infrared (SWIR) 2, Panchromatic, Cirrus, Thermal Infrared (TIRS) 1, Thermal Infrared (TIRS) 2. We won’t use the last four layers and we will learn how to remove those in following sections.

Visualize single and multi-band imagery

You can plot individual layers of a RasterStack of a multi-spectral image.

par(mfrow = c(2,2))
plot(b2, main = "Blue", col = gray(0:100 / 100))
plot(b3, main = "Green", col = gray(0:100 / 100))
plot(b4, main = "Red", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))


Check the legends. They represent the range of values in each layer, and they range from 0 and 1. Notice the difference in shading and range of legends between the different bands. This is because different Earth surface features reflect the incident solar radiation differently. Each layer represent how much incident solar radiation is reflected for that particular wavelength. For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter in NIR wavelength. However water absorbs most of the incident energy in NIR regions and appears dark.

However we did not get much information from these grey-scale plots, but these are often combined to create more interesting plots. To combine three bands, we can use plotRGB. To make a “true/natural color” image (that is, something that looks like a normal photograph), we need to select the bands that we want to render in the red, green and blue regions. For this Landsat image, r = 3 (red), g = 2(green), b = 1(blue) will plot the true color composite (vegetation in green, water blue etc). You can also supply additional arguments to plotRGB to improve the visualization (e.g. a linear stretch of the values, using strecth = "lin").

plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")


This RGB-plot has lot more information than the previous ones, although it has been created using the same bands. Another popular image visualization method in remote sensing is known “false color” image where red, green and blue bands are replaced by other bands. Selecting r = NIR, g = red, b = green will plot a “false color”” composite. This representation is popular as it makes it easy to see the vegetation (in red).

par(mfrow = c(1,2))
plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat True Color Composite")
plotRGB(landsatFCC, r=1, g=2, b=3, axes=TRUE, stretch="lin", main="Landsat False Color Composite")


Note: Always check for package documentation (help(plotRGB)) for other arguments that can be added (like scale) to improve or modify the image.

Exercise 1 Use the RasterStack landsat to create a true and false color composite (hint remember the position of the bands in the stack).

Subset and rename spectral bands

You can select specific layers (bands) using subset function, or via indexing.

# select first 3 bands only
landsatsub1 <- subset(landsat, 1:3)
# same
landsatsub2 <- landsat[[1:3]]
# Number of bands in orginal and new data
## [1] 11
## [1] 3
## [1] 3

As mentioned above, we have no use here for the last four bands inlandsat. You can remove those using

landsat <- subset(landsat, 1:7)

Set the names of the bands using the following:

## [1] "LC08_044034_20170614_B1" "LC08_044034_20170614_B2"
## [3] "LC08_044034_20170614_B3" "LC08_044034_20170614_B4"
## [5] "LC08_044034_20170614_B5" "LC08_044034_20170614_B6"
## [7] "LC08_044034_20170614_B7"
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
## [1] "ultra.blue" "blue"       "green"      "red"        "NIR"
## [6] "SWIR1"      "SWIR2"

Spatial subset or crop

Spatial subsetting can be used to limit analysis to a geographic subset of the image. Spatial subsets can be created with the crop function, using an extent object, or another spatial object from which an Extent can be extracted/created,.

# Using extent
## class      : Extent
## xmin       : 594090
## xmax       : 639000
## ymin       : 4190190
## ymax       : 4227540
e <- extent(624387, 635752, 4200047, 4210939)
# crop landsat by the extent
landsatcrop <- crop(landsat, e)

Exercise 2 Interactive selection from the image is also possible. Use drawExtent and drawPoly to select an area of interest. Note: drawing will not work for plots generated within the markdown document. Please run the plot and draw commands from the RStudio console.

Exercise 3 Use the RasterStack landsatcrop to create a true and false color composite (hint remember the position of the bands in the stack).

Saving results to disk

At this stage we may want to save the raster to disk using the function writeRaster. Multiple file types are supported. We will use the commonly used GeoTiff format. While the layer order is preserved, layer names are unfortunately lost in the GeoTiff format.

writeRaster(landsatcrop,filename = "cropped-landsat.tif", overwrite = TRUE)

To keep the layer names, you can used the ‘raster-grd’ format:

writeRaster(landsatcrop, filename = "cropped-landsat.grd", overwrite = TRUE)

The disadvantage of this format is that not many other programs can read the data, in contrast to the GeoTiff format.

Note: Check for package documentation (help(writeRaster)) for additional helpful arguments that can be added.

Relation between bands

A scatterplot matrix can be helpful in exploring relationships between raster layers. This can be done with the pairs() function of the raster package.

Plot of reflection in the ultra-blue wavelength against reflection in the blue wavelength.

pairs(landsatcrop[[1:2]], main = "Ultra-blue versus Blue")


Plot of reflection in the red wavelength against reflection in the NIR wavelength.

pairs(landsatcrop[[4:5]], main = "Red versus NIR")


The first plot reveals high correlations between the blue wavelength regions. Because of the high correlation, we can just use one of the Blue bands without losing much information.

This distribution of points in second plot (between NIR and Red) is unique due to its triangular shape. Vegetation reflects very highly in the NIR range than red and creates the upper corner close to NIR (y) axis. Water absorbs energy from all the bands and occupies the location close to origin. The furthest corner is created due to highly reflecting surface features like bright soil or concrete.

Extract raster values

Often we want to get the values of raster cells for specific geographic locations or area. The extract function is used to get raster values at the locations of other spatial data. You can use points, lines, polygons or an Extent (rectangle) object. You can also use cell numbers to extract values. When using points, extract returns the values of a Raster* object for the cells in which a set of points fall.

# load the polygons with land use land cover information
samp <- readRDS('data/rs/samples.rds')
# generate 300 point samples from the polygons;
ptsamp <- spsample(samp, 300, type = 'random')
# add the class information to the point samples from polygons
ptsamp$class <- over(ptsamp, samp)$class
# extract values with points
df <- extract(landsat, ptsamp)
# To see some of the reflectance values
##      ultra.blue      blue      green        red        NIR      SWIR1
## [1,]  0.1379691 0.1290126 0.12775482 0.15560021 0.22853172 0.27767318
## [2,]  0.1101888 0.0879602 0.08026151 0.04788366 0.56883532 0.11875495
## [3,]  0.1352366 0.1306825 0.14317387 0.19207680 0.38185486 0.33379769
## [4,]  0.2017490 0.2012068 0.21020667 0.23414850 0.27179617 0.28823447
## [5,]  0.1333065 0.1141357 0.09225412 0.07171705 0.05042097 0.03623804
## [6,]  0.1105358 0.0874831 0.08321087 0.04701620 0.59245187 0.15438578
##           SWIR2
## [1,] 0.19363822
## [2,] 0.04150784
## [3,] 0.18745759
## [4,] 0.24750735
## [5,] 0.02932006
## [6,] 0.05449802

Spectral profiles

A plot of the spectrum (all bands) for pixels representing a certain earth surface features (e.g. water) is known as a spectral profile. Such profiles demonstrate the differences in spectral properties of various earth surface features and constitute the basis for image analysis. Spectral values can be extracted from any multispectral data set using extract function. In the above example, we extracted values of Landsat data for the samples. These samples include: cropland, water, fallow, built and open. First we compute the mean reflectance values for each class and each band.

ms <- aggregate(df, list(ptsamp$class), mean)
# instead of the first column, we use rownames
rownames(ms) <- ms[,1]
ms <- ms[,-1]
##          ultra.blue       blue      green        red        NIR      SWIR1
## built     0.1787321 0.17204419 0.17406668 0.19112260 0.24847852 0.24790147
## cropland  0.1124401 0.09043855 0.08602197 0.05467287 0.48582778 0.15340649
## fallow    0.1331698 0.11651040 0.10211767 0.10796453 0.16892134 0.22391816
## open      0.1391130 0.13724872 0.15152051 0.20347356 0.33900667 0.35123016
## water     0.1340993 0.11716628 0.09994984 0.07940147 0.04941895 0.03397931
##               SWIR2
## built    0.20483975
## cropland 0.06696569
## fallow   0.19192028
## open     0.21178264
## water    0.02748395

Now we will plot the mean spectra of these features.

# Create a vector of color for the land cover classes for use in plotting
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
#transform ms from a data.frame to a matrix
ms <- as.matrix(ms)
# First create an empty plot
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
# add the different classes
for (i in 1:nrow(ms)){
  lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
# Title
title(main="Spectral Profile from Landsat", font.main = 2)
# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")


The spectral profile shows (dis)similarity in the reflectance of different features on the earth’s surface (or above it). Unsurprisingly, the spectral signatures of ‘crop’ and ‘vegetation’ are similar. ‘Water’ shows relatively low reflection in all wavelengths, and ‘built’, ‘fallow’ and ‘open’ have relatively high reflectance in the longer wavelengts.