快手10000赞只需0.5元(快手1万赞多少钱)

R语言栅格时间序列回归分析——整体和逐像元计算，并行计算

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

使用的R语言程序包

terra 包，栅格计算，支持并行计算

tidyverse 包，R语言数据处理可视化必备，内含 ggplot2

ggpmisc 包，作为 ggplot2 的补充，时间序列图中标注回归方程用

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方有点小，解释性略差。

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

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

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

ggplot2绘制时间序列变化图

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

图文编辑：王梦娇

审编：黄莘绒

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

猜你喜欢

1. 地学招聘 | 临沂大学资源环境学院2021年人才招聘启事

2. 会议通知 | 亚洲地理学会第一届青年地理学家研讨会

3. 期刊目录 | 《中国生态旅游》第11卷第4期

4. 专题征稿 | 《生态学报》“青藏高原生物多样性保护与生态安全”专刊征稿通知

扫描二维码，关注我们

都看到这里了，点个【在 看】 再走呗~