The first "n" points near a polygon? (spatialite or postgis query)

Question

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.

Answer 1

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;

Answer 2

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