其他
使用rgee进行Climate Data Record NDVI月度最大值合成
使用rgee进行Climate Data Record NDVI月度最大值合成
前面的推文给大家介绍了1982-2020年逐月NDVI数据产品的获取：
今天讲讲这个数据如何制作的。
本数据集基于NOAA CDR AVHRR NDVI V5数据，使用R语言rgee包调用Google Earth Engine服务处理，在本地使用R语言terra包进行波段融合、裁剪等处理得到。
rgee月度最大值合成
启动 rgee设置区域范围，由于中国国境矢量节点数过多，不能直接在GEE中处理，需要先外扩一个四至范围
#引用所需的R包，启动GEE服务
library(rgee)
library(sf)
ee_Initialize(drive = T)
数据预览
在GEE中获取数据之前，最好先进行一下预览，看看数据计算结果是不是想要的。
# 设置范围和数据集
ee_roi <- read_sf("./SHP/ChinaSizhi.SHP") %>% #这里需要用一个节点少的中国四至范围
sf_as_ee()
cdr_ndvi <- ee$ImageCollection("NOAA/CDR/AVHRR/NDVI/V5")$select("NDVI")
#数据预览
ndvi_composite = cdr_ndvi$
filter(ee$Filter$date('2020-07-01', '2020-07-31'))$
max()
ndviParams <- list(palette = c(
"#d73027", "#f46d43", "#fdae61",
"#fee08b", "#d9ef8b", "#a6d96a",
"#66bd63", "#1a9850"
))
Map$addLayer(ndvi_composite, ndviParams, "CDR_NDVI")
MVC合成输出
CDR NDVI数据覆盖1981-2021年，为保证整年，我选择了1982-2020年的数据。需要注意的是，下面使用了两个循环，分别获取1-11月还有12月的数据。
for (i in 1:11) {
for (j in 1982:2020) {
date_start = ee$Date$fromYMD(j, i, 1)
date_end = ee$Date$fromYMD(j,i+1,1)
ndvi_composite = cdr_ndvi$
filter(ee$Filter$date(date_start, date_end))$
max()
tiffnames= as.character.Date(eedate_to_rdate(date_start))
cdr_ndvi_cmr <- ee_as_raster(
image = ndvi_composite,
region = ee_roi$geometry(),
scale = 5000,
dsn = paste0('./CDR_NDVI_MVC/', tiffnames),
via = "drive"
)
}
}
for (j in 1982:2020) {
date_start = ee$Date$fromYMD(j, 12, 1)
date_end = ee$Date$fromYMD(j+1,1,1)
ndvi_composite = cdr_ndvi$
filter(ee$Filter$date(date_start, date_end))$
max()
tiffnames= as.character.Date(eedate_to_rdate(date_start))
cdr_ndvi_cmr <- ee_as_raster(
image = ndvi_composite,
region = ee_roi$geometry(),
scale = 5000,
dsn = paste0('./CDR_NDVI_MVC/', tiffnames),
via = "drive"
)
}
数据后处理
前面是在rgee中获取NDVI数据，数据下载完成后是单波段的逐月数据，由于使用的是四至范围，还需要进行一些裁剪，以及波段合成等操作，让数据成为多波段数据，并且只保留中国境内的值，范围外的值设为Nodata。
library(terra)
#读取CDR NDVI月最大值数据
cdrnames= list.files(path = 'D:/R/CMRrgee/CDR_NDVI_MVC')
cdr = paste0("D:/R/CMRrgee/CDR_NDVI_MVC/", cdrnames)
cdr_mon = rast(cdr)
#读取中国范围
maskdata <- vect("./SHP/China.shp")
cdr_mon <- trim(mask(cdr_mon, maskdata))
#提取时间信息修改图层名称
#CDR时间
CDR_temp = seq(as.Date("1982-01-01"), as.Date("2020-12-01"), "month")
#CDR图层改名
names(cdr_mon) = CDR_temp
#值域修改波段融合数据输出
cdr_mon2 = cdr_mon/10000
#CDR裁剪
writeRaster(cdr_mon2, filename = "./CDRNDVI/CDRMonthlyMax.tif", names=cdr_mon2@ptr[["names"]])
参考文献
徐洋, 杨雅萍. 1982–2020年中国5 km分辨率逐月NDVI数据集[J/OL]. 中国科学数据, 2021. (2021-7-11). DOI: 10.11922/11-6035.csd.2021.0041.zh.
更多rgee相关文章请查看rgee专题