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

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

```

```
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

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
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/
```
Associated Topics:
College Higher-Dimensional Geometry
College Trigonometry
High School Higher-Dimensional Geometry
High School Trigonometry

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