Find the Cities

Date: 06/19/2001 at 12:03:18
From: Rod Doornbosch
Subject: Find All Cities (Using latitude and Longitude), x miles off a 
straight line between two cities

First of all I'd like to say WOW, I've never seen so many answers 
about calculating distance using latitude and longitude. I'll try to 
keep this simple. 

I have a database with thousands of cities (latitude and longitude). 
Given the latitude and longitude of Detroit, Michigan and the latitude 
and longitude of Miami, Florida, draw a straight line between the two 
and find every city (using latitude and longitude) within 10 miles of 
that line all the way from Detroit to Miami.

To keep it simple, let's assume the earth is flat. I only need to be 
able to do this for North America, and I'm hoping that I can use a 
computer SQL statement to get the cities I need. 

Looking forward to your answer.

Date: 06/21/2001 at 08:56:30
From: Doctor Rick
Subject: Re: Find All Cities (Using latitude and Longitude), x miles 
off a straight line between two cities

Hi, Rod. We keep getting new variations on the theme of latitude-
longitude calculations - it's a popular topic. Thanks for your 
addition.

Assume the earth is flat? Over the 10 miles from the "straight line" 
we can safely make this assumption, but in order to find the "straight 
line" we must take the curvature of the earth into account. The 
"straight line" is really a great circle. I can't think how to do such 
a hybrid calculation. What follows is fully based on a spherical 
earth.

I don't know what is available to you in SQL, but the straightforward 
solution to your problem is to calculate, for each city in the 
database, how far it falls from the great circle connecting Detroit 
and Miami. Then you could find the cities whose distances from the 
line are less than 10 miles.

Here is a way to find the distance from the line. You have spherical 
coordinates of points A and B. Convert these to cartesian coordinates, 
regard them as vectors, and take the cross product, vector D. Every 
point on the great circle makes a 90-degree (pi/2 radians) angle with 
D. 

For a given city C, take the dot product of vector D with vector C, 
and from this obtain the angle between C and D (in radians). Take the 
difference from pi/2 radians, multiply by the radius of the earth to 
convert it to miles. This is the distance from D to the great circle 
between A and B.

Here's how it works out in detail. We have coordinates

  A: lat1, lon1
  B: lat2, lon2
  C: lat3, lon3

I will work with cartesian coordinates such that the x-y plane passes 
through the equator and the x-z plane passes through A. I will also 
take the radius of the earth to be 1 for convenience, until the last 
step. Then

  A = (cos(lat1), 0, sin(lat1))
  B = (cos(lat2)*cos(dlon2), cos(lat2)*sin(dlon2), sin(lat2))
  C = (cos(lat3)*cos(dlon3), cos(lat3)*sin(dlon3), sin(lat3))

where

  dlon2 = lon2 - lon1
  dlon3 = lon3 - lon1

Take the cross product of A and B:

  D = A x B = (-sin(lat1)*cos(lat2)*cos(dlon2),
               sin(lat1)*cos(lat2)*cos(dlon2)-cos(lat1)*sin(lat2),
               cos(lat1)*cos(lat2)*sin(dlon2))

Then we take the dot product of C with D:

  C . D = -sin(lat1)*cos(lat2)*cos(dlon2)*cos(lat3)*cos(dlon3) +
          (sin(lat1)*cos(lat2)*cos(dlon2)-cos(lat1)*sin(lat2))*
          cos(lat3)*sin(dlon3) +
          cos(lat1)*cos(lat2)*sin(dlon2)*sin(lat3)

The dot product is also

  C . D = |C|*|D|*cos(theta)

where theta is the angle between C and D. Since |C| = |D| = 1, we have

  theta = arccos(C . D)

and the distance of C from the great circle is

  d = R*|arcsin(C . D)|

The dot product looks pretty messy, but it isn't as bad as it looks. 
You can pre-compute the parts of the expression that will not change 
from city to city. Thus

  a = -sin(lat1)*cos(lat2)*cos(dlon2)
  b = sin(lat1)*cos(lat2)*cos(dlon2)-cos(lat1)*sin(lat2)
  c = cos(lat1)*cos(lat2)*sin(dlon2)

  dotprod = (a*cos(dlon3) + b*sin(dlon3))*cos(lat3) + c*sin(lat3)

Only the last calculation must be recomputed for each city in the 
database. Furthermore, rather than compute d for each city and 
compare d to 10 miles, you can compute (just once)

  dotmax = sin(m/R)
  dotmin = -dotmax

where m is the maximum distance (10 miles in your case). For each 
city, check whether dotprod is between dotmin and dotmax.

Just be sure to convert the latitudes and longitudes from degrees to 
radians if your trig functions expect angles in radians: multiply each 
by pi/180. If the trig functions take angles in degrees, you must 
convert (pi/2 +or- d/R) to degrees (multiply by 180/pi) before taking 
the cosine.

I haven't tested this, so let me know if it works.

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