Introduction

This is an R Markdown Notebook which illustrates and extends instructions from the tutorial on remote sensing image analysis written by Aniruddha Ghosh and Robert J. Hijmans.

In the original tutorial Ghosh and Hijmans describe how to access and explore satellite remote sensing data with R. They also show how to use them to make maps. We illustrate the results for every command. We also complement the original instructions using additional explanations and code.

We will primarily use a spatial subset of a 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 each band.

Study Area

We will start by exploring and visualizing the study area using the leaflet package. Using the Landsat Adquisition Tool, we find that the Landsat Path 44 Row 34 scene center correspond to latitude 37.910 and longitude -122.091.

library(leaflet)
m <- leaflet() %>% #setView(lng = -122.091, lat = 37.910, zoom = 10)
addTiles() %>%
addMarkers(lng = -122.091, lat = 37.910, popup = "The Landsat scene center")

m

Let’s assume that we have done some field work to collect land cover samples. These samples usually are stored using a vector geospatial format such as shapefile, geopackage, or keyhole markup language. After reading such data from R, we can save every dataset as a .rds object. Later, in another R session, we can use the readRDS function to load the saved object and to assign its contents to a variable.

# load the polygons with land use land cover information
lc_samples <- readRDS('./rs/samples.rds')

We can check now what is the class for the loaded object.

class(lc_samples)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

As leaflet uses the EPSG 4326 coordinate reference system, we need to check what is the CRS for the landcover samples.

library(sp)
proj4string(lc_samples)
## [1] "+proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Note that lc_samples uses a UTM projection. We need to project our data to longlat using spTransform from either rgdal or sp libraries.

lc_rep <- spTransform(lc_samples, CRS("+init=epsg:4326"))

Lets check the CRS for the new landcover object.

proj4string(lc_rep)
## [1] "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Leaflet allows to use different base map layers. We can take advantage of the tab-completion feature to select the preferred base map just scrolling through the list of 110 providers!

We can use the popup argument to display the landcover class associated to every polygon.

Also, we need to provide a SpatialPolygonsDataFrame in the addPolygons() call.

leaflet() %>%
  addProviderTiles(providers$Esri.WorldImagery, options= providerTileOptions(opacity = 0.99)) %>%
  addPolygons(data = lc_rep, popup=lc_rep$class,
    stroke = FALSE, fillOpacity = 0.5, smoothFactor = 0.5
  )

Image properties

Create RasterLayer objects for single Landsat layers (bands)

library(raster)
# Blue
(b2 <- raster('./rs/LC08_044034_20170614_B2.tif'))
## 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     : /Users/ials/Documents/unal/remotesensing/rs/LC08_044034_20170614_B2.tif 
## names      : LC08_044034_20170614_B2 
## values     : 0.0748399, 0.7177562  (min, max)
# Green
(b3 <- raster('./rs/LC08_044034_20170614_B3.tif'))
## 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     : /Users/ials/Documents/unal/remotesensing/rs/LC08_044034_20170614_B3.tif 
## names      : LC08_044034_20170614_B3 
## values     : 0.04259216, 0.6924697  (min, max)
# Red
(b4 <- raster('./rs/LC08_044034_20170614_B4.tif'))
## 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     : /Users/ials/Documents/unal/remotesensing/rs/LC08_044034_20170614_B4.tif 
## names      : LC08_044034_20170614_B4 
## values     : 0.02084067, 0.7861769  (min, max)
# Near Infrared (NIR)
(b5 <- raster('./rs/LC08_044034_20170614_B5.tif'))
## 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     : /Users/ials/Documents/unal/remotesensing/rs/LC08_044034_20170614_B5.tif 
## names      : LC08_044034_20170614_B5 
## values     : 0.0008457669, 1.012432  (min, max)

Check data values stored in each Landsat layer (band).

head(values(b2))
## [1] 0.09889016 0.10969001 0.11717184 0.11268274 0.11754051 0.11931879
head(values(b5))
## [1] 0.2429532 0.2839189 0.2794515 0.2641191 0.2868899 0.2949356

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(b2)
## CRS arguments:
##  +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# Number of cells, rows, columns
ncell(b2)
## [1] 1863765
# Dimensions  nrow, ncol, nbands
dim(b2)
## [1] 1245 1497    1
# spatial resolution
res(b2)
## [1] 30 30
# Number of bands
nlayers(b2)
## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareRaster(b2,b3)
## [1] TRUE

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

(s <- stack(b5, b4, b3))
## 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_B5, LC08_044034_20170614_B4, LC08_044034_20170614_B3 
## min values :            0.0008457669,            0.0208406653,            0.0425921641 
## max values :               1.0124315,               0.7861769,               0.6924697

We can also create the RasterStack using the filenames.

# first create a list of raster layers to use
(filenames <- paste0('./rs/LC08_044034_20170614_B', 1:11, ".tif"))
##  [1] "./rs/LC08_044034_20170614_B1.tif" 
##  [2] "./rs/LC08_044034_20170614_B2.tif" 
##  [3] "./rs/LC08_044034_20170614_B3.tif" 
##  [4] "./rs/LC08_044034_20170614_B4.tif" 
##  [5] "./rs/LC08_044034_20170614_B5.tif" 
##  [6] "./rs/LC08_044034_20170614_B6.tif" 
##  [7] "./rs/LC08_044034_20170614_B7.tif" 
##  [8] "./rs/LC08_044034_20170614_B8.tif" 
##  [9] "./rs/LC08_044034_20170614_B9.tif" 
## [10] "./rs/LC08_044034_20170614_B10.tif"
## [11] "./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 you will see how to remove those in following sections.

Single band and composite maps

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

plot(b2, main = "Blue", col = gray(0:100 / 100))

In case, we cannot see anything in the displayed band, we can use the stretch function of the raster library. See this link to understand what the function parameters.

plot(stretch(b2, minq=0.1, maxq=0.9), main = "Blue", col = gray(0:100 / 100))

plot(stretch(b3, minq=0.1, maxq=0.9), main = "Green", col = gray(0:100 / 100))

plot(stretch(b4, minq=0.1, maxq=0.9), main = "Red", col = gray(0:100 / 100))

plot(stretch(b5, minq=0.1, maxq=0.9), main = "NIR", col = gray(0:100 / 100))

We do not gain that much information from these grey-scale plots; they are often combined into a “composite” to create more interesting plots. You can learn more about color composites in remote sensing here and also in the section below.

To make a “true (or natural) color” image, that is, something that looks like a normal photograph (vegetation in green, water blue etc), we need bands in the red, green and blue regions. For this Landsat image, band 4 (red), 3 (green), and 2 (blue) can be used. The plotRGB method can be used to combine them into a single composite. You can also supply additional arguments to plotRGB to improve the visualization (e.g. a linear stretch of the values, using strecth = “lin”).

landsatRGB <- stack(b4, b3, b2)
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Landsat True Color Composite")

The true-color composite reveals much more about the landscape than the earlier gray images. Another popular image visualization method in remote sensing is known “false color” image in which NIR, red, and green bands are combined. This representation is popular as it makes it easy to see the vegetation (in red).

landsatFCC <- stack(b5, b4, b3)
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main="Landsat False Color Composite")

Subset and rename bands

We 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 the original and new data
nlayers(landsat)
## [1] 11
nlayers(landsatsub2)
## [1] 3

We won’t use the last four bands in landsat. You can remove those using

landsat <- subset(landsat, 1:7)

For clarity, it is useful to set the names of the bands.

names(landsat)
## [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')
names(landsat)
## [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.

# Using extent
extent(landsat)
## 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)
extent(landsatcrop)
## class      : Extent 
## xmin       : 624390 
## xmax       : 635760 
## ymin       : 4200060 
## ymax       : 4210950

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.

x1 <- writeRaster(landsatcrop, filename="./rs/cropped-landsat.tif", overwrite=TRUE)

Alternatively you can used the ‘raster-grd’ format.

x2 <- writeRaster(landsatcrop, filename="./rs/cropped-landsat.grd", overwrite=TRUE)

An advantage of this format is that it saves the layer names. 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.

Check number of bands.

nlayers(landsatcrop)
## [1] 7

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 pixel 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(lc_samples, 300, type='regular')

Check the raster library help to understand how the over function works.

# add the land cover class to the points
ptsamp$class <- over(ptsamp, lc_samples)$class

Check the raster library help to understand how the over function works.

# extract values with points
df <- extract(landsat, ptsamp)
# To see some of the reflectance values
head(df)
##      ultra.blue      blue      green        red       NIR     SWIR1
## [1,]  0.1317018 0.1236995 0.12671387 0.16566272 0.2926585 0.3119377
## [2,]  0.1462317 0.1423932 0.15499298 0.20029597 0.3212629 0.4008088
## [3,]  0.1367113 0.1198393 0.10719607 0.10695753 0.1846817 0.2182523
## [4,]  0.1331114 0.1152201 0.09830463 0.09617936 0.1685036 0.2047851
## [5,]  0.1381426 0.1418510 0.16323383 0.23143770 0.4074448 0.3517324
## [6,]  0.1421980 0.1428486 0.16037123 0.21816559 0.3528601 0.3732887
##          SWIR2
## [1,] 0.1954599
## [2,] 0.2716660
## [3,] 0.1815806
## [4,] 0.1684386
## [5,] 0.1980840
## [6,] 0.2296811

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 row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms

Now we 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). ‘Water’ shows relatively low reflection in all wavelengths, and ‘built’, ‘fallow’ and ‘open’ have relatively high reflectance in the longer wavelengts.