Associated Topics || Dr. Math Home || Search Dr. Math

### Points within a Given Distance on a Sphere

```Date: 07/25/2004 at 13:09:45
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

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
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
"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.

```

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

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
Subject: Points within a given distance on a sphere

Dr Rick,

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

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.

```
Associated Topics:
College Coordinate Plane Geometry

Search the Dr. Math Library:

 Find items containing (put spaces between keywords):   Click only once for faster results: [ Choose "whole words" when searching for a word like age.] all keywords, in any order at least one, that exact phrase parts of words whole words

Submit your own question to Dr. Math
Math Forum Home || Math Library || Quick Reference || Math Forum Search