1985～2022 年中国各省市区县土地覆盖类型占比面板数据

今天给大家分享一份 1985～2022 年中国各省市区县土地覆盖类型面板数据。原始数据为栅格数据：「The 30 m annual land cover datasets and its dynamics in China from 1990 to 2022」，下载链接为：https://zenodo.org/record/8176941

每年的数据都非常大，最近有小伙伴想要把这份数据裁剪汇总成各省市区县的面板数据使用，因此我就处理了下。

直接从上面的链接即可下载到 1985 年和 1990～2022 年共 34 年的栅格数据，例如 2022 年的栅格数据是这样的：

该图是使用 R 语言绘制，绘图代码为：

library(terra)
library(tidyverse)
mycrs <- "+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
rast("CLCD/CLCD_v01_2022_albert.tif") -> rst
vect("九段线.geojson") %>%
  terra::project(mycrs) -> jdx
vect("海岸线/海岸线.shp") %>%
  terra::project(mycrs) -> hax
rst[rst == 0] <- NA

# png
png("2022年中国土地覆盖类型分布.png",
    width = 4000, height = 4000, res = 600)
terra::plot(jdx, axes = F)
title(main = "2022 年中国土地覆盖类型分布", cex.main = 1.5)
terra::plot(rst, col = c("#FAE39C", "#446F33", "#33A02C", "#ABD37B", "#1E69B4", "#A6CEE3", "#CFBDA3", "#E24290", "#289BE8"),
            type = "classes", legend = "bottomleft", axes = F,
            add = T, plg = list(legend = c(classdf$class)))
terra::plot(jdx, add = T, lwd = 1.5, axes = F)
terra::plot(hax, add = T, lwd = 1, col = "#709ae1", axes = F)
dev.off()

在之前的课程中我们讲解过如何把栅格数据分区域裁剪汇总成面板数据，感兴趣的小伙伴可以学习这两个课程：

不过这次的数据处理方法和课程里面讲解的有所差异，因为这个数据并不是要分区域求均值之类的统计量，而是分区域统计每种植被类型的像元数。该数据一共包含了如下几种类型：

library(tidyverse)
readxl::read_xlsx("CLCD_classificationsystem.xlsx") %>%
  knitr::kable(align = "c")

IDClassColor1Cropland250,227,1562Forest68,111,513Shrub51,160,444Grassland171,211,1235Water30,105,1806Sonw/Ice166,206,2277Barren207,189,1638Impervious226,66,1449Wetland40,155,232

所以这里数据处理的过程和之前课程中讲解的差异是我们需要编写一个 countfun 函数：

countfun <- function(x){
  x <- x[!is.na(x)]
  x %>%
    dplyr::as_tibble() %>%
    dplyr::count(value)
}

该函数会统计每个区域中不同取值的像元数量。然后就可以在 parLapply() 循环里面使用这个函数进行汇总栅格数据了：

parLapply(cl, 1:nrow(df), function(i){
  x <- df$x[i]
  y <- df$y[i]
  x %>%
    stringr::str_match("CLCD/CLCD_v01_(.*)_albert\\.tif") %>%
    .[,2] -> flname
  terra::vect(county) %>%
    terra::project("+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs") -> countyvect
  terra::rast(x) %>%
    terra::extract(countyvect[y],
                   fun = "countfun") %>%
    readr::write_rds(paste0("rds/", flname, "_", y, ".rds"))
})