這個方法可用於類別變數的插值。實際上是利用點生產沃諾尼多邊形,每個點所在的沃諾尼多邊形的值就等於點值,這樣就實現了由點到面的變化,完成了插值。**裡的p是乙個由很多點組成的空間物件,boundry即是我們的研究區域。p裡的所有點應該在boundry的範圍之內。在製作柵格、插值的過程中,都是基於p的外廓邊界(bbox,這個中文是我隨便起的,不知道對不對),即由所有座標的極值圍合形成矩形,所以為獲得我們需要的研究邊界範圍內的插值資料還需要用boundry來進行裁剪。
1. sp物件版本
library(dismo)
pv ba
#boundry如果內部包含了許多更小的多邊形,需要先進行合併
pv_b
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包裡的函式
不得不說,由於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是個很靈活的函式,蒙板的物件是向量或者柵格都可以。
r語言插補法 R語言 缺失值處理之多重插補
筆者寄語 缺失值是資料清洗過程中非常重要的問題 其他方法可見 r語言 異常值檢驗 離群點分析 異常值處理 筆者在進行mice包的多重插補過程中遇到相當多的問題。大致的步驟簡介如下 缺失資料集 mcmc估計插補成幾個資料集 每個資料集進行插補建模 glm lm模型 將這些模型整合到一起 pool 評價...
R語言中的空間最鄰近點問題
在空間分析中,nearest point problem是乙個十分基本而重要的問題。對於小樣本量的點集而言,通過計算距離矩陣並進行排序即可暴力解決。而對於上萬甚至幾十萬以上的點集而言,通過計算距離矩陣顯然需要耗費太多的資源,一般採用kd樹等方法進行搜尋。背後的演算法我也不太懂,這裡僅把r語言中面向s...
做空間插值 matlab
計算各站點冬 夏兩季降水值 對降水值插值,轉換到格點資料上,便於做空間相關 clear clc 1 讀取所有站點資料 file uigetfile multiselect on 得到有資料的站點名 info load e info.txt data l length file summer wint...