其他
R语言获取和处理MCD12Q1土地覆被数据
R语言处理MCD12Q1土地覆被数据
MCD12Q1数据介绍
MCD12Q1是基于MODIS的一套年度土地覆被数据，其中的IGBP分类较为常用。目前MCD12Q1数据为2001年至2019年，每年一期。具体IGBP分类如下：
IGBP分类
1常绿针叶林，树冠高于2m，覆盖度>60% 2常绿阔叶林 3落叶针叶林 4落叶阔叶林 5混合林，针叶阔叶在40-60%之间，树冠高>2米，覆盖度>60% 6茂盛的灌木，覆盖度>60% 7稀疏的灌木，覆盖度10-60% 8有林草原，具有草本和其它林下系统的土地，树覆盖度30-60%，树冠高大于2m 9稀树草原，具有草本和其它林下系统的土地，树木覆盖度10-30%，树冠高大于2m 10草地，草高<2米 11永久性湿地，30-60%的水域，植被>10% 12农田，60%以上耕地 13城市和建成区，至少30%的不透水面，包含建筑材料、沥青路面、车辆等 14农田/自然植被混合区域，较小块的耕地，有40-60%的部分是树木和草本植物 15永久性积雪，60%以上区域，一年10个月以上被积雪覆盖 16荒地，60%以上区域无植被覆盖，植被覆盖度低于10% 17水域，60%以上的水域
MCD12Q1数据获取
常规方法
从网上下载MODIS数据，然后使用HEG或者MRT等软件进行坐标系转换和拼接。具体可以看以前的两篇推文，不再多说，太麻烦了，反正我是不想用了，尤其是MODIS正弦投影的转换，太麻烦，HEG或者MRT软件都很难用。
Google Earth Engine
自从用了这个，我就再也不用常规方法了，太方便了！强烈推荐！更多rgee文章可以从公众号回复rgee，查看rgee专题。
library(rgee)
library(sf)
ee_Initialize(drive = T)
# Define a region of interest with sf
ee_roi <- read_sf("./SHP/ChinaSizhi.SHP") %>%
sf_as_ee()
cmr_lucc = ee$ImageCollection("MODIS/006/MCD12Q1")$select('LC_Type1')$
filterBounds(ee_roi)
cmr_lucc1 = cmr_lucc$mosaic()$clip(ee_roi)
igbpLandCoverVis = list(
min = 1.0,
max = 17.0,
palette= c('05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
'69fff8', 'f9ffa4', '1c0dff')
)
Map$addLayer(cmr_lucc1, igbpLandCoverVis, 'cmrlucc')+
Map$addLayer(ee_roi, visParams = list(palette = "red"), "roi")
Map$centerObject(ee_roi)
igbp_lucc <- ee_as_raster(
image = cmr_lucc1,
region = ee_roi$geometry(),
scale = 500,
dsn = './LUCC/lucc2015_500.tif',
via = "drive"
)
#批量下载
for (j in 2001:2019) {
date_start = ee$Date$fromYMD(j, 1, 1)
date_end = ee$Date$fromYMD(j,12,31)
cmrlucc_composite = cmr_lucc$
filter(ee$Filter$date(date_start, date_end))$
mosaic()$clip(ee_roi)
tiffnames= as.character.Date(eedate_to_rdate(date_start))
cdr_ndvi_cmr <- ee_as_raster(
image = cmrlucc_composite,
region = ee_roi$geometry(),
scale = 500,
dsn = paste0('./CNlucc/', tiffnames),
via = "drive"
)
}
数据批量裁剪
由于中国范围矢量节点过多，GEE不能直接上传裁剪，需要下载下来数据后在本地进行裁剪。使用R对数据进行循环批量裁剪。
library(terra)
lucclist = list.files(path = "./CNlucc/", pattern = ".tif$")
luccdir = paste0("./CNlucc/", lucclist)
China = vect("D:/R/Yunnan/SHP/YRB.shp")
t1=proc.time()
for(i in 1:length(luccdir)){
lucc = rast(luccdir[i])
lucc_cn = trim(mask(lucc, China))
writeRaster(lucc_cn, filename = paste0("./lucc_YRB/", lucclist[i]))
print(paste("finish", i))
}
t2=proc.time()
t=t2-t1
print(paste0('执行时间：',t[3][[1]],'秒'))
栅格重分类
IGBP分类过于复杂，不便使用，为了方便研究，我对IGBP分类根据我的需要进行了合并。
重分类森林、草地、农田、其他
1森林：1、2、3、4、5 2草地：9、10 3农田：12、14 4其他
栅格批量重分类
重分类我使用了raster包，terra包重分类函数有BUG，计算结果和重分类矩阵指定的值不一致。同样还是使用了循环并行计算，同时对多年的土地覆被数据进行计算，提高计算速度。
library(doParallel)
lclist = list.files(path = "./YearlyLUCC/", pattern = ".tif$")
lcdir = paste0("./YearlyLUCC/", lclist)
#重分类矩阵，默认左开右闭
m = c(0, 5, 1,
5,8,4,
8, 10, 2,
10,11,4,
11,12,3,
12,13,4,
13,14,3,
14,18,4)
rclmat = matrix(m, ncol = 3, byrow = T)
cl = makeCluster(4) #指定核心数
registerDoParallel(cl)
foreach(i=1:length(lclist), .packages = 'raster')%dopar%{
lucc = raster(lcdir[i])
lucc_rcl = reclassify(lucc, rclmat) #默认左开右闭
writeRaster(lucc_rcl, filename = paste0("./RclYearLUCC/", lclist[i]))
}
stopCluster(cl)
回复“IGBP中国”有套路获取2001-2019年逐年中国IGBP数据和重分类后数据，以及QGIS制图工程文件。
lucc_cn为2001-2019逐年IGBP分类数据 lucc_cn_rcl为按上文重分类后的 IGBP中国.qgz为上图制图工程文件
点击阅读原文，查看更多土地覆被相关视频课程