Exploration

In this chapter we describe how to explore satellite remote sensing data with R. We also show how to use them to make maps.

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 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 separate files for each band, combined into a single zip file.

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

Image properties

Create SpatRaster objects for single Landsat layers (bands)

library(terra)
# Blue
b2 <- rast('data/rs/LC08_044034_20170614_B2.tif')
# Green
b3 <- rast('data/rs/LC08_044034_20170614_B3.tif')
# Red
b4 <- rast('data/rs/LC08_044034_20170614_B4.tif')
# Near Infrared (NIR)
b5 <- rast('data/rs/LC08_044034_20170614_B5.tif')

Print the variables to check. E.g.

b2
## class       : SpatRaster
## dimensions  : 1245, 1497, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## source      : LC08_044034_20170614_B2.tif
## name        : LC08_044034_20170614_B2
## min value   :               0.0748399
## max value   :               0.7177562

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 SpatRaster object.

# coordinate reference system (CRS)
crs(b2)
## [1] "PROJCRS[\"WGS 84 / UTM zone 10N\",\n    BASEGEOGCRS[\"WGS 84\",\n        ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n            MEMBER[\"World Geodetic System 1984 (Transit)\"],\n            MEMBER[\"World Geodetic System 1984 (G730)\"],\n            MEMBER[\"World Geodetic System 1984 (G873)\"],\n            MEMBER[\"World Geodetic System 1984 (G1150)\"],\n            MEMBER[\"World Geodetic System 1984 (G1674)\"],\n            MEMBER[\"World Geodetic System 1984 (G1762)\"],\n            MEMBER[\"World Geodetic System 1984 (G2139)\"],\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]],\n            ENSEMBLEACCURACY[2.0]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 10N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-123,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Engineering survey, topographic mapping.\"],\n        AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],\n        BBOX[0,-126,84,-120]],\n    ID[\"EPSG\",32610]]"
# Number of cells, rows, columns
ncell(b2)
## [1] 1863765
dim(b2)
## [1] 1245 1497    1
# spatial resolution
res(b2)
## [1] 30 30
# Number of layers (bands in remote sensing jargon)
nlyr(b2)
## [1] 1
# Do the bands have the same extent, number of rows and columns, projection, resolution, and origin
compareGeom(b2,b3)
## [1] TRUE

You can create a SpatRaster with multiple layers from the existing SpatRaster (single layer) objects.

s <- c(b5, b4, b3)
# Check the properties of the multi-band image
s
## class       : SpatRaster
## dimensions  : 1245, 1497, 3  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## sources     : LC08_044034_20170614_B5.tif
##               LC08_044034_20170614_B4.tif
##               LC08_044034_20170614_B3.tif
## names       : LC08_04403~0170614_B5, LC08_04403~0170614_B4, LC08_04403~0170614_B3
## min values  :          0.0008457669,            0.02084067,            0.04259216
## max values  :          1.0124315023,            0.78617686,            0.69246972

You can also create the multi-layer SpatRaster using the filenames.

# first create a list of raster layers to use
filenames <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
filenames
##  [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 <- rast(filenames)
landsat
## class       : SpatRaster
## dimensions  : 1245, 1497, 11  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610)
## sources     : LC08_044034_20170614_B1.tif
##               LC08_044034_20170614_B2.tif
##               LC08_044034_20170614_B3.tif
##               ... and 8 more source(s)
## names       : LC08_~14_B1, LC08_~14_B2, LC08_~14_B3, LC08_~14_B4,  LC08_~14_B5,  LC08_~14_B6, ...
## min values  :  0.09641791,   0.0748399,  0.04259216,  0.02084067, 0.0008457669, -0.007872183, ...
## max values  :  0.73462820,   0.7177562,  0.69246972,  0.78617686, 1.0124315023,  1.043204546, ...

Above we created a SpatRaster 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.

Single band and composite maps

You can plot individual layers 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))

image0

The legends of the maps created above can range between 0 and 1. Notice the difference in shading and range of legends between the different bands. This is because different surface features reflect the incident solar radiation differently. Each layer represent how much incident solar radiation is reflected for a particular wavelength range. For example, vegetation reflects more energy in NIR than other wavelengths and thus appears brighter. In contrast, water absorbs most of the energy in the NIR wavelength and it appears dark.

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. With plotRGB we can combine them into a single composite image. Note that use of strecth = "lin" (otherwise the image will be pitch-dark).

landsatRGB <- c(b4, b3, b2)
plotRGB(landsatRGB, stretch = "lin")

image1

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 <- c(b5, b4, b3)
plotRGB(landsatFCC, stretch="lin")

image2

Question 1: Use the plotRGB function with the multi-band (11 bands) ``landsat`` image to create a true and false color composite (hint remember the position of the bands in the image).

Subset and rename 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 the original and new data
nlyr(landsat)
## [1] 11
nlyr(landsatsub1)
## [1] 3
nlyr(landsatsub2)
## [1] 3

We won’t use the last four bands in landsat. You can remove those by selecting the ones we want.

landsat <- subset(landsat, 1:7)

For clarity, it is useful to set the names of the bands. (source)[https://www.usgs.gov/faqs/what-are-band-designations-landsat-satellites?qt-news_science_products=0#qt-news_science_products]

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 a SpatExtent object, or another spatial object from which an Extent can be extracted.

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

Question 2: Use the ``landsatcrop`` image to plot a true and false color composite

Saving results to disk

At this stage we may want to save the raster to disk with writeRaster. Multiple file types are supported. We will use the commonly used GeoTiff format.

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

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.

A plot of reflection in the ultra-blue wavelength against reflection in the blue wavelength.

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

image3

A plot of reflection in the red wavelength against reflection in the NIR wavelength.

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

image4

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 (see Baret et al)[http://www.ipgp.fr/~jacquemoud/publications/baret1993a.pdf].

Extract cell values

Often we want to get the values of raster cells (pixels in remote sensing jargon) 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 SpatRaster object for the cells in which a set of points fall.

# load the polygons with land use land cover information
samp <- readRDS('data/rs/lcsamples.rds')
# generate 50 point samples from the polygons
set.seed(555)
ptsamp <- spatSample(samp, 50, 'regular')
# We use the x-y coordinates to extract the spectral values for the locations
df <- extract(landsat, ptsamp)
# To see some of the reflectance values
head(df)
##   ID ultra-blue       blue      green        red       NIR     SWIR1      SWIR2
## 1  1  0.1152851 0.09450950 0.09054089 0.06332441 0.4896148 0.1536918 0.07551219
## 2  2  0.1158706 0.09381554 0.08429519 0.05584258 0.4386517 0.1711277 0.08882765
## 3  3  0.1151116 0.09260110 0.09442276 0.05775099 0.5586210 0.1705638 0.07067610
## 4  4  0.1100153 0.08713611 0.08986861 0.05057278 0.4605116 0.1765059 0.07436280
## 5  5  0.1109695 0.08767828 0.08394821 0.04983544 0.6520245 0.1484437 0.04799209
## 6  6  0.1070009 0.08368797 0.06970022 0.04556321 0.3409542 0.1594387 0.07312667

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[,-1], list(ptsamp$class), mean)
# instead of the first column, we use row names
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms
##          ultra-blue       blue      green        red        NIR      SWIR1
## built     0.1715777 0.16020858 0.15984533 0.17228251 0.22624379 0.21830656
## cropland  0.1123755 0.08990475 0.08546264 0.05381490 0.49006296 0.16329528
## fallow    0.1348246 0.11801761 0.10233108 0.10520815 0.15550623 0.23542800
## open      0.1396925 0.13829443 0.15342000 0.20654889 0.34112771 0.35808941
## water     0.1342292 0.11744489 0.10053341 0.07981103 0.04949535 0.03381014
##               SWIR2
## built    0.18381969
## cropland 0.07174958
## fallow   0.21624273
## open     0.21559214
## water    0.02739490

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 Signatures", font.main = 2)
# Legend
legend("topleft", rownames(ms),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

image5

The spectral signatures (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.