1

I am having trouble determining the distance (in miles) between a point and 2 geographic coordinates using lat/lon.

# line segment point 1
x1 = -75.9667128 # longitude
y1 =  41.66279222 # latitude 

# line segment point 2
x2 = -75.96248381 # longitude
y2 = 41.65800548 # latitude 

# point to measure the orthogonal distance to
x3 = -75.96017288 # latitude 
y3 = 41.67049662 # latitude 

d = abs((x2-x1)*(y1-y3) - (x1-x3)*(y2-y1)) / np.sqrt(np.square(x2-x1) + np.square(y2-y1))

This yields 0.010002193890447786

When I convert this from degrees to miles using a scaling factor of 69.2 (I found this online), I get a different answer than when using a ruler, measuring between the line and the point, on a GIS site to double check myself.

Am I missing something here? When I plot out something like:

plt.plot([x1,x1+0.010002193890447786],[y1,y1])

the length of that line matches up with the correct distance.

someds
  • 57
  • 6
  • You cannot use lat, lon like that. First of all, at ~40⁰ the distance between two meridians 1⁰ apart is ~70% of the distance between two parallels 1⁰ apart. [This answer](https://stackoverflow.com/a/5538413/2749397) could help you to clarify your problem. – gboffi Apr 29 '21 at 19:06
  • Thank you, gboffi. Where do I find the appropriate equation for the shortest distance between a point and a great circle? Some of my points will not be orthogonal with the line segments, so the great circle must extend out, as in the equation I posted above. Thanks! Edit: I just noticed your link. – someds Apr 29 '21 at 19:31
  • If your points are so close, you can simply use a rectangular projection to have the Cartesian coordinates of the points, then you can use *correctly* the cross product formula (that you have used in the question but using latitudes and longitudes) to get the requested distance. See the answer I'll post in a few minutes … – gboffi Apr 29 '21 at 19:36
  • Have you considered using Haversin's formula? – Petr Fořt Fru-Fru Apr 29 '21 at 19:38
  • gboffi, yes, all of my points will be within a mile of each other. – someds Apr 29 '21 at 19:39
  • Petr Fořt Fru-Fru, I use the haversine formula for points at greater distances. My real hang-up is finding the distance of the point closest to the great circle line. – someds Apr 29 '21 at 19:41

1 Answers1

0

Here 40000 is the circumference of the Earth, in km, and 40000/2π is the Earth's radius — I don't remember the measurement in miles by the heart,

In [221]: # your points, in degrees
     ...: deg1, deg2, deg3 = (-75.96671280, 41.66279222), (-75.96248381, 41.65800548), (-75.96017288, 41.67049662)
     ...: # in radians
     ...: rad1, rad2, rad3 = ((lon*np.pi/180, lat*np.pi/180) for lon, lat in (deg1, deg2, deg3))
     ...: # an approximation of the location of the centre of the three points
     ...: lon_ctr, lat_ctr = (sum(rads)/3 for rads in zip(rad1, rad2, rad3))
     ...: # the position of the points wrt the centre
     ...: drad1, drad2, drad3 = ((lon-lon_ctr, lat-lat_ctr) for lon, lat in (rad1, rad2, rad3))
     ...: # an approximation to the Earth radius, and the radius of the parallel
     ...: rkm_lat = 40000/2/np.pi ; rkm_lon = rkm_lat*np.cos(lat_ctr)
     ...: # the position of the points, in km , wrt the centre
     ...: km1, km2, km3 = ((rkm_lon*dlon, rkm_lat*dlat) for dlon, dlat in (drad1, drad2, drad3))
     ...: #
     ...: print(";  ".join("p%d = %9.5f km, %9.5f km"%(i, x, y) for i, (x, y) in enumerate((km1, km2,km3), 1)))
p1 =  -0.29796 km,  -0.10806 km;  p2 =   0.05307 km,  -0.63992 km;  p3 =   0.24489 km,   0.74798 km
gboffi
  • 22,939
  • 8
  • 54
  • 85
  • 1
    Thanks for the reply, gboffi. I'm working through this snippet. I get a name error: NameError: name 'rad_lonc' is not defined – someds Apr 29 '21 at 20:37
  • 1
    I figured it out. I changed the variable lon_ctr to rad_lonc. – someds Apr 29 '21 at 20:41
  • I changed the name of a variable but I didn't edit all the occurrences, now I have checked again and it should work as it's now posted, after my last edit. – gboffi Apr 29 '21 at 20:42