3

I am trying to check the error that is introduced when you compute the distance of two points on earth with the euclidean distance instead of using the great circle distance (gcd). I have two points that are defined by their lattitude and longtitude. I used the python geopy framework for the great circle distance. Here the code for the gcd:

def measure(self, a, b):
        a, b = Point(a), Point(b)

        lat1, lng1 = radians(degrees=a.latitude), radians(degrees=a.longitude)
        lat2, lng2 = radians(degrees=b.latitude), radians(degrees=b.longitude)

        sin_lat1, cos_lat1 = sin(lat1), cos(lat1)
        sin_lat2, cos_lat2 = sin(lat2), cos(lat2)

        delta_lng = lng2 - lng1
        cos_delta_lng, sin_delta_lng = cos(delta_lng), sin(delta_lng)

        d = atan2(sqrt((cos_lat2 * sin_delta_lng) ** 2 +
                       (cos_lat1 * sin_lat2 -
                        sin_lat1 * cos_lat2 * cos_delta_lng) ** 2),
                  sin_lat1 * sin_lat2 + cos_lat1 * cos_lat2 * cos_delta_lng)

        return self.RADIUS * d

So or two points:

p1=[39.8616,-75.0748], p2=[-7.30933,112.76]

the

gcd = 78.8433004543197 klm

using the great_circle(p1,p2).kilometers function from geopy

I then transformed these two points in cartesian coordinates using this formula:

  def spherical_to_cartesian(r,la,lo):
       x=r*np.sin(90-la)*np.cos(lo)
       y=r*np.sin(90-la)*np.sin(lo)
       z=r*np.cos(90-la)
       return (x,y,z)

where r=6372.795, which results in the following cartesians coordinates

p1=[ -765.81579368,  -256.69640558,  6321.40405587], 
p2=[480.8302149,-168.64726394,-6352.39140142]

Then by typing: np.linalg.norm(p2-p1) i am getting 1103.4963114787836 as their euclidean norm which doesn't seem reasonable compared with ~78klm from the gcd. Am i inffering sth wrong?

curious
  • 1,524
  • 6
  • 21
  • 45
  • Are you sure you're comparing in the correct units? – jonrsharpe Mar 21 '14 at 15:22
  • 1
    not sure if this is your only error but np.sin and np.cos operate in radians... – John Greenall Mar 21 '14 at 15:25
  • The gdc seems very small. You have two points many degrees apart. I think 1 degree difference in position would be 100 km. – Salix alba Mar 21 '14 at 15:33
  • you are right!The problem was the way i am computing the cartesians. I should first convert to radians.Now i am getting 78.842797622886096 which is very close to gcd – curious Mar 21 '14 at 15:35
  • From a [NJ turnpike tollbooth](https://www.google.com/maps/@39.861,-75.074,435m/data=!3m1!1e3) to somewhere in [Surabaya, Indonesia](https://www.google.com/maps/@-7.30933,112.76,435m/data=!3m1!1e3) if you were wondering... – Nick T Mar 21 '14 at 17:37

1 Answers1

1

Python includes two functions in the math package; radians converts degrees to radians, and degrees converts radians to degrees.

The method sin() returns the sine of x, in radians.

import math
def spherical_to_cartesian(r,la,lo):
   rlo = math.radians(lo)
   rla = math.radians(90-la)
   x=r*np.sin(rla)*np.cos(rlo)
   y=r*np.sin(rla)*np.sin(rlo)
   z=r*np.cos(rla)
   return (x,y,z)
sabbahillel
  • 4,357
  • 1
  • 19
  • 36