R语言空间插值的几种方法及案例应用

作者简介

勾蒙蒙，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)