From some millions of points I am trying to select those points that are within a polygon and save them into a new table. Both datasets with spatial index and both on projected coordinates.

The query I am trying with is:

SELECT t1.point_id, t1.geom
FROM points t1
INNER JOIN boundaries t2
ON ST_WITHIN(t1.geom, t2.geom)

The query has been running for hours.

On the same machine I have done exactly the same thing using ArcMap in 25 minutes. I really don't understand this. I was under the impression that doing geoprocessing tasks directly on a RDBS (PostgreSQL or SQL Server) was the fastest option (provided the query is the right one).

Is there anything wrong with the query? Is there any way to optimise it? Am I missing something here?

PS: I tried with ST_Intersect instead of ST_WITHIN with the same result.

As suggested in 'Comments', I have tried with FROM being the smaller table (table with 3 polygons).

SELECT t1.point_id, t1.geom
FROM boundaries t2
INNER JOIN points t1
ON ST_WITHIN(t1.geom, t2.geom)

and also

SELECT t1.point_id, t1.geom
FROM boundaries t2
INNER JOIN points t1
ON ST_Contains(t2.geom, t1.geom)

Unfortunately, these two options seem to be pretty slow too. 2 hours running and still processing. I must admit that the polygon geometry is pretty complex (England boundary). But still, as indicated in the descriptions, other GIS tools take 25 minutes in doing the same job.

Explain analysis: (from last query written above, the one with "ST_Contains")

"Nested Loop  (cost=111.22..26768.81 rows=3007 width=44)"
"  ->  Seq Scan on gb_bdry t2  (cost=0.00..17.88 rows=3 width=32)"
"        Filter: ((country)::text ~~ 'ENGLAND'::text)"
"  ->  Bitmap Heap Scan on final_def_sgm_centroids t1  (cost=111.22..8906.96 rows=1002 width=44)"
"        Recheck Cond: (t2.geom ~ geom)"
"        Filter: _st_contains(t2.geom, geom)"
"        ->  Bitmap Index Scan on final_def_sgm_centroids_geom_idx  (cost=0.00..110.96 rows=3007 width=0)"
"              Index Cond: (t2.geom ~ geom)"

Different approach following suggestion from losbaltica:

First we limit the number of points with &&:

SELECT t1.point_id, t1.geom
FROM points t1
WHERE t1.geom && (SELECT geom FROM boundaries WHERE country = 'ENGLAND')

Then we do the actual intersection:

SELECT t1.point_id, t1.geom
FROM points t1,
    (SELECT geom FROM boundaries WHERE country = 'ENGLAND') AS t2
WHERE ST_Intersects(t1.geom, t2.geom)
  • please add the output of explain SELECT t1.point_id, t1.geom FROM points t1 INNER JOIN boundaries t2 ON ST_WITHIN(t1.geom, t2.geom) WHERE t2.country LIKE "England"
    – Ian Turton
    Oct 10, 2019 at 17:01
  • 5
    When comparing the relationships with many features on one side and few on the other, be sure to make your FROM be the smaller table. This is especially true when doing a points in poly query -- always use the poly to find points.
    – Vince
    Oct 10, 2019 at 18:02
  • I second what @Vince says. The order of your elements is relevant in that you should use the smaller number of features (polygons) as the "FROM" element.
    – Kartograaf
    Oct 10, 2019 at 18:07
  • switch table order; run VACUUM ANALYZE <tables>; try adding ... ON t1.geom && t2.geom AND ST_Intersects(...) ... to force index lookup
    – geozelot
    Oct 10, 2019 at 19:55
  • What your explain analyse is saying? Are your indexes are hit? What indexes are you using?
    – Losbaltica
    Oct 11, 2019 at 10:32

2 Answers 2


Looks like your query is a bit overcomplicated. To find all the points that are in "England you can simple run following query.

SELECT * FROM points
WHERE geom && (SELECT geom FROM boundaries WHERE country = "England");

Then if indexes are created properly you should have the following query plan: enter image description here

Now the problems with your current solution:

  1. Use "Like" only when you want to check if the column contains a string. Like its more expensive than normal equal.
  2. Don't use "ST_WITHIN" or "ST_Contains" if you looking only for an intersection. "ST_WITHIN" is more for distance searches and "ST_Contains" is not recommended for point intersections but for more complex geometry checks. Both commands are more expensive than "st_intersects".
  3. When you intersecting use "ST_INTERSECTS" or "&&" both got high index hit rate and should give you fast results.
  4. Don't do JOINs unless you want to get some information from the other table. Joins are great things but use them wisely. Sometimes simple "where" condition is all you need.
  5. You are missing a b-tree index on boundaries. Add it using the following query:

    CREATE INDEX ON boundaries (country);

I hope it helps you a bit to sort your issue out.

  • Thanks for all the tips. However, I have run the query as you suggested and I am not getting any output. 0 points selected (I have mapped both data sets to make sure there are points actually within the polygons; in fact hundreds of thousands). What does "table1.geom && table2.geom" exactly do? Oct 11, 2019 at 11:02
  • Is doing the same thing as "ST_INTERSECTS" but using bounding box. More details: postgis.net/docs/overlaps_geometry_box2df.html Try to replace && with "ST_INTERSECTS" and see if it make any changes. If not then you might have some geometry issues. Some more info about these 2 methods: dba.stackexchange.com/a/191669/138151
    – Losbaltica
    Oct 11, 2019 at 11:24
  • You may also try if SELECT geom FROM boundaries WHERE country = "England" is returning anything. Maybe your country filter is not equal to "England" and you do need LIKE.
    – Losbaltica
    Oct 11, 2019 at 11:29
  • I made a silly mistake. It's 'ENGLAND', not 'England'. That's why it was not selecting any. I have run the query using && and it takes seconds. SPECTACULAR. Just out of curiosity I have run the same thing but using ST_Intersects. I have stopped it after half an hour. How come there's such a huge difference? Oct 11, 2019 at 12:06
  • 2
    It is a lot different! && checks for overlaps of bounding boxes, ST_Intersects does also check if any vertice of a geometry overlaps with the other; for && it's comparing two vectors against two others, for ST_Intersects it also need to traverse all vertices of both geometries! The difference in results: consider POLYGON((0 0, 3 3, 6 0, 0 0)) and POINT(1 2): && will return TRUE, ST_Intersects will return FALSE.
    – geozelot
    Oct 11, 2019 at 12:38

You should definitely try to subdivide your boundary polygons! The spatial index serves only to select candidate geometries based on their bounding boxes. Then real geometry has to be rechecked for each candidate to produce real results (as you can see in the execution plan). So... if your boundary polygon is complex, that recheck is very expensive! That can be avoided by subdividing large polygons to smaller pieces and then perform spatial operations on that smaller polygons.

Firs create a new table with subdivided polygons:

create table boundaries_subdiv as select country, st_subdivide(geom, 50) geom from boundaries;

The second parameter on st_subdivide defines the maximum number of vertices that resulting polygons can have. You can try other numbers, but 50 to 100 seems to be a good option (by my observations...)

Then create indexes on that table

create index boundaries_subdiv_geom on boundaries_subdiv using gist (geom);
create index boundaries_subdiv_country on boundaries_subdiv (country);

Analyze table:

analyze boundaries_subdiv;

Then perform your query, but with subdivided boundaries:

SELECT t1.point_id, t1.geom
FROM points t1
INNER JOIN boundaries_subdiv t2
ON ST_WITHIN(t1.geom, t2.geom)

Or more refined one with st_intersects:

create table result AS
  select t1.point_id, t1.geom
    from points t1
    join boundaries_subdiv t2 ON st_intersects(t1.geom, t2.geom)
   where t2.country = "ENGLAND"

That should help a lot! Further reading here.

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.