新版本在这里：https://rstata.duanshu.com/#/brief/course/890e946ea9454a66ad5215516588cc94
之前给大家分享过中国各省市区县的土地覆盖数据：
最近有培训班的小伙伴想要乡镇的土地覆盖数据，考虑到原始数据非常大，于是我就帮忙处理了下。
原始数据为栅格数据：「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 年的栅格数据。
由于这份数据对我的电脑来说也是非常巨大的，所以我这里实际上是一个个乡镇进行裁剪汇总的（所以即使用了多线程还是挺耗时的）。经过处理即可得到下面的文件：
1985～2022 年中国各乡镇土地覆被类型占比（百分比）.xlsx
1985～2022 年中国各乡镇土地覆被类型占比（百分比）.dta
其中 dta 格式的是供 Stata 读取的。数据包含了乡镇名称和行政区划代码、各种覆被类型占比、行政区划面积（直接根据地理矢量数据计算的），预览如下：
图表展示效果更好些。下图展示了 2022 年各乡镇耕地占比：
处理方法 在之前的课程中我们讲解过如何把栅格数据分区域裁剪汇总成面板数据，感兴趣的小伙伴可以学习这两个课程：
不过这次的数据处理方法和课程里面讲解的有所差异，因为这个数据并不是要分区域求均值之类的统计量，而是分区域统计每种植被类型的像元数。该数据一共包含了如下几种类型：
library( tidyverse) readxl:: read_xlsx( "CLCD_classificationsystem.xlsx" ) %>% knitr:: kable( align = "c" )
所以这里数据处理的过程和之前课程中讲解的差异是我们需要编写一个 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" ) ) } )
土地覆盖类型分布 我还绘制了一幅图来展示土地覆盖类型分布，以 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( "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( )
对绘图感兴趣的小伙伴可以参考这个课程学习：
点击这里跳转到 RStata 短书平台获取附件：旧版本｜1985～2022 年中国各乡镇土地覆盖类型面板数据
×
评论