4

I am using PostGIS and have two tables. One containing points of start-locations in EPSG:4326 (about 3 million records in total) and another one containing county-polygons in UTM EPSG:25832 (about 30 records). I wanted to calculate the amount of start points per county. So I used this query:

SELECT
count(f.ride_uid) as amount,
l.name,
l.u_key
FROM schema_a.rides as f

LEFT JOIN
schema_b.counties as l
ON ST_Intersects(f.start_geom,ST_Transform(l.geometry,4326)) -- f.start_geom are points in 4326 and l.geometry are polygons in 25832 (UTM)

WHERE
f.start_date_year = 2021
and f.start_date_month in (1,2,3)

GROUP BY
l.name,
l.u_key

But it ran endlessly and I cancelled it after 2 hours, wondering whats wrong. Tried it again after checking the code a second time, but again no result after 15 minutes. Then trying

SELECT
count(f.ride_uid) as amount,
l.name,
l.u_key
FROM schema_a.rides as f

LEFT JOIN
schema_b.counties as l
ON ST_Intersects(ST_Transform(f.start_geom,25832),l.geometry) -- f.start_geom are points in 4326 and l.geometry are polygons in 25832 (UTM)

WHERE
f.start_date_year = 2021
and f.start_date_month in (1,2,3)

GROUP BY
l.name,
l.u_key

gives the correct result in about 1 second.

Why does it matter which geometry I transform? Why does ST_Intersects(ST_Transform(f.start_geom,25832),l.geometry) work so fast and ST_Intersects(f.start_geom,ST_Transform(l.geometry,4326)), well, not at all?

1
  • 2
    because you lose the spatial index when you transform (unless you add it specially). Run Explain on the queries to be sure
    – Ian Turton
    Commented May 3, 2021 at 9:36

2 Answers 2

7

Because the GEOMETRY returned from ST_Transform is not covered by any index.

The planner chooses to run the points table against the polygon table and its GIST index on geometry, which is built with GEOMETRY(<TYPE>, 4326); using a transformed GEOMETRY, or any other derived GEOMETRY, does not match the index definition, and a sequential scan between both tables is performed, including the calculation of intersection between every point and every polygon.


You could add a functional index on the transformed GEOMETRY:

CREATE INDEX ON <polygon_table>
  USING GIST ( (ST_Transform(geom, 4326) )

;
4
  • 1
    You may wonder 'why not using the index on the point table'? That's a delicate matter. Is there an index on the point table geometries? What are the table counts (many polygons vs few points; index sizes)? What are the geometric extents (how packed is the index; how much positive matches are expected)? Are the statistics updated (VACUUM ANALYZE <both_tables>;)? Does removing the GROUP BY statement make a difference? Indexes are most performant when you look for the few among the many; if table a has 1M rows, and b 1k, run a query against the index on b.
    – geozelot
    Commented May 3, 2021 at 10:11
  • @MrXsquared and is there an index on that points table? Are both tables vacuum analyzed? The planner should absolutely choose to run against the point table index!
    – geozelot
    Commented May 3, 2021 at 10:30
  • ..except it decides to utilize indexes on those other filter columns.
    – geozelot
    Commented May 3, 2021 at 10:33
  • 1
    Ah, @JGH is right about the 'why not the index on points'; I completely ignored the fact that you are using a LEFT JOIN - it forces a seq scan on the points table by design.
    – geozelot
    Commented May 3, 2021 at 12:58
5

@Geozelot has addressed why the polygon layer index is not usable.

Now there are more issues here.

The index on the point layer is also not usable because the query is starting from the point layer and then is doing a left join to the polygon layer. A full table scan of the point layer is to be done anyways, so the spatial index is not used. (and, as answered by @Geozelot, the polygon index is not usable because its geometry is transformed)

The left join makes sense only if you want to count points that don't intersect with any polygon, but such points can be found differently.

With such point vs polygon counts, the query should start from the polygon layer and make a regular join with the point layer. Let's phrase the two joins:

  • For 30 polygons, check which points are contained (and make use of the index on the 3M points)
  • For 3 M points, check in which polygon they are (and make use of the index on the 30 polygons)

it would be something like

SELECT
  count(f.ride_uid) as amount,
  l.name,
  l.u_key
FROM 
  schema_b.counties as l
JOIN schema_a.rides as f
  ON ST_Intersects(f.start_geom,ST_Transform(l.geometry,4326)) --ok to transform the polygon, we want to use the point index
WHERE
  f.start_date_year = 2021
  and f.start_date_month in (1,2,3)
GROUP BY
  l.name,
  l.u_key

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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