問題引入
2000-2010年逐年全國降水柵格數(shù)據(jù)浴鸿，共11張柵格圖像擎场。對每個柵格點進(jìn)行線性回歸，可以得到在該點處的回歸系數(shù)，最終的得到全國11年降水線性回歸系數(shù)分布圖厌蔽。
R包加載
raster、rgdal包用于柵格數(shù)據(jù)的輸入輸出夺姑，broom包可以方便的得到線性回歸的各種參數(shù)
library(raster)
library(broom)
library(rgdal)
數(shù)據(jù)導(dǎo)入
setwd("D:/1data/rain") #數(shù)據(jù)所在文件夾
raster1<-stack(list.files(pattern='*.tif$')) #堆棧所有柵格數(shù)據(jù)
time<-1:nlayers(raster1)
編寫函數(shù)
線性回歸系數(shù)可以直接提取忍坷，而顯著性p和r^2借助broom包的glance（）函數(shù)才能提取，缺點是效率很低深胳，運行時間可能是直接提取的5~10倍绰疤。
#提取線性回歸系數(shù)
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（）函數(shù)可以對柵格每個像素執(zhí)行上述函數(shù)非常方便
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)
小結(jié)
目前r語言處理柵格數(shù)據(jù)的國內(nèi)攻略還比較少，希望能給大家?guī)韼椭?/p>
