CODEWITHSUNDEEP
rspatialr-sf
Diego Uriarte
Diego Uriarte

Reputation: 13

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

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

Upvotes: 1

Views: 6122

Answers (1)

Spacedman
Spacedman

Reputation: 94182

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).

Upvotes: 3

Related Questions