# Optimizing a nearest neighbor calculation using PostGIS

I'm using PostGIS to calculate the nearest neighbors of polygons. What I want to calculate is the minimum distance from each polygon, to the nearest polygon.

Thus far I have got great help from Mike Toews' answer (which I quote with a minor change) here:

``````SELECT
a.hgt AS a_hgt,
b.hgt AS b_hgt,
ST_Distance(a.the_geom, b.the_geom) AS distance_between_a_and_b
FROM
public."TestArea" AS a, public."TestArea" AS b
WHERE
a.hgt !=  b.hgt AND ST_Distance(a.the_geom, b.the_geom) < 400
``````

Then I calculated the minimum:

``````SELECT a_hgt, MIN(distance_between_a_and_b)
FROM public."lon_TestArea"
GROUP BY a_hgt
``````

However, my challenge is to calculate this for a large number of polygons (1,000,000). As the above calculation compares each polygon to every other polygon, I wondered how I could improve the calculation so that I do not have to perform 10^12 calculations.

One thought I had was to buffer each polygon, and then to calculate the nearest neighbors of all values within the buffer for that polygon, and record the minimum. I'm not sure if that is the best approach, or whether there is a function in PostGIS that I should be using.

EDIT: Using one of Nicklas' suggestions, I'm experimenting with `ST_Dwithin()`:

``````CREATE TABLE mytable_withinRange AS SELECT
a.hgt AS a_hgt,
b.hgt AS b_hgt,
ST_DWithin(a.the_geom, b.the_geom, 400)
FROM
public."lon_TestArea" AS a, public."lon_TestArea" AS b
`````` This returns a table of the ID of each polygon, and whether it is within a certain distance or not. Is it possible to construct an `IF/ELSE` type statement using SQL? (I read about using the `CASE` condition) Or should I try joining the table I produce to the original table and then running the query again using ST_Distance?

• take a look at the second example in the boston gis link in my answer. you should use st_dwithin in the where part of the query. Jul 9 '11 at 9:53

There's a big "Nearest Neighbor" section on the BostonGIS page.

EDIT:

``````CREATE TABLE mytable_withinRange AS SELECT
a.hgt AS a_hgt,
b.hgt AS b_hgt
FROM
public."lon_TestArea" AS a, public."lon_TestArea" AS b
WHERE
ST_DWithin(a.the_geom, b.the_geom, 400)
``````

Concerning the CASE statement:

``````SELECT a,
CASE WHEN a=1 THEN 'one'
WHEN a=2 THEN 'two'
ELSE 'other'
END
FROM test;
``````
• Do you know if the line `WHERE ST_DWithin(a.the_geom, b.the_geom, 400)` will prevent distances greater then `400` to be calculated or just recorded? Also, can a case statement be used for numerical calculations? for example: `CASE WHEN ST_DWithin(a.the_geom, b.the_geom, 400) == TRUE THEN ST_DWithin(a.the_geom, b.the_geom)`
– djq
Jul 10 '11 at 12:59
• @celenius If distance is more than 400 m nothing in the select part will be calculated. I don't get why you want to put case in the mix. Jul 10 '11 at 17:21
• @Nicklas ok - I understand. I thought it might have meant that only distances less than 400 were stored; this makes it a lot easier than I though. Thanks!
– djq
Jul 10 '11 at 17:46

Hallo

There is some things do consider to make things move faster, and some things that might be possible in the future.

First, you mentioned that you are considering using a buffer to find polygons in some minimum range to avoid calculating all combinations.

As discussed in another link from Boston gis the right way to do that in PostGIS is using ST_Dwithin. ST_Dwithin uses the index to find the neighbors in a certain range.

It depends on the dataset of course if it is enough to just use a fixed value for st_DWithin for all polygons or if you need to do something like underdark and wildintellect are discussing.

A second thing is to use PostGIS 1.5+ here. That is because polygon to polygon calculations is much faster since 1.5 if their bounding boxes doesn't intersect. You can read more about that here..

The third thing to mention is the future.

In PostgreSQL 9.1 there will be something called knn-gist. That is an index that can not just answer yes or no but also return the result ordered directly from the index. You can read about that here.

But there will still be a lot of work to do on the PostGIS side before knn gist will help for things like this. There is a ticket for that here.

Regards

Nicklas

• Thanks for the suggestions Nicklas; as I've found it tricky to get pgAdmin/PostGIS up and running I think I will avoid using 1.5 at the moment. It seems like ST_Dwithin() is a way to solve this.
– djq
Jul 8 '11 at 21:35
• installing 1.5 will not affect the relation between postgresql and pgadmin. you can have more than one version of postgis in the database server and then you load one of them in the database. so you can have one 1.4 and one 1.5 database one the same database server. Jul 9 '11 at 9:55

The following pages related to Nathan Kerr's masters work provide some good insight on this direct issue. My coworker tried the Bostongis method here and here, but had some issues getting it to work right.

Another approach to think about that's similar to the buffer is to do an expanding/contracting rectangle. Basically pass 1 do a bounding box (it's a straight +x units to your original polygon's bbox) intersect that you think will catch at least one intersect. For the data that got an intersect perform a sub query that tests those matches for the closest. For the data the failed to match expand the bounding box and repeat.

It's clearly a recursive programming problem, and might be better done in Python with Shapely than 100% directly in postgis.