Skip to main content
Geographic Information Systems

Return to Answer

deleted 8 characters in body
Source Link
John Powell
John Powell
  • 13.7k
  • 5
  • 48
  • 62
  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 ON TRUE;pts;

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;

  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 ON TRUE;

  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;

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;
Source Link
John Powell
John Powell
  • 13.7k
  • 5
  • 48
  • 62

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 ON TRUE;

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.