Basic mathematical operations

The terra package supports many mathematical operations. Math operations are generally performed per pixel (grid cell). First we will do some basic arithmetic operations to combine bands. In the first example we write a custom math function to calculate the Normalized Difference Vegetation Index (NDVI). Learn more about vegetation indices here and NDVI.

We use the same Landsat data as in Chapter 2.

library(terra)
## terra 1.6.48
rfiles <- paste0('data/rs/LC08_044034_20170614_B', 1:11, ".tif")
landsat <- rast(rfiles)
landsatRGB <- landsat[[c(4,3,2)]]
landsatFCC <- landsat[[c(5,4,3)]]

Vegetation indices

Let’s define a general function for a ratio based (vegetation) index. In the function below, img is a muti-layer SpatRaster object and i and k are the indices of the layers (layer numbers) used to compute the vegetation index.

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}
# For Landsat NIR = 5, red = 4.
ndvi <- vi(landsat, 5, 4)
plot(ndvi, col=rev(terrain.colors(10)), main = "NDVI")

image0

You can see the variation in greenness from the plot.

Below is an alternative way to accomplish this. First write a general function that can compute 2-layer NDVI type indices.

vi2 <- function(x, y) {
    (x - y) / (x + y)
}

And use that function as an argument in lapp

nir <- landsat[[5]]
red <- landsat[[4]]
ndvi2 <- lapp(c(nir, red), fun = vi2)
# or in one line
#ndvi2 <- lapp(landsat[[5:4]], fun=vi2)
plot(ndvi2, col=rev(terrain.colors(10)), main="Landsat-NDVI")

image1

Question 1: Adapt the code shown above to compute indices to identify i) water and ii) built-up. Hint: Use the spectral profile plot to find the bands having maximum and minimum reflectance for these two classes. Or read about ithere..

Histogram

We can explore the distribution of values contained within our raster using hist to produces a histogram. Histograms are often useful in identifying outliers and bad data values in our raster data.

hist(ndvi, main = "NDVI values", xlab = "NDVI", ylab= "Frequency",
    col = "wheat", xlim = c(-0.5, 1),  breaks = 30, xaxt = "n")
## Warning: [hist] a sample of54% of the cells was used
axis(side=1, at = seq(-0.6, 1, 0.2), labels = seq(-0.6, 1, 0.2))

image2

We will refer to this histogram for the following sub-section on thresholding.

Question 2: Make histograms of the values the vegetation indices developed in question 1.

Thresholding

We can apply basic rules to get an estimate of spatial extent of different Earth surface features. Note that NDVI values are standardized and ranges between -1 to +1. Higher values indicate more green cover.

Cells with NDVI values greater than 0.4 are definitely vegetation. The following operation masks all cells that are perhaps not vegetation (NDVI < 0.4).

veg <- clamp(ndvi, 0.4, values=FALSE)
plot(veg, main='Vegetation')

image3

Let’s map the area that corresponds to the peak between 0.25 and 0.3 in the NDVI histogram.

m <- c(-Inf, 0.25, NA,  0.25, 0.3, 1,  0.3, Inf, NA)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
land <- classify(ndvi, rcl)
plot(land, main = 'What is it?')

image4

You can plot land on top of the landsatRGB raster to find out.

plotRGB(landsatRGB, r=1, g=2, b=3, axes=TRUE, stretch="lin")
plot(land, add=TRUE, legend=FALSE)

image5

You can also create classes for different intensity of vegetation.

m <- c(-1,0.25, 0.3, 0.4, 0.5, 1)
vegc <- classify(ndvi, m)
plot(vegc, col = rev(terrain.colors(4)), main = 'NDVI based thresholding')

image6

Question 3: Is it possible to find water using thresholding of NDVI or any other indices?

Principal component analysis

Multi-spectral data are sometimes transformed to helps to reduce the dimensionality and noise in the data. The principal components transform is a generic data reduction method that can be used to create a few uncorrelated bands from a larger set of correlated bands.

You can calculate the same number of principal components as the number of input bands. The first principal component (PC) explains the largest percentage of variance and other PCs explain additional the variance in decreasing order.

set.seed(1)
sr <- spatSample(landsat, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")

image7

This is known as vegetation and soil-line plot (Same as the scatter plot in earlier section).

pca <- prcomp(sr, scale = TRUE)
pca
## Standard deviations (1, .., p=11):
##  [1] 2.53668771 1.40078059 1.08362915 0.92501695 0.54958227 0.41473655
##  [7] 0.27030406 0.12220817 0.08661844 0.04763013 0.03609876
##
## Rotation (n x k) = (11 x 11):
##                                PC1        PC2         PC3         PC4
## LC08_044034_20170614_B1  0.2937880  0.3642123 -0.28481332 -0.06676012
## LC08_044034_20170614_B2  0.3357803  0.3382035 -0.15857739 -0.03635911
## LC08_044034_20170614_B3  0.3612145  0.2650391  0.07275822  0.04141597
## LC08_044034_20170614_B4  0.3676650  0.1641879  0.10447910  0.03578747
## LC08_044034_20170614_B5  0.1591361 -0.1852694  0.71300561  0.32594984
## LC08_044034_20170614_B6  0.3472817 -0.2134275  0.22853576  0.15831890
## LC08_044034_20170614_B7  0.3499019 -0.2367876 -0.11725588  0.06147126
## LC08_044034_20170614_B8  0.3496278  0.1964745  0.08489385  0.02433066
## LC08_044034_20170614_B9  0.1325628 -0.1117693  0.33293128 -0.92447725
## LC08_044034_20170614_B10 0.2534434 -0.4831768 -0.30276292 -0.02013211
## LC08_044034_20170614_B11 0.2507597 -0.4850418 -0.30570381 -0.02197088
##                                  PC5         PC6         PC7          PC8
## LC08_044034_20170614_B1   0.49633695  0.17556592 -0.23610914  0.219084848
## LC08_044034_20170614_B2   0.22593198  0.09672106  0.06516109  0.191800572
## LC08_044034_20170614_B3   0.06775249 -0.01134076  0.29589321 -0.502516608
## LC08_044034_20170614_B4  -0.33203379 -0.06966965  0.60728435  0.005032259
## LC08_044034_20170614_B5   0.51515049 -0.06796061 -0.07734807 -0.094468368
## LC08_044034_20170614_B6  -0.28840709  0.33603872 -0.01251499  0.639691364
## LC08_044034_20170614_B7  -0.24635064  0.53401272 -0.39446741 -0.489769588
## LC08_044034_20170614_B8  -0.33312920 -0.65599341 -0.53808102  0.022845490
## LC08_044034_20170614_B9   0.04429935  0.04329372 -0.02425017 -0.001980250
## LC08_044034_20170614_B10  0.16429483 -0.24521079  0.13682816  0.057604000
## LC08_044034_20170614_B11  0.19644322 -0.24453814  0.11434312 -0.028609160
##                                    PC9          PC10         PC11
## LC08_044034_20170614_B1  -0.1023702943  0.5365796233  0.127065695
## LC08_044034_20170614_B2  -0.1218742739 -0.7691469170 -0.196283512
## LC08_044034_20170614_B3   0.6640135096  0.0551533999  0.059318851
## LC08_044034_20170614_B4  -0.5347992170  0.2350928333  0.021665820
## LC08_044034_20170614_B5  -0.1995974299 -0.0311387747  0.004029261
## LC08_044034_20170614_B6   0.3828069102  0.0639560393 -0.021372642
## LC08_044034_20170614_B7  -0.2438645765 -0.0551104634  0.011338146
## LC08_044034_20170614_B8   0.0005337710  0.0005066492 -0.001773935
## LC08_044034_20170614_B9  -0.0009025984 -0.0006276713  0.001823281
## LC08_044034_20170614_B10  0.0039273450 -0.1709760476  0.686896213
## LC08_044034_20170614_B11  0.0433142125  0.1576520508 -0.684766019
screeplot(pca)

image8

We use a function to restrict prediction to the first two principal components

pca_predict2 <- function(model, data, ...) {
  predict(model, data, ...)[,1:2]
}
pci <- predict(landsat, pca, fun=pca_predict2)
plot(pci)

image9

The first principal component highlights the boundaries between different land use classes. it is difficult to understand what the second principal component is highlighting. Lets try thresholding again:

# quick look at the histogram of second component
hist <- pci[[2]]
m <- c(-Inf,-3,NA, -3,-2,0, -2,-1,1, -1,0,2, 0,1,3, 1,2,4, 2,6,5, 6,Inf,NA)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
rcl
##      [,1] [,2] [,3]
## [1,] -Inf   -3   NA
## [2,]   -3   -2    0
## [3,]   -2   -1    1
## [4,]   -1    0    2
## [5,]    0    1    3
## [6,]    1    2    4
## [7,]    2    6    5
## [8,]    6  Inf   NA
pcClass <- classify(pci[[2]], rcl)

Now plot the results

par(mfrow=c(1,2))
plotRGB(landsatFCC, stretch = "lin", main="False Color", mar=c(3.1, 3.1, 2.1, 2.1))
plot(pcClass, main="PCA")

image10

To learn more about the information contained in the vegetation and soil line plots read this paper by Gitelson et al. Details about PCA and an extension of PCA in remote sensing, Tasseled-cap Transformation.