QGIS+R语言栅格土地覆被转移矩阵计算与土地覆被面积变化柱状图、和弦图绘制，计算多期土地覆被不变区域

以前给大家介绍过土地覆被数据处理和转移矩阵制作方法，主要是依靠ArcGIS，操作比较麻烦，数据更改后得重新点很多次，不如编程方法一次编程后面可以稍微改一下输入重复运行方便，而且R语言可以直接统计出图，所以这次我使用QGIS+R语言来完成土地覆被数据转移矩阵计算和统计绘图工作。

• 使用ArcGIS完成土地覆被制图与转移矩阵计算，栅格矢量两种方法：
• 栅格数据土地覆被变化制图与土地利用转移矩阵
• 土地覆被变化制图与土地利用转移矩阵

terra包的重分类函数重分类的结果看着不太对，而且裁剪大数据量的栅格速度很慢，所以ESA CCI土地覆被的裁剪和重分类工作我选择QGIS来完成。

QGIS土地覆被数据重分类

本研究使用的是ESA CCI土地覆被数据，数据来源于欧空局，获取方式见往期推文：

• 土地覆被数据获取和处理方法合集

土地覆被分类体系选择

ESA CCI使用的是FAO的土地覆被分类体系，有22种类型，比较复杂，同时给出了和IPCC土地覆被变化分类的对照表，这里我选用了发表于《Applied Geography》的9种类型的分类（农田、森林、草地、湿地、居民地、稀疏植被、裸地、水域、永久冰雪），详见： 对历史上和新出现的全球土地覆盖变化的洞察力：ESA-CCI-LC数据集的案例

QGIS重分类

使用Reclassify by table工具进行重分类，如下图所示，输入输出设置比较简单，主要是重分类的表比较复杂，好在QGIS的工具可以复制设置信息为JSON，在后续的重分类中只需要粘贴设置即可。

重分类JSON，分类方法为发表于《Applied Geography》的9种类型的分类：

{
  "area_units": "m2",
  "distance_units": "meters",
  "ellipsoid": "EPSG:7030",
  "inputs": {
    "DATA_TYPE": 1,
    "INPUT_RASTER": "E:/R/Mongolia/ESALUCC/ESALUCC2020.tif",
    "NODATA_FOR_MISSING": false,
    "NO_DATA": -999.0,
    "OUTPUT": "TEMPORARY_OUTPUT",
    "RANGE_BOUNDARIES": 0,
    "RASTER_BAND": 1,
    "TABLE": [
      "9",
      "40",
      "1",
      "49",
      "100",
      "2",
      "159",
      "170",
      "2",
      "109",
      "111",
      "3",
      "129",
      "131",
      "3",
      "179",
      "181",
      "4",
      "189",
      "191",
      "5",
      "119",
      "123",
      "6",
      "139",
      "154",
      "6",
      "199",
      "203",
      "7",
      "209",
      "211",
      "8",
      "219",
      "221",
      "9"
    ]
  }
}

QGIS创建栅格属性表

为了方便查看土地覆被类型，可以创建栅格属性表，QGIS3.30以上版本默认有栅格属性表功能，低版本需要借助Raster Attribute Table插件。

安装好插件之后，先在QGIS图层属性的符号化选项卡中选择唯一值的可视化方式，进行Classify，完成栅格的像元值统计及调色工作，只有在完成像元统计后才能建立属性表

在栅格图层上右击创建栅格属性表，建议选择DBF格式，方便后续使用。之后打开就可以查看属性表，由于标注使用的汉字，属性表会出现乱码。

通过查看栅格文件我们可以发现，有一个cpg文件，这个cpg文件中存储有属性表的编码信息，我们可以使用记事本打开cpg文件，输入UTF-8

QGIS裁剪重分类工作流

QGIS提供了Graphical Modeler工具可以将多个处理过程组合成一个工作流，和ArcGIS的模型构建器有点类似，但是又不一样，经过实验我发现这个的组合结果是组合成了一个新的工具，在新的工具中可以指定输入和输出参数，这样就可以把裁剪和重分类工作结合到一起来，同时可以使用Batch processing批处理，还是比较方便的。

R语言计算土地覆被转移矩阵

使用R语言计算土地覆被转移矩阵，思路就是使用 terra 包读取栅格数据，然后将栅格转为数据框， tidyverse 分类汇总不同时期不同土地覆被像元值的组合，然后将土地覆被类别汉语名称挂接到对应的像元值，通过长数据转宽数据，宽数据转长数据等多次转换得到土地覆被转移矩阵，将多期转移矩阵组合为一个LIST，使用 openxlsx 包输出到一个EXCEL中，稍加整饰即可得到论文中的土地覆被转移矩阵。

读取土地覆被栅格数据并预览

library(terra)
library(tidyverse)
library(circlize)
library(ggsci)
library(readxl)

#读取三期土地覆被数据
ipcc2001 = rast("./ipcc2001.tif")
ipcc2010 = rast("./ipcc2010.tif")
ipcc2020 = rast("./ipcc2020.tif")
#将三期土地覆被合成为多波段SpatRaster：
ipcclc = c(ipcc2001, ipcc2010, ipcc2020)
names(ipcclc) = c("2001", "2010", "2020")
plot(ipcclc, nc = 3)

计算土地覆被转移矩阵并输出EXCEL

#读取土地覆被属性表
code_ipcc = foreign::read.dbf("./ipcc2001.tif.vat.dbf") %>% 
  select(c("Value", "Class"))

#计算土地覆被转移矩阵
#研究区2001-2010年土地覆被转移矩阵
ipccTM0110 = values(ipcclc) %>%     #将三期土地覆被数据转化为矩阵
  as.data.frame() %>% 
  na.omit() %>% 
  set_names("lc2001", "lc2010", "lc2020") %>% 
  select(lc2001, lc2010) %>%                 #选择2001，2010年数据
  filter(lc2001 >0) %>%                      #去除空值0的部分
  group_by(lc2001, lc2010) %>%              #分类汇总，类似于ArcGIS中的Combine
  summarise(., count = n()) %>% 
  left_join(code_ipcc, by = c("lc2001" = "Value")) %>% 
  left_join(code_ipcc, by = c("lc2010" = "Value"), suffix = c("_2001", "_2010")) %>% 
  select(Class_2001, Class_2010, count) %>% 
  mutate(Area = count*300*300/1000/1000) %>%     #Area字段将像元值计算为平方公里
  pivot_wider(names_from = Class_2010, values_from = Area, values_fill = 0) %>%
  pivot_longer(cols = -Class_2001, names_to = "Class_2010", values_to = "Area") %>%
  group_by(Class_2001, Class_2010) %>%
  summarise(Area_sum = sum(Area))%>%               #统计地类面积
  pivot_wider(names_from = Class_2010, values_from = Area_sum, names_prefix = "Class_2010_")   #将长数据转宽数据，转移矩阵形式

#研究区2010-2020年土地覆被转移矩阵
ipccTM1020 = values(ipcclc) %>%     #将三期土地覆被数据转化为矩阵
  as.data.frame() %>% 
  na.omit() %>% 
  set_names("lc2001", "lc2010", "lc2020") %>% 
  select(lc2010, lc2020) %>%                 #选择2010，2020年数据
  filter(lc2010 >0) %>%                      #去除空值0的部分
  group_by(lc2010, lc2020) %>%              #分类汇总，类似于ArcGIS中的Combine
  summarise(., count = n()) %>% 
  left_join(code_ipcc, by = c("lc2010" = "Value")) %>%      #给像元值挂接地类名称
  left_join(code_ipcc, by = c("lc2020" = "Value"), suffix = c("_2010", "_2020")) %>% 
  select(Class_2010, Class_2020, count) %>% 
  mutate(Area = count*300*300/1000/1000) %>%     #Area字段将像元值计算为平方公里
  pivot_wider(names_from = Class_2020, values_from = Area, values_fill = 0) %>%
  pivot_longer(cols = -Class_2010, names_to = "Class_2020", values_to = "Area") %>%
  group_by(Class_2010, Class_2020) %>%
  summarise(Area_sum = sum(Area))%>%               #统计地类面积
  pivot_wider(names_from = Class_2020, values_from = Area_sum, names_prefix = 




    
"Class_2020_")   #将长数据转宽数据，转移矩阵形式

#研究区2001-2020年土地覆被转移矩阵
ipccTM0120 = values(ipcclc) %>%     #将三期土地覆被数据转化为矩阵
  as.data.frame() %>% 
  na.omit() %>% 
  set_names("lc2001", "lc2010", "lc2020") %>% 
  select(lc2001, lc2020) %>%                 #选择2001，2020年数据
  filter(lc2001 >0) %>%                      #去除空值0的部分
  group_by(lc2001, lc2020) %>%              #分类汇总，类似于ArcGIS中的Combine
  summarise(., count = n()) %>% 
  left_join(code_ipcc, by = c("lc2001" = "Value")) %>% 
  left_join(code_ipcc, by = c("lc2020" = "Value"), suffix = c("_2001", "_2020")) %>% 
  select(Class_2001, Class_2020, count) %>% 
  mutate(Area = count*300*300/1000/1000) %>%     #Area字段将像元值计算为平方公里
  pivot_wider(names_from = Class_2020, values_from = Area, values_fill = 0) %>%
  pivot_longer(cols = -Class_2001, names_to = "Class_2020", values_to = "Area") %>%
  group_by(Class_2001, Class_2020) %>%
  summarise(Area_sum = sum(Area))%>%               #统计地类面积
  pivot_wider(names_from = Class_2020, values_from = Area_sum, names_prefix = "Class_2020_")   #将长数据转宽数据，转移矩阵形式

#输出土地覆被转移矩阵至EXCEL
ipccxlsx = list(ipccTM0110, ipccTM1020, ipccTM0120)
names(ipccxlsx) = c("ipccTM0110", "ipccTM1020", "ipccTM0120")  #修改list中的表名，对应EXCEL中sheet名
openxlsx::write.xlsx(ipccxlsx, "./lcTrMat20012020.xlsx")

R语言绘制土地覆被面积变化柱状图

#计算不同时期土地覆被面积绘制柱状图
ipcclcdf = values(ipcclc) %>% 
  as.data.frame() %>% 
  na.omit() %>% 
  pivot_longer(cols = c("2001", "2010", "2020"), names_to = "year", values_to = "value") %>%  #
  group_by(value, year) %>% 
  summarise(np = n()) %>%              #计算不同土地覆被类型对应的像元数
  left_join(code_ipcc, by = c("value" = "Value")) %>%    #将土地覆被汉语类型和Value中数字对应起来
  mutate(area = np*300*300/1000/1000) %>% 
  na.omit()

openxlsx::write.xlsx(ipcclcdf, "./ipcclc_Russiadf.xlsx")
# 绘制柱状图
p1 = ggplot(ipcclcdf, aes(x = Class, y = area, fill = as.factor(year))) +
  geom_col(position = "dodge") +
  theme_bw()+
  labs(x = "土地覆被类型", y = "面积(km2)", fill = "Year")+
  theme(legend.title=element_blank(),
        legend.position = c(0.9, 0.85),
        legend.text = element_text(size = 8),
        legend.background = element_blank(),
        axis.text = element_text(size = 8),
        axis.title = element_text(size = 10),
        text=element_text(family="SimSun"))+
  scale_color_npg()+
  scale_fill_npg()
p1
ggsave("./土地覆被面积.jpg", p1, width = 14, height = 10, units = "cm", dpi = 600)

R语言绘制土地覆被转移矩阵和弦图

绘制和弦图，输入数据，我是单独又存了一个EXCEL，这个EXCEL是整理后的，如下图所示，第一行、第一列都是地类名称，其他是转移矩阵的数值

和弦图的绘制我使用了 circlize 包来完成，这个包的细节设置比较复杂，我也不太熟悉，我选择输出PDF文件在AI里可以手动修饰一下。

a1 = read_excel(path = "./ESALUCC/plot/test.xlsx", sheet = 1, col_names = T) %>% 
  as.data.frame()
rownames(a1) = a1[,1]
a1[,1] 
a1 = as.matrix(a1)

a2 = read_excel(path = "./ESALUCC/plot/test.xlsx", sheet = 2, col_names = T) %>% 
  as.data.frame()
rownames(a2) = a2[,1]
a2[,1] 
a2 = as.matrix(a2)

a3 = read_excel(path = "./ESALUCC/plot/test.xlsx", sheet = 3, col_names = T) %>% 
  as.data.frame()
rownames(a3) = a3[,1]
a3[,1] 
a3 = as.matrix(a3)

grid_col = c(农田 = "#13cf55", 森林 = "#17aa03", 草地 = "#8fc91c", 湿地 = "#08dca0", 居民地 = "#de1171",
             灌丛 = "#455dca", 稀疏植被 = "#cfe145", 裸地 = "#c0cecb", 水域 = "#007df2", 永久积雪 = "grey")
chordDiagram(a1, grid.col = grid_col)

par(mfrow = c(1,3))
chordDiagram(a1, grid.col = grid_col, annotationTrack = c("grid"))
circos.track(track.index = 1, panel.fun = function(x, y) {
  circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, 
              facing = "clockwise", niceFacing = TRUE, adj = c(0.2, 0.5))
}, bg.border = NA)
chordDiagram(a2, grid.col = grid_col, annotationTrack = c("grid"))
circos.track(track.index = 1, panel.fun = function(x, y) {
  circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, 
              facing = "clockwise", niceFacing = TRUE, adj = c(0.2, 0.5))
}, bg.border = NA)
chordDiagram(a3, grid.col = grid_col, annotationTrack = c("grid"))
circos.track(track.index = 1, panel.fun = function(x, y) {
  circos.text(CELL_META$xcenter, CELL_META$ylim[1], CELL_META$sector.index, 
              facing = "clockwise", niceFacing = TRUE, adj = c(0.2, 0.5))
}, bg.border = NA)
circos.clear()

library(showtext)
#输出中文pdf
showtext_auto(enable = TRUE)
font_add('SimSun', 'simsun.ttc')   #添加中文字体
pdf('GFtest.pdf',width = 11.69,height = 11.69)   #输出PDF，指定长宽   
par(mfrow = c(1,3))
chordDiagram(a1, grid.col = grid_col)
chordDiagram(a2, grid.col = grid_col)
chordDiagram(a3, grid.col = grid_col)
circos.clear()
dev.off()

R语言计算土地覆被多年不变的部分

当一组数标准差为0时，我们可以认为这一组数都相等，根据这个，我们可以通过判断不同时期的土地覆被像元值标准差是否为0，获得土地覆被多年不变的区域。

代码的主要内容：

• 将多个时期的土地覆被读取为一个多波段的SpatRaster
• 判断每个像元位置对应的一组像元值标准差是否为0，为0的返回第一个像元值，不为零的也就是变化的设为空值
• 输出结果，输出的结果是不变区域的土地覆被类型，变化区域为空值

#计算土地覆被不发生变化的区域
lclist = fs::dir_ls(path = "./ESALUCCV2/", regexp = ".tif$")
lc = rast(lclist)
names(lc) = c("2001", "2005", "2010", "2015", "2020")
plot(lc)
#利用标准差等于0计算土地覆被不发生改变的部分
fun_sdeq = function(x){
  ifelse(sd(x)==0, x[1], NA)
}
lceq = app(lc, fun_sdeq, cores = 28)
writeRaster(lceq, "./ESALUCCV2/plot/lceq.tif")

点击 阅读原文 学习视频课程