2017-09-14 95 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 
相关问题