R语言地理空间分析(四)空间插值

泰森多边形

#读入数据
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 插值为什么要建立格网
  • 格网的分辨率就是栅格的分辨率?
  • 通过格网来创建地图?
原文地址:https://www.cnblogs.com/cqy-wt1124/p/15126983.html