鄰近多邊形插值

這個方法可用於類別變數的插值。實際上是利用點生產沃諾尼多邊形，每個點所在的沃諾尼多邊形的值就等於點值，這樣就實現了由點到面的變化，完成了插值。程式碼裡的p是一個由很多點組成的空間物件，boundry即是我們的研究區域。p裡的所有點應該在boundry的範圍之內。在製作柵格、插值的過程中，都是基於p的外廓邊界（bbox，這個中文是我隨便起的，不知道對不對），即由所有座標的極值圍合形成矩形，所以為獲得我們需要的研究邊界範圍內的插值資料還需要用boundry來進行裁剪。
1. sp物件版本

library(dismo)
pv <- voronoi(p)
ba <- aggregate(boundry)
 
#boundry如果內部包含了許多更小的多邊形，需要先進行合併
pv_b <- intersect(v, ba)

2.sf物件版本

sf物件完成同樣的事情似乎要複雜的多，首先是對sfc物件的型別有著嚴格的要求，必須是MULTIPOINT型別,在st_voronoi函式裡的物件也必須是sfc物件，返回的物件型別是GOMETRYCOLLECTION，完全不方便，所以其實更簡單的方法應該是用sp的方法最後結果轉成sf型別。純sf型別的方法就比較繁瑣了，以下是我摸索出來的。

pv_single=st_coordinates(p)  %>% st_multipoint %>%  st_sfc(crs=st_crs(dta)) %>% 
st_voronoi(st_as_sfc(st_bbox(p)))
#v仍然還是一個sfc物件，而且是隻有一個元素的sfc物件，而我們需要的是每一行即是一個多邊形物件的資料框，顯然不是我們需要的型別，需要進一步處理。
pv_sfc=pv[[1]]
 
 %>%unclass %>% st_sfc(crs=st_crs(p))
這個時候，就可以和原來的資料框組合生成我們需要的沃諾尼多邊形的sf物件了。
pv_sf = st_sf(p,pv)
ba=st_union(boundry)#我們這裡設定boundry是sf物件，那麼合併內部多邊形也需要用sf包裡的函式
pv_b = st_intersection(pv,ba)

不得不說，由於sf物件是新近開發出來的，而R語言中用於地統計分析的軟體包像gstat,dismo 開發得都比較早，基本上還都只是適配sp物件的。

最鄰近插值

前面利用voronoi多邊形的方法更加適用於類別變數的插值，而對於連續型變數，顯然柵格更加適合。這裡需要用到gstat

library(gstat)
#製作柵格蒙板
p_r=raster(p,res)#製作一個空白的柵格，其空間範圍是p的外廓邊界（bbox,即由所有座標的極值圍合形成矩形）,解析度是res
ba=aggregate(boundry)
gs= gstat(formula=var~1,location=p,nmax=n,set=list(idp=0))#nmax即是計算插值時最多把多少個鄰近點計算在內，idp=0即是不考慮距離的影響，不設定即預設採用反距離權重進行插值。
p_interpolate=interpolate(p_r,gs)#本質上就是利用gs物件裡的模型給空白的柵格物件的每一個單元格進行賦值
p_ba = mask(p_interpolate,ba)#mask是個很靈活的函式，蒙板的物件是向量或者柵格都可以。