其他
R语言月均温栅格转年均温栅格
R语言月均温转年均温
前面给大家推荐过彭守璋老师发布的1km分辨率全国的月度气候数据。本套数据是逐月的，对于年际变化研究来说使用不太方便，需要把月度数据合称为年度数据，在这举个栗子，教大家如何进行年度数据的合成。
月均温转年均温
本套数据多数为3年一个nc文件，我采用的处理思路如下：
使用 terra包读取nc文件提取每年12个月的数据 对12个月的数据进行合成，月均温转年均温采用求均值的方法，栅格计算 app函数对每年的数据进行波段合成 c导出处理好的数据 writeRaster
数据处理中遇到的问题
实际处理过程中发现，由于2018年及以后的数据和2018年以前的数据范围不一致，导致波段合成失败。对此，我使用了crop和resample两个函数对栅格范围进行统一。以下是合成2001-2020年年度均值数据的代码：
library(terra)
# 2000-2002
tmp_00_02 = rast("./tmp/tmp_2000_2002.nc")
temp_00_02 = seq(as.Date("2000-01-01"), as.Date("2002-12-01"), "month")
names(tmp_00_02) = temp_00_02
tmp_03_05 = rast("./tmp/tmp_2003_2005.nc")
temp_03_05 = seq(as.Date("2003-01-01"), as.Date("2005-12-01"), "month")
names(tmp_03_05) = temp_03_05
tmp_06_08 = rast("./tmp/tmp_2006_2008.nc")
temp_06_08 = seq(as.Date("2006-01-01"), as.Date("2008-12-01"), "month")
names(tmp_06_08) = temp_06_08
tmp_09_11 = rast("./tmp/tmp_2009_2011.nc")
temp_09_11 = seq(as.Date("2009-01-01"), as.Date("2011-12-01"), "month")
names(tmp_09_11) = temp_09_11
tmp_12_14 = rast("./tmp/tmp_2012_2014.nc")
temp_12_14 = seq(as.Date("2012-01-01"), as.Date("2014-12-01"), "month")
names(tmp_12_14) = temp_12_14
tmp_15_17 = rast("./tmp/tmp_2015_2017.nc")
temp_15_17 = seq(as.Date("2015-01-01"), as.Date("2017-12-01"), "month")
names(tmp_15_17) = temp_15_17
#计算2001-2002年气温
tmp01 = tmp_00_02[[13:24]]
plot(tmp01)
tmpmean01 = app(tmp01, "mean", cores=4)
tmp02 = tmp_00_02[[25:36]]
plot(tmp02)
tmpmean02 = app(tmp02, "mean", cores=4)
#计算2003-2005年气温
tmp03 = tmp_03_05[[1:12]]
tmp04 = tmp_03_05[[13:24]]
tmp05 = tmp_03_05[[25:36]]
tmpmean03 = app(tmp03, "mean", cores=4)
tmpmean04 = app(tmp04, "mean", cores=4)
tmpmean05 = app(tmp05, "mean", cores=4)
#计算2006-2008年气温
tmp06 = tmp_06_08[[1:12]]
tmp07 = tmp_06_08[[13:24]]
tmp08 = tmp_06_08[[25:36]]
tmpmean06 = app(tmp06, "mean", cores=4)
tmpmean07 = app(tmp07, "mean", cores=4)
tmpmean08 = app(tmp08, "mean", cores=4)
#计算2009-2011年气温
tmp09 = tmp_09_11[[1:12]]
tmp10 = tmp_09_11[[13:24]]
tmp11 = tmp_09_11[[25:36]]
tmpmean09 = app(tmp09, "mean", cores=4)
tmpmean10 = app(tmp10, "mean", cores=4)
tmpmean11 = app(tmp11, "mean", cores=4)
#计算2012-2014年气温
tmp12 = tmp_12_14[[1:12]]
tmp13 = tmp_12_14[[13:24]]
tmp14 = tmp_12_14[[25:36]]
tmpmean12 = app(tmp12, "mean", cores=4)
tmpmean13 = app(tmp13, "mean", cores=4)
tmpmean14 = app(tmp14, "mean", cores=4)
#计算2015-2017年气温
tmp15 = tmp_15_17[[1:12]]
tmp16 = tmp_15_17[[13:24]]
tmp17 = tmp_15_17[[25:36]]
tmpmean15 = app(tmp15, "mean", cores=4)
tmpmean16 = app(tmp16, "mean", cores=4)
tmpmean17 = app(tmp17, "mean", cores=4)
#计算2018年气温
tmp18 = rast("./tmp/tmp_2018.nc")
tmpmean18 = app(tmp18, "mean", cores=4)
#计算2019年气温
tmp19 = rast("./tmp/tmp_2019.nc")
tmpmean19 = app(tmp19, "mean", cores=4)
#计算2020年气温
tmp20 = rast("./tmp/tmp_2020.nc")
tmpmean20 = app(tmp20, "mean", cores=4)
tmpmean01_10 = c(tmpmean01, tmpmean02, tmpmean03, tmpmean04, tmpmean05, tmpmean06, tmpmean07, tmpmean08, tmpmean09, tmpmean10)
names(tmpmean01_10) = seq(2001,2010,1)
writeRaster(tmpmean01_10, "./tmpmean/Yeartmp2001_2010.tif")
tmpmean11_17 = c(tmpmean11, tmpmean12, tmpmean13, tmpmean14, tmpmean15, tmpmean16, tmpmean17)
names(tmpmean11_17) = seq(2011,2017,1)
writeRaster(tmpmean11_17, "./tmpmean/Yeartmp2011_2017.tif")
china = vect("./SHP/China.shp")
tmpmean18_20 = c(tmpmean18, tmpmean19, tmpmean20)
tmpmean18_20cn = trim(mask(tmpmean18_20, china))
writeRaster(tmpmean18_20cn, "./tmpmean/Yeartmp2018_2020.tif")
#合成2001-2020数据
tmpmean11_17 = rast("./tmpmean/Yeartmp2011_2017.tif")
tmpmean18_20 = rast("./tmpmean/Yeartmp2018_2020.tif")
tmpmean01_10 = rast("./tmpmean/Yeartmp2001_2010.tif")
tmpmean11_17 = crop(tmpmean11_17, tmpmean18_20)
tmpmean01_10 = crop(tmpmean01_10, tmpmean18_20)
tmpmean01_10 = resample(tmpmean01_10, tmpmean18_20)
tmpmean11_17 = resample(tmpmean11_17, tmpmean18_20)
tmpmean01_20 = c(tmpmean01_10, tmpmean11_17, tmpmean18_20)
temp = seq(2001,2020,1)
names(tmpmean01_20) = temp
writeRaster(tmpmean01_20, "./tmpmean/Yeartmp2001_2020.tif")
合成结果检查
合成之后，检查一下计算结果，我计算了一个整体的时间序列，生成了每年的合成数据的预览图。R语言代码如下：
#查看合成后的数据结果
library(tidyverse)
tmpseries = global(tmpmean01_20, "mean", cores=6, na.rm=T)
tmpseries = rownames_to_column(tmpseries, "temp")
tmpseries$temp = as.numeric(tmpseries$temp)
p1 = ggplot(data = tmpseries, aes(temp, mean))+
labs(x= "Year", y= "TMP")+
geom_point()+
geom_line()+
theme(title=element_text(size=22),
axis.title=element_text(size=20),#调整标题大小
axis.text.x=element_text(size=18),#x轴标签大小
axis.text.y=element_text(size=18),
legend.position="bottom",
legend.text=element_text(size=18),
legend.title =element_blank())+#y轴标签大小
theme_bw()
p1
参考文献
点击阅读原文，查看更多遥感数据处理视频课程