Filtrado geográfico de Python por distancia.

Necesito filtrar geocodificaciones para estar cerca de una ubicación. Por ejemplo, quiero filtrar una lista de códigos geográficos de restaurantes para identificar esos restaurantes a 10 millas de mi ubicación actual.

¿Puede alguien indicarme una función que convierta una distancia en deltas de latitud y longitud? Por ejemplo:

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 

Nota: Estoy usando la norma suprema para la distancia.

Parece que no ha habido una buena implementación de Python. Afortunadamente, la barra lateral de “Artículos relacionados” es nuestra amiga. Este artículo de SO apunta a un excelente artículo que da las matemáticas y una implementación de Java. La función real que necesita es bastante corta y está incrustada en mi código Python a continuación. Probado hasta el punto mostrado. Lea las advertencias en los comentarios.

 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's 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 

Así es como se calculan las distancias entre pares lat / largos utilizando la fórmula de haversine:

 import math R = 6371 # km dLat = (lat2-lat1) # Make sure it's 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; 

Ahora es trivial probar “d” (también en km) contra su umbral. Si quieres algo más que km, ajusta el radio.

Lo siento, no puedo darte una solución inmediata, pero no entiendo el código de tu código (ver comentario).

También tenga en cuenta que en estos días probablemente desee utilizar la ley esférica de los cosenos en lugar de Haversine. Las ventajas en la estabilidad numérica ya no valen la pena, y es mucho más fácil de entender, codificar y usar.

Si almacena datos en MongoDB, realiza búsquedas de geolocalización bien indexadas para usted, y es superior a las soluciones puras de Python anteriores porque manejará la optimización para usted.

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

La respuesta de John Machin me ayudó mucho. Solo hay un pequeño error: las latitudes y longitudes se intercambian en el boundigbox de boundigbox :

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

Esto resuelve el problema. La razón es que las longitudes no cambian su distancia por grado, pero sí las latitudes. Su distancia depende de la longitud.