问题引入
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语言处理栅格数据的国内攻略还比较少，希望能给大家带来帮助
