2013-05-08 77 views
7

我目前在一個項目中工作,我有一個點要素 - 點要素包括142個點 - 和多個多邊形(約10個)。我想計算每個點與R中最近的多邊形特徵之間的距離。點特徵到最近多邊形的距離R

我目前的方法很乏味,有點冗長。我目前正計劃計算每個單點和每個多邊形之間的距離。例如,我會計算142個點與多邊形A之間的距離,142個點與多邊形B之間的距離,142個點與多邊形C之間的距離等。下面是這些距離計算之一的示例代碼:

dist_cen_polya <- dist2Line(centroids_coor, polygonA_shp) 

做完這些計算後,我會寫一段代碼來選擇每個點和最近的多邊形之間的最小/最近距離。問題在於這個過程很乏味。

有沒有人知道一個程序包/代碼會盡量減少計算的工作量/計算時間?我真的想用一個包比較單個指向最近的多邊形特徵或計算一個點和所有感興趣的多邊形之間的距離?

謝謝。

+0

從你的最後一段來看,你似乎有一個數學問題:找到一個更好的算法比做一套比較,對吧?這可能更適合數學SE。 – Frank 2013-05-08 19:34:39

+0

'spatstat'包可能可以做你想做的。我不是那個工具集的專家,所以無法確定。 – 2013-05-08 19:39:48

+0

在這裏涉及的數字中,10個多邊形和142個點(1420個距離!)的蠻力應該不成問題。如果你不喜歡for循環,'plyr'軟件包會幫助你。 – Gregor 2013-05-08 19:43:18

回答

13

這裏我使用rgeos拓撲庫中的gDistance函數。我正在使用一個強力雙循環,但它是驚人的快。它需要不到2秒,142點和10個多邊形。我相信有一種更優雅的方式來執行循環。

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 
    [email protected] <- data.frame([email protected], PolyID=min.dist, PDist=PolyDist) 

    # PLOT RESULTS 
    require(classInt) 
    (cuts <- classIntervals([email protected]$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矢量包含要素投影單元中的實際笛卡爾最小距離。

+0

@傑弗裏埃文斯我真的非常感謝你的解決方案。我已經使用我的數據和代碼,它真的很好!謝謝!! – user1738753 2013-05-09 14:09:09

+5

@ Jeffrey Evans,似乎gDistance中的byid參數給出了與雙循環相同的結果? 'gDistance(pts,polys,byid = T)' - 一個10 x 142的矩陣,然後您可以操作它來獲得最小距離等。 – CCID 2013-10-10 20:48:06

1

在一個多邊形中,你有很多線條。如果點位於多邊形內或位於邊上,多邊形之間的距離爲零。

所以,其實你正在尋找兩種情況:

  1. 檢查點位於任意多邊形的內部(記住它可能會更然後一個多邊形
  2. 獲取所有邊的集合,計算每個從邊緣點的距離。 最近的距離,讓您邊屬於太多邊形的距離。

所以這是一個簡單的算法,如果我們認爲每個多邊形10層的邊緣需要O(10 * 10 )* 142你所有的觀點。這使得100 * 142 = 14200操作。 => O(m * deltaE)* n(m是多邊形的數量,deltaE是每個多邊形的平均邊數,n是點的數量)。

現在我們檢查一下我們是否可以加快速度。首先想到的是,我們可以爲每個多邊形使用邊界框檢查或邊界圓。

另一個想法是爲一組角度準備每個多邊形的最近邊緣。例如,如果您有8個角度(每45°),則可以從列表中刪除被另一個邊替代的所有邊(因此,刪除的邊的任意點總是會產生比同一邊的任何點都大的距離多邊形

這種方式通常可以大幅降低複雜度如果您考慮矩形,則只有一個或兩個邊而不是4個。如果您想到常規8邊多邊形,最終可能會出現通常是一個或兩個以及每個多邊形最多3個邊線

如果將法向量添加到每條邊,則可以計算點是否可能在裏面,並且必須執行掃描線或任何檢查(或者您知道它的konvex )確定。

也有映射指數可能像在等距離的甘露中用x和y分隔二維空間。這裏你只需要測試9個扇區中的polygones。

下一個版本可能使用R樹,每個節點的每個邊界框(圓)必須檢查最小預期距離和最大預期距離。因此,不需要檢查導致比另一節點的最大距離大得多的最小距離的節點的多邊形。

另一件事是如果你有像給定地圖數據一樣的給定樹。在街道地圖中,您始終擁有世界 - >地區 - >國家 - >縣 - >城市 - >城市區域 - > ...

此路徑可以搜索全球地圖中最近的位置,多邊形在合理時間內大部分爲< 10ms。

可以這麼說,你在這裏有很多選擇。並預處理您的多邊形列表,並通過使用多邊形的二進制空間分區樹或使用有角度的方法提取相關的邊緣,或甚至更有趣的事物。隨你便。我希望你最終可以像O(log(n)* log(deltaE))那樣在對數範圍內做一些事情,使其變成O(log(n))作爲平均複雜度。

相關問題