Points within a Given Distance on a Sphere

Date: 07/25/2004 at 13:09:45
From: Adam
Subject: Points within a given distance on a sphere

Suppose I have a table of US cities (about 80,000 entries) with 
latitudes and longitudes, a latitude and longitude of City1, and a 
Distance1 in miles.  Now what I'd like to do is find all cities within 
Distance1 miles of City1.  My first pass at the solution was to 
compute Distance2 for all cities from City1 using the Great Circle 
algorithm.  If Distance2 < then Distance1 then I add the city to my 
list of cities within Distance1.  This takes a certain amount of 
computational time which I'd like to decrease as much as possible.

My thought was to do the following:

1) Compute the maximum longitude and latitude from Distance1.  That 
is, keep latitude constant and find the longitude that is Distance1 
miles from City1.  Then keep longitude constant and compute the 
latitude that is Distance1 miles from City1.  This should give me a 
square area (okay slightly curved since it's on a sphere). 

2) Scan my cities tables for all cities with latitudes and longitudes 
that fall within the maximum latitude and longitude computed in step 
1 and store those cities in a second table.  This search should be a 
simple comparison without the need to compute the Great Circle 
distance for every city.  Thus this search should be fairly quick.

3) Then I'd perform a search of my second smaller city list using the 
Great Circle algorithm to retain any cities that fall within 
Distance1, and discard any cities that are at the corners of the 
square.

How do you compute the latitudes and longitude in #1?  Are there any
pitfalls that are not evident?

Date: 07/25/2004 at 21:49:47
From: Doctor Rick
Subject: Re: Points within a given distance on a sphere

Hi, Adam.

This seems like a reasonable approach.  There are items in the 
archives that deal more directly with the questions you are asking. 
Check these out.  This first one gives the formulas for the maximum 
and minimum latitude and longitude of points a given distance from a 
given point:

  Maximum Difference, Longitude and Latitude
    http://mathforum.org/library/drmath/view/51804.html 

Then I recommend using the haversine formula for distance:

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

The formula given there is

  dlon = lon2 - lon1
  dlat = lat2 - lat1
  a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
  c = 2 * atan2(sqrt(a), sqrt(1-a)) 
  d = R * c

Your goal is to determine whether (lat2,lon2) is within distance D of 
(lat1,lon1).  Let's manipulate this formula a bit.  You want to know 
whether

  R*c < D
  c < D/R
  atan2(sqrt(a), sqrt(1-a)) < D/(2R)
  sqrt(a)/sqrt(1-a) < tan(D/(2R))
  a/(1-a) < tan(D/(2R))^2
  (1-a)/a < 1/tan(D/(2R))^2
  1/a - 1 < 1/tan(D/(2R))^2
  1/a < 1 + 1/tan(D/(2R))^2
  a > 1/(1 + 1/tan(D/(2R))^2

The right side can be computed once; I'll call this value A.  You just 
need to determine for each point 2 whether

  a > A
  (sin(dlat/2))^2 + cos(lat1)*cos(lat2)*(sin(dlon/2))^2 > A

In this formula, cos(lat1) can be computed just once.  For each point, 
you need to compute

  sin((lat2-lat1)/2)
  sin((lon2-lon1)/2)
  cos(lat2)

and put them in the formula above.  That's three trig function 
computations per point, plus some multiplications and additions.  Will 
this suffice?

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

Date: 07/26/2004 at 13:12:54
From: Adam
Subject: Points within a given distance on a sphere

Dr. Rick,

Thank you for the reply.  I understand how you manipulated the
Haversine formula to reduce the amount of computations needed to
determine whether Lat2, Lon2 is within or beyond distance D from
Lat1, Lon1, and therefore I believe I can convert this into the
software I'm designing.

On the other hand I got lost in the math in the article titled:
"Maximum Difference, Longitude and Latitude" which you referenced in
your reply.  Basically, the article mentions something about the
"maximum" longitude at the equator.  I'm not sure whether that gives 
me the lon and lat values I'm seeking.  I'll repeat my problem with
examples to illustrate. 

Suppose I have the Lat1 and Lon1 for Boston and a distance of 10
miles.  I'd like to compute a latitude which is 10 miles North of
Boston, the latitude which is 10 miles South of Boston, the longitude
which is 10 miles West of Boston, and the longitude 10 miles East of
Boston.  My guess is that I can then use those four lat and lon values
to compare against my table of ZipCodes, Lats, and Lons.

So, if I compare the Lon for 10 miles West of Boston against the Lon
for San Francisco I should find SF is beyond the Lon for 10 miles West
of Boston.  That means I can discard SF for consideration and this 
also means that instead of performing all the trig functions in any
variation of the distance formula (Great Circle, Haversine, etc) the
test for whether the city is within 10 miles of Boston reduces to a
simple comparison of integer values. 

The Lon for Cambridge would fall between the 10 miles West and 10
miles East of Boston lons.  Of course, I'd have to write my software 
to check whether Cambridge is within the 10 mile range North and South
along the latitude.  But this winds up that I'd keep Cambridge on the
list. 

Once I have compiled this second smaller list of zipcodes/cities that
are within a square area 10 miles wide and 10 miles high I would then
use the Haversine formula as you manipulated it to refine the list
further by discarding any zipcodes which are at the corners of the
square area. 

Again, the "maximum" article you referenced probably does provide the
answer I need but I just don't see it at the moment.

Adam

Date: 07/26/2004 at 14:26:57
From: Doctor Rick
Subject: Re: Points within a given distance on a sphere

Hi, Adam.

This is what I understood you to be seeking.  Yes, the article I 
referenced gives you the information you want.  You can just skip the 
explanation that confuses you, and use the formula for max_dlon.  This 
is the maximum difference in longitude between point A (at latitude 
lat1) and any point 1000 units away from point A.  Replace 1000 with 
the distance you want, and R with the radius of the earth in the same 
units.  Thus

  max_lon = lon1 + arcsin(sin(D/R)/cos(lat1))
  min_lon = lon1 - arcsin(sin(D/R)/cos(lat1))

where D is the distance you are interested in.  For latitude, it's 
simple:

  max_lat = lat1 + (180/pi)(D/R)
  min_lat = lat1 - (180/pi)(D/R)

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

Date: 07/26/2004 at 18:49:49
From: Adam
Subject: Points within a given distance on a sphere

Dr Rick,

Thank you for your kind reply and confirming that the Maximum/Minimum
Distance article contained the algorithms I needed.

Initially, I forgot to convert radians to degrees before doing the
min/max computations, but once I realized what the problem was and
added the degrees -> radians conversion the formula worked as expected.

With these latest changes my software performance has increased by a
factor of about 6.  The original time to cull out the zipcodes within 
a specified distance was about 0.2 seconds.  Now the time is 0.0156288
(Distance = 5 miles) to 0.0468736 (Distance = 40 miles).  Given these
changes, I can't think of anything else I can do to get more speed out
of the process.

Thank you for your help.

Adam