Updating attribute values of points outside area using PostGIS

Question

I have two datasets in PostGIS 1.5 (on Postgres 9.1) one of polygons and one of points.

I am hoping to create a SQL command that will select all of the points in a dataset that are not within of the polygon features and assign a value of -9999 to one of the columns. In this case the boundary is a continent so -9999 defines a feature off the mainland.

I have several hundred thousand features and have tried to achieve this more or less manually in QGIS. However once the features are selected I cannot interact with QGIS without it crashing.

---

Having a couple of problems with performance that has this function taking 96+ hours and counting. I tried stopping previously and set a LIMIT of 100 on the SELECT and it was still taking a few hours.

This is with the GiST and BTree indexing.

I tried running SELECT on point_id only without the UPDATE function and got similar results.

Previously in QGIS it was taking about 16 hours to select all points within the polygons.

Will look at some performance tuning and retry a bit later

---

I am 99% sure Index are on but there is always a little room for doubt =D

Result of the first query is

QUERY PLAN
Nested Loop Left Join (cost=0.00..12955070.02 rows=1 width =4)
Join Filter _st_intersects(point_table.the_geom, polygon_table.the_geom)
Filter: (polygon_table.polygon_id IS NULL)
-> Seq Scan on point_table (cost=0.00..209494.58 rows=1562658 width=104)
-> Index Scan using poly_table_geom_idx on polygon_table (cost=0.00..7.89 rows=1 width=22520) Index Cond: (point_table.the_geom && the_geom)
(6 Rows)

Result of second query is

sum
59257303
(1 row)

Result of comparing the two queries with LIMIT 1000 is that the first takes about 15 minutes and the second is still running after an hour.

I have split the features using ST_Dump and at present I am setting up the Indexing, Keys etc for that.

---

Have complete ST_Dump

Nested Loop cost increased to 13142577, Seq Scan stayed the same, Index Scan width increased to 1162940

Point layer contains 1728921 points

Answer 1

The timings you have had seems strange. Something is not as expected.

To check that the indexes really are used, what is the text result from:

EXPAIN
SELECT point_id FROM 
point_table LEFT JOIN polygon_table 
ON ST_Intersects(point_table.geom, polygon_table.geom)
WHERE polygon_id IS NULL;

How many vertex points does your polygon data set hold? What is the result of:

SELECT SUM(ST_NPoints(the_geom)) FROM polygon_table;

To analyse what is happening I think the easiest is to run parts of the query to see how that behaves, like you have done with "limit 100".

So, if you run the following queries and study the timing it should give a hint about what the problem might be:

SELECT count(*) 
FROM point_table INNER JOIN polygon_table 
ON point_table.geom && polygon_table.geom LIMIT 1000;

compare that timing and result to running the actual ST_Intersects function:

SELECT count(*) 
FROM point_table INNER JOIN polygon_table 
ON ST_Intersects(point_table.geom,polygon_table.geom) 
LIMIT 1000;

If the result from query 1 is a much higher number than query 2 and query 2 takes a lot more time than query 1, the problem is that the index selectivity. If it is MultiPolygons in the data set the first thing then is to split them to Polygons with ST_Dump and build a new index on the result.

There is more things that can be done like splitting the bigger polygons with a grid, but that is more complex.

Start with that to see if that is the slow part. If not it is time to continue testing with the LEFT JOIN and the "polygon_id is not null" part.

---

Yes, this is a PostGIS task.

I have not tested the below and I am not sure that it is the most efficient way to do the update, but the subquery should be quite fast. The idea is to find all points that ARE intersecting with the polygons and then ask for the points that is not found in the polygons. Something like this should work:

UPDATE point_table SET p_id=-9999 WHERE point_id IN
(SELECT point_id FROM point_table LEFT JOIN polygon_table 
ON ST_Intersects(point_table.geom, polygon_table.geom)
WHERE polygon_id IS NULL);

Put gist indexes on the geoemtry columns and btree index on the point id.