terra包进行遥感数据处理（3）

欢迎查看terra包进行遥感数据处理（1）（2）

library(terra)

## terra version 0.9.11

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

定义植被指数函数

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

ndvi直方图

NDVI值是标准化的，并且范围在-1到+1之间。较高的值表示更多的绿色覆盖。

#从直方图可看到峰值在0.25-0.3之间
hist(ndvi,
     main = "Distribution of NDVI values",
     xlab = "NDVI",
     ylab= "Frequency",
     col = "wheat",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05),
     labels = seq(-0.5,1, 0.05))

#NDVI值大于0.4的栅格认为是植被,
#可通过classify函数提取植被栅格。
veg <- classify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main='Vegetation')

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 = 'NDVI值为0.25-0.3')

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

主成分

为了降低数据的维数和噪声，有时需要对多光谱数据进行变换。主成分变换是一种通用的数据降维方法，它可以从一组较大的相关波段中创建几个不相关波段。

set.seed(1)
sr <- spatSample(landsat, 10000)
#植被和土壤线图
plot(sr[,c(4,5)], main = "NIR-Red plot")

pca <- prcomp(sr, scale = TRUE)
pca

screeplot(pca)

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

# 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)
par(mfrow=c(1,2))
plotRGB(landsatFCC, r = 1, g = 2, b = 3, stretch = "lin", axes=TRUE, main = "False Color")
plot(pcClass, main="PCA")