2

I've two spatial dataset, the fist one is a polygons layers (My_polygons) and the second one is a points layer (My_points).

I need to find, for each polygons, the "nearest" "N" points (N=10) using a spatial query in spatialite or postgis and store this information somewhere in my database for future queries.

9

One nice way of doing this is using the new LATERAL JOIN syntax in Postgres in conjunction with the <-> operator. A lateral join essentially runs the right hand query for each row in the left hand one, and can avoid some quite nasty array logic and sub-queries, which you would otherwise need to use to find k nearest neighbours.

CREATE TABLE sometable AS
SELECT 
    poly.id, array_agg(pts.id) 
  FROM
     (SELECT id, geom FROM polygons) poly
  CROSS JOIN LATERAL
     (SELECT id, geom 
        FROM points pt
        ORDER BY pt.geom <-> poly.geom LIMIT 10
  ) pts;

There is a good blog by Paul Ramsey explaining this technique.

There are a few gotchas though. The <-> operator operates on the spatial index, but generally requires that one of the two geometries is a constant. For small tables with LIMIT 1, this may not matter, but, as always EXPLAIN is your friend. I have found, in practice, that for very large tables (I have done this for tables of 50 million points, to find each points 6 nearest neighbours, which in the naive case gives you 50 x 50 million possible combinations) you might need to add

 ST_DWithin(pt.geom, poly.geom, distance)

to the query, which will use the spatial index, but requires you to know beforehand a value for distance that is guaranteed to find at least 10 points for each polygon.

  SELECT 
    poly.id, array_agg(pts.id) 
  FROM
     (SELECT id, geom FROM polygons) poly
  CROSS JOIN LATERAL
     (SELECT id, geom 
        FROM points pt
        WHERE ST_DWithin(pt.geom, poly.geom, some_distance)
        ORDER BY pt.geom <-> poly.geom LIMIT 10
  ) pts;

I have found, in practice, that for fairly large tables where you are looking for several nearest neighbours, this can be done in a loop, where you start with a small search distaince in ST_DWithin, and gradually increase it for those polygons you haven't round the nearest 10 points for (which obviously requires a where clause in the initial select). Sadly, while it is syntactically very nice LATERAL JOINS and the <-> operator by themselves do not guarantee fast execution automatically, when you are looking for the k nearest neighbours for multiple geometries.

Note. It is possible to use LEFT JOIN LATERAL syntax too, but in this instance, you have to add ON TRUE, after the right-hand table, eg,

SELECT 
  poly.id, array_agg(pts.id) 
FROM
  (SELECT id, geom FROM polygons) poly
LEFT JOIN LATERAL
  (SELECT id, geom 
    FROM points pt
    WHERE ST_DWithin(pt.geom, poly.geom, some_distance)
    ORDER BY pt.geom <-> poly.geom LIMIT 10
) pts ON TRUE;
3

Something like this in postgres

SELECT PG.id,

       (SELECT array_agg(id) FROM
          (SELECT PT.id FROM My_points ORDER BY PG.geom <-> PT.geom LIMIT 10) FOO
       ) as points_id

INTO new_table
FROM My_polygons as PG

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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