Как я могу получить правильные расстояния (в метрах) от st_distance (sp package) в R, используя проекцию? - PullRequest
Новая
       39

Как я могу получить правильные расстояния (в метрах) от st_distance (sp package) в R, используя проекцию?

1 голос
/ 13 апреля 2019

Я хочу рассчитать расстояние между точками.Я знаю, что есть несколько способов сделать это в R ( см. Здесь для одного примера ), я подумал, что было бы лучше использовать функцию st_distance из пакета sf, но когда я использую проекцию, отличную от WGS84(crs = 4326), я получаю расстояния в десятичных градусах, а не в метрах.

Однако, когда я устанавливаю проекцию crs = 32718, я получаю расстояние в десятичных градусах.Есть ли способ преобразовать это в метры (или получить метры в первую очередь).Я не понимаю, почему, когда я устанавливаю проекцию crs = 4326, я получаю расстояние в метрах.

Я включил воспроизводимый пример: 

library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(tidyverse)
library(maptools)
#> Loading required package: sp
#> Checking rgeos availability: TRUE


crs <- CRS("+init=epsg:32718")

df <- tibble::tribble(
  ~documento, ~cod_mod, ~nlat_ie, ~nlong_ie,
  "00004612", 238840, -8.37661, -74.53749,
  "00027439", 238758, -8.47195, -74.80497,
  "00074909", 502518, -8.83271, -75.21418,
  "00074909", 612663, -8.82781, -75.05055,
  "00074909", 612812, -8.64173, -74.96442,
  "00102408", 237255, -13.4924, -72.9337,
  "00102408", 283341, -13.5317, -73.6769,
  "00109023", 238717, -9.03639, -75.50947,
  "00109023", 238840, -8.37661, -74.53749,
  "00109023", 1122464, -8.37855, -74.57039,
  "00124708", 238717, -9.03639, -75.50947,
  "00124708", 238840, -8.37661, -74.53749,
  "00124708", 1122464, -8.37855, -74.57039,
  "00186987", 612663, -8.82781, -75.05055,
  "00186987", 1121383, -8.36195, -74.57805,
  "00237970", 327379, -3.55858, -80.45579,
  "00238125", 1137678, -3.6532, -80.4266,
  "00238125", 1143577, -3.50163, -80.27616,
  "00239334", 1143577, -3.50163, -80.27616,
  "00239334", 1372333, -3.6914, -80.2521
)

df_spatial <- df

coordinates(df_spatial) <- c("nlong_ie", "nlat_ie")

proj4string(df_spatial) <- crs

# Now we create a spatial dataframe with coordinates in the average location of each documento

df_mean_location <- df %>%
  group_by(documento) %>%
  summarize(
    mean_long = mean(nlong_ie),
    mean_lat = mean(nlat_ie)
  )

df_mean_location_spatial <- df_mean_location

coordinates(df_mean_location_spatial) <- c("mean_long", "mean_lat")
proj4string(df_mean_location_spatial) <- crs

df_spatial_st <- st_as_sf(df_spatial)
df_mean_location_spatial_st <- st_as_sf(df_mean_location_spatial)

distancias1 <- st_distance(df_spatial_st, df_mean_location_spatial_st, by_element = TRUE)

distancias1
#> Units: [m]
#>  [1] 0.00000000 0.00000000 0.15248325 4.99880005 0.10219044 5.26515886
#>  [7] 5.06614947 7.38054767 7.53880558 7.43549151 1.17475732 0.28396349
#> [13] 0.63815871 4.99880005 0.37683694 7.52071866 7.47784143 0.18844161
#> [19] 0.10677741 0.09564457

Когда я меняю crs <- CRS ("+ init = epsg: 4326"), я получаю правильные результаты (в метрах): </p> 

 [1]      0.00      0.00  16792.18 552085.93  11258.44 581428.01 560043.61 816269.42 834131.40 822686.13 129481.67  31286.98  70373.13 552085.93
[15]  41565.46 832000.85 827230.50  20928.56  11835.41  10577.04

1 Ответ

1 голос
/ 13 апреля 2019

EPSG 32718 - декартова система координат в метрах.Присваивая этот CRS набору данных, вы говорите, что «эти числа являются метрами, а начало координат не в (0,0) градусах (где экватор встречается по гринвичскому меридиану), а в начале зоны 18 системы UTM».Таким образом, вы получите расстояние в метрах.

EPSG 4326 - это система координат длинной широты с определенной формой эллипсоидальной земли.Координаты - это широта в градусах.st_distance определяет это и вычисляет расстояние по большому кругу между точками на основе эллипсоида.Если вы хотите получить расстояние в десятичных градусах, тогда назначьте NA CRS, и вы получите безразмерные расстояния, которые являются пифагорейскими расстояниями в широте (например, очень неверно в реальном выражении около полюсов).

...