自然邻点插值。计算多边形相交面积时出错

我正在尝试写this http://en.wikipedia.org/wiki/Natural_neighborR 中的算法。它已经存在于任何包中吗？！？

这就是我所做的（在SO和各种博客文章的帮助下）：

library(rgdal)
library(ggmap)
require("maptools")
require("plyr")

locations<- unique(cbind(data22[,1], data22[,2]))
      [,1]    [,2]
  [1,] 24.9317 60.1657
  [2,] 24.9415 60.1608
  [3,] 24.9331 60.1577
  [4,] 24.9228 60.1477
  [5,] 24.9370 60.1545
  [6,] 24.9491 60.1559
  [7,] 24.9468 60.1591
  [8,] 24.9494 60.1675
  [9,] 24.9561 60.1609
 [10,] 24.9218 60.1632
 [11,] 24.9213 60.1605
 [12,] 24.9219 60.1557
 [13,] 24.9208 60.1704
 [14,] 24.9233 60.1714
 [15,] 24.9469 60.1737
 [16,] 24.9440 60.1738
 [17,] 24.9531 60.1714
 [18,] 24.9601 60.1736
 [19,] 24.9304 60.1687
 [20,] 24.9312 60.1659
 [21,] 24.9313 60.1658
 [22,] 24.9418 60.1608
 [23,] 24.9336 60.1577
 [24,] 24.9213 60.1494
 [25,] 24.9415 60.1538
 [26,] 24.9560 60.1620
 [27,] 24.9610 60.1587
 [28,] 24.9142 60.1635
 [29,] 24.9072 60.1636
 [30,] 24.9132 60.1582
 [31,] 24.9166 60.1668
 [32,] 24.9146 60.1742
 [33,] 24.9259 60.1751
 [34,] 24.9308 60.1742
 [35,] 24.9524 60.1690
 [36,] 24.9601 60.1709
 [37,] 24.9570 60.1742
 [38,] 24.9324 60.1655
 [39,] 24.9426 60.1610
 [40,] 24.9332 60.1581
 [41,] 24.9274 60.1480
 [42,] 24.9393 60.1539
 [43,] 24.9466 60.1550
 [44,] 24.9478 60.1593
 [45,] 24.9431 60.1670
 [46,] 24.9559 60.1615
 [47,] 24.9623 60.1581
 [48,] 24.9144 60.1632
 [49,] 24.9077 60.1634
 [50,] 24.9110 60.1575
 [51,] 24.9212 60.1685
 [52,] 24.9193 60.1739
 [53,] 24.9270 60.1752
 [54,] 24.9305 60.1746
 [55,] 24.9517 60.1700
 [56,] 24.9598 60.1710
 [57,] 24.9565 60.1737
 [58,] 24.9306 60.1686
 [59,] 24.9361 60.1621
 [60,] 24.9415 60.1580
 [61,] 24.9312 60.1561
 [62,] 24.9253 60.1528
 [63,] 24.9501 60.1589
 [64,] 24.9467 60.1591
 [65,] 24.9458 60.1630
 [66,] 24.9374 60.1715
 [67,] 24.9438 60.1707
 [68,] 24.9527 60.1674
 [69,] 24.9556 60.1604
 [70,] 24.9205 60.1698
 [71,] 24.9141 60.1633
 [72,] 24.9082 60.1633
 [73,] 24.9118 60.1569
 [74,] 24.9220 60.1683
 [75,] 24.9231 60.1630
 [76,] 24.9475 60.1735
 [77,] 24.9434 60.1735
 [78,] 24.9535 60.1713
 [79,] 24.9605 60.1739
 [80,] 24.9307 60.1685
 [81,] 24.9373 60.1618
 [82,] 24.9402 60.1582
 [83,] 24.9311 60.1560
 [84,] 24.9257 60.1527
 [85,] 24.9485 60.1589
 [86,] 24.9460 60.1635
 [87,] 24.9374 60.1709
 [88,] 24.9519 60.1673
 [89,] 24.9554 60.1595
 [90,] 24.9228 60.1629
 [91,] 24.9215 60.1602
 [92,] 24.9217 60.1556
 [93,] 24.9212 60.1706
 [94,] 24.9239 60.1715
 [95,] 24.9466 60.1735
 [96,] 24.9436 60.1740
 [97,] 24.9532 60.1715
 [98,] 24.9609 60.1738
 [99,] 24.9354 60.1626
[100,] 24.9351 60.1626
[101,] 24.9374 60.1579
[102,] 24.9300 60.1542
[103,] 24.9263 60.1529
[104,] 24.9522 60.1589
[105,] 24.9435 60.1622
[106,] 24.9369 60.1721
[107,] 24.9580 60.1615
[108,] 24.9620 60.1586
[109,] 24.9545 60.1545

# Carson's Voronoi polygons function
# http://stackoverflow.com/a/9405831/489704
# http://www.carsonfarmer.com/2009/09/voronoi-polygons-with-r/
voronoipolygons <- function(x) {
  require(deldir)
  require(sp)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  z <- deldir(crds[, 1], crds[, 2])
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1, ])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)
  voronoi <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[, 1],
    y=crds[,2], row.names=sapply(slot(SP, 'polygons'), 
    function(x) slot(x, 'ID'))))
}


v2 <- voronoipolygons(locations)
a=fortify(v2)

bbox<-c(24.90, 60.14, 
        24.97, 60.18)
predgrid <- expand.grid(lon=seq(from=bbox[1], to=bbox[3], length.out=10), 
                        lat=seq(from=bbox[2], to=bbox[4], length.out=10))

N <- 100
loc <- as.matrix(cbind(predgrid[,1:2]))
proj4string(v2) <- CRS("+proj=longlat +ellps=WGS84 +no_defs+ellps=WGS84 
                       +towgs84=0,0,0 ") 
int<-matrix(nrow=N, ncol=109)
for (i in 1:N){
  loc.temp<-rbind(locations, loc[i,])
  vor.temp<-voronoipolygons(loc.temp)
  proj4string(vor.temp) <- CRS("+proj=longlat +ellps=WGS84 +no_defs+ellps=WGS84 
                               +towgs84=0,0,0 ") 
  for (j in 1:109){
    s<-gIntersection(SpatialPolygons(vor.temp@polygons[110]), 
                     SpatialPolygons(vor.temp@polygons[i]))
    if(is.null(s)) {
      int[i,j]=0
    } else {
      # my.area <- vor.temp@polygons[[i]]@Polygons[[1]]@area 
      int[i,j] <- 
        gArea(gIntersection(SpatialPolygons(vor.temp@polygons[i]), 
                            SpatialPolygons(vor.temp@polygons[overlaid.poly])))
    }
  }
}



max(int)

所以基本上我画了一次添加一个点的voronoi缨子，并尝试计算多边形之间的交集，问题是：max(int)给出0，就好像没有插值一样，为什么会这样呢？是计算面积gIntersection正确的？它们的值非常小，我不知道度量单位。

EDIT: Tassellation before new location and the one after it . I am not sure about the areas given back by the tassellation, that is where I would think the error is, I am not sure though.

这段代码的问题在于第二个 for 循环，其中函数gIntersection已经用过。在此函数中多边形的交集：

SpatialPolygons(vor.temp@polygons[110])

并且 vor.temp 的所有多边形都已计算出来，显然这些多边形之间不会有交集，因为它们都属于一个 Voronoi。在里面gIntersection第二个参数应更改为：

SpatialPolygons(vor.temp@polygons[i])

to:

SpatialPolygons(v2@polygons[i])

Now the gIntersection函数将获得交集多边形SpatialPolygons(vor.temp@polygons[110])以及来自第一个 Voronoi 的所有多边形 (v2)

另一个问题是使用gArea功能。我没有明白使用的目的overlaid.poly。函数的输入gArea是一个多边形，所以我们可以将该函数的输入设置为s因为已经使用获得了相交多边形gIntersection作为对象的函数s。所以在 if 语句中我们会有：

gArea(s)

代替

gArea(gIntersection(SpatialPolygons(vor.temp@polygons[i]),SpatialPolygons(vor.temp@polygons[overlaid.poly])))

由此产生的对象int将是可用于执行插值的权重矩阵。