Расстояние между набором точек и многоугольником с sf в R

У меня есть фрейм данных точек на карте и интересующая область, описанная как многоугольник точек. Я хочу рассчитать расстояние между каждой из точек многоугольника, в идеале используя пакет sf.

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. просто использовать st_distance веселье - не работает, неправильный ответ

st_distance(pnts_sf, area)

2. В вызове изменения - все неправильные ответы

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, поэтому мне интересно, что здесь не так. Насколько я могу судить, первые 2 подхода так или иначе заканчиваются рассмотрением точек «в целом» (одна из точек находится внутри многоугольника области, поэтому я предполагаю, что один из неправильных ответов - 0). Третий подход предусматривает рассмотрение момента за раз, что является моим намерением.
Есть идеи, как я могу заставить работать вызов mutate?

Я на R 3.4.1 с

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

r sf spatial
person davidski    schedule 14.09.2017    source источник
comment
Вам нужно установить by_element = TRUE в вызове st_distance()? Также +1 для воспроизводимого примера   -  person Phil    schedule 14.09.2017
comment
@Phil aghh Я понял - это очень-очень-очень ошибка новичка! Я получил _1 _ / _ 2_ баллов из другого источника, который использует epsg 4326 crs. Эта часть ускользнула от моего удивления. Итак, после создания фрейма данных, преобразование в объект sf должно быть сначала pnts_sf <- st_as_sf(pnts, crs = 4326L, coords = c("long", "lat")) (потому что это была моя начальная система координат!), А затем еще один вызов для преобразования в тот же crs, что и у area - pnts_sf <- st_transform(crs = 7415L). Тогда вызов st_distance() дает правильные результаты. Мораль истории = всегда следите за своим crs!   -  person davidski    schedule 14.09.2017
comment
Очень полезная мораль. Хотите написать свое решение в виде ответа, чтобы люди увидели, что оно решено?   -  person Phil    schedule 14.09.2017
comment
да, ты прав, я упакую это в ответ. спасибо за помощь, ваш комментарий побудил меня обнаружить ошибку.   -  person davidski    schedule 14.09.2017

Ответы (1)


arrow_upward
5
arrow_downward

Получается, что вся путаница возникла из-за моей глупой оплошности. Вот разбивка:

  • Фрейм данных 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 transformation to another crs
pnts_sf <- st_transform(pnts_sf, crs = 7415L)

st_distance(pnts_sf, area)

--------------------------
Units: m
          [,1]
[1,] 3998.5701
[2,]    0.0000
[3,]  751.8097
person davidski    schedule 14.09.2017
comment
Когда я слежу за кодом, как я его вижу, я не вижу, что точка 2 находится внутри многоугольника. Возможно ли, что отсутствует ступенька? - person gregmacfarlane; 07.06.2019