1

By this website http://www.movable-type.co.uk/scripts/latlong.html,

I would like to calculate distances between latitude 37.919 and 37.465, longitude 15.358 and 14.601, the results are 50.48 km and 66.81 km, respectively.

e.g. Point1: 37.919, 14.601 and Point2: 37.465, 14.601, the website gives 50.48 km,

Point1: 37.465, 15.358 and Point2: 37.465, 14.601, I got the distance among them, 66.81 km.

However, I use the pyproj.Proj to calculate the distances, which differs greatly from the results above.

Below is my python code:

import pyproj.Proj
p=Proj(proj='eqc',lat_0=37.692, lon_0=14.980)
p(14.601,37.465) #return (-42190.08701065082, -25269.524410073795)
p(15.358,37.919) #return (42078.76751985713, 25269.52441007167)

The distance between 37.919 and 37.465 is 25269.52441007167-(-25269.524410073795), roughly 50539 m = 50.539 km, which matches approximately the value by the website: 50.44 km

But, the distance between 15.358 and 14.601 is roughly 84269 m, 84.269 km, which differs greatly from the one by the website: 66.81 km.

I don't know why?

TomazicM
  • 25,601
  • 22
  • 29
  • 39
Makoto
  • 11
  • 1
  • Why do you think that an Equirectangular projection will give you correct distances? – Devdatta Tengshe Aug 18 '14 at 09:48
  • @Devdatta Tenshe: Because the google earth uses Equirectangular projection at local scale, so I think the pyproj should also give the same result. What am I missing? – Makoto Aug 19 '14 at 01:07
  • For example, I got the same result for latitude 37.919 and 37.465, but for longitude 15.358 and 14.601 the pyproj gives 84.269 km while the website gives 66.81 km. What is wrong with them? – Makoto Aug 19 '14 at 01:17

1 Answers1

1

According to the website, the first box is calculating the distance with the Great-circle distance between two point. This is the geodesic calculation in PROJ (https://proj.org/geodesic.html).

To perform this calculation in pyproj, you can use pyproj.Geod.inv(). The distances are in meters.

>>> from pyproj import Geod
>>> gg = Geod(ellps="WGS84")
>>> gg.inv(14.601,37.919, 14.601,37.465)[-1]
50389.758721051636  # website: 50.44 km
>>> gg.inv(15.358,37.465, 14.601,37.465)[-1]
66969.11203130838  # website: 66.81 km

Depending on the ellipsoid being used by the website, you can adjust the parameters used to create the pyproj.Geod class.

snowman2
  • 7,321
  • 12
  • 29
  • 54