干货分享 | R语言栅格时间序列回归分析——整体和逐像元计算，并行计算，快到飞起！

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值 参考文献

1. 彼得·布鲁斯, 安德鲁·布鲁斯. 面向数据科学家的实用统计学[M]. 盖磊, 译. 人民邮电出版社, 2018.

2. 地理学数学方法SPSS与R语言应用

3. ggplot2绘制时间序列变化图

4. R语言并行计算Sen+MK趋势分析

图文编辑：王梦娇

审编：黄莘绒

终审： 颜子明 黄宗财 鲁嘉颐

最后编辑：2023-06-16

作者：萌小白

一个热爱网络的青年！

站内专栏 邮箱