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/
Search the Dr. Math Library:
Ask Dr. MathTM
© 1994- The Math Forum at NCTM. All rights reserved.