1

I need to filter geocodes for near-ness to a location. For example, I want to filter a list of restaurant geocodes to identify those restaurants within 10 miles of my current location.

Can someone point me to a function that will convert a distance into latitude & longitude deltas? 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 am using the supremum norm for distance.

Andrew B.
  • 1,010
  • 6
  • 19
  • I'm sorry, I don't understand. Do you want inverse_haversine to return a callable that takes the "other" parameter and returns True or False? Or are you planning to pass "other" in some other way? – drxzcl Jul 05 '10 at 21:37
  • (1) "can someone point me": someone == google (2) "provides dlat, dlng such that" stuff that doesn't mention dlat, dlng -- please edit your question. (3) What is "the supremum norm for distance distance metric"? – John Machin Jul 05 '10 at 21:50
  • @ John Machin. Of course google is your friend, too, for learning about the supremum norm. – Andrew B. Jul 05 '10 at 21:57
  • @Ranieri Sorry, edited the function to clarify dlat, dlng. Inverse haversine should give the max deltas in lat/lng, relative to the reference geocode, to define a max_distance square region around the reference geocode. – Andrew B. Jul 05 '10 at 21:59
  • @AndrewB: I know what "supremum norm" is. I asked "What is "the supremum norm for distance distance metric". Also "inverse haversine" is a function that is the inverse of the haversine function (http://mathworld.wolfram.com/InverseHaversine.html); you are misusing the term to describe your problem which is in some sense the inverse of the easier (point1, point2) -> distance problem (which may be calculated using haversines OR other methods). – John Machin Jul 05 '10 at 22:56
  • @John Machin. Good point on "inverse haversine" -- you are correct that what I am seeking is different (though it may rely on) the inverse haversine function. I'll work up an appropriate edit for the question. Btw, it would be nice if you directed more of your commenting effort on this question toward an actual solution.... – Andrew B. Jul 06 '10 at 01:35
  • @AndrewB: I've put a bit of effort into an actual solution. Over to you. – John Machin Jul 07 '10 at 21:01
  • @AndrewB: The general idea is that a questioner should accept an answer or give feedback as to why a serious answer attempt is not acceptable or ask for clarification ... not just maintain silence. I have a few improvments and further suggestions but am not interested in typing into a black hole. – John Machin Jul 09 '10 at 23:20

4 Answers4

6

There seems not to have been a good Python implementation. Fortunately the SO "Related articles" sidebar is our friend. This SO article points to an excellent article that gives the maths and a Java implementation. The actual function that you require is rather short and is embedded in my Python code below. Tested to extent shown. Read warnings in 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'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
Community
  • 1
  • 1
John Machin
  • 81,303
  • 11
  • 141
  • 189
1

This is how you calculate distances between lat/long pairs using the haversine formula:

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;

It is now trivial to test "d" (also in km) against your threshold. If you want something else than km, adjust the radius.

I'm sorry I can't give you a drop-in solution, but I do not 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 a hell of a lot simple to understand, code and use.

drxzcl
  • 2,952
  • 1
  • 26
  • 28
  • Yes, I have the haversine formula. I am looking for the implementation of the inversion. – Andrew B. Jul 05 '10 at 22:00
  • Fixed the code snippet. Basically I want to compute the max/min lat/lng values that are within a max_distance range of a geocode. So inverse_haversine should return a dlat, dlng tuple that are the +/- deltas that would still be within range. Does that help? – Andrew B. Jul 05 '10 at 22:02
0

If you store data in MongoDB, it does nicely indexed geolocation searches for you, and is superior to the pure-Python solutions above because it will handle optimization for you.

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

A. Jesse Jiryu Davis
  • 71
  • 3
  • 1
    Chris Columbus wouldn't have "discovered" America with this gear: "The current implementation assumes an idealized model of a flat earth, meaning that an arcdegree of latitude (y) and longitude (x) represent the same distance everywhere." – John Machin Jul 06 '10 at 06:21
  • Yes, that *will* be a limitation near the poles. A location N or S of the center point will be considered closer than it should be, and a location E or W will be considered farther. – A. Jesse Jiryu Davis Jul 14 '10 at 15:37
0

John Machin's answer helped me much. There is just a small mistake: latitudes and longitudes are swapped in 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 don't changes their distance per degree - but latitudes do. Their distance is depending on the longitude.

nlaan
  • 46
  • 3