2

Background

I want to add a model manager function that filters a queryset based on the proximity to coordinates. I found this blog posting with code that is doing precisely what I want.

Code

The snippet below seems to make use of geopy functions that have since been removed. It coarsely narrows down the queryset by limiting the range of latitude and longitude.

    # Prune down the set of all locations to something we can quickly check precisely
    rough_distance = geopy.distance.arc_degrees(arcminutes=geopy.distance.nm(miles=distance)) * 2
    queryset = queryset.filter(
        latitude__range=(latitude - rough_distance, latitude + rough_distance), 
        longitude__range=(longitude - rough_distance, longitude + rough_distance)
    )

Problem

Since some of the used geopy functions have been removed/moved, I'm trying to rewrite this stanza. However, I do not understand the calculations---barely passed geometry and my research has confused me more than actually helped me.

Can anyone help? I would greatly appreciate it.

Belmin Fernandez
  • 8,307
  • 9
  • 51
  • 75

6 Answers6

2

In case anybody else is looking at this now, since I tried to use geopy and just hit up against it, the modern equivalent of the rough_distance snippet above is:

import geopy
rough_distance = geopy.units.degrees(arcminutes=geopy.units.nautical(miles=1))
Tim Perry
  • 11,766
  • 1
  • 57
  • 85
1

It looks like distance in miles is being converted to nautical miles, which are each equal to a minute of arc, which are 1/60th of an arc degree each. That value is then doubled, and then added and subtracted from a given latitude and longitude. These four values can be used to form a bounding box around the coordinates.

You can lookup any needed conversion factors on Wikipedia. There's also a relevant article there titled Horizontal position representation which discusses pros and cons of alternatives to longitude and latitude positioning which avoid some of their complexities. In other words, about the considerations involved with replacing latitude and longitude with another horizontal position representation in calculations.

martineau
  • 119,623
  • 25
  • 170
  • 301
  • fantastic! This section of the Wikipedia article helps tremendously: http://en.wikipedia.org/wiki/Minute_of_arc#Cartography – Belmin Fernandez Dec 01 '10 at 17:04
  • 1
    After peeking at the latest geopy-0.94 sources, it looks like the only thing missing that the code in the snippet uses is the `geopy.distance.arc_degrees(arcminutes=xxx)` function, which you ought to be able to easily replace. Interestingly, comments in the `distance.py` code indicate that it is already using the WGS-84 ellipsoid internally, so there may be no need for pyproj and PROJ.4 (although it maybe be useful to you for you application). – martineau Dec 01 '10 at 20:16
1

The Earth is not a sphere, only approximately so. If you need a more accurate calculation, use pyproj. Then you can calculate the location based a reference ellipsoid (e.g. WGS84).

Andy Barbour
  • 8,783
  • 1
  • 24
  • 35
1

martineau's answer is right on, in terms of what the snippet actually does, but it is important to note that 1 minute of arc represents very different distances depending on location. At the equator, the query covers the least axis aligned bounding box enclosing a circle of diameter distance, but off the equator, the bounding box does not completely contain that circle.

SingleNegationElimination
  • 151,563
  • 33
  • 264
  • 304
  • Thanks. I guess the code is trying to limit it very roughly. Is this your evaluation as well or do you have a suggestion on how it could be more precise? – Belmin Fernandez Dec 01 '10 at 18:17
1

This code from the blog is sloppy:

  def near(self, latitude=None, longitude=None, distance=None):
    if not (latitude and longitude and distance):
      return []

If latitude == 0 (equator) or longitude == 0 (Greenwich meridian), it returns immediately. Should be if latitude is None or longitude is None .......

@TokenMacGuy's answer is an improvement, but:

(a) The whole idea of the "bounding box" is to avoid an SQL or similar query calculating a distance to all points that otherwise satisfy the query. With appropriate indexes, the query will execute much faster. It does this at the cost of leaving the client to (1) calculate the coordinates of the bounding box (2) calculate and check the precise distance for each result returned by the query.

If step 2 is omitted, you get errors, even at the equator. For example "find all pizza shops in a 5-mile radius" means you get answers up to 7.07 miles (that's sqrt(5*2 + 5*2)) away in the corners of the box.

Note that the code that you show seems to be arbitrarily doubling the radius. This would mean you get points 14.1 miles away.

(b) As @TokenMacGuy said, away from the equator, it gets worse. The bounding box so calculated does not include all points that you are interested in -- unless of course you are overkilling by doubling the radius.

(c) If the circle of interest includes either the North or South Pole, the calculation is horribly inexact, and needs adjusting. If the circle of interest is crossed by the 180-degree meridian (i.e. the International Date Line without the zigzags), the results are a nonsense; you need to detect this case and apply a 2-part query (one part for each side of the meridian).

For solutions for problems (b) and (c), see this article.

John Machin
  • 81,303
  • 11
  • 141
  • 189
  • Very extensive answer +1. Thank you! Researching your points right now and might have some additional questions soon. – Belmin Fernandez Dec 02 '10 at 00:23
  • @John Machin: So, summarizing your points, the code might find pizza joints too far away near the equator and would be no good at all finding them near either Pole. Somehow I suspect that neither of those would be big problems for the OP's application. ;-) – martineau Dec 07 '10 at 10:55
  • @martineau: Wrong. The bounding box produced at the equator is twice as high and twice as wide as it should be. That's **FOUR** times as much work as necessary. The "twice" is the first silliness. The second is that the correct box is not square once you get away from the equator, so at e.g. Seattle (lat: 47 deg) it's doing only 2 * 1.37 == 2.73 times the work it should. Now think what happens if someone notices and fixes the obvious(?) first silliness -- a bounding box of +/- 0.1447 deg for a 10-mile distance, when in Seattle it should be +/- 0.2122 deg longitude and you are MISSING points. – John Machin Dec 08 '10 at 03:55
  • @martineau: Summary: Two wrongs don't make a right. We know nothing about the OP's application, nor about the application of people who read and believe accepted answers uncritically. – John Machin Dec 08 '10 at 03:59
  • @John Machin: It was only meant to be a little comic relief. It's true, I don't know the purpose of the OP's or other applications of everyone who might read these posts, but I'm not worried because we have people like yourself around to set us straight. For the record, when I said "might find pizza joints too far away near the equator" I was alluding to the fact that the bounding box was too big near there. – martineau Dec 08 '10 at 04:41
  • @martineau: OK ... I'll chill out; speaking of which, I do hope that that blogger is not the navigator on a nuclear submarine under the polar icecap :-O – John Machin Dec 08 '10 at 05:08
0

If the coordinates on the earth are known, you can use geopy to get a good estimate of the decimal degrees to miles (or any distance units) scale at that point:

SCALE_VAL = 0.1
lat_scale_point = (cur_lat + SCALE_VAL, cur_long)
long_scale_point = (cur_lat, cur_long + SCALE_VAL)
cur_point = (cur_lat, cur_long)

lat_point_miles = distance.distance(cur_point, lat_scale_point).miles
long_point_miles = distance.distance(cur_point, long_scale_point).miles

# Assumes that 'radius_miles` is the range around the point you want to look for
lat_rough_distance = (radius_miles / lat_point_miles) * SCALE_VAL
long_rough_distance = (radius_miles / long_point_miles) * SCALE_VAL

Some caveats:

  • Special-case handling for the the scale points is needed around polls or prime meridean
  • Depending on how large or small you want your radius to be, you could pick a more appropriate SCALE_VAL
crchurchey
  • 188
  • 6