在之前的课程「使用 R 语言处理 Merra2 数据获取各省市区县的比湿、降水量、风速和气压数据」中我们介绍了普通的栅格数据如何分省市区县进行裁剪和汇总。今天我们继续介绍分区域汇总计算栅格数据的面积。
附件中的IGBPLUCC_CN文件夹存放了 2001-2021 年中国500m分辨率IGBP土地覆被数据，该数据是徐老师处理的：中国500m分辨率IGBP土地覆被数据（2001-2021年） 今天的课程中我们将以这个数据为例讲解。
对于面积统计量，省份和城市的可以通过加总区县的得到，所以我们仅仅计算区县的就好了。
加载相应的 R 包：
library(tidyverse)
library(sf)
library(terra)
|
读取区县的矢量数据：
read_sf("2021行政区划/县.shp") %>%
select(-contains("类型")) -> county
county
|
#> Simple feature collection with 2877 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension: XY
#> Bounding box: xmin: 73.50235 ymin: 3.83703 xmax: 135.0957 ymax: 53.56362
#> Geodetic CRS: WGS 84
#> # A tibble: 2,877 × 7
#> 省 省代码 市 市代码 县 县代码 geometry
#> <chr> <dbl> <chr> <dbl> <chr> <dbl> <MULTIPOLYGON [°]>
#> 1 安徽省 340000 安庆市 340800 大观区 340803 (((117.0457 30.53298, 117.0528 30.…
#> 2 安徽省 340000 安庆市 340800 怀宁县 340822 (((116.9897 30.74026, 116.9893 30.…
#> 3 安徽省 340000 安庆市 340800 潜山市 340882 (((116.7036 31.04533, 116.7046 31.…
#> 4 安徽省 340000 安庆市 340800 宿松县 340826 (((116.4555 30.12877, 116.4619 30.…
#> 5 安徽省 340000 安庆市 340800 太湖县 340825 (((115.839 30.75738, 115.8464 30.7…
#> 6 安徽省 340000 安庆市 340800 桐城市 340881 (((117.2153 30.69533, 117.2175 30.…
#> 7 安徽省 340000 安庆市 340800 望江县 340827 (((116.8906 30.35758, 116.8984 30.…
#> 8 安徽省 340000 安庆市 340800 宜秀区 340811 (((117.0716 30.77248, 117.0955 30.…
#> 9 安徽省 340000 安庆市 340800 迎江区 340802 (((117.1982 30.69515, 117.1986 30.…
#> 10 安徽省 340000 安庆市 340800 岳西县 340828 (((116.5255 31.02535, 116.5256 31.…
#> # ℹ 2,867 more rows
|
索引所有的栅格数据：
fs::dir_ls("IGBPLUCC_CN", regexp = "tif$") -> ls
ls
|
#> IGBPLUCC_CN/2001_01_01.tif IGBPLUCC_CN/2002_01_01.tif
#> IGBPLUCC_CN/2003_01_01.tif IGBPLUCC_CN/2004_01_01.tif
#> IGBPLUCC_CN/2005_01_01.tif IGBPLUCC_CN/2006_01_01.tif
#> IGBPLUCC_CN/2007_01_01.tif IGBPLUCC_CN/2008_01_01.tif
#> IGBPLUCC_CN/2009_01_01.tif IGBPLUCC_CN/2010_01_01.tif
#> IGBPLUCC_CN/2011_01_01.tif IGBPLUCC_CN/2012_01_01.tif
#> IGBPLUCC_CN/2013_01_01.tif IGBPLUCC_CN/2014_01_01.tif
#> IGBPLUCC_CN/2015_01_01.tif IGBPLUCC_CN/2016_01_01.tif
#> IGBPLUCC_CN/2017_01_01.tif IGBPLUCC_CN/2018_01_01.tif
#> IGBPLUCC_CN/2019_01_01.tif IGBPLUCC_CN/2020_01_01.tif
#> IGBPLUCC_CN/2021_01_01.tif
|
读取第一个栅格数据：
#> class : SpatRaster
#> dimensions : 10518, 13712, 1 (nrow, ncol, nlyr)
#> resolution : 0.004491576, 0.004491576 (x, y)
#> extent : 73.50016, 135.0887, 6.319648, 53.56205 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source : 2001_01_01.tif
#> name : LC_Type1
#> min value : 1
#> max value : 17
|
上面读取的 county 是个 sf 对象，在参与 terra 包的栅格数据的运算时，需要使用 vect() 函数转换成 SpatVector 对象。虽然 vect() 其实也可以直接读取 shp 数据，不过使用 read_sf() 读取后的更方便处理。
vect(county)[1] 就表示第一个区县：
#> class : SpatVector
#> geometry : polygons
#> dimensions : 1, 6 (geometries, attributes)
#> extent : 116.745, 117.0528, 30.35715, 30.58828 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> names : 省 省代码 市 市代码 县 县代码
#> type : <chr> <num> <chr> <num> <chr> <num>
#> values : 安徽省 3.4e+05 安庆市 3.408e+05 大观区 3.408e+05
|
从栅格数据中裁剪出来这个区的部分：
rst %>%
mask(vect(county)[1]) %>%
crop(vect(county)[1]) -> rstn
plot(rstn)
|
使用 cellSize() 函数可以计算栅格数据的每个斑块的面积：
cellSize(rstn, unit = "km") -> arearst
|
在 crs = 4326 坐标系下，不同纬度的斑块面积是不同的。
plot(arearst)
plot(rstn, add = T)
|
这里可以使用 zonal() 函数分区统计不同覆被类型斑块的总面积：
zonal(arearst, rstn, fun = sum, na.rm = T)
|
#> LC_Type1 area
#> 1 5 0.4290403
#> 2 8 1.5020209
#> 3 9 89.9452921
#> 4 10 3.2201476
#> 5 11 17.8218621
#> 6 12 91.0787940
#> 7 13 28.9742177
#> 8 14 27.9211682
#> 9 16 1.5040513
#> 10 17 41.2358192
|
这里得到的结果就是大观区各种土地覆被类型的面积统计结果了。
然后我们就可以循环统计各个区县、各个年份的了。为了更快的计算，我们这里采用多线程，首先根据自己电脑的实际情况创建多线程，并把所需的对象和 R 包分发给各个线程：
library(parallel)
detectCores()
|
makeCluster(16) -> cl
clusterExport(cl, "county")
clusterExport(cl, "ls")
dir.create("res")
clusterEvalQ(cl, ({
library(tidyverse)
library(sf)
library(terra)
})) -> tempres
|
我们需要循环各年、各区县，为了充分发挥多线程的功能，下面创建一个交叉数据集循环：
as.character(ls) %>%
crossing(county$县代码) %>%
set_names("rstname", "ctyindex") -> cdf
cdf
|
#> # A tibble: 60,417 × 2
#> rstname ctyindex
#> <chr> <dbl>
#> 1 IGBPLUCC_CN/2001_01_01.tif 110101
#> 2 IGBPLUCC_CN/2001_01_01.tif 110102
#> 3 IGBPLUCC_CN/2001_01_01.tif 110105
#> 4 IGBPLUCC_CN/2001_01_01.tif 110106
#> 5 IGBPLUCC_CN/2001_01_01.tif 110107
#> 6 IGBPLUCC_CN/2001_01_01.tif 110108
#> 7 IGBPLUCC_CN/2001_01_01.tif 110109
#> 8 IGBPLUCC_CN/2001_01_01.tif 110111
#> 9 IGBPLUCC_CN/2001_01_01.tif 110112
#> 10 IGBPLUCC_CN/2001_01_01.tif 110113
#> # ℹ 60,407 more rows
|
循环过程中可能会因为各种原因中断，为了重新运行的时候跳过之前运行成功的，可以使用下面的程序：
dir.create("res")
cdf %>%
mutate(value = paste0("res/",
str_extract(rstname, "\\d{4}"),
"_", ctyindex, ".rds")) %>%
anti_join(
fs::dir_ls("res") %>%
as.character() %>%
as_tibble()
) -> rstfl
rstfl
|
#> # A tibble: 27 × 3
#> rstname ctyindex value
#> <chr> <dbl> <chr>
#> 1 IGBPLUCC_CN/2001_01_01.tif 110101 res/2001_110101.rds
#> 2 IGBPLUCC_CN/2001_01_01.tif 110102 res/2001_110102.rds
#> 3 IGBPLUCC_CN/2001_01_01.tif 110105 res/2001_110105.rds
#> 4 IGBPLUCC_CN/2001_01_01.tif 110106 res/2001_110106.rds
#> 5 IGBPLUCC_CN/2001_01_01.tif 110107 res/2001_110107.rds
#> 6 IGBPLUCC_CN/2001_01_01.tif 110108 res/2001_110108.rds
#> 7 IGBPLUCC_CN/2001_01_01.tif 110109 res/2001_110109.rds
#> 8 IGBPLUCC_CN/2001_01_01.tif 110111 res/2001_110111.rds
#> 9 IGBPLUCC_CN/2001_01_01.tif 110112 res/2001_110112.rds
#> 10 IGBPLUCC_CN/2001_01_01.tif 110113 res/2001_110113.rds
#> # ℹ 17 more rows
|
clusterExport(cl, "rstfl")
|
rstfl 中的 value 变量就是每次运行要保存的文件名称。
然后就可以使用 parLapply() 多线程循环了：
parLapply(cl, 1:nrow(rstfl), function(x){
rast(rstfl$rstname[x]) -> rst
county %>%
dplyr::filter(县代码 == rstfl$ctyindex[x]) %>%
vect() -> countytemp
try({
rst %>%
terra::crop(countytemp) %>%
terra::mask(countytemp) -> rsttemp
cellSize(rsttemp, unit = "km") -> newrst
zonal(newrst, rsttemp, fun = sum, na.rm = T) %>%
set_names("class", "area_km2") %>%
mutate(file = rstfl$value[x],
县代码 = countytemp$县代码) %>%
write_rds(rstfl$value[x])
rm(rst)
rm(countytemp)
rm(rsttemp)
rm(newrst)
})
}) -> res
|
这里我没有使用 if (!file.exists(rstfl$value[x])){} 跳过，是因为这种判断其实也会浪费时间，在要处理的文件特别多的情况下，不如这里用的 anti_join() 方法更高效。
运行完之后就可以合并所有的 rds 文件了：
parLapply(cl, fs::dir_ls("res"), readr::read_rds) %>%
bind_rows() %>%
as_tibble() -> df
df
|
#> # A tibble: 486,534 × 4
#> class area_km2 file 县代码
#> <dbl> <dbl> <chr> <dbl>
#> 1 13 50.9 res/2001_110101.rds 110101
#> 2 13 58.8 res/2001_110102.rds 110102
#> 3 17 0.191 res/2001_110102.rds 110102
#> 4 12 7.65 res/2001_110105.rds 110105
#> 5 13 496. res/2001_110105.rds 110105
#> 6 8 0.767 res/2001_110106.rds 110106
#> 7 9 1.34 res/2001_110106.rds 110106
#> 8 10 4.41 res/2001_110106.rds 110106
#> 9 11 0.958 res/2001_110106.rds 110106
#> 10 12 1.15 res/2001_110106.rds 110106
#> # ℹ 486,524 more rows
|
再处理下：
df %>%
as_tibble() %>%
tidyr::extract(file, "year", "/(\\d{4})_", convert = T) -> df
df
|
#> # A tibble: 486,534 × 4
#> class area_km2 year 县代码
#> <dbl> <dbl> <int> <dbl>
#> 1 13 50.9 2001 110101
#> 2 13 58.8 2001 110102
#> 3 17 0.191 2001 110102
#> 4 12 7.65 2001 110105
#> 5 13 496. 2001 110105
#> 6 8 0.767 2001 110106
#> 7 9 1.34 2001 110106
#> 8 10 4.41 2001 110106
#> 9 11 0.958 2001 110106
#> 10 12 1.15 2001 110106
#> # ℹ 486,524 more rows
|
df %>%
left_join(
county %>%
st_drop_geometry() %>%
select(-contains("类型"))
) %>%
rename(年份 = year, 该覆被类型区域面积_km2 = area_km2,
覆被类型 = class) %>%
select(5:9, 4, 3, everything()) %>%
filter(覆被类型 != 0) -> countydf
countydf
|
#> # A tibble: 486,534 × 9
#> 省 省代码 市 市代码 县 县代码 年份 覆被类型 该覆被类型区域面积_km2
#> <chr> <dbl> <chr> <dbl> <chr> <dbl> <int> <dbl> <dbl>
#> 1 北京市 110000 北京… 110000 东城… 110101 2001 13 50.9
#> 2 北京市 110000 北京… 110000 西城… 110102 2001 13 58.8
#> 3 北京市 110000 北京… 110000 西城… 110102 2001 17 0.191
#> 4 北京市 110000 北京… 110000 朝阳… 110105 2001 12 7.65
#> 5 北京市 110000 北京… 110000 朝阳… 110105 2001 13 496.
#> 6 北京市 110000 北京… 110000 丰台… 110106 2001 8 0.767
#> 7 北京市 110000 北京… 110000 丰台… 110106 2001 9 1.34
#> 8 北京市 110000 北京… 110000 丰台… 110106 2001 10 4.41
#> 9 北京市 110000 北京… 110000 丰台… 110106 2001 11 0.958
#> 10 北京市 110000 北京… 110000 丰台… 110106 2001 12 1.15
#> # ℹ 486,524 more rows
|
countydf %>%
haven::write_dta("2001~2021年各区县不同作物覆被类型的土地面积.dta", label = "数据处理：微信公众号 RStata")
|
城市的通过加总区县即可得到：
countydf %>%
group_by(省, 省代码, 市, 市代码, 年份, 覆被类型) %>%
summarise(该覆被类型区域面积_km2 = sum(该覆被类型区域面积_km2, na.rm = T)) %>%
ungroup() -> citydf
citydf
|
#> # A tibble: 84,449 × 7
#> 省 省代码 市 市代码 年份 覆被类型 该覆被类型区域面积_km2
#> <chr> <dbl> <chr> <dbl> <int> <dbl> <dbl>
#> 1 上海市 310000 上海市 310000 2001 9 518.
#> 2 上海市 310000 上海市 310000 2001 10 119.
#> 3 上海市 310000 上海市 310000 2001 11 469.
#> 4 上海市 310000 上海市 310000 2001 12 2641.
#> 5 上海市 310000 上海市 310000 2001 13 2727.
#> 6 上海市 310000 上海市 310000 2001 14 575.
#> 7 上海市 310000 上海市 310000 2001 16 40.5
#> 8 上海市 310000 上海市 310000 2001 17 874.
#> 9 上海市 310000 上海市 310000 2002 9 784.
#> 10 上海市 310000 上海市 310000 2002 10 134.
#> # ℹ 84,439 more rows
|
citydf %>%
haven::write_dta("2001~2021年各城市不同作物覆被类型的土地面积.dta", label = "数据处理：微信公众号 RStata")
|
省份的通过加总城市的即可得到：
citydf %>%
group_by(省, 省代码, 年份, 覆被类型) %>%
summarise(该覆被类型区域面积_km2 = sum(该覆被类型区域面积_km2, na.rm = T)) %>%
ungroup() -> provdf
provdf %>%
haven::write_dta("2001~2021年各省份不同作物覆被类型的土地面积.dta", label = "数据处理：微信公众号 RStata")
|
这样就完成了栅格数据的分区域面积汇总～
点击这里跳转到 RStata 短书平台获取附件：栅格数据的裁剪与面积汇总——基于 R 语言的方法
评论