如何使用投影从 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 分配给数据集，您表示 "these numbers are metres, and the origin is not at (0,0) degrees (where equator meets Greenwich meridian) but at the origin of zone 18 of the UTM system"。所以你得到以米为单位的距离。

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

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

How can I get correct distances (in meters) from st_distance (sp package) in R using a projection?

r

spatial

sf