如何使用投影从R中的st_distance（sp包）获得正确的距离（以米为单位）？

Question

我想计算点到点之间的距离。我知道在R（see here for one example）中有几种方法可以做到这一点，我认为最好使用sf包中的st_distance函数，但是当我使用与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”）时，我确实得到了正确的结果（以米为单位）：

 [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

Answer 1

EPSG 32718是一个以米为单位的笛卡尔坐标参考系统。通过将CRS分配给数据集，您说“这些数字是米，而原点不是（0,0）度（赤道遇到格林威治子午线），而是在UTM系统的区域18的原点”。所以你得到一米的距离。

EPSG 4326是一种长距离参考系统，具有特殊形状的椭球形。坐标是纬度。 st_distance发现了这一点，并根据椭球计算了点之间的大圆距离。如果你想要十进制度数的距离，那么分配一个NA CRS，你将获得无单位距离，这是lat-long中的毕达哥拉斯距离（例如，在极点附近的实际上非常错误）。