1

I try to calculate the distance between two geogprahical points, using st_distance. The points are stored in a data frame, but for simplicity, I take one pair of points to show the problem.

Here are the two points (an example of a pair of points):

a <-st_sfc(st_point(c(70.321117, 22.896517)))

b <-st_sfc(st_point(c(68.743738549303, 17.0129893190695)))

st_distance(a,b)

Which returns 6.09.

I understand this should be a distance in degrees, and to calculate the same in km I need to multiply by 111 km. But the distance seems to be too large (ca 666 km), compared to online calculators of great circle distance (I got 293 km).

Why is the difference so large? I do I calculate the distances in a wrong way?

  • How many km "fit" in one degree varies with the latitude. At the equator your 111 km are about correct, but at the polar circle one degree is about 38 km. You could convert your data to EPSG 3857, which is based on m. – Erik Aug 11 '20 at 15:00
  • @Erik No, the last thing the OP wants is to compute distances in Web Mercator. That's actually worse than Cartesian degrees, because it has units of meters, but they are known to be highly inaccurate away from the Equator. – Vince Aug 11 '20 at 15:16
  • @Erik, thanks for the comment. Is there any reference explaining how many meters in 1 degrees are in sub-polar areas? – marina_esp Aug 11 '20 at 15:57
1

sf package is tightly coupled with unitspackage, together they do really neat things, as giving you the units of the output of a function; in your case, the objects lack CRS, once set, you'll get:

a <- st_sfc(st_point(c(70.321117, 22.896517))) %>% 
  st_set_crs(4326)

b <- st_sfc(st_point(c(68.743738549303, 17.0129893190695))) %>% 
  st_set_crs(4326)

st_distance(a,b)

Units: [m]
         [,1]
[1,] 671910.7

From the function help: "Compute Euclidian or great circle distance between pairs of geometries; compute, the area or the length of a set of geometries."; on the parameter which :" character; for Cartesian coordinates only: one of Euclidean, Hausdorff or Frechet; for geodetic coordinates, great circle distances are computed; see details"

1

I think your coordinates might be reversed?

When I initially copied your coordinates, they put me out in the Indian ocean, and the distance was almost 700KM. But I thought that might be wrong, so I reversed them, and I think it works.

Here's the SQL I used (referencing a blog post below as well):

;with cte_geog as (
select 

st_makepoint(22.896517, 70.321117)::geography as geog_a
, st_makepoint(17.0129893190695, 68.743738549303)::geography as geog_b

)

select 

st_distance(d.geog_a, d.geog_b)/1000 as distance_km
, st_makeline(d.geog_a::geometry, d.geog_b::geometry)::geography as geog_line
, *
from cte_geog as d

The distance returned is 289.15KM

Here's the result in the DBeaver spatial viewer:

enter image description here

Referencing this blog post from Paul Ramsey (https://info.crunchydata.com/blog/postgis-and-the-geography-type), as long as your working with Geography the ST_Distance() function will do all the work of calculating the distances on the sphere.

  • Great, I see that you got exactly the distance I need after swapping the coordinates. But I am still confused, if the distance should be 290 km or 670 km. @Elio Diaz in his approach below got 671 km. Was it a wrong distance? – marina_esp Aug 12 '20 at 10:17
  • @marina_esp I think the distance of 290KM is correct - the distance of 670Km is when the points are reversed and placed in the Indian Ocean... – DPSSpatial Aug 12 '20 at 14:57
  • Thank you, @DPSSpatiall, now everything is clear! – marina_esp Aug 14 '20 at 6:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.