R語言空間插值的幾種方法及案例應用

R語言空間插值的幾種方法及案例應用

來自專欄 R語言中文社區6 人贊了文章

作者簡介

勾蒙蒙，R語言資深愛好者。

個人公眾號： R語言及生態系統服務。

##載入程序包

library(raster)

library(sp)

library(rgdal)

library(gstat)

library(raster)

library(maptools)

##設置工作空間

setwd("C:/Users/lx/Desktop/sun")

數據為環京津冀地區153個站點2002年7月降雨數據

##讀取數據

Data<-na.omit(read.csv("Data.csv",header=T))

head(Data)

Month Plot Year rain X Y

1 7 54365 2002 139.8 125.3500 41.28333

2 7 54259 2002 197.5 124.9167 42.01000

3 7 54493 2002 317.6 124.7833 40.71667

4 7 54497 2002 179.4 124.3333 40.05000

5 7 54351 2002 159.8 124.0833 41.91667

6 7 54254 2002 88.90 124.0500 42.53333

製圖顯示數據

plot(sort(Data$rain), ylab=" 年降雨量(mm)", las=1, xlab=站點)

##導入京津冀地區矢量地圖

bound<-readOGR("bound.shp")

plot(bound,col="grey")

##設定降雨數據的投影為WGS84

dsp <- SpatialPoints(Data[,5:6], proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

dsp <- SpatialPointsDataFrame(dsp,Data)

展示一個直觀的地圖

cuts<-c(0,30,50,100,150,300)#設置間距

blues <- colorRampPalette(c("red","blue"))(5)#設置顏色梯度，即漸變色。c("red","blue")(5)代表顏色從黃色漸變到橘色，再漸變到藍色，再到深藍色,5則代表長度為5。例子：plot(20:1, bg = blues[rank(5:1)], cex = 2, pch = 22)

pols <- list("sp.polygons",bound, fill = "lightgray")#構建京津冀的SpatialPolygons對象

spplot(dsp,"rain", cuts=cuts, col.regions=blues, sp.layout=pols, pch=20, cex=2)

將經緯度轉成平面坐標，使插值結果與數據保持一致，這裡用到的坐標系也是WGS84，+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

WGS84<- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")#設置參考系WGS84

dsp1<-spTransform(dsp,WGS84)#將經緯度轉成平面坐標，使用WGS參考系

bound1<-spTransform(bound,WGS84)

使用鄰域多邊形插值

library(dismo)

v <- voronoi(dsp1)

plot(v)

這個圖看起來怪怪的，接下來，我們將其範圍限定在京津冀地區

bound2<-aggregate(bound1)#聚合，降低解析度

v1<-intersect(v,bound2)#將兩個圖層相交

spplot(v1, "rain", col.regions=rev(get_col_regions()))#繪製多邊形圖

這個圖看起來則要舒適的多，但其輸出的結果是多邊形，接下來用」rasterize」將結果柵格化，為詳細介紹shape格式轉柵格的過程，這裡增加一個小專題：Shape轉柵格

在shape轉柵格之前，首先需要建議一個新的空白的柵格，並指定控制柵格解析度的行列，用extent制定空間範圍

blank_raster<-raster(nrow=100,ncol=100,extent(bound))

接下來給柵格賦值

values(blank_raster)<-1

plot(blank_raster)

因為給柵格的賦值都為1，因此上圖顯示的也只有一個顏色，你也可以給柵格賦值多個，即nrow*ncol

values(blank_raster)<-1:(100*100)

plot(blank_raster,col=rainbow(100))

下圖則看起來要好的多了，下面將用到rasterize進行正式轉換了，這裡需要注意的是柵格是有解析度的，因此上面設置的100*100的解析度有可能不一定合適，我們用循環看一下

layout(matrix(1:4, ncol=2, byrow=TRUE))

res<-c(20,100,500,1000)

for(r in res){

blank_raster<-raster(nrow=r,ncol=r,extent(bound))

values(blank_raster)<-1

bound_raster<-rasterize(bound,blank_raster)

bound_raster[!(is.na(bound_raster))] <- 1

plot(bound_raster,main=paste("Res: ",r,"*",r))

plot(bound,add=T)

}

毫無疑問，1000*1000的解析度吻合效果要更好，但也不一定越高越好，會影響處理速度，這裡我們選擇1000*1000的解析度

言歸正傳，將結果柵格化：

vr <- rasterize(v1,bound_raster,"rain")

plot(vr)

使用最近鄰點插值

gs<-gstat(formula=rain~1,location=dsp1,nmax=5,set=list(idp=0))

nn<-interpolate(bound_raster,gs)

nnmask<-mask(nn,vr)##掩膜提取

plot(nnmask)

反距離加權插值

gs <- gstat(formula=rain~1, locations=dsp1)

idw <- interpolate(bound_raster, gs)

idwmask<-mask(idw,vr)

plot(idwmask)

普通克里金插值

克里金法是通過一組具有 z 值的分散點生成估計表面的高級地統計過程。與插值工具集中的其他插值方法不同，選擇用於生成輸出表面的最佳估算方法之前，有效使用克里金法工具涉及 z 值表示的現象的空間行為的交互研究。

求變異函數,首先繪製半變異圖

v<- variogram(log(rain) ~ 1, data =dsp)

plot(v,plot.number=T)

根據半變異圖可知，已知點的自相關關係semivariance隨著距離distance的增加而增加，通過其分布，可初步確定用線性函數來擬合：確定哪些函數可以用show.vgms()實現。

v.fit<-fit.variogram(v,model=vgm(1,"Lin",0))

plot(v,v.fit)

Grid<-as(bound_raster,"SpatialGridDataFrame")#首先現將邊界柵格轉成空間網格

kri<-krige(formula=rain~1,model=v.fit,locations=dsp,newdata=Grid,nmax=12, nmin=10)#location為已知點的坐標；newdata為需要插值的點的位置；nmax和nmin分別代表最多和最少搜索點的個數

spplot(kri["var1.pred"])

泛克里金插值

用泛克里金法需謹慎，因其假定數據中存在覆蓋趨勢，應該僅在了解數據中存在某種趨勢並能夠提供科學判斷描述泛克里金法時，才可使用該方法

v<- variogram(log(rain) ~ X+Y, data =dsp)

plot(v,plot.number=T)

v.fit<-fit.variogram(v,model=vgm(1,"Lin",0))

plot(v,v.fit)

uk <- krige(ANN_PREC ~ X + Y, dsp, Grid, v.fit)

Grid<-as(bound_raster,"SpatialGridDataFrame")

kri<-krige(formula=rain~1,model=v.fit,locations=dsp,newdata=Grid,nmax=12, nmin=10)

spplot(kri["var1.pred"])

薄盤樣條函數

install.packages("fields")

library(fields)

m <- Tps(coordinates(dsp), dsp$rain)

tps <- interpolate(bound_raster, m)

tps <- mask(tps, bound_raster)

plot(tps)

推薦閱讀：