2012-03-20 93 views
10

我試圖找到一組位置座標和一組線(道路或河流)之間的正交距離。點的集合以緯度/經度對的形式出現,並且線條位於shapefile(.shp)中。使用maptoolsPBSmapping,將它們繪製在地圖上不成問題。但我的基本問題是要找到人們從一個地方到達一條公路或一條河的最小距離。在R中有沒有辦法做到這一點?R/GIS:找到位置和最近線之間的正交距離

回答

19

如果我理解正確,那麼可以使用rgeos包中的gDistance來做到這一點。

讀入以線條爲SpatialLines/DataFrame和點作爲SpatialPoints/DataFrame然後循環在每個點的每個時間計算的距離:

require(rgeos) 
## untested code 
shortest.dists <- numeric(nrow(sp.pts)) 
for (i in seq_len(nrow(sp.pts)) { 
    shortest.dists[i] <- gDistance(sp.pts[i,], sp.lns) 
} 

這裏sp.pts是空間點對象,並且sp.lns是空間線對象。

您必須循環,以便您只將sp.pts中的單個座標與sp.lns中所有線幾何的整體進行比較,否則您將獲得跨所有點的聚合值的距離。

由於您的數據處於經度/緯度範圍,因此gDistance函數假定笛卡爾距離,因此您應將這兩行和點轉換爲合適的投影。

更多討論和實例(編輯)

這將是整齊就要上線的最近點/秒,而不是僅僅的距離,但是這將打開另一個選項是你是否需要就近協調沿着一條線或實際的交點具有比任何現有頂點更近的線段。如果您的頂點足夠密集以至於差異無關緊要,請在sp包中使用spDistsN1。你必須從集合中的每一行提取所有座標(不難,但有點難看),然後遍歷每個感興趣的點計算距離線頂點的距離 - 然後你可以找到哪個是最短的,並選擇從頂點集合中協調,所以你可以很容易地獲得距離和座標。因爲該函數可以使用橢圓距離和longlat = TRUE的參數,所以不需要進行投影。

library(maptools) 

## simple global data set, which we coerce to Lines 
data(wrld_simpl) 

wrld_lines <- as(wrld_simpl, "SpatialLinesDataFrame") 

## get every coordinate as a simple matrix (scary but quick) 
wrld_coords <- do.call("rbind", lapply([email protected], function(x1) do.call("rbind", lapply([email protected], function(x2) [email protected][-nrow([email protected]), ])))) 

以交互方式檢查出來,您必須修改它以保存座標或最小距離。這將繪製線條並等待您點擊圖中的任何位置,然後它會從您的點擊到最接近的頂點上劃一條線。

## no out of bounds clicking . . . 
par(mar = c(0, 0, 0, 0), xaxs = "i", yaxs = "i") 

plot(wrld_lines, asp = "") 

n <- 5 

for (i in seq_len(n)) { 
xy <- matrix(unlist(locator(1)), ncol = 2) 
    all.dists <- spDistsN1(wrld_coords, xy, longlat = TRUE) 
    min.index <- which.min(all.dists) 
    points(xy, pch = "X") 
lines(rbind(xy, wrld_coords[min.index, , drop = FALSE]), col = "green", lwd = 2) 
} 
+0

尼斯和翔實的答案。感謝那。它看起來像'rgeos'沒有任何函數返回(或一個)最近的位置。你碰巧有一個功能的建議嗎?它看起來像'spatstat'中的'project2segment()'這樣做,但不幸的是它不使用'sp'類對象,這將更加容易... – 2012-03-20 06:57:41

+0

我會更新,我正在考慮更多,一些微妙之處。 – mdsumner 2012-03-20 07:08:10

+0

在spatstat和sp之間轉換並不難,在maptools中有強制工具,如spst < - as.psp(wrld_lines) – mdsumner 2012-03-20 12:02:05

3

geosphere封裝具有dist2line功能,這是否爲經度/緯度數據。它可以使用Spatial *對象或矩陣。

line <- rbind(c(-180,-20), c(-150,-10), c(-140,55), c(10, 0), c(-140,-60)) 
pnts <- rbind(c(-170,0), c(-75,0), c(-70,-10), c(-80,20), c(-100,-50), 
     c(-100,-60), c(-100,-40), c(-100,-20), c(-100,-10), c(-100,0)) 
d <- dist2Line(pnts, line) 
d 

插圖結果

plot(makeLine(line), type='l') 
points(line) 
points(pnts, col='blue', pch=20) 
points(d[,2], d[,3], col='red', pch='x') 
for (i in 1:nrow(d)) lines(gcIntermediate(pnts[i,], d[i,2:3], 10), lwd=2) 
相關問題