问题引入

2000-2010年逐年全国降水栅格数据，共11张栅格图像。对每个栅格点进行线性回归，可以得到在该点处的回归系数，最终的得到全国11年降水线性回归系数分布图。

R包加载

raster、rgdal包用于栅格数据的输入输出，broom包可以方便的得到线性回归的各种参数

library(raster)
library(broom)
library(rgdal)

数据导入

setwd("D:/1data/rain")                                    #数据所在文件夹
raster1<-stack(list.files(pattern='*.tif$'))         #堆栈所有栅格数据
time<-1:nlayers(raster1)

编写函数

线性回归系数可以直接提取，而显著性p和r^2借助broom包的glance（）函数才能提取，缺点是效率很低，运行时间可能是直接提取的5~10倍。

#提取线性回归系数
fun1 <- function(x) { if (is.na(x[1])){ NA } else lm(x ~ time)$coefficients[2] }      
##提取线性回归的显著性
fun2 <- function(x) { if (is.na(x[1])){ NA } else glance(lm(x ~ time))$p.value }
##提取线性趋势r^2
fun3<-function(x) { if (is.na(x[1])){ NA } else glance(lm(x ~ time))$r.squared }

逐栅格计算

raster包中的calc（）函数可以对栅格每个像素执行上述函数非常方便

rain.b<-calc(raster1,fun1)
rain.p<-calc(raster1,fun2)
rain.r2<-calc(raster1,fun3)

栅格输出

writeRaster(rain.b,filename = file.path( "E:/1data/rain","rain.b.tif"),format="GTiff", overwrite=TRUE)
writeRaster(rain.p,filename = file.path( "E:/1data/rain","rain.p.tif"),format="GTiff", overwrite=TRUE)
writeRaster(rain.r2,filename= file.path("E:/1data/rain","rain.r2.tif"),format="GTiff", overwrite=TRUE)

小结

目前r语言处理栅格数据的国内攻略还比较少，希望能给大家带来帮助