Filter zipcodes by proximity in Django with the Spherical Law of Cosines

Question

I'm trying to handle proximity search for a basic store locater in Django. Rather than haul PostGIS around with my app just so I can use GeoDjango's distance filter, I'd like to use the Spherical Law of Cosines distance formula in a model query. I'd like all of the calculations to be done in the database in one query, for efficiency.

An example MySQL query from The Internet implementing the Spherical Law of Cosines like this:

SELECT id, ( 
    3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * 
    cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * 
    sin( radians( lat ) ) ) 
) 
AS distance FROM stores HAVING distance < 25 ORDER BY distance LIMIT 0 , 20;

The query needs to reference the Zipcode ForeignKey for each store's lat/lng values. How can I make all of this work in a Django model query?

Revulsive

(a) That is NOT the haversine formula; it's the spherical law of cosines formula; see (for example) http://www.movable-type.co.uk/scripts/latlong.html and look at the respective Wikipedia articles. (b) I do trust that you'll replace the hard-coded user coordinates with variables :-) (c) The gentle reader should be warned that your unit of distance is somewhat archaic (something to do with 1000 * (length of Roman legion's standard pace), I believe) :-) – Jethro

Jeez, what's so hard about GeoDjango? Just install it :) – Freeland

Answer 1

It's possible the execute raw SQL queries in Django.

My suggestion is, write the query to pull a list of IDs (which it looks like you're doing now), then use the IDs to pull the associated models (in a regular, non-raw-SQL Django query). Try to keep your SQL as dialect-independent as possible, so that you won't have to worry about one more thing if you ever have to switch databases.

To clarify, here's an example of how to do it:

def get_models_within_25 (self):
    from django.db import connection, transaction
    cursor = connection.cursor()

    cursor.execute("""SELECT id, ( 
        3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * 
        cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * 
        sin( radians( lat ) ) ) )
        AS distance FROM stores HAVING distance < 25
        ORDER BY distance LIMIT 0 , 20;""")
    ids = [row[0] for row in cursor.fetchall()]

    return MyModel.filter(id__in=ids)

As a disclaimer, I can't vouch for this code, as it's been a few months since I've written any Django, but it should be along the right lines.

Natala

It works just fine, only needs triple-quoting (or to be turned into a single string). – Metallography

Just a follow up to this. The original query returns a "distance" field (which shows the distance between the 2 sets of log/lat). How would i perform the second part, but with this extra 'distance' field? – Rochette

I'm asuming that "ids = [row[0] for row in cursor.fetchall()] return MyModel.filter(id__in=ids)" gets all the IDS from the first query, then using djangos ORM to select all rows from a model? – Rochette

Answer 2

To follow up on Tom's answer, it won't work in SQLite by default because of SQLite's lack of math functions by default. No problem, it's pretty simple to add:

class LocationManager(models.Manager):
    def nearby_locations(self, latitude, longitude, radius, max_results=100, use_miles=True):
        if use_miles:
            distance_unit = 3959
        else:
            distance_unit = 6371

        from django.db import connection, transaction
        from mysite import settings
        cursor = connection.cursor()
        if settings.DATABASE_ENGINE == 'sqlite3':
            connection.connection.create_function('acos', 1, math.acos)
            connection.connection.create_function('cos', 1, math.cos)
            connection.connection.create_function('radians', 1, math.radians)
            connection.connection.create_function('sin', 1, math.sin)

        sql = """SELECT id, (%f * acos( cos( radians(%f) ) * cos( radians( latitude ) ) *
        cos( radians( longitude ) - radians(%f) ) + sin( radians(%f) ) * sin( radians( latitude ) ) ) )
        AS distance FROM location_location WHERE distance < %d
        ORDER BY distance LIMIT 0 , %d;""" % (distance_unit, latitude, longitude, latitude, int(radius), max_results)
        cursor.execute(sql)
        ids = [row[0] for row in cursor.fetchall()]

        return self.filter(id__in=ids)

Carmella

Seems like magic, how does this work? I guess it's pretty inefficient? – Rhodium

The 'from mysite import settings' can be made more generic by using 'from django.conf import settings' – Metaphysic

Answer 3

Just to follow up on jboxer's answer, here's the whole thing as part of a custom manager with some of the hard-coded stuff turned into variables:

class LocationManager(models.Manager):
    def nearby_locations(self, latitude, longitude, radius, max_results=100, use_miles=True):
        if use_miles:
            distance_unit = 3959
        else:
            distance_unit = 6371

        from django.db import connection, transaction
        cursor = connection.cursor()

        sql = """SELECT id, (%f * acos( cos( radians(%f) ) * cos( radians( latitude ) ) *
        cos( radians( longitude ) - radians(%f) ) + sin( radians(%f) ) * sin( radians( latitude ) ) ) )
        AS distance FROM locations_location HAVING distance < %d
        ORDER BY distance LIMIT 0 , %d;""" % (distance_unit, latitude, longitude, latitude, int(radius), max_results)
        cursor.execute(sql)
        ids = [row[0] for row in cursor.fetchall()]

        return self.filter(id__in=ids)

Metallography

Answer 4

To follow up on Tom, if you want to have a query that also works in postgresql, you can not use AS because you will get an error saying 'distance' does not exist.

You should put the whole spherical law expresion in the WHERE clause, like this (It also works in mysql):

import math
from django.db import connection, transaction
from django.conf import settings

from django .db import models

class LocationManager(models.Manager):
    def nearby_locations(self, latitude, longitude, radius, use_miles=False):
        if use_miles:
            distance_unit = 3959
        else:
            distance_unit = 6371

        cursor = connection.cursor()

        sql = """SELECT id, latitude, longitude FROM locations_location WHERE (%f * acos( cos( radians(%f) ) * cos( radians( latitude ) ) *
            cos( radians( longitude ) - radians(%f) ) + sin( radians(%f) ) * sin( radians( latitude ) ) ) ) < %d
            """ % (distance_unit, latitude, longitude, latitude, int(radius))
        cursor.execute(sql)
        ids = [row[0] for row in cursor.fetchall()]

        return self.filter(id__in=ids)

Please note that you have to select the latitude and longitude, otherwise you can not use it in the WHERE clause.

Bischoff

Answer 5

Following jboxer's response

def find_cars_within_miles_from_postcode(request, miles, postcode=0):

    # create cursor for RAW query
    cursor = connection.cursor()

    # Get lat and lon from google
    lat, lon = getLonLatFromPostcode(postcode)

    # Gen query
    query = "SELECT id, ((ACOS(SIN("+lat+" * PI() / 180) * SIN(lat * PI() / 180) + COS("+lat+" * PI() / 180) * COS(lat * PI() / 180) * COS(("+lon+" - lon) * PI() / 180)) * 180 / PI()) * 60 * 1.1515) AS distance FROM app_car HAVING distance<='"+miles+"' ORDER BY distance ASC"

    # execute the query
    cursor.execute(query)

    # grab all the IDS form the sql result
    ids = [row[0] for row in cursor.fetchall()]

    # find cars from ids
    cars = Car.objects.filter(id__in=ids)

    # return the Cars with these IDS
    return HttpResponse( cars )

This returns my cars from x amount of miles, this works well. However the raw query returned how far they were from a certain location, i think the fieldname was 'distance'.

How can i return this field 'distance' with my car objects?

Rochette

Answer 6

Using some of the proposed answers above, I was getting incosistent results so I decided to check the equation again using [this link]http://www.movable-type.co.uk/scripts/latlong.html as a reference, the equation is d = acos(sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1) ) * 6371 where d is the distance to be calculated,

lat1,lon1 is the coordinate of the base point and lat2,lon2 is the coordinate of the other points which in our case are points in the database.

From the above answers, the LocationManager class looks like this

class LocationManager(models.Manager):
def nearby_locations(self, latitude, longitude, radius, max_results=100, use_miles=True):
    if use_miles:
        distance_unit = 3959
    else:
        distance_unit = 6371

    from django.db import connection, transaction
    from mysite import settings
    cursor = connection.cursor()
    if settings.DATABASE_ENGINE == 'sqlite3':
        connection.connection.create_function('acos', 1, math.acos)
        connection.connection.create_function('cos', 1, math.cos)
        connection.connection.create_function('radians', 1, math.radians)
        connection.connection.create_function('sin', 1, math.sin)

    sql = """SELECT id, (acos(sin(radians(%f)) * sin(radians(latitude)) + cos(radians(%f))
          * cos(radians(latitude)) * cos(radians(%f-longitude))) * %d)
    AS distance FROM skills_coveragearea WHERE distance < %f
    ORDER BY distance LIMIT 0 , %d;""" % (latitude, latitude, longitude,distance_unit, radius, max_results)
    cursor.execute(sql)
    ids = [row[0] for row in cursor.fetchall()]

    return self.filter(id__in=ids)

Using the site [link]http://www.movable-type.co.uk/scripts/latlong.html as check, my results where consistent.

Serpigo

Answer 7

It's also possible to do this using Django's database functions, which means you can add a distance_miles column using a .annotate() call and then sort by it. Here's an example:

from django.db.models import F
from django.db.models.functions import ACos, Cos, Radians, Sin

locations = Location.objects.annotate(
    distance_miles = ACos(
        Cos(
            Radians(input_latitude)
        ) * Cos(
            Radians(F('latitude'))
        ) * Cos(
            Radians(F('longitude')) - Radians(input_longitude)
        ) + Sin(
            Radians(input_latitude)
        ) * Sin(Radians(F('latitude')))
    ) * 3959
).order_by('distance_miles')[:10]

Wall

Answer 8

@classmethod def nearby_locations(cls, latitude, longitude, radius, max_results=1000, use_miles=False): if use_miles: distance_unit = 3959 else: distance_unit = 6371000

from django.db import connection, transaction
from django.conf import settings
cursor = connection.cursor()

sql = """SELECT id, (%f * acos( cos( radians(%f) ) * cos( radians( latitude ) ) *
cos( radians( longitude ) - radians(%f) ) + sin( radians(%f) ) * sin( radians( latitude ) ) ) )
AS distance FROM yourapp_yourmodel 
GROUP BY id, latitude, longitude
HAVING (%f * acos( cos( radians(%f) ) * cos( radians( latitude ) ) *
cos( radians( longitude ) - radians(%f) ) + sin( radians(%f) ) * sin( radians( latitude ) ) ) ) < %d
ORDER BY distance OFFSET 0 LIMIT %d;""" % (distance_unit, latitude, longitude, latitude, distance_unit, latitude, longitude, latitude, int(radius), max_results)

cursor.execute(sql)
ids = [row[0] for row in cursor.fetchall()]

return cls.objects.filter(id__in=ids)

Cosmopolis