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

Question

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

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

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

Answer 1

您可以使用rgeos package和gDistance方法。这将需要你准备几何图形，根据你拥有的数据创建对象（我假设它是一个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

否则：

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 

Answer 3

我决定回，写了一个理论上的解决方案，只是为后人。这不是最简洁的例子，但对于那些想要知道如何去解决这样的问题的人来说，它是完全透明的。

理论算法

首先，我们的假设。

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

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

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

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

代码

为了使这项工作正常，我们要避免任何“的”循环和希望只使用矢量功能看多边形顶点的整个列表时。幸运的是，这在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 

Answer 4

我在geosphere封装中使用distm()功能时被呈现点和顶点坐标系来计算distence。另外，您可以通过​​替代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