必知必会（二）（下）

四川大学在读研究生

接上文

7.2.2 泛克里金插值

用泛克里金法需谨慎，因其假定数据中存在覆盖趋势，

应该仅在了解数据中存在某种趋势并能够提供科学判断描述泛克里金法时，才可使用该方法

这在地质统计领域用得比较多，如矿藏的分布。

我这儿没有相关数据，仅仅用气温数据，可能不太准确，但是思路和流程是对的。

首先建立趋势模型，根据将大部分点，绘制样本实验变异函数图。

如下图，通过调节cutoff和width将大多数点显示在范围内，

可以使用plot.number = TRUE显示点的数量。

 1library(gstat) 

 2
 3### 添加变量X和Y
 4TEM_sp2  5TEM_sp2$X  6TEM_sp2$Y  7trend_1  8
 9TEM_v 10 cutoff = 30, # cutoff为对角线长度调整，其与width会相互影响
11 width = 2) # width表示相邻2个点之间的distance，width越小，点越多
12plot(TEM_v) 

然后根据点的走势，对照表，选择合适的拟合模型，不同的模型，vgm()内参数不一样。

这里我们选择spherical模型。 其3个参数分别表示如下意思：

1TEM_v_fit 2 fit.ranges = FALSE, fit.sills = FALSE,
3 model = vgm(psill = 18, model = "Sph", range = 28, nugget = 2.5))
4
5plot(TEM_v, TEM_v_fit) 
6

运算量非常大，经常溢出，这里不进行运算。

 1library(gstat)
 2library(raster)
 3
 4# 根据上面的拟合模型进行克里金插值的计算
 5TEM_krg  6 locations = TEM_sp2, # 数据点坐标
 7 newdata = grd_TEM, # 需要插值点的位置
 8 model = TEM_v_fit
 9 ) 
10
11TEM_raster 12TEM_mask 13
14rm(TEM_v, TEM_v_fit, TEM_krg, TEM_raster, TEM_sp2)

1library(tmap) 

2library(ggplot2)
3
4# tmap绘图
5tmap_TEM()
6
7## ggplot2绘图
8TEM_mask_df 9ggplot_TEM()

7.3

akima插值

akima插值不支持sp数据对象的插值，只支持dataframe和matrix对象插值。

插值结果为dataframe对象, 只能形成SpatialPixelsDataFrame栅格类型，

与前面的sp对象插值不同，sp对象插值结果为SpatialGridDataFrame栅格类型。

目前dataframe对象插值更简单，

但是SpatialPixelsDataFrame对象不支持多个多边形边界进行筛选。 over()不支持多个多边形边界进行筛选 ，mask()不支持SpatialPixelsDataFrame对象。

所以只能一个个边界进行筛选，然后索引内部元素进行合并。

7.3.1 插值运算

 1library(akima)
 2library(sp)
 3library(raster)
 4
 5# 对整个TEM_1th进行插值 

 6TEM_interp  7 xo = seq(60, 140, by = 0.1), # 指定插值经度范围
 8 yo = seq(10, 60, 0.1), # 指定插值纬度范围
 9 linear = FALSE, # 表示是线性插值还是样条插值
10 extrap = TRUE) # 表示是否接受外插，有的栅格只能外插才能得到，
11
12# 生成网格数据
13TEM_grd 14TEM_grd 15 length(TEM_interp$z))))
16# 分离地图边界
17df_as_sp 18 map_subset 19 Sr1 20 Srs1 21 SpP 22 partmap_sp 23 Sr = SpP,
24 data = data.frame(Names = "coords", row.names = row.names(SpP)))
25 return(partmap_sp)
26}
27
28Mainland_sp 29Hainan_sp 30
31# 筛选各个边界内的栅格数据
32Mainland_overcheck 33Hainan_overcheck 34Mainland_grd 35Hainan_grd 36
37# 栅格合并
38Mainland_grd 39Hainan_grd 40
41grd_bind_noTaiwan 42
43rm(TEM_interp, TEM_grd, df_as_sp,
44 Mainland_overcheck, Hainan_overcheck, 
45 Mainland_grd, Hainan_grd )# 移除中途数据
46
47# 生成等值线数据
48breaks_lines 

7.3.2 ggplot2绘图

 1library(ggplot2)
 2
 3ggplot(data = grd_bind_noTaiwan) + 
 4 # 所有栅格
 5 geom_raster(aes(x=x,y=y,fill=kde)) + 
 6 # 增加等值线
 7 geom_contour(aes(x=x,y=y,z=kde),
 8 color ="white",breaks=breaks_lines) + 
 9 # 中国省级地图轮廓 

10 geom_polygon(data = Chinaprovinces_df,
11 aes(x = long, y = lat, group = group), 
12 color = "black", fill = "transparent", size = 0.5) + 
13 # 台湾地图
14 geom_polygon(data = Taiwan_df,
15 aes(x = long, y = lat, group = group), 
16 color = NA, fill = "grey", size = 0.5) + 
17 # 各省名字，坐标为省会所在位置
18 geom_text(data = provinces,
19 aes(x = x, y = y, label = shortname), 
20 color = "green", size = 2) + 
21 # 增加南海九段线
22 geom_line(data = Nine_lines, 
23 aes(x=long, y=lat, group=ID), 
24 colour="black", size=1) + 
25 coord_cartesian() + # geom_raster只能搭配笛卡尔坐标系
26 # 栅格颜色填充标度
27 scale_fill_gradient2(low = "blue", mid = "white", midpoint = 0, 
28 high= "red",na.value = "grey", # 设定缺失值为灰色
29 name = "气温(℃)") + #
30 labs(title = "中国平均气温分布图\n(2017年1月1日)", 
31 caption = "注：图中数据不含台湾地区") + 
32 theme_void()+
33 theme(
34 legend.position=c(0.2,0.2),
35 legend.background=element_blank(),
36 plot.title = element_text(colour = "magenta", size = 13, 
37 face = "bold", hjust = 0.5),
38 legend.title = element_text(face = "bold", colour = "deeppink")
39 )

8、插入小地图

tmap中只能插入与leaflet中类似的小地图，即在线小地图。

实际上还是用PS截图拼图更加方便，甚至用PPT也行。

小地图会增加上百行代码。

参

考来源

gstat插值参数确定

https://mgimond.github.io/Spatial/spatial-interpolation.html

tmap空间插值

https://mgimond.github.io/Spatial/interpolation-in-r.html

R中的点插值

https://www.cdrc.ac.uk/wp-content/uploads/2016/11/Practical_11.html

tmap实例

https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html

tmap绘制人口调查地图

http://zevross.com/blog/2018/10/02/creating-beautiful-demographic-maps-in-r-with-the-tidycensus-and-tmap-packages/

tmap分面及布局

https://geocompr.robinlovelace.net/adv-map.html

sf介绍

https://cran.r-project.org/web/packages/sf/vignettes/sf1.html

SpatialLinesDataFrame创建

https://gis.stackexchange.com/questions/163286/how-do-i-create-a-spatiallinesdataframe-from-a-dataframe

SpatialPolygonsDataFrame创建

https://www.rdocumentation.org/packages/sp/versions/1.3-1/topics/SpatialPolygonsDataFrame-class

SPDF转化为SLDF

https://gis.stackexchange.com/questions/200679/convert-spatialpolygonsdf-boundaries-to-spatiallinesdf-keeping-information-on-po

sp对象及栅格数据介绍(最全)

https://rspatial.org/spatial/8-rastermanip.html

合并多个SPDF

https://gis.stackexchange.com/questions/155328/merging-multiple-spatialpolygondataframes-into-1-spdf-in-r

Kriging插值

https://rpubs.com/nabilabd/118172

IDW插值

https://nceas.github.io/oss-lessons/spatial-data-gis-law/4-tues-spatial-analysis-in-r.html

akima插值与ggplot2绘图

https://stackoverflow.com/questions/50533738/best-method-of-spatial-interpolation-for-geographic-heat-contour-maps

IDW插值与ggplot2

http://aasa.ut.ee/LOOM02331/R_idw_interpolation.html

kknn插值与ggplot2

https://timogrossenbacher.ch/2018/03/categorical-spatial-interpolation-with-r/

反距离加权插值

https://www.jianshu.com/p/b38c5e464d16

将rasterLayer对象转化为SpatialGrid/SpatialPixels对象

https://gis.stackexchange.com/questions/43707/how-to-produce-spatial-grid-from-raster

KNN插值原理

http://www.kuqin.com/algorithm/20120817/329048.html

往期精彩：

閱讀更多 天善智能 的文章

關鍵字: 线性插值 html 克里