目录
- 0.问题导入
- 1.随机生成示例数据
- 2.生成插值目标精度的数据框（本实例以0.5°×0.5°为目标）
- 3.将样本点转化为坐标点并对降水NA值进行筛除
- 4.基于样本点进行克里金插值
- 5.插值结果可视化并出图（图1）
- 6.总结
- 7.本篇所使用的软件包（没有的需要通过install.packages('')进行安装）
- 8.致谢
- 9.写在最后
0. 问题导入
当我们想把点状数据扩展成面状数据并进行可视化的时候，第一步操作就是插值。通常，我们会在ArcGIS中完成插值操作，今天这篇给大家带来R版的的插值算法（算法已打包成函数，只需要输入同示例数据相同结构的data.frame，即可实现插值与出图，文末附获取方式～）。
1. 随机生成示例数据
1.1 根据extent 生成示例区域shp
如何根据extent批量生成shp请参考R-如何基于并行计算根据extent生成shp，并将一个list的shp文件塌缩为一个shp？
setwd('L:\\JianShu\\20200101')
ex = extent(105,112,42,47)
shp = as(ex,'SpatialPolygons')
crs(shp) = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
1.1 在示例区域边界范围内随机生成坐标点及对应位置的降水值
long = runif(20,105,112)
lat = runif(20,42,47)
pr = runif(20,0,15)
df = data.frame(long = long,lat = lat,pr = pr)
head(df)
long lat pr
1 107.0188 45.48182 7.6899736
2 111.6168 44.08171 3.1733193
3 111.7310 46.96119 3.1563925
4 107.1516 46.35613 0.1777086
5 105.0868 42.40558 2.7821381
6 108.8859 44.23798 6.0884563
2. 生成插值目标精度的数据框（本实例以0.5°×0.5°为目标）
long_shp = extent(shp)[1:2]
lat_shp = extent(shp)[3:4]
long_p = seq(min(long_shp),max(long_shp),by =resolution )
lat_p = seq(min(lat_shp),max(lat_shp),by = resolution)
s_p = expand.grid(long_p,lat_p)
inusa = map.where(x = s_p[,1],y = s_p[,2])
inusa = !is.na(inusa)
s_p = s_p[inusa,]
colnames(s_p) = c('long','lat')
coordinates(s_p) = ~ long + lat
gridded(s_p) = T
fullgrid(s_p) = F
proj4string(s_p) = crs(shp)
3. 将样本点转化为坐标点并对降水NA值进行筛除
本篇示例数据不含NA值，但如果你的数据中有含有NA值，这一步将避免因NA值而报错的出现。
pr = as.numeric(df[,3])
na_index = which(is.na(pr))
df = data.frame(long=df[,1],lat=df[,2],pr=as.numeric(df[,3]))
df2 = df
if(length(na_index)>0){
df = df[-na_index,]
}
coordinates(df) = ~ long+lat
proj4string(df) = crs(shp)
print("Pass Phase 2")
4. 基于样本点进行克里金插值
在这一步，我们将完成克里金插值，并将插值结果输出为名为raster.tif 的栅格文件，并储存在（setwd('L:\JianShu\20200101')）路径下。
f1 = as.formula(pr ~ long+lat)
var.sampl = variogram(f1,df,cloud =F,
cutoff = 100000,
width = 89900)
dat.fit = fit.variogram(var.sampl,fit.ranges = F,
fit.sills = F,
vgm(psill = 14,model = 'Sph',
range = 590000,
nugget = 0))
dat.krig = krige(f1,df,s_p,dat.fit)
dat.r = raster(dat.krig)
writeRaster(dat.r,'raster.tif',format = 'GTiff',overwrite = T)
5.插值结果可视化并出图（图1）
在这一步，我们将完成插值结果的可视化与出图。图件存储路径与raster.tif 存储路径一致，同步骤4.
krg.res = as.data.frame(dat.r,xy = T)
colnames(krg.res) = c('long','lat','pr')
mycolor = colorRampPalette(brewer.pal(11,'Spectral'))(12)
p = ggplot()+
geom_tile(data = krg.res,aes(x = long,y = lat,fill = pr),
color = 'transparent')+
geom_point(data = df2,aes(x = long,y = lat),color = 'black',
shape =1, size = 5)+
scale_fill_gradientn(colours = mycolor)+
theme_bw()+
theme(
axis.text = element_text(face = 'bold',color = 'black',size = 12),
axis.title = element_text(face = 'bold',color = 'black',size = 14, hjust = .5),
legend.text = element_text(face = 'bold',color = 'black',size = 12),
legend.title = element_blank(),
legend.position = 'bottom',
legend.direction = 'horizontal',
legend.key.width = unit(2,'cm')
)+
xlab('Longitude')+
ylab('Latitude')
png('krige.png',
height = 20,
width = 20,
units = 'cm',
res = 800)
print(p)
dev.off()
图1 克里金插值结果可视化
6. 总结
本篇主要解决了以下几个问题：
- 如何在R中根据点数据进行克里金插值？
- 如何将克里金插值结果在R中基于ggplot2进行可视化？
7. 本篇所使用的软件包（没有的需要通过install.packages('')进行安装）
library(rgdal)
library(tmap)
library(raster)
library(ggplot2)
library(maptools)
library(maps)
library(gstat)
library(RColorBrewer)
8. 致谢
首先，祝大家元旦快乐！感谢大家的持续关注，小编会继续努力，持续更新下去的！
大家如果觉得有用，还麻烦大家关注点赞，也可以扩散到朋友圈，多谢大家啦～
大家如果在使用本代码的过程有遇到问题的，可以留言评论，也可以私信我哈~~
写在最后
本篇已完成克里金插值R函数的打包，只需要输入
1. 同篇头df相同的数据；
2. 插值边界的shp文件；
3. 设置好插值精度（默认为0.5°）。
即可完成插值，并输出
1. 插值结果的Tiff文件；
2. 插值结果的可视化图件。
使用举例:
setwd('L:\\JianShu\\20200101')
kriging(df,shp,0.5)
是不是瞬间简单多了，需要函数的朋友需要转发朋友圈并集够20个赞哈，然后微信私信我即可获取克里金插值打包函数文件。
小编联系方式
网友评论