自然邻居插值.多边形相交面积计算中的错误 [英] natural neighbour interpolation. error in calculating polygon intersection area
问题描述
我正在尝试在R中编写此算法.它已经存在于任何软件包中吗?！?
I am trying to write this algorithm in R. Does it exist in any package already?!?
这就是我所做的(在SO和各种博客文章的帮助下):
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)
所以基本上我画一次将voronoi镶嵌添加一个点，然后尝试计算多边形之间的交点，问题是:max(int)给出0，好像没有插值，为什么会这样呢? gIntersection计算的面积正确吗?它们的价值很小，我不知道措施的统一性.
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.
在新位置之前进行细化处理及其后的内容. 我不确定流失所给的领域，也就是我认为错误所在，但我不确定.
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.
推荐答案
代码中的小而重要的错误:内部循环在 j 上，因此第二个参数必须为SpatialPolygons(v2$polygons[j]).现在，代码对我而言非常完美！
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!
这篇关于自然邻居插值.多边形相交面积计算中的错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
