7

我目前正在做一个项目，我有一个点特征——点特征包括一个 142 个点——和多个多边形（大约 10 个）。我想计算 R 中每个点与最近的多边形特征之间的距离。

我目前的方法很乏味，而且有点啰嗦。我目前正计划计算每个点和每个多边形之间的距离。例如，我将计算 142 点与多边形 A 之间的距离、142 点与多边形 B 之间的距离、142 点与多边形 C 之间的距离等。以下是其中一种距离计算的示例代码：

dist_cen_polya <- dist2Line(centroids_coor, polygonA_shp)

完成这些计算后，我会编写一个代码来选择每个点和最近的多边形之间的最小/最近距离。问题是这个过程很乏味。

有谁知道可以最大限度地减少计算的工作量/计算时间的包/代码？我真的很想使用一个包来比较一个指向最近的多边形特征或计算一个点和所有感兴趣的多边形之间的距离？

谢谢你。

4

2 回答 2

14

这里我使用的是 rgeos 拓扑库中的gDistance函数。我正在使用蛮力双循环，但速度惊人。142 个点和 10 个多边形只需不到 2 秒。我确信有一种更优雅的方式来执行循环。 

   require(rgeos)

    # CREATE SOME DATA USING meuse DATASET
    data(meuse)
      coordinates(meuse) <- ~x+y
        pts <- meuse[sample(1:dim(meuse)[1],142),]  
    data(meuse.grid) 
      coordinates(meuse.grid) = c("x", "y") 
        gridded(meuse.grid) <- TRUE 
          meuse.grid[["idist"]] = 1 - meuse.grid[["dist"]]    
        polys <- as(meuse.grid, "SpatialPolygonsDataFrame")
          polys <- polys[sample(1:dim(polys)[1],10),]   
    plot(polys)
      plot(pts,pch=19,cex=1.25,add=TRUE)      

    # LOOP USING gDistance, DISTANCES STORED IN LIST OBJECT
    Fdist <- list()
      for(i in 1:dim(pts)[1]) {
        pDist <- vector()
          for(j in 1:dim(polys)[1]) { 
            pDist <- append(pDist, gDistance(pts[i,],polys[j,])) 
          }
        Fdist[[i]] <- pDist
      } 

    # RETURN POLYGON (NUMBER) WITH THE SMALLEST DISTANCE FOR EACH POINT  
    ( min.dist <- unlist(lapply(Fdist, FUN=function(x) which(x == min(x))[1])) ) 

    # RETURN DISTANCE TO NEAREST POLYGON
    ( PolyDist <- unlist(lapply(Fdist, FUN=function(x) min(x)[1])) ) 

    # CREATE POLYGON-ID AND MINIMUM DISTANCE COLUMNS IN POINT FEATURE CLASS
    pts@data <- data.frame(pts@data, PolyID=min.dist, PDist=PolyDist)

    # PLOT RESULTS
    require(classInt)
    ( cuts <- classIntervals(pts@data$PDist, 10, style="quantile") )
       plotclr <- colorRampPalette(c("cyan", "yellow", "red"))( 20 )
         colcode <- findColours(cuts, plotclr)
    plot(polys,col="black")
      plot(pts, col=colcode, pch=19, add=TRUE)

min.dist 向量表示多边形的行号。例如，您可以通过使用此向量来对最近的多边形进行子集化。

near.polys <- polys[unique(min.dist),]

PolyDist 向量包含要素投影单位中的实际笛卡尔最小距离。

于 2013-05-09T00:12:46.257 回答
1

在多边形中，您有很多线条。如果点位于多边形内或位于边上，则多边形之间的距离为零。

所以实际上你正在寻找两种情况：

  1. 检查该点是否位于任何多边形内（记住它可能比单个多边形更多
  2. 获取所有边缘的集合并计算点到边缘的每个距离。最近的距离也为您提供了边缘所属的多边形的距离。

所以这是一个简单的算法，如果我们考虑每个多边形有 10 条边，那么所有点都需要 O(10 * 10) * 142。这使得 100 * 142 = 14200 次操作。=> O(m * deltaE) * n （m 是多边形数，deltaE 是每个多边形的平均边数，n 是点数）。

现在我们检查是否可以加快速度。首先想到的是我们可以为每个多边形使用边界框检查或边界圆。

另一个想法是为每个多边形准备一组角度的最近边。例如，如果您有 8 个角度（每 45°），您可以从列表中删除被另一条边取代的所有边（因此，被删除边的任何点将始终产生比同一条另一边的任何点更大的距离多边形。

这种方式通常可以大大降低复杂性。如果您想到一个矩形，那么您只有一个或两个边而不是 4 个。如果您想到一个规则的 8 边多边形，您最终可能会得到一个或两个边，每个多边形最多有 3 个边。

如果将法线向量添加到每个边缘，则可以计算该点是否在内部，并且您必须执行扫描线或任何检查（或者您知道它的 konvex）才能确定。

也有可能的映射索引，例如以等距方式将二维空间通过 x 和 y 分开。这样，您只需测试位于 9 个扇区中的多边形。

下一个版本可能会使用 R-tree，每个节点的每个边界框（圆）都必须检查最小预期距离和最大预期距离。因此，不需要检查导致最小距离远大于另一个节点的最大距离的节点的多边形。

另一件事是，如果你有一个给定的树，就像你对地图数据的顺序一样。在街道地图中，您总是有世界 -> 地区 -> 国家 -> 县 -> 城市 -> 城市部门 ->...

通过这种方式，您可以在包含数百万个多边形的全球地图中搜索最近的位置，并且在合理的时间内（通常小于 10 毫秒）。

可以这么说，您在这里有很多选择。并预处理您的多边形列表并通过使用多边形的二进制空间分区树或使用角度方法或什至使用更花哨的方法来提取相关边缘。由你决定。我希望你最终会在对数范围内做一些事情，比如 O(log(n) * log(deltaE)) 变成 O(log(n)) 作为平均复杂度。

于 2014-06-13T10:00:41.850 回答