查看原文
其他

R语言栅格时间序列回归分析——整体和逐像元计算,并行计算,快到飞起!

走天涯徐小洋 走天涯徐小洋地理数据科学 2022-07-17

R语言栅格时间序列回归分析——整体和逐像元计算,并行计算,快到飞起!

前面给大家分享了2001-2020年1km分辨率中国降水量数据,有同学问如何做回归分析,接下来以一元线性回归举例,从整体上和逐像元两种方法来计算。

使用的R语言程序包

  • terra包,栅格计算,支持并行计算
  • tidyverse包,R语言数据处理可视化必备,内含ggplot2
  • ggpmisc包,作为ggplot2的补充,时间序列图中标注回归方程用

整体回归

为了减小点计算量,我把数据进行了裁剪,保留京津冀地区。降水量为年降水量,单位为0.1mm。

京津冀降水分布,0.1mm

对整个区域的计算和时间序列分析较为简单,只需对每个栅格进行统计,计算整体均值,生成时间序列,使用ggplot2包进行制图,还使用了ggpmisc包进行了统计结果的标注。

library(terra)
library(tidyverse)
library(ggpmisc)
pre_cn = rast("D:/R/CnNDVI/Pre/YearPre2001_2020.tif")
jjj = vect("./SHP/JJJ.shp")
pre_jjj = trim(mask(pre_cn, jjj))
writeRaster(pre_jjj, "./Pre/YearPreJJJ2001_2020.tif")

#整体的线性回归
pre_ts = global(pre_jjj, "mean", na.rm=T, cores=4) #像元统计
year = seq(2001, 2020, 1)                  
pre_ts = cbind(year, pre_ts)                        #生成时间序列
Pre.formula <- y~x
p1 <- ggplot(data = pre_ts, aes(year, mean))+
  labs(x="Year", y="Pre")+
  theme_bw()+
  theme(title = element_text(size = 20),
        axis.title = element_text(size = 16),    #调整标题大小
        axis.text.x = element_text(size = 14),    #x轴标签大小
        axis.text.y = element_text(size = 14), 
        legend.position = "bottom",
        legend.text = element_text(size = 14),
        legend.title = element_blank())+     #y轴标签大小
  geom_point(size=1)+
  geom_line(size=0.5)+
  geom_smooth(method = "lm", formula = Pre.formula, size=0.5)+
  stat_poly_eq(aes(label =  paste(stat(eq.label), stat(rr.label), stat(p.value.label),sep = "*\", \"*")),
               formula = Pre.formula, parse = TRUE, 
               size=8)  
p1
时间序列图

上面是线性回归时间序列制图的结果,下面是R语言回归分析结果:

lmdt = lm(mean~year, data = pre_ts)   #对降水和年份进行一元线性回归
summary(lmdt)  #查看回归结果
线性回归结果

解释分析

  • Intercept:截距,即当X=0的预测值
  • 回归系数:回归线的斜率,52.60是计算的斜率结果,表明京津冀地区2001-2020年降水量呈增长态势,平均每年增长5.26mm。
  • Residuals:残差,观测值和拟合值之间的差异
  • p-value:p值,在这是0.03041 <0.05,我们可以认为通过了假设检验,也就是说,2001-2020年京津冀地区降水增长是可信的。
  • Multiple R-squared:R方,决定系数,可以被模型解释的变异的比例,值介于0到1之间,在这里为0.2347,R方有点小,解释性略差。

逐像元回归计算

逐像元回归计算,还是使用的terra包,重点是想办法把回归的结果提取出来,这个需要对回归结果比较了解。代码如下:

fun_linear = function(x){
  if(length(na.omit(x))<20) return(c(NA, NA, NA))  #删除数据不连续含有NA的像元
  year = seq(2001, 2020, 1) #自变量
  lmdt = lm(x~year)     #回归
  a1 = summary(lmdt)
  slope = a1$coefficients[2]   #斜率
  rsquared = a1$r.squared      #R2
  pvalue = a1$coefficients[8]  #p值
  return(c(slope, rsquared, pvalue))
}

pre_pixellinear = app(pre_jjj, fun_linear, cores=4)
names(pre_pixellinear) = c("slope""r2""p-value")
plot(pre_pixellinear)
writeRaster(pre_pixellinear, "./Pre/pre_pixellinear.tif")
逐像元回归结果,斜率,R方,P值

答疑部分

这种情况下可以把数据转为数据框,对这两幅影像进行回归计算。

  • values函数可以将栅格数据转为vector向量,然后转数据框
  • na.omit删除空值
  • 对生成的数据框的两列进行回归计算,代码如下:
pre_ts2 = pre_jjj[[1:2]]
Vpre_ts2 = as.data.frame(values(pre_ts2))
Vpre_ts2 = na.omit(Vpre_ts2)
names(Vpre_ts2) = c("x1""x2")
lm2 = lm(x1~x2, Vpre_ts2)
summary(lm2)
回归计算结果

参考文献

  1. 彼得·布鲁斯, 安德鲁·布鲁斯. 面向数据科学家的实用统计学[M]. 盖磊, 译. 人民邮电出版社, 2018.
  2. 地理学数学方法SPSS与R语言应用
  3. ggplot2绘制时间序列变化图
  4. R语言并行计算Sen+MK趋势分析

推荐这本书:

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

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存