Planar Approximation: Latitude and Longitude

Date: 04/18/2003 at 11:01:48
From: Jessica
Subject: When can two points (lat and long) be treated as planar

How can I calculate the distance between two points below which they 
can be treated as if they were in a plane rather than on a sphere?  

I am trying to calculate the midpoint between cases of legionella and 
their nearest neighbor.  The cases are usually within 15 miles, but 
how can I calculate the maximum distance at which I need to treat the 
two cases as if they were on a sphere (i.e. account for the curvature 
of the earth)? 

I've seen many people mention that below 20 km or below 10 miles the 
earth's curvature doesn't matter - how do they calculate this?

Date: 04/18/2003 at 11:38:18
From: Doctor Rick
Subject: Re: When can two points (lat and long) be treated as planar

Hi, Jessica.

It's a complicated matter. How much error the planar approximation 
introduces in the distance between points depends not only on the 
distance but on the latitude and on the direction between the points.

For instance, if one point is on the equator, distances north and 
south, and distances east and west, are both exact regardless of the 
distance (ignoring effects of nonsphericity). Errors do enter when 
the direction is at an angle.

At the poles, the error is extreme. In fact, if one point is at a 
pole, the correct formula for distance is (90-lat)*pi*R/180 where lat 
is the latitude of the second point and R is the radius of the earth. 
This is very different from the Pythagorean formula!

If some distances will be great enough to require the spherical 
formula, why not just use it all the time? I know it takes more trig 
calculations, and if one point is fixed, the one trig function in the 
planar calculation can be pre-calculated, so there may be very good 
reasons for using the planar calculation whenever possible. But is 
that the case here?

I would use the Haversine formula for distance when distances may be 
quite small:

   Deriving the Haversine Formula
   http://mathforum.org/library/drmath/view/51879.html 

If you'd like, I could generate a table showing the errors for 
various distances and directions from a given point. Let me know if 
this will help, and tell me what latitude to use for the reference 
point.

- Doctor Rick, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 04/18/2003 at 12:41:12
From: Jessica
Subject: When can two points (lat & long) be treated as planar

Thanks for your quick response Dr. Rick!  You've convinced me to use 
the Haversine formula. 

I am programming this in SAS and need to calculate the distance 1.3 
billion times in the program (this simulation is run every day too), 
so efficiency in calculations will factor in.  All of the cases are 
in NYC (Lat range is 40.5 to ~41 dd and Long is -73.7 to -74.3 dd).  
I think an error table will be very useful; thank you for the 
offer.   

I have been using the law of cosines for spherical trig formula to 
calculate the distance between two points, but will replace it with 
the Haversine formula for all distance calculations.  

After I've found each case's two closest neighbors using the Haversine 
formula, I need to calculate the midpoint.  If the distance is small 
(i.e. less than 20 km) can I use a simple average of the two lats or 
longs (see below)?  

    midptlong = ((x1 + x2)/2)
	midptlat = ((y1 + y2)/2)  

How would I calculate the midpoint if the distance were greater than 
20 km?  

Many thanks again!
Jessica

Date: 04/18/2003 at 20:36:17
From: Doctor Rick
Subject: Re: When can two points (lat and long) be treated as planar

Hi, Jessica.

See:

   Latitude and Longitude of a Point Halfway between Two Points 
   http://mathforum.org/library/drmath/view/51822.html

You just need the formulas for dlon3 and lat3; to find the longitude 
of the midpoint, use

  lon3 = lon1 + dlon3

Since you need an efficient algorithm and you are (in the first 
phase) only interested in *comparing* distances, I suggest that you 
only calculate 'a' in the Haversine algorithm. It is the square of 
half the chord length between the two points, so of two choices for 
point B, the one that has a smaller 'a' value is closer to point A. 
This saves you two square roots and an atan2 for each point in your 
search. Of course, you can pre-calculate cos(lat1) for the given 
point A, leaving only three trig calculations per point in the search.

For the error table: I set up a spreadsheet to compute a table of 
relative errors due to the planar approximation for points at various 
distances and bearings. The table has a column for each distance (in 
kilometers) and a row for each bearing (in degrees from north), from a 
point at a given latitude (41 degrees). For each distance and bearing, 
I computed the latitude and longitude of the second point using the 
exact spherical formula:

LAT, LON GIVEN DISTANCE (d*R) AND AZIMUTH (tc)
lat2=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
IF (cos(lat)=0)
  lon2=lon1      // endpoint a pole
ELSE
  lon2=mod(lon1-asin(sin(tc)*sin(d)/cos(lat1))+pi,2*pi)-pi
ENDIF

Then I computed the distance using the approximate formula:

PLANAR APPROXIMATION
x = (lon2-lon1)*cos(lat1)*pi/180
y = (lat2-lat1)*pi/180
d = R * sqrt(x^2 + y^2)

The distances I found were:

CALCULATED DISTANCE BY PLANAR APPROXIMATION
     distance (km)    
     1            2            5            10           50
  0  1            2            5            10           50
 15  1.000004416  2.000017668  5.000110532  10.00044283  50.01121206
 30  1.000014774  2.000059113  5.000369752  10.00148097  50.03741972
 45  1.000024125  2.000096516  5.000603557  10.00241644  50.06085552
 60  1.000025585  2.000102347  5.000639801  10.00256008  50.06417495
 75  1.000016472  2.000065881  5.000411598  10.00164534  50.0409209
 90  0.999999996  1.999999969  4.999999515   9.999996121 49.99951521
105  0.999983524  1.999934086  4.999587879   9.998350475 49.95855614
120  0.999974418  1.999897681  4.999360636   9.997443424 49.9362629
135  0.999975884  1.999903555  4.999397548   9.997592399 49.94024953
150  0.999985233  1.999940949  4.999631225   9.998526849 49.96355763
165  0.999995587  1.999982354  4.999889816   9.99955996  49.98913625
180  1            2            5            10           50

Finally I computed the relative error, that is, 

  (approx. distance - actual distance)/(actual distance)

RELATIVE ERROR
     distance (km)    
     1            2            5           10           50
  0  1.40554E-13   5.3535E-13   1.4051E-13  -1.74083E-14  1.42109E-14
 15  4.41569E-06   8.83417E-06  2.21064E-05  4.42826E-05  0.000224241
 30  1.47744E-05   2.95567E-05  7.39505E-05  0.000148097  0.000748394
 45  2.41246E-05   4.8258E-05   0.000120711  0.000241644  0.00121711
 60  2.55851E-05   5.11736E-05  0.00012796   0.000256008  0.001283499
 75  1.64723E-05   3.29404E-05  8.23196E-05  0.000164534  0.000818418
 90 -3.87847E-09  -1.55143E-08 -9.6964E-08  -3.87854E-07 -9.69588E-06
105 -1.64765E-05  -3.29571E-05 -8.24242E-05 -0.000164952 -0.000828877
120 -2.55816E-05  -5.11596E-05 -0.000127873 -0.000255658 -0.001274742
135 -2.41158E-05  -4.82227E-05 -0.00012049  -0.00024076  -0.001195009
150 -1.47666E-05  -2.95254E-05 -7.3755E-05  -0.000147315 -0.000728847
165 -4.41291E-06  -8.82303E-06 -2.20367E-05 -4.4004E-05  -0.000217275
180  1.40554E-13   1.40554E-13  1.4051E-13  -1.74083E-14  1.42109E-14

You can see that, at latitude 41 degrees, the greatest error at 10 km 
is 0.0256%, or 2.56 m. Not bad, I'd say. At 50 km the greatest error 
is 0.128%, or 64 m. That's better than I expected, frankly. At 100 km 
the error is still under 258 meters (0.26%).

- Doctor Rick, The Math Forum
  http://mathforum.org/dr.math/ 

Date: 04/23/2003 at 10:13:18
From: Jessica
Subject: Thank you (When can two points (lat and long) be treated as 
planar)

Many thanks to you!  That was brilliant and very useful!

One final question - why is the cos(lat1) term added to the x term?  

PLANAR APPROXIMATION
x = (lon2-lon1)*cos(lat1)*pi/180
y = (lat2-lat1)*pi/180
d = R * sqrt(x^2 + y^2)

Thanks again,
Jessica

Date: 04/23/2003 at 10:43:08
From: Doctor Rick
Subject: Re: When can two points (lat and long) be treated as planar

Hi, Jessica.

The lines of longitude get closer together as they head toward the 
poles. Therefore a degree of longitude corresponds to a shorter 
distance as the latitude increases. 

To make this quantitative, consider a line of latitude, which is a 
circle. Looking at a side view of the earth (a cross-section along a 
line of longitude), we see that the radius of the latitude circle, r, 
and the radius of the earth, R, form a right triangle with an angle 
that is the latitude:

                       *********
                 ******    |    ******
              ***          +----------***
            **             |     r      /**
          **               |           / | **
        **                 |          /  |   **
       *                   |         /   |     *
      *                    |        /    |      *
     *                     |       /     |       *
    *                      |      /R     |        *
   *                       |     /       |         *
  *                        |    /        |          *
  *                        |   /         |          *
 *                         |  /          |           *
 *                         | /           |           *
 *                         |/ lat        |           *
 *-------------------------*-------------+-----------*
                                  r

Thus we see that cos(lat) = r/R, so r = R*cos(lat).

Now, the difference in longitude between two points at the same 
latitude is a central angle of the latitude circle. The distance 
between the points is the radius of the circle times the angle in 
radians. To convert the longitude difference to radians, we multiply 
it by pi/180. Then we multiply by r = R*cos(lat) to get the distance 
x:

  x = (lon2-lon1)*(pi/180)*R*cos(lat)

- Doctor Rick, The Math Forum
  http://mathforum.org/dr.math/