Python geocode filtering by distance

Question

I need to filter geocodes for proximity to a location. For example, I want to filter a restaurant’s geocode list to identify these restaurants within 10 miles of my current location.

Can someone point me to a function that converts distance to latitude and longitude? For example:

class GeoCode(object): """Simple class to store geocode as lat, lng attributes.""" def __init__(self, lat=0, lng=0, tag=None): self.lat = lat self.lng = lng self.tag = None def distance_to_deltas(geocode, max_distance): """Given a geocode and a distance, provides dlat, dlng such that |geocode.lat - dlat| <= max_distance |geocode.lng - dlng| <= max_distance """ # implementation # uses inverse Haversine, or other function? return dlat, dlng 

Note. I use the supremum rate for distance.

Answer 1

It seems that there was no good Python implementation. Fortunately, the SO "Related Articles" sidebar is our friend. This SO article points to a great article that gives math and Java implementation. The actual function you need is quite short and is built into my Python code below. Tested to the extent. Read the warnings in the comments.

 from math import sin, cos, asin, sqrt, degrees, radians Earth_radius_km = 6371.0 RADIUS = Earth_radius_km def haversine(angle_radians): return sin(angle_radians / 2.0) ** 2 def inverse_haversine(h): return 2 * asin(sqrt(h)) # radians def distance_between_points(lat1, lon1, lat2, lon2): # all args are in degrees # WARNING: loss of absolute precision when points are near-antipodal lat1 = radians(lat1) lat2 = radians(lat2) dlat = lat2 - lat1 dlon = radians(lon2 - lon1) h = haversine(dlat) + cos(lat1) * cos(lat2) * haversine(dlon) return RADIUS * inverse_haversine(h) def bounding_box(lat, lon, distance): # Input and output lats/longs are in degrees. # Distance arg must be in same units as RADIUS. # Returns (dlat, dlon) such that # no points outside lat +/- dlat or outside lon +/- dlon # are <= "distance" from the (lat, lon) point. # Derived from: http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates # WARNING: problems if North/South Pole is in circle of interest # WARNING: problems if longitude meridian +/-180 degrees intersects circle of interest # See quoted article for how to detect and overcome the above problems. # Note: the result is independent of the longitude of the central point, so the # "lon" arg is not used. dlat = distance / RADIUS dlon = asin(sin(dlat) / cos(radians(lat))) return degrees(dlat), degrees(dlon) if __name__ == "__main__": # Examples from Jan Matuschek article def test(lat, lon, dist): print "test bounding box", lat, lon, dist dlat, dlon = bounding_box(lat, lon, dist) print "dlat, dlon degrees", dlat, dlon print "lat min/max rads", map(radians, (lat - dlat, lat + dlat)) print "lon min/max rads", map(radians, (lon - dlon, lon + dlon)) print "liberty to eiffel" print distance_between_points(40.6892, -74.0444, 48.8583, 2.2945) # about 5837 km print print "calc min/max lat/lon" degs = map(degrees, (1.3963, -0.6981)) test(*degs, dist=1000) print degs = map(degrees, (1.3963, -0.6981, 1.4618, -1.6021)) print degs, "distance", distance_between_points(*degs) # 872 km 

Answer 2

So you calculate the distance between lat / long pairs using the haversine formula:

 import math R = 6371 # km dLat = (lat2-lat1) # Make sure it in radians, not degrees dLon = (lon2-lon1) # Idem a = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(lat1) * math.cos(lat2) * math.sin(dLon/2) * math.sin(dLon/2) c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) d = R * c; 

Now it’s trivial to check "d" (also in km) at your doorstep. If you want something more than km, adjust the radius.

Sorry, I can’t give you an insert solution, but I don’t understand your code skeleton (see comment).

Also note that these days you probably want to use the spherical law of cosines rather than Haversine. The advantages in numerical stability are no longer worth it, and it's damn easy to understand, encode and use.

Answer 3

If you store data in MongoDB, it tracks geolocation searches well and surpasses the pure-Python solutions above because it will handle you.

http://www.mongodb.org/display/DOCS/Geospatial+Indexing

Answer 4

The answer to John Machin's question helped me a lot. There is only a small mistake: latitude and longitude are replaced by boundigbox :

 dlon = distance / RADIUS dlat = asin(sin(dlon) / cos(radians(lon))) return degrees(dlat), degrees(dlon) 

this solves the problem. The reason is that longitudes do not change their distance by a degree, but by latitude. Their distance depends on longitude.