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