R中点特征到最近多边形的距离

Posted 2023-02-19

【中文标题】R中点特征到最近多边形的距离

【英文标题】：Distance of point feature to nearest polygon in R

【发布时间】：2013-05-03 03:20:58

【问题描述】：

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

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

dist_cen_polya <- dist2Line(centroids_coor, polygonA_shp)

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

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

谢谢。

【问题讨论】：

从你的最后一段来看，你似乎有一个数学问题：找到一个比进行以下比较更好的算法，对吧？这可能更适合数学 SE。

spatstat 包也许可以做你想做的事。我不是该工具集的专家，因此无法确定。

有了这里涉及的数字，10 个多边形和 142 个点（1420 个距离！）蛮力应该不是问题。如果你不喜欢 for 循环，plyr 包应该可以帮助你。

如果您的多边形足够小可以近似为点（即从点到多边形中心的距离足够接近从点到多边形边缘的距离），您可以使用 Voronoi 镶嵌（又名 Dirichlet镶嵌）在多边形的中心（使用deldir包），然后一个点所在的Voronoi瓦片将对应于它最近的多边形。

@shujaa 我必须对另外二十组点和多边形执行相同的程序。因此，我将接近 20,000 距离。我希望让它不那么啰嗦。

【参考方案1】：

这里我使用的是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 向量包含要素投影单位中的实际笛卡尔最小距离。

【讨论】：

@ Jeffrey Evans 我非常感谢您提供的解决方案。我已经将我的数据与代码一起使用了，它真的很好用！谢谢！！

@ Jeffrey Evans，似乎 gDistance 中的 byid 参数给出的结果与您的双循环相同？ gDistance(pts,polys, byid=T) - 一个 10 x 142 矩阵，然后您可以对其进行操作以获得最小距离等。

【参考方案2】：

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

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

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

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

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

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

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

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

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

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

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

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

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