泰森多邊形
#讀入數據
point<-st_read("dian.shp")
point<-st_transform(point,32619)
point<-as(point,"Spatial")
shunyi<-st_read("shunyi.shp")
#創建泰森多邊形，dirichlet()，來自spatstat包
ts<-dirichlet(as.ppp(point))
#將泰森多邊形轉換為空間數據
ts<-as(ts,"SpatialPolygons")
#給泰森多邊形添加坐標系
proj4string(ts)<-proj4string(point)
#over()來自sp包，提取ts對應point位置的屬性值，返回數據框
ts_dataframe<-over(ts,point,fn=mean)
#將屬性數據加入泰森多邊形
ts<-SpatialPolygonsDataFrame(ts,ts_dataframe)
shunyi<-st_transform(shunyi,32619)
shunyi<-as(shunyi,"Spatial")
shunyi<-aggregate(shunyi)
ts<-intersect(shunyi,ts)
反距離加權插值
反距離加權插值是一種常用的插值方法，gstat包有反距離加權的模型函數idw()
idw(formula, locations, data, newdata, nmax = Inf, nmin = 0,
omax = 0, maxdist = Inf, block = numeric(0),
na.action = na.pass, idp = 2.0, debug.level = 1)
參數解釋：
- formula: 插值表達式，一般為
var~1； - locations: 已知點對象的空間數據；
- data: 已知點對象的屬性數據，若
locations中包含了相關屬性數據應省略； - newdata: 插值點對象；
- nmax、nmin: 插值任意位置的屬性值至多和至少所需的已知點數；
- maxdist: 距離閾值；與插值點超過這個距離的已知點不參與該插值點屬性值的計算；
- idp: 距離衰減參數，默認取值為2。
創建插值所需的網格：
#創建網格,sf包
shunyigrid<-st_make_grid(shunyi,n=c(100,100))
#插值，生成結果為面數據
idw1<-gstat::idw(Cd ~ 1,point,shunyigrid,idp=2)
#將插值結果存為文本
save(idw1,file = "idw1.Rdata")
#將格網數據裁剪為順義區
idw1<-st_intersection(idw1,st_union(shunyi))
#可視化
ggplot()+geom_sf(data=idw1,aes(fill=var1.pred))
tmap包
tm_shape(inter_idw)+tm_polygons("var1.pred",border.alpha = 0)+tm_shape(shunyibj)+tm_lines()
不知道為什么不可以用tmap包，出現Error: Shape contains invalid polygons. Please fix it or set tmap_options(check.and.fix = TRUE) and rerun the plot錯誤
已解決。設置tmap_options(check.and.fix=TRUE)就可以進行正常畫圖。
問題
1 模型如何映射到整個地圖
- 將地圖格網化？
2 插值為什么要建立格網 - 格網的分辨率就是柵格的分辨率？
- 通過格網來創建地圖？