Spatial index/query (finding k nearest points)

Question

I have +10k points (latitude, longitude) and I'm building an app that shows you the k nearest points to a user's location.

I think this is a very common problem and I don't want to reinvent the wheel. I'm learning about Quadtrees. It seems to be a good approach to solve this spatial problem.

I'm using these tools:

• Python 2.5
• MySQL
• MongoDb

Building the Quadtree is not that hard: http://donar.umiacs.umd.edu/quadtree/points/pointquad.html But once I've created the tree and saved it to a db (MySQL or MongoDb), how I run the query?

I need to run queries like these:

1. Find all points within 10 km of the user's location.
2. Find the 6 (or at least 6) nearest points to the user's location.

What's the standard and common approach to do it?

EDIT 1:

I've loaded the +10k points into MongoDB (Geospatial indexing) and it works fine at first glance. Anyway I've found PostGis:

PostGIS is an extension to the PostgreSQL object-relational database system which allows GIS (Geographic Information Systems) objects to be stored in the database.

So I think I'll give PostGis a try.

I've also found SimpleGeo. You can store points/places in the cloud and then query them via an API: https://simplegeo.com/docs/tutorials/python#how-do-radial-nearby-query

Answer 1

MongoDB has support for spatial indexes built-in, so all you'd need to do is load your points using the correct format, create the spatial index, and then run your queries.

For a quick example, I loaded the center points for all 50 states in the mongo shell:

> db.places.ensureIndex({loc: "2d"})
> db.places.save({name: "AK", loc: {long: -152.2683, lat: 61.3850}})
> db.places.save({name: "AL", loc: {long: -86.8073, lat: 32.7990}})
> db.places.save({name: "AR", loc: {long: -92.3809, lat: 34.9513}})
> db.places.save({name: "AS", loc: {long: -170.7197, lat: 14.2417}})
> ...

Next, to query for the 6 nearest points to a given location:

> db.places.find({loc: { $near: {long: -90, lat: 50}}}).limit(6)
{"name" : "WI", "loc" : { "long" : -89.6385, "lat" : 44.2563 } }
{"name" : "MN", "loc" : { "long" : -93.9196, "lat" : 45.7326 } }
{"name" : "MI", "loc" : { "long" : -84.5603, "lat" : 43.3504 } }
{"name" : "IA", "loc" : { "long" : -93.214, "lat" : 42.0046 } }
{"name" : "IL", "loc" : { "long" : -89.0022, "lat" : 40.3363 } }
{"name" : "ND", "loc" : { "long" : -99.793, "lat" : 47.5362 } }

Next, to query for all points within 10km of a given location. Since I'm calculating the nearest states, I'll use 888km (which is approximately 8 degrees of latitude):

> db.places.find({loc: { $near: {long: -90, lat: 50}, $maxDistance: 8}})
{"name" : "WI", "loc" : { "long" : -89.6385, "lat" : 44.2563 } }
{"name" : "MN", "loc" : { "long" : -93.9196, "lat" : 45.7326 } }

Since one degree of latitude is approximately 111.12km, you'd use a $maxDistance: 0.08999 to represent 10km for your application.

Updated By default MongoDB assumes an "idealized flat earth model" but this results in inaccuracies since longitude lines converge at the poles. MongoDB versions 1.7+ support spherical distance calculations, which provides the increased precision.

Here is an example of running the above query using spherical distance. the maxDistance is in radians, so we need to divide by the earth's average radius:

> db.runCommand({geoNear: "places", near: [-90, 50], spherical: true, 
                 maxDistance: 800/6378});
(summarizing results as they're too verbose to include)
"MN"  dis: 0.087..
"WI"  dis: 0.100..
"ND"  dis: 0.120..

Answer 2

You may want to look at kdtree entry in wikipedia. This would be useful when you have more than two dimensions too (unlike quadtrees). I suggest the kd-tree because the entry has python code for creating and querying the tree.

Answer 3

If you want to use MongoDB, then read their docs carefully. The default model is flat earth. It assumes that a degree of longitude has the same length as a degree of latitude.

I quote: """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. This is only true at the equator where they are both about equal to 69 miles or 111km. However, at the 10gen offices at { x : -74 , y : 40.74 } one arcdegree of longitude is about 52 miles or 83 km (latitude is unchanged). This means that something 1 mile to the north would seem closer than something 1 mile to the east."""

You need their "new spherical model". Be warned: you need to use (longtitude, latitude) in that order -- again, read their docs carefully.