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

Question

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

我的第一个想法是计算多边形的所有线的最小距离(无限扩展,因此它们是线,而不是线条),然后使用线条的起点和毕达哥拉斯计算从点到每条单独线的距离.

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

Answer 1

您可以使用rgeos包和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库,是几何计算中事实上的标准之一.如果您不想重新发明轮子,这是一个很好的方法.

Answer 2

我决定回归并撰写理论解决方案,仅供后人使用.这不是最简洁的例子,但对于那些想要知道如何手动解决这类问题的人来说,它是完全透明的.

理论算法

首先,我们的假设.

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

在编码之前,我们应该用基本的术语写出我们想要做的事情.我们可以假设多边形和多边形外部点之间的最短距离将始终是两个事物之一:多边形的顶点或两个顶点之间的线上的点.考虑到这一点,我们执行以下步骤:

1. 计算所有顶点和目标点之间的距离.
2. 找到最接近目标点的两个顶点.
3. 如果:(a)两个最接近的顶点不相邻或(b)两个顶点中的任一个的内角大于或等于90度,则最近的顶点是最近的点.计算最近点和目标点之间的距离.
4. 否则,计算两点之间形成的三角形的高度.

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

代码

为了使其正常工作,我们希望避免任何"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