計算多邊形和點之間的距離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