#读入数据
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 插值为什么要建立格网
• 格网的分辨率就是栅格的分辨率？
• 通过格网来创建地图？