计算R中多边形和点之间的距离答案

【问题标题】：Calculating the distance between polygon and point in R计算R中多边形和点之间的距离
【发布时间】：2023-03-25 12:32:01
【问题描述】：

我有一个不一定是凸面的没有交点的多边形和这个多边形外的一个点。我想知道如何在二维空间中最有效地计算欧几里得距离。 R有标准方法吗？

我的第一个想法是计算多边形所有线的最小距离（无限延伸，因此它们是线，而不是线段），然后使用线段的起点计算从点到每条线的距离和毕达哥拉斯。

你知道实现高效算法的包吗？

【问题讨论】：

我不知道有什么预包装的东西可以满足您的要求。你能给我们更多的信息吗？在计算距离之前，您是否对多边形有所了解，或者它可以是任何形状？既然你提到了毕达哥拉斯，我假设这是在笛卡尔空间中。
根据建议的方法，听起来多边形不允许内部相交，这表明顶点具有旋转顺序（顺时针或逆时针）。这是真的吗？
你看过spatstat 包吗？我知道他们至少有一些用于线段和点集的距离算法。

标签： r gis

【解决方案1】：

您可以使用rgeos package 和gDistance 方法。这将要求您准备几何图形，从您拥有的数据创建 spgeom 对象（我假设它是一个 data.frame 或类似的东西）。 rgeos 文档非常详细（请参阅 CRAN 页面中软件包的 PDF 手册），这是gDistance 文档中的一个相关示例：

pt1 = readWKT("POINT(0.5 0.5)")
pt2 = readWKT("POINT(2 2)")
p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
p2 = readWKT("POLYGON((2 0,3 1,4 0,2 0))")
gDistance(pt1,pt2)
gDistance(p1,pt1)
gDistance(p1,pt2)
gDistance(p1,p2)

readWKT 也包含在 rgeos 中。

Rgeos 基于 GEOS 库，它是几何计算中事实上的标准之一。如果您不想重新发明轮子，这是一个不错的选择。

【讨论】：

【解决方案2】：

我决定返回并编写一个理论解决方案，以供后代使用。这不是最简洁的示例，但对于那些想知道如何手动解决此类问题的人来说，它是完全透明的。

理论算法

首先，我们的假设。

我们假设多边形的顶点以顺时针或逆时针的旋转顺序指定多边形的点，并且多边形的线不能相交。这意味着我们有一个正常的几何多边形，而不是一些奇怪定义的矢量图形。
我们假设这是一组笛卡尔坐标，使用代表二维平面上的位置的“x”和“y”值。
我们假设该点必须在多边形的内部区域之外。
最后，我们假设所需距离是该点与多边形周长上所有无限个点之间的最小距离。

现在在编码之前，我们应该用基本的术语写出我们想要做什么。我们可以假设多边形与多边形外点之间的最短距离始终是以下两种情况之一：多边形的一个顶点或两个顶点之间的一条线上的一个点。考虑到这一点，我们执行以下步骤：

计算所有顶点到目标点的距离。
找到离目标点最近的两个顶点。
如果是： (a) 两个最近的顶点不相邻或 (b) 两个顶点中任一个的内角大于或等于 90 度，那么最近的顶点就是最近的点。计算最近点与目标点之间的距离。
否则，计算两点之间形成的三角形的高度。

我们基本上只是查看一个顶点是否最接近该点，或者一条线上的一个点是否最接近该点。我们必须使用一些三角函数来完成这项工作。

代码

为了使其正常工作，我们希望避免任何“for”循环，并希望在查看整个多边形顶点列表时只使用矢量化函数。幸运的是，这在 R 中非常简单。我们接受一个带有“x”和“y”列的数据框作为多边形的顶点，我们接受一个带有一个“x”和“y”值的向量作为点的位置。

get_Point_Dist_from_Polygon <- function(.polygon, .point){

    # Calculate all vertex distances from the target point.
    vertex_Distance <- sqrt((.point[1] - .polygon$x)^2 + (.point[2] - .polygon$y)^2)

    # Select two closest vertices.
    min_1_Index <- which.min(vertex_Distance)
    min_2_Index <- which.min(vertex_Distance[-min_1_Index])

    # Calculate lengths of triangle sides made of
    # the target point and two closest points.
    a <- vertex_Distance[min_1_Index]
    b <- vertex_Distance[min_2_Index]
    c <- sqrt(diff(.polygon$x[c(min_1_Index, min_2_Index)])^2 + diff(.polygon$y[c(min_1_Index, min_2_Index)])^2)

    if(abs(min_1_Index - min_2_Index) != 1 |
        acos((b^2 + c^2 - a^2)/(2*b*c)) >= pi/2 | 
        acos((a^2 + c^2 - b^2)/(2*a*c)) >= pi/2
        ){
        # Step 3 of algorithm.
        return(vertex_Distance[min_1_Index])
    } else {
        # Step 4 of algorithm.
        # Here we are using the law of cosines.
        return(sqrt((a+b-c) * (a-b+c) * (-a+b+c) * (a+b+c)) / (2 * c))
    }
}

演示

polygon <- read.table(text="
x,  y
0,  1
1,  0.8
2,  1.3
3,  1.4
2.5,0.3
1.5,0.5
0.5,0.1", header=TRUE, sep=",")

point <- c(3.2, 4.1)

get_Point_Dist_from_Polygon(polygon, point)
# 2.707397

【讨论】：

也许我误解了这个解决方案，但是距离该点最近的两个顶点形成的多边形的面不一定是多边形的最接近该点的面。
我完全同意理查德的观点。 Dinre，对于一个简单的示例，您的解决方案失败了，其中查询点位于底部有很长边的四边形下方，其中离查询点最近的两个顶点位于边缘的正上方（不在同一侧边缘作为查询点）。

【解决方案3】：

否则：

p2poly <- function(pt, poly){
    # Closing the polygon
    if(!identical(poly[1,],poly[nrow(poly),])){poly<-rbind(poly,poly[1,])}
    # A simple distance function
    dis <- function(x0,x1,y0,y1){sqrt((x0-x1)^2 +(y0-y1)^2)}
    d <- c()   # Your distance vector
    for(i in 1:(nrow(poly)-1)){
        ba <- c((pt[1]-poly[i,1]),(pt[2]-poly[i,2])) #Vector BA
        bc <- c((poly[i+1,1]-poly[i,1]),(poly[i+1,2]-poly[i,2])) #Vector BC
        dbc <- dis(poly[i+1,1],poly[i,1],poly[i+1,2],poly[i,2]) #Distance BC
        dp <- (ba[1]*bc[1]+ba[2]*bc[2])/dbc          #Projection of A on BC
        if(dp<=0){ #If projection is outside of BC on B side
            d[i] <- dis(pt[1],poly[i,1],pt[2],poly[i,2])
            }else if(dp>=dbc){ #If projection is outside of BC on C side
                d[i] <- dis(poly[i+1,1],pt[1],poly[i+1,2],pt[2])
                }else{ #If projection is inside of BC
                    d[i] <- sqrt(abs((ba[1]^2 +ba[2]^2)-dp^2))
                    }
        }
    min(d)
    }

例子：

pt <- c(3,2)
triangle <- matrix(c(1,3,2,3,4,2),byrow=T, nrow=3)
p2poly(pt,triangle)
[1] 0.3162278

【讨论】：

底部有一个例子。 Poly可以是矩阵、数组或data.frame，只要每一行都是一个顶点坐标即可。这是我能想到的最简单的算法。
我喜欢你的努力，但接受的答案更容易实施。

【解决方案4】：

我使用geosphere 包中的distm() 函数来计算点和顶点在坐标系中出现时的距离。此外，您可以通过替换 dis <- function(x0,x1,y0,y1){sqrt((x0-x1)^2 +(y0-y1)^2)} 轻松地进行一些替换 distm() 。

algo.p2poly <- function(pt, poly){
  if(!identical(poly[1,],poly[nrow(poly),])){poly<-rbind(poly,poly[1,])}
  library(geosphere)
  n <- nrow(poly) - 1
  pa <- distm(pt, poly[1:n, ])
  pb <- distm(pt, poly[2:(n+1), ])
  ab <- diag(distm(poly[1:n, ], poly[2:(n+1), ]))
  p <- (pa + pb + ab) / 2
  d <- 2 * sqrt(p * (p - pa) * (p - pb) * (p - ab)) / ab
  cosa <- (pa^2 + ab^2 - pb^2) / (2 * pa * ab)
  cosb <- (pb^2 + ab^2 - pa^2) / (2 * pb * ab)
  d[which(cosa <= 0)] <- pa[which(cosa <= 0)]
  d[which(cosb <= 0)] <- pb[which(cosb <= 0)]
  return(min(d))
}

例子：

poly <- matrix(c(114.33508, 114.33616,
                 114.33551, 114.33824,
                 114.34629, 114.35053,
                 114.35592, 114.35951, 
                 114.36275, 114.35340, 
                 114.35391, 114.34715,
                 114.34385, 114.34349,
                 114.33896, 114.33917,
                 30.48271, 30.47791,
                 30.47567, 30.47356, 
                 30.46876, 30.46851,
                 30.46882, 30.46770, 
                 30.47219, 30.47356,
                 30.47499, 30.47673,
                 30.47405, 30.47723, 
                 30.47872, 30.48320),
               byrow = F, nrow = 16)
pt1 <- c(114.33508, 30.48271)
pt2 <- c(114.6351, 30.98271)
algo.p2poly(pt1, poly)
algo.p2poly(pt2, poly)

结果：

> algo.p2poly(pt1, poly)
[1] 0
> algo.p2poly(pt2, poly)
[1] 62399.81

【讨论】：