必知必會（二）（下）

四川大學在讀研究生

接上文

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 克里