一文搞定多年平均月度气温、降水计算

在前面的推文可以白嫖的数据资源！土地覆被、气象、GDP、人口、土壤等等各种空间数据中，给大家介绍过CRU气候数据，这个数据是逐月的，如果我们想知道近60年的1-12月月均温是多少，如何计算呢？其它参数比如降水，也可以参照这个方法进行计算。

已经计算好的数据在这，需要的自取：

数据准备

下载CRU平均气温逐月数据。

• 链接：https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.04/
• tmp平均气温
• 0.5度分辨率

在这里我下载了1901-2019年的nc数据

平均气温1901-2019

计算思路

• 读取nc
• terra包，rast
• 研究区范围裁剪
• 矢量范围读取vect
• 裁剪mask
• 掩膜外的值都设为NA并调整数据范围trim
• 提取所需年份对应的波段subset
• 导出栅格TIFF

以上是使用的terra包，terra包支持并行运算，效率高，下面的多年数据合成需要使用gimms包，是基于raster包写的，因此需要再使用raster包把前面裁剪好的栅格读取为RasterStack，然后使用gimms包进行多年平均计算

• 读取TIFF作为RasterStack
• 生成切片
• 月度合成
• 导出数据
• cellStats像元统计
• ggplot2绘制多年平均月度气温变化曲线

代码实现

数据裁剪转换TIFF部分

library(terra)
library(tidyverse)
library(raster)
library(gimms)
#读取温度
tmp19012019 = rast("cru_ts4.04.1901.2019.tmp.dat.nc", "tmp")
#中国矢量
cmr = vect("D:/R/CMR/SHP/CHina.shp")
#提取19600101-20191201逐月温度
tmpcmrts = terra::subset(tmp19012019, 709:1428)
# 近60年时间序列
temp <- seq.Date(from = as.Date("1960/01/01",format = "%Y/%m/%d"), by = "month", length.out = 720)
#数据裁剪范围外部分去除，设为NA
TempCMR = trim(mask(tmpcmrts, cmr))
#给每个图层赋予月份时间作为名称
names(TempCMR) = temp
#将数据导出为TIFF
writeRaster(TempCMR, filename ="./TMPCMR/TempChina.tiff", names=TempCMR@ptr[["names"]])

处理完成后，我们可以使用plot进行一下预览

plot(TempCMR)

plot(TempCMR[["2019-12-01"]])

多年平均计算

#使用GIMMS包进行逐月多年平均计算
#读取前面导出的TIFF
tempstack = stack("./TMPCMR/TempChina.tiff")
#生成切片
yearIndices = rep(1:12, 60)
ymc = monthlyComposite(tempstack, fun=mean, indices=yearIndices, cores=4)
writeRaster(ymc, filename = "./TMPCMR/TempChinaMon.tiff", bylayer=F, format="GTiff" )

近60年中国各月平均气温

生成曲线图

#月均温变化
#月均温计算
TempMon = cellStats(ymc, mean)
Mon = rep(1:12)
TempSeries = as.data.frame(cbind(TempMon, Mon))
#制图
p1 <- ggplot(data = TempSeries, aes(Mon, TempMon))+
  labs(x="Month", y="Temperature (℃)")+
  theme_bw()+
  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轴标签大小
  scale_x_continuous(breaks = c(seq(1,12, by =1)))+
  geom_point()+
  geom_line()
p1

逐月平均气温变化曲线

参考文献

点击阅读原文加入培训班查看视频课程讲解