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

Question

I want to calculate the distance between to points. I know there are several ways to do it in R (see here for one example), I thought it would be best to use the st_distance function from the sf package, but when I use a projection different to WGS84 (crs = 4326), I get the distances in decimal degrees and not in meters.

However, when I set the projection to crs = 32718, I get the distance in decimal degrees. Is there a way to convert this to meters (or to get meters in the first place). What I don't understand is why when I set the projection to crs = 4326, I do get the distance in meters.

I included a reproducible example:

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

When I change the crs <- CRS("+init=epsg:4326"), I do get the correct results (in meters):

 [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

Is there any difference in your expression for distancias1 and distancias2? They look identical. — Spacedman
Also st_distance is from the sf package, and most of the other spatial functions you use are coming from the sp package which is attached with maptools - you probably don't need maptools here (or tidyverse) to generate a minimal example. — Spacedman
@Spacedman you are right, I left a variable (distancias2) that was not needed, so I edited the question. I'll take into account your comments next time to only use the packages strictly needed — Diego Uriarte
Cool, also there's the GIS stack site where a lot of spatial Qs get asked: gis.stackexchange.com/questions/tagged/r — Spacedman

Spacedman Spacedman · Accepted Answer · 2019-04-13T15:18:49

EPSG 32718 is a cartesian coordinate reference system in metres. By assigning that CRS to a data set, you are saying "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". So you get a distance in metres.

EPSG 4326 is a lat-long reference system with a particular shape of ellipsoid earth. The coordinates are lat-long degrees. st_distance spots this and works out the great circle distance between points based on the ellipsoid. If you want the distance in decimal degrees then assign an NA CRS and you'll get unitless distances, which are the pythagorean distances in lat-long (and so very wrong in real terms near the poles, for example).

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

1 Answers