Optimal query in GeoDjango & Postgis for locations within a distance in meters

Question

I have geometries in a Postgis database with GeoDjango's default SRID, WGS84, and have found lookups directly in degrees to be much faster than in kilometres, because the database can skip the projections I would think.

Basically, Place.objects.filter(location__distance__lte=(point, D(km=10))) is several orders of magnitude slower than Place.objects.filter(location__dwithin=(point, 10)) as the first query produces a full scan of the table. But sometimes I need to lookup places with a distance threshold in kilometres.

Is there a somewhat precise way to convert the 10 km to degrees for the query? Maybe another equivalent lookup with the same performance that I should be using instead?

Cesya

Answer 1

You have several approaches to deal with your problem, here are two of them:

If you do not care much about precision you could use dwithin and use a naive meter to degree conversion degree(x meters) -> x / 40000000 * 360. You would have nearly exact results nearby the equator, but as you go north or south your distance would get shrink (shit we are living on a sphere). Imagine a region that is a circle in the beginning and shrinks to a infinite narrow elipse approaching one of the poles.

If you care about precision you can use:

max_distance = 10000 # distance in meter
buffer_width = max_distance / 40000000. * 360. / math.cos(point.y / 360. * math.pi)
buffered_point = point.buffer(buffer_width)
Place.objects.filter(
    location__distance__lte=(point, D(m=max_distance)),
    location__overlaps=buffered_point
)

The basic idea is to query for a all points that are within a circle around your point in degree. This part is very performant as the circle is in degreee and the geo index can be used. But the circle is sometimes a bit too big, so we left the filter in meters to filter out places that may be a bit farer away than the allowed max_distance.

Gonyea

btw, point.y can be -90 or 90, so you might have a division by zero exception – Gonyea

Answer 2

A small update to the answer of frankV.

    max_distance = 10000 # distance in meter
buffer_width = max_distance / 40000000. * 360. / math.cos(point.y / 360. * math.pi)
buffered_point = point.buffer(buffer_width)
Place.objects.filter(
    location__distance__lte=(point, D(m=max_distance)),
    location__intersects=buffered_point
)

I found the __overlaps doesn't work with postgresql and a point, but __intersects does.

To be sure it helps your query to speed-up, check the explain plan of the query (queryset.query to get the used query.)

Matron