2017-09-14 57 views
1

我有一個地圖上的點的數據框和感興趣的區域描述爲點的多邊形。我想計算每個點與多邊形之間的距離,理想情況下使用sf包。一組點與一個多邊形之間的距離R在R

library("tidyverse") 
library("sf") 

# area of interest 
area <- 
    "POLYGON ((121863.900623145 486546.136633659, 121830.369032584 486624.24942906, 121742.202408334 486680.476675484, 121626.493982203 486692.384434804, 121415.359596921 486693.816446951, 121116.219703244 486773.748535465, 120965.69439283 486674.642759986, 121168.798757601 486495.217550029, 121542.879304342 486414.780364836, 121870.487595417 486512.71203006, 121863.900623145 486546.136633659))" 

# convert to sf and project on a projected coord system 
area <- st_as_sfc(area, crs = 7415L) 

# points with long/lat coords 
pnts <- 
    data.frame(
    id = 1:3, 
    long = c(4.85558, 4.89904, 4.91073), 
    lat = c(52.39707, 52.36612, 52.36255) 
    ) 

# convert to sf with the same crs 
pnts_sf <- st_as_sf(pnts, crs = 7415L, coords = c("long", "lat")) 

# check if crs are equal 
all.equal(st_crs(pnts_sf),st_crs(area)) 

我想知道爲什麼下面的方法不給我正確的答案。

1.Simply使用st_distance樂趣doesn't工作,錯誤的答案

st_distance(pnts_sf, area) 

2.In一個發生變異電話 - 所有錯誤的答案

pnts_sf %>% 
    mutate(
    distance = st_distance(area, by_element = TRUE), 
    distance2 = st_distance(area, by_element = FALSE), 
    distance3 = st_distance(geometry, area, by_element = TRUE) 
) 

但是這種方法似乎工作和給出正確的距離。

3. map在長期/緯度 - 正常工作

pnts_geoms <- 
    map2(
    pnts$long, 
    pnts$lat, 
    ~ st_sfc(st_point(c(.x, .y)) , crs = 4326L) 
) %>% 
    map(st_transform, crs = 7415L) 

map_dbl(pnts_geoms, st_distance, y = area) 

我是新來的空間數據,我努力學習sf封裝所以我不知道是怎麼回事錯在這裏。據我所知,前兩種方法以某種方式最終考慮點「整體」(其中一點是在區域多邊形內,所以我想這就是爲什麼一個錯誤答案是0)。第三種方法是在我的意圖下考慮一個時間點。
任何想法我怎樣才能讓mutate電話工作以及?

我在R 3.4.1與

> packageVersion("dplyr") 
[1] ‘0.7.3’ 
> packageVersion("sf") 
[1] ‘0.5.5’ 
+1

你需要設置'by_element = TRUE;在'st_distance()'調用?另外,+1可重現的例子 – Phil

+0

@Phil aghh我想通了 - 這是一個非常非常新秀的錯誤!我通過使用「epsg 4326」crs的另一個來源獲得了「長」/「拉特」。這部分逃脫了我的想法。所以在創建數據框之後,對'sf'對象的轉換應該是'pnts_sf < - st_as_sf(pnts,crs = 4326L,coords = c(「long」,「lat」))'first(因爲那是我的初始座標系!),然後另一個調用轉換爲與'area'相同的crs具有''pnts_sf < - st_transform(crs = 7415L)'。然後'st_distance()'調用產生正確的結果。故事的道德= *總是*保持跟蹤你的crs! – davidski

+1

一個非常有用的道德。注意將解決方案寫成答案,以便人們可以看到解決方案? – Phil

回答

1

因此,原來整個混亂是由我的一個小疏忽愚蠢造成的。這裏的故障:

  • points數據幀來自不同的來源比area多邊形(!)。
  • 監督我一直試圖將它們設置爲crs 7415這是一個合法但不正確的舉動,並最終導致錯誤的答案。
  • 正確的做法是將它們轉換爲sf對象,它們源自crs,將它們轉換爲area對象所在的對象,然後繼續計算距離。

全部放在一起:

# this part was wrong, crs was supposed to be the one they were 
# originally coded in 
pnts_sf <- st_as_sf(pnts, crs = 4326L, coords = c("long", "lat")) 

# then apply the transfromation to another crs 
pnts_sf <- st_transform(crs = 7415L) 

st_distance(pnts_sf, area) 

-------------------------- 
Units: m 
      [,1] 
[1,] 3998.5701 
[2,] 0.0000 
[3,] 751.8097 
相關問題