列出属于另一个SpatialPoint的gBuffer的所有SpatialPoints

问题描述 投票:0回答:2

我想创建一个data.frame,其中每一行都列出了位于另一个SpatialPoint缓冲区内的SpatialPoints的Set {}。我正在使用rgeos :: gBuffer来创建我的每个点周围的缓冲区。

以下是一些点位置:

head(x)
        Lon            Lat
1       839171.2    3861540
2       838852.4    3861143
3       838945.9    3861240
4       824506.8    3865499
5       838851.8    3861160
6       827834.7    3878655
7       888196.5    3929905
8       508308.4    4031569
9       838750.5    3864169
10      983995.6    3993308

制作一个* Spatial对象:

coordinates(x) <- ~Lon + Lat

围绕点创建一个缓冲区:

xbuff <- gBuffer(x, width=1000, byid=TRUE)

现在,我如何找到落在xbuff中创建的10个缓冲区中的每个缓冲区中的点列表(如果有的话)?

r geospatial sp
2个回答
4
投票

如何找到属于10个缓冲区中每个缓冲区的点列表(如果有)

你可以做到

x<-read.table(header=T, text="
Lon            Lat
1       839171.2    3861540
2       838852.4    3861143
3       838945.9    3861240
4       824506.8    3865499
5       838851.8    3861160
6       827834.7    3878655
7       888196.5    3929905
8       508308.4    4031569
9       838750.5    3864169
10      983995.6    3993308")
library(sp)
library(rgeos)
coordinates(x) <- ~Lon + Lat
xbuff <- gBuffer(x, width=1000, byid=TRUE)
over(xbuff, x[1:5,], T)
# $`1`
# [1] 1 2 3 5
# 
# $`2`
# [1] 1 2 3 5
# 
# $`3`
# [1] 1 2 3 5
# 
# $`4`
# [1] 4
# 
# $`5`
# [1] 1 2 3 5
# 
# $`6`
# integer(0)
# 
# $`7`
# integer(0)
# 
# $`8`
# integer(0)
# 
# $`9`
# integer(0)
# 
# $`10`
# integer(0)

2
投票

您不需要对此问题使用缓冲区和交叉操作。相反,您可以计算距离矩阵。

x <- matrix(c(839171.2, 838852.4, 838945.9, 824506.8, 838851.8, 827834.7, 888196.5, 508308.4, 838750.5, 983995.6, 3861540, 3861143, 3861240, 3865499, 3861160, 3878655, 3929905, 4031569, 3864169, 3993308), ncol=2)

library(raster)
d <- pointDistance(x, lonlat=FALSE)
diag(d) <- NA
r <- apply(d, 1, function(i) which(i < 1000))
r

这应该在计算上更有效。但是,如果你有很多点,距离矩阵可能变得非常(太大)。在这种情况下,您可以循环数据块。

chunksize <- 5
nr <- nrow(x)
s  <- seq(1, nr, chunksize)
r  <- vector(length=nr, mode='list')
for (i in 1:length(s)) {
    start <- s[i]
    end  <- min(nr, start + chunksize)
    y <- x[start:end, ,drop=FALSE]
    d <- pointDistance(y, x, lonlat=FALSE)
##  d[cbind(1:nrow(y), start:end)] <- NA
    r[start:end] <- apply(d, 1, function(i) which(i < 1000))
}

这包括联络点。您可以将“对角线”设置为NA以避免这种情况,但如果范围内没有点,则会导致错误,因此我评论了帽子。

© www.soinside.com 2019 - 2024. All rights reserved.