问题描述
我正在尝试在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!
这篇关于自然邻居插值.多边形相交面积计算中的错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！