将多边形分配给 R 数据框中的数据点
Assign polygon to data point in R dataframe
我有两个数据框:
points
包含一系列具有 x, y
坐标的点。
poly
包含两个多边形的坐标(我实际上有100多个,但这里保持简单)。
我希望能够向数据框 points
添加一个名为 Area
的额外列,其中包含该点所在的多边形的名称。
poly <- data.frame(
pol= c("P1", "P1","P1","P1","P1","P2","P2","P2","P2", "P2"),
x=c(4360, 7273, 7759, 4440, 4360, 8720,11959, 11440,8200, 8720),
y=c(1009, 9900,28559,28430,1009,9870,9740,28500,28040,9870))
points <- data.frame(
object = c("P1", "P1","P1","P2","P2","P2"),
timestamp= c(1485670023468,1485670023970, 1485670024565, 1485670025756,1485670045062, 1485670047366),
x=c(6000, 6000, 6050, 10000, 10300, 8000),
y=c(10000, 20000,2000,5000,20000,2000))
plot(poly$x, poly$y, type = 'l')
text(points$x, points$y, labels=points$object )
所以基本上在这个例子中,前两行应该有 Area= "P1"
而最后一个点应该是空白的,因为该点不包含在任何多边形中。
我已尝试使用函数 in.out
,但无法按照我的描述构建我的数据框。
非常感谢任何帮助!
您可以使用包 sp
中的函数 point.in.polygon
:
points$Area = apply(points, 1, function(p)ifelse(point.in.polygon(p[3], p[4], poly$x[which(poly$pol==p[1])], poly$y[which(poly$pol==p[1])]), p[1], NA))
给你
object timestamp x y Area
1 P1 1.48567e+12 6000 10000 P1
2 P1 1.48567e+12 6000 20000 P1
3 P1 1.48567e+12 6050 2000 <NA>
4 P2 1.48567e+12 10000 5000 <NA>
5 P2 1.48567e+12 10300 20000 P2
6 P2 1.48567e+12 8000 2000 <NA>
虽然这是使用 for
循环,但实际上速度非常快。
library(mgcv)
x <- split(poly$x, poly$pol)
y <- split(poly$y, poly$pol)
todo <- 1:nrow(points)
Area <- rep.int("", nrow(points))
pol <- names(x)
# loop through polygons
for (i in 1:length(x)) {
# the vertices of i-th polygon
bnd <- cbind(x[[i]], y[[i]])
# points to allocate
xy <- with(points, cbind(x[todo], y[todo]))
inbnd <- in.out(bnd, xy)
# allocation
Area[todo[inbnd]] <- pol[i]
# update 'todo'
todo <- todo[!inbnd]
}
points$Area <- Area
其效率的两个原因:
for
循环是通过多边形,而不是点。因此,如果您有 100 个多边形和 100000 个点要分配,则循环只有 100 次迭代而不是 100000 次。在每次迭代中,C 函数 in.out
的矢量化能力被利用;
- 它以渐进的方式工作。一旦一个点被分配,它将被排除在以后的分配之外。
todo
变量控制通过循环分配的点数。随着时间的推移,工作集正在减少。
我有两个数据框:
points
包含一系列具有x, y
坐标的点。poly
包含两个多边形的坐标(我实际上有100多个,但这里保持简单)。
我希望能够向数据框 points
添加一个名为 Area
的额外列,其中包含该点所在的多边形的名称。
poly <- data.frame(
pol= c("P1", "P1","P1","P1","P1","P2","P2","P2","P2", "P2"),
x=c(4360, 7273, 7759, 4440, 4360, 8720,11959, 11440,8200, 8720),
y=c(1009, 9900,28559,28430,1009,9870,9740,28500,28040,9870))
points <- data.frame(
object = c("P1", "P1","P1","P2","P2","P2"),
timestamp= c(1485670023468,1485670023970, 1485670024565, 1485670025756,1485670045062, 1485670047366),
x=c(6000, 6000, 6050, 10000, 10300, 8000),
y=c(10000, 20000,2000,5000,20000,2000))
plot(poly$x, poly$y, type = 'l')
text(points$x, points$y, labels=points$object )
所以基本上在这个例子中,前两行应该有 Area= "P1"
而最后一个点应该是空白的,因为该点不包含在任何多边形中。
我已尝试使用函数 in.out
,但无法按照我的描述构建我的数据框。
非常感谢任何帮助!
您可以使用包 sp
中的函数 point.in.polygon
:
points$Area = apply(points, 1, function(p)ifelse(point.in.polygon(p[3], p[4], poly$x[which(poly$pol==p[1])], poly$y[which(poly$pol==p[1])]), p[1], NA))
给你
object timestamp x y Area
1 P1 1.48567e+12 6000 10000 P1
2 P1 1.48567e+12 6000 20000 P1
3 P1 1.48567e+12 6050 2000 <NA>
4 P2 1.48567e+12 10000 5000 <NA>
5 P2 1.48567e+12 10300 20000 P2
6 P2 1.48567e+12 8000 2000 <NA>
虽然这是使用 for
循环,但实际上速度非常快。
library(mgcv)
x <- split(poly$x, poly$pol)
y <- split(poly$y, poly$pol)
todo <- 1:nrow(points)
Area <- rep.int("", nrow(points))
pol <- names(x)
# loop through polygons
for (i in 1:length(x)) {
# the vertices of i-th polygon
bnd <- cbind(x[[i]], y[[i]])
# points to allocate
xy <- with(points, cbind(x[todo], y[todo]))
inbnd <- in.out(bnd, xy)
# allocation
Area[todo[inbnd]] <- pol[i]
# update 'todo'
todo <- todo[!inbnd]
}
points$Area <- Area
其效率的两个原因:
for
循环是通过多边形,而不是点。因此,如果您有 100 个多边形和 100000 个点要分配,则循环只有 100 次迭代而不是 100000 次。在每次迭代中,C 函数in.out
的矢量化能力被利用;- 它以渐进的方式工作。一旦一个点被分配,它将被排除在以后的分配之外。
todo
变量控制通过循环分配的点数。随着时间的推移,工作集正在减少。