R语言使用Delta降尺度算法对降雨数据降尺度

为了防止走失，请大家记得星标公众号

前面给大家介绍过彭守璋老师的这套数据：【数据分享】中国1km分辨率逐月降水、最低最高平均气温数据1901-2019，目前这套数据已经更新到2022年。在彭老师的论文中，介绍了这套数据的制作方法：ESSD|中国1km分辨率气温降水数据集（1901-2017年）

看着里面的算法还是比较简单的，就是计算量比较大，拿R语言复现一下这个制作过程，做一些自己所需的数据。

Delta降尺度降水数据

以京津冀区域为例进行试验，这个降尺度的计算量很大，全球的数据数据量非常巨大，尝试了一下全球的，说需要1600GB的存储空间，太占地方了，还是弄个小块的试试。

本复现研究的降尺度的数据时间范围为1982-2021年。

这里涉及到多年逐月降水量的计算，在公众号以前的推文中使用gimms包试验过，但是后来发现gimms包对于分辨率较高的数据有时候计算不出结果，而且gimms包是基于raster包制作的，raster包目前已经退休，因此本复现研究完全使用terra包完成，terra包是我目前使用的R包能运算数据量最大的，目前还没有遇到过软件性能不足的问题（只要电脑内存够用，我目前没发现terra包有啥算不动的）。

• gimms包计算多年平均月度气温降水见以前推文：
• 一文搞定多年平均月度气温、降水计算
• 本次使用了terra包tapp函数计算，用法和gimms包monthlyComposite函数类似

#使用DELTA降尺度算法
library(terra)
library(tidyverse)

#读取CRU降水数据
cru_pre = rast("G:/Geodata/CRUTS4.06/cru_ts4.06.1901.2021.pre.dat.nc", "pre")    #读取降水数据1901-2021年
plot(cru_pre)
summary(cru_pre)
cru_ts = seq(as.Date("1901-01-01"), as.Date("2021-12-31"), "month")      #生成1901-2021时间序列
cru_tsdf = as.data.frame(cru_ts)
cru_pre19822021 = cru_pre[[973:1452]]   #提取1982-2021年降水数据
plot(cru_pre19822021)

#获取1970-2000年和WorldClim对应的多年气象数据
cru_pre19702000 = cru_pre[[829:1200]]
#计算多年逐月降水CRUCLIM_PRE
yearIndices = rep(1:12, 31)
cru_pre_WC = tapp(cru_pre19702000, index=yearIndices, fun=mean, cores=20)

#计算1982-2021年CRUCLIM_PRE
cru_pre_wc40 = rep(cru_pre_WC,40)
plot(cru_pre_wc40)

#各个年月的降水异常值计算
An_PRE = cru_pre19822021/cru_pre_wc40
#异常值重采样，要和WorldClim一致

#读取WORLDCLIM降水数据
wcpreclist = list.files(path = "G:/Geodata/WorldCLIMB/wc2.1_30s_prec", pattern = ".tif$")
wcprecdir = paste0("G:/Geodata/WorldCLIMB/wc2.1_30s_prec/", wcpreclist)
wcprec = rast(wcprecdir)

#读取京津冀矢量数据
JJJSizhi = vect("./SHP/JJJ.shp")
plot(JJJSizhi)
#获取京津冀范围WorldClim降水数据
wcprec_JJJ = trim(mask(wcprec, JJJSizhi))
#获取京津冀范围异常值
An_PRE_JJJ = trim(mask(An_PRE, JJJSizhi))

#重采样后的逐年月降水异常值
An_PRE_resample = resample(An_PRE_JJJ, wcprec_JJJ)  #默认双线性重采样方法

#重采样后的异常值乘以WorldClim降水数据，获取降尺度结果
#worldclim只有12个月份，相乘之前需要将波段复制成和异常值一致的波段数，对应起来之后才能相乘
wcprec_JJJ480 = rep(wcprec_JJJ, 40)
#计算降尺度结果
PREdelta19822021 = An_PRE_resample*wcprec_JJJ480
cru_ts = seq(as.Date("1982-01-01"), as.Date("2021-12-31"), "month")  
names(PREdelta19822021) = cru_ts
plot(PREdelta19822021)
writeRaster(PREdelta19822021, filename = paste0("./DeltaClim/pre/", PREdelta19822021@ptr[["names"]], ".tif"))

R语言计算结果和彭老师数据对比验证

在这我主要是进行了2021年数据的相关性检验。

逐像元相关性

#对比降尺度计算结果
library(terra)
library(tidyverse)

#读取彭老师计算结果
pre_2021 = rast("./pre_2021.nc")
plot(pre_2021)

#读取矢量数据
JJJSizhi = vect("./SHP/JJJ.shp")
plot(JJJSizhi)

#彭老师裁剪至京津冀
pre_2021_JJJ = trim(mask(pre_2021, JJJSizhi))
plot(pre_2021_JJJ)

# 彭老师降水为0.1mm，单位换算为mm需要除以10
pre_2021_JJJ = pre_2021_JJJ/10
plot(pre_2021_JJJ)

#读取降尺度后的数据
Delta_PRElist = list.files(path = "./DeltaClim2/pre/", pattern = ".tif$")
Delta_PREdir = paste0("./DeltaClim2/pre/", Delta_PRElist)
Delta_PRE = rast(Delta_PREdir)
Delta_PRE2021 = Delta_PRE[[469:480]]
plot(Delta_PRE2021)

#将两套数据对齐裁剪
pre_2021_JJJ = extend(pre_2021_JJJ, Delta_PRE2021)
pre_2021_JJJ = crop(pre_2021_JJJ, Delta_PRE2021)
pre_2021_JJJ = resample(pre_2021_JJJ, Delta_PRE2021)

#自做降尺度结果和彭老师结果相关
z = c(Delta_PRE2021, pre_2021_JJJ)
#逐像元计算r和P
fun_cor =  function(x) {
  Rs = Hmisc::rcorr(x[1:12], x[13:24], type = "spearman")
  Rx = Rs$r[2]
  Px = Rs$P[2]
  return(c(Rx, Px))
}

#计算自做降尺度结果和彭老师结果的相关和P值
r_NDVI_PRE = app(z, fun_cor, cores=20)
names(r_NDVI_PRE) = c("r", "P-value")
plot(r_NDVI_PRE)
summary(r_NDVI_PRE)
writeRaster(r_NDVI_PRE, filename = "r_Me_Peng_PRE.tif", names=r_NDVI_PRE@ptr[["names"]])

散点图相关性

• summary()查看数据的值域，用于确定散点图坐标轴范围
• 本例中，降水数据最大值339,因此在这坐标轴设置区间0~360，间隔40

# 制作散点图
Vdelta = as.vector(values(Delta_PRE2021))
Vpeng = as.vector(values(pre_2021_JJJ))
V = as.data.frame(cbind(Vdelta,Vpeng))
summary(V)
#制作散点图
library(tidyverse)
library(ggpmisc)
NDVI.formula <- y~x
ggplot(V, aes(x=Vdelta, y=Vpeng))+
  geom_hex(bins = 100) +
  scale_fill_continuous(type = "viridis") +
  geom_smooth(method = "lm", formula = NDVI.formula) +
  stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), stat(p.value.label), sep = "*\", \"*")),
               formula = NDVI.formula, parse = TRUE)+
  geom_abline(slope = 1, intercept = 0, linetype="dotted", size=1.2)+
  coord_equal(ratio = 1)+
  scale_x_continuous(breaks = seq(0,360,40),limits = c(0,360))+  #x轴设置区间0~360，间隔40
  scale_y_continuous(breaks = seq(0,360,40),limits = c(0,360))+
  labs(x="Me", y= "Peng")+    #X,Y轴标题
  theme_bw()+
  theme(axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),  #y轴标题字号
        axis.text.x = element_text(size = 10),   #坐标轴刻度字号
        axis.text.y = element_text(size = 10),
        legend.position = "none")+  #隐藏图例
  theme(text=element_text(size=16,  family="serif"))

更多遥感数据处理与R语言数据分析课程请加入地理数据科学培训班：

点击阅读原文查看视频课程讲解