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’
你需要設置'by_element = TRUE;在'st_distance()'調用?另外,+1可重現的例子 – Phil
@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
一個非常有用的道德。注意將解決方案寫成答案,以便人們可以看到解決方案? – Phil