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

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

非监督分类

利用K-means算法进行影像分类

library(terra)

## terra version 0.9.11

landsat5 <- rast('data/rs/centralvalley-2011LT5.tif')
names(landsat5) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
ndvi <- (landsat5[['NIR']] - landsat5[['red']]) / (landsat5[['NIR']] + landsat5[['red']])
# 设定输出范围
e <- ext(-121.807, -121.725, 38.004, 38.072)
#裁剪
ndvi <- crop(ndvi, e)
ndvi

## class       : SpatRaster
## dimensions  : 252, 304, 1  (nrow, ncol, nlyr)
## resolution  : 0.0002694946, 0.0002694946  (x, y)
## extent      : -121.807, -121.725, 38.00413, 38.07204  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## data source : memory
## names       :        NIR
## min values  : -0.3360085
## max values  :  0.7756007

set.seed(99)
#栅格数据转为R中的vector格式
nr <- values(ndvi)
str(nr)

##  num [1:76608, 1] 0.245 0.236 0.272 0.277 0.277 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "NIR"

set.seed(99)
kmncluster <- kmeans(na.omit(nr), centers = 10, iter.max = 500,nstart = 5, algorithm="Lloyd")
str(kmncluster)

## List of 9
##  $ cluster     : int [1:76608] 4 4 3 3 3 3 3 4 4 4 ...
##  $ centers     : num [1:10, 1] 0.55425 0.00498 0.29997 0.20892 -0.20902 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:10] "1" "2" "3" "4" ...
##   .. ..$ : chr "NIR"
##  $ totss       : num 6459
##  $ withinss    : num [1:10] 5.69 6.13 4.91 4.9 5.75 ...
##  $ tot.withinss: num 55.8
##  $ betweenss   : num 6403
##  $ size        : int [1:10] 8932 4550 7156 6807 11672 8624 8736 5040 9893 5198
##  $ iter        : int 108
##  $ ifault      : NULL
##  - attr(*, "class")= chr "kmeans"

# 聚类结果变成栅格属性
knr <- ndvi
values(knr) <- kmncluster$cluster
knr

## class       : SpatRaster
## dimensions  : 252, 304, 1  (nrow, ncol, nlyr)
## resolution  : 0.0002694946, 0.0002694946  (x, y)
## extent      : -121.807, -121.725, 38.00413, 38.07204  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs
## data source : memory
## names       : NIR
## min values  :   1
## max values  :  10

#定义每类颜色显示
mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff",
             "#0000ff","#00ff00","#cbbeb5",
             "#c3ff5b", "#ff7373", "#00ff00", "#808080")
par(mfrow = c(1,2))
plot(ndvi, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
plot(knr, main = '分监督分类', col = mycolor )

监督分类法

有各种监督分类算法存在，算法的选择会影响分类结果。这里将探索CART算法。在监督分类中，可通过实地考察、参考空间数据或高分辨率图像(如谷歌地图上可用的图像)的解译，获得一些关于土地覆盖类型的先验知识。在研究区域内确定已知土地覆盖类型的地点。这些区域通常被称为训练站点，这些点的光谱属性被用来训练分类算法。 下面的例子使用了分类和回归树(CART)分类器，步骤如下:

• 创建用于分类的样点
• 从样本点的Landsat数据中提取单元值
• 使用训练样本训练分类器
• 使用训练好的模型对Landsat数据进行分类
• 评估模型的准确性

library(terra)
raslist <- paste0('data/rs/LC08_044034_20170614_B', 2:7, ".tif")
landsat <- rast(raslist)
names(landsat) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')

已知土地利用类型

library(sp)
# 已知面状土地利用类型数据
samp <- readRDS('data/rs/samples.rds')
plot(samp)
#显示类型
text(samp, samp$class)

set.seed(1)
#从samp随机产生1000个样本点
ptsamp <- spsample(samp, 1000, type='random')
# 给点加上土地利用类型信息
ptsamp$class <- over(ptsamp, samp)$class
# 把`ptsamp` 转换为 `SpatVector`类型
ptsamp <- vect(ptsamp)

提取样本点光谱信息

xy <- as.matrix(geom(ptsamp)[,c('x','y')])
df <- extract(landsat, xy)[,-1]

训练分类器

library(rpart)
sampdata <- data.frame(class = ptsamp$class, df)
cartmodel <- rpart(as.factor(class)~., data = sampdata,
                   method = 'class', minsplit = 5)
print(cartmodel)

## n= 1000
##
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
##
##  1) root 1000 612 water (0.07 0.099 0.16 0.28 0.39) 
##    2) NIR>=0.08310244 612 328 open (0.11 0.16 0.26 0.46 0) 
##      4) SWIR1< 0.2725443 303 159 fallow (0.19 0.33 0.48 0.0066 0) 
##        8) SWIR2< 0.1270392 100   1 cropland (0.01 0.99 0 0 0) *
##        9) SWIR2>=0.1270392 203  59 fallow (0.28 0 0.71 0.0099 0) 
##         18) blue>=0.1296632 58   1 built (0.98 0 0 0.017 0) *
##         19) blue< 0.1296632 145   1 fallow (0 0 0.99 0.0069 0) *
##      5) SWIR1>=0.2725443 309  27 open (0.039 0 0.049 0.91 0) *
##    3) NIR< 0.08310244 388   0 water (0 0 0 0 1) *

plot(cartmodel, uniform=TRUE,
     main="分类树")
text(cartmodel, cex = 1)

对影像分类

classified <- predict(landsat, cartmodel,
                      na.rm = TRUE)
classified

## class       : SpatRaster
## dimensions  : 1245, 1497, 5  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=10 +datum=WGS84 +units=m +no_defs
## data source : memory
## names       :     built,  cropland,    fallow,      open,     water
## min values  :         0,         0,         0,         0,         0
## max values  : 0.9827586, 0.9900000, 0.9931034, 0.9126214, 1.0000000

plot(classified)

class <- c("built","cropland","fallow","open","water")
mycolor <- c('darkred', 'yellow', 'burlywood', 'cyan', 'blue')
classdf <- data.frame(classvalue = c(1,2,3,4,5),
                      classnames = class,
                      color = mycolor)
lulc <- app(classified, fun = which.max)
lulcc <- as.factor(lulc)
levels(lulcc) <- as.character(classdf$classnames)
plot(lulcc)

模型评估

利用kappa评价模型精度。在该例中样本数据的土地利用类型有5类，依次以其中一类作为模型测试数据，另4类作为模型训练数据。

set.seed(99)
k<-5
#对每类数据随机进行采样
j<-sample(rep(1:k,each=round(nrow(sampdata))/k))
table(j)

## j
##   1   2   3   4   5
## 200 200 200 200 200

x <- list()
for (k in 1:5) {
    train <- sampdata[j!= k, ]
    test <- sampdata[j == k, ]
    cart <- rpart(as.factor(class)~., data=train,
                  method = 'class',minsplit = 5)
    #用训练数据的分类器对测试数据分类
    pclass <- predict(cart, test, na.rm = TRUE)
    # 概率最大的类
    pclass <- apply(pclass, 1, which.max)
    # 将测试数据已知类别与模拟生成的类别合并在一起
    x[[k]] <- cbind(test$class, as.integer(pclass))
}

计算混淆矩阵

y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
#混淆矩阵
conmat <- table(y)
# 行列名修改为土地利用类型
colnames(conmat) <- classdf$classnames
rownames(conmat) <- classdf$classnames
print(conmat)

##           predicted
## observed   built cropland fallow open water
##   built       56        1      1   12     0
##   cropland     0       98      1    0     0
##   fallow       0        1    143   15     0
##   open         1        0      2  281     0
##   water        0        0      0    0   388

总体精度

n <- sum(conmat)
n

## [1] 1000

# 对角线上为分类正确的个数
diag <- diag(conmat)
# 总体精度
OA <- sum(diag) / n
OA

## [1] 0.966

Kappa

rowsums <- apply(conmat, 1, sum)
#已知类型的每类个数
p <- rowsums / n
# predicted cases per class
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
kappa

## [1] 0.9530854

Producer and user accuracy

#Producer accuracy
PA <- diag / colsums
# User accuracy
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
knitr::kable(outAcc)

	producerAccuracy	userAccuracy
built	0.9824561	0.8000000
cropland	0.9800000	0.9898990
fallow	0.9727891	0.8993711
open	0.9123377	0.9894366
water	1.0000000	1.0000000