Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

enter image description here

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?

share|improve this question
    
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. –  Nicklas Avén Jul 9 '11 at 9:53
add comment

3 Answers 3

up vote 6 down vote accepted

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


EDIT:

How about

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;
share|improve this answer
    
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
1  
@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. –  Nicklas Avén 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
add comment

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

share|improve this answer
    
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
2  
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. –  Nicklas Avén Jul 9 '11 at 9:55
    
Thanks for the encouragement! Successfully using v1.5.2 –  djq Jul 10 '11 at 14:56
add comment

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.