ST_Intersection fills up LIMIT with empty geometry points

Question

I have managed to get the intersection between a circle around a point and objects in a geometry using this PostGIS query:

SELECT ST_Intersection(sea_level.geom, ST_Buffer(ST_Transform(ST_SetSRID(ST_MakePoint(13.01667,55.71667), 4326), 3006), 10000))
FROM sea_level
LIMIT 20000;

This works fine if I set the limit high enough to get "valuable" objects, but most of the returned rows are just "empty" objects (as GeoJSON):

{"type":"Polygon","crs":{"type":"name","properties":{"name":"EPSG:3006"}},"coordinates":[]}

If I set the limit too low, I will only get these "empty" rows, and none of the valuable rows with the actual intersecting objects. What are these "empty" rows and why are they returned?

Answer 1

Because my comment about proximity searches seems to have caused confusion: the correct way to implement a high performance proximity filter to help with this specific scenario would be:

SELECT ST_Intersection(sl.geom, ST_Buffer(poi, 10000)) AS geom
FROM   sea_level AS sl
CROSS JOIN LATERAL
       ST_Transform(ST_SetSRID(ST_MakePoint(13.01667, 55.71667), 4326), 3006) AS poi
WHERE  ST_DWithin(sl.geom, poi, 10000)
;

Here

• ST_DWithin quickly filters for all sl.geom that would not fall inside the ST_Buffer, and thus
  • reduces the actual (and costly) ST_Intersection operation, on (less but still costly) ST_Buffer'ed geometries, to only those sl.geom that are relevant
  • in effect does not include NULL values for non-intersecting features
• the CROSS JOIN LATERAL moves the point-of-interest creation (poi) to a cache-able function value expression, mainly to avoid repeated calls

Without the need to intersect an ST_Buffer, the filter would e.g. use a call to ST_Intersects. In general, for where it's applicable, and functions interchangeable, the benefit of using ST_DWithin is its algorithmic complexity, which may (dependent on geometry types) be lesser than that of a spatial overlap check except for the worst-case.

---

If you see any significant performance boost is questionable: the planner may neglect any indexes completely for a table that small, and a few thousand spatial overlap calculations on simple point buffers are fairly fast these days. But add two orders of magnitude to the row count, or more complex geometries, and the difference in execution times will be minutes to seconds.

Answer 2

To avoid returning empty intersection results, run the intersection only on geometries which actually intersect by using a ST_Intersects condition:

SELECT ST_Intersection(sea_level.geom, poi.geom)
FROM sea_level JOIN
  (SELECT ST_Buffer(ST_Transform(ST_SetSRID(
              ST_MakePoint(13.01667,55.71667), 4326), 3006), 10000) AS geom) poi
ON ST_Intersects(sea_level.geom, poi.geom)
LIMIT 20000;

This also improves performance, because:

• ST_Intersects takes advantage of a spatial index on sea_level
• the ST_Intersects calculation is faster than ST_Intersection, so it's better to use it first