natural neighbour interpolation. error in calculat

I am trying to write this algorithm in R. Does it exist in any package already?!?

This is what I did (with help from SO and various blog posts):

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)

So basically I draw the voronoi tassellation adding one point at a time and try to calculate the intersection between polygons, problem is: max(int) gives 0, as if there was no interpolation, why is this so? is the area calculated with gIntersection correct? they are super small value and I have no idea of the measure unity.

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.

回答1:

The problem with this code is in the second for loop where the function gIntersection has been used. In this function the intersection of the polygon:

SpatialPolygons(vor.temp@polygons[110])

and all the polygons of vor.temp has been calculated and it is obvious that there will be no intersection between those polygons because they all are belong to one Voronoi. In the gIntersection the second argument should be changed from:

SpatialPolygons(vor.temp@polygons[i])

to:

SpatialPolygons(v2@polygons[i])

Now the gIntersection function will obtain the intersection polygon of the SpatialPolygons(vor.temp@polygons[110]) and all the polygons from the first Voronoi (v2)

The other problem is using gArea function. I did not get the purpose of using overlaid.poly. The input of function gArea is a polygon, so we can set input of this function as s since the intersection polygon has been obtain already using gIntersection function as object s. So in the if statement we would have:

gArea(s)

instead of

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

The resultant object of int would be the weight matrix which can be used to perform interpolation.

回答2:

Small but important error in the code: the inner loop is over j thus the second argument must be SpatialPolygons(v2$polygons[j]). Now the code does it perfectly for me!