[GIS] PostGIS Intermittent INDEX Performance

postgispostgresqlspatial-index

I have a table containing approximately 55 million data points ( point is a geometry with SRID 4326 ) and for my query I need to join this to an area table ( currently have ~1800 areas ) which contains a variety of different ranging from large polygons ( 2000 km square ) to fairly small ( small being about 100 km square ).

The initial query selected by the user narrows down the initial 55 million points to around ~300,000 points dependent on date range etc they select. Then the join is done and dependent on what area set they have selected to use once the query is complete this normally narrows it down to a ~150,000.

The problem I'm having is that some times the query just grinds to a halt and instead of taking the expected ~25 seconds it can take up to ~18 minutes. At this point have to usually do a VACUUM ANALYSE and then run a few queries before it starts behaving itself again. No data has been added, updated or removed from the data or areas tables at this point.

I've played around with everything I can think of and this still seems to keep happening with no constancy too it. Both the data.point column and the area.polygon column have GIST INDEXES on them.

I found removing the INDEX from the data.point column seemed to make things a bit more stable but is slower ~35 seconds normally. However removing an INDEX seems to be a very bad choice as should it not be helping not hindering?

I'm using PostgreSQL 9.1.4 with PostGIS 1.5

Here's the query I'm running

    select * FROM data, area WHERE st_intersects (data.point, area.polygon) AND 
(readingdatetime BETWEEN '1948-01-01' AND '2012-11-19') AND datasetid IN(3) AND
 "polysetID" = 1 AND area.id IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11)

EXPLAIN

Nested Loop  (cost=312.28..336.59 rows=5 width=2246) (actual time=1445.973..11557.824 rows=12723 loops=1)
  Join Filter: _st_intersects(data.point, area.polygon)
  ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.017..0.229 rows=35 loops=1)
        Index Cond: ("polysetID" = 1)
        Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))
  ->  Bitmap Heap Scan on data  (cost=312.28..316.29 rows=1 width=297) (actual time=328.771..329.136 rows=641 loops=35)
        Recheck Cond: ((point && area.polygon) AND (datasetid = 3))"
        Filter: ((readingdatetime >= '1948-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-11-19 00:00:00'::timestamp without time zone))
        ->  BitmapAnd  (cost=312.28..312.28 rows=1 width=0) (actual time=328.472..328.472 rows=0 loops=35)
              ->  Bitmap Index Scan on data_point_index  (cost=0.00..24.47 rows=276 width=0) (actual time=307.115..307.115 rows=1365770 loops=35)
                    Index Cond: (point && area.polygon)
              ->  Bitmap Index Scan on data_datasetid_index  (cost=0.00..284.37 rows=12856 width=0) (actual time=1.522..1.522 rows=19486 loops=35)
                    Index Cond: (datasetid = 3)
Total runtime: 11560.879 ms

My create tables

CREATE TABLE data
(
  id bigserial NOT NULL,
  datasetid integer NOT NULL,
  readingdatetime timestamp without time zone NOT NULL,
  value double precision NOT NULL,
  description character varying(255),
  point geometry,
  CONSTRAINT "DATAPRIMARYKEY" PRIMARY KEY (id ),
  CONSTRAINT enforce_dims_point CHECK (st_ndims(point) = 2),
  CONSTRAINT enforce_geotype_point CHECK (geometrytype(point) = 'POINT'::text OR point IS NULL),
  CONSTRAINT enforce_srid_point CHECK (st_srid(point) = 4326)
);

CREATE INDEX data_datasetid_index ON data USING btree (datasetid);
ALTER TABLE data CLUSTER ON data_datasetid_index;

CREATE INDEX "data_datasetid_readingDatetime_index" ON data USING btree (datasetid , readingdatetime );
CREATE INDEX data_point_index ON data USING gist (point);

CREATE INDEX "data_readingDatetime_index" ON data USING btree (readingdatetime );

CREATE TABLE area
(
  id serial NOT NULL,
  polygon geometry,
  "polysetID" integer NOT NULL,
  CONSTRAINT area_primary_key PRIMARY KEY (id )
)

CREATE INDEX area_polygon_index ON area USING gist (polygon);
CREATE INDEX "area_polysetID_index" ON area USING btree ("polysetID");
ALTER TABLE area CLUSTER ON "area_polysetID_index";

Hope that all makes a degree of sense if need to know anything else please ask.

Short summary is really that the INDEXES seem to work some times but not others.

Could anyone suggest anything I could try to work out what's happening?

Thanks in advance.

EDIT:

Another example

select * FROM data, area WHERE st_intersects ( data.point, area.polygon) AND 
(readingdatetime BETWEEN '2009-01-01' AND '2012-01-19') AND datasetid IN(1,3) AND
 "polysetID" = 1 AND area.id IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11) 

Run on copy of table with point index

Nested Loop  (cost=0.00..1153.60 rows=35 width=2246) (actual time=86835.883..803363.979 rows=767 loops=1)
  Join Filter: _st_intersects(data.point, area.polygon)
  ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.021..16.287 rows=35 loops=1)
        Index Cond: ("polysetID" = 1)
        Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))
  ->  Index Scan using data_point_index on data  (cost=0.00..1133.30 rows=1 width=297) (actual time=17202.126..22952.706 rows=33 loops=35)
        Index Cond: (point && area.polygon)
        Filter: ((readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone) AND (datasetid = ANY ('{1,3}'::integer[])))
Total runtime: 803364.120 ms

Run on copy of table without point index

Nested Loop  (cost=2576.91..284972.54 rows=34 width=2246) (actual time=181.478..235.608 rows=767 loops=1)
  Join Filter: ((data_new2.point && area.polygon) AND _st_intersects(data_new2.point, area.polygon))
  ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.149..0.196 rows=35 loops=1)
        Index Cond: ("polysetID" = 1)
        Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))
  ->  Bitmap Heap Scan on data_new2  (cost=2576.91..261072.36 rows=90972 width=297) (actual time=4.808..5.599 rows=2247 loops=35)
        Recheck Cond: ((datasetid = ANY ('{1,3}'::integer[])) AND (readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone))
        ->  Bitmap Index Scan on "data_new2_datasetid_readingDatetime_index"  (cost=0.00..2554.16 rows=90972 width=0) (actual time=4.605..4.605 rows=2247 loops=35)
              Index Cond: ((datasetid = ANY ('{1,3}'::integer[])) AND (readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone))
Total runtime: 235.723 ms

As you can see the query is significantly slower when the point index is being used.

EDIT 2 ( Pauls Suggested Query ):

WITH polys AS (
  SELECT * FROM area
  WHERE "polysetID" = 1 AND area.id IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11)
)
SELECT * 
FROM polys JOIN data ON ST_Intersects(data.point, polys.polygon)
WHERE datasetid IN(1,3) 
AND (readingdatetime BETWEEN '2009-01-01' AND '2012-01-19');

EXPLAIN

Nested Loop  (cost=20.04..1155.43 rows=1 width=899) (actual time=16691.374..279065.402 rows=767 loops=1)
  Join Filter: _st_intersects(data.point, polys.polygon)
  CTE polys
    ->  Index Scan using "area_polysetID_index" on area  (cost=0.00..20.04 rows=1 width=1949) (actual time=0.016..0.182 rows=35 loops=1)
          Index Cond: ("polysetID" = 1)
          Filter: (id = ANY ('{28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11}'::integer[]))
  ->  CTE Scan on polys  (cost=0.00..0.02 rows=1 width=602) (actual time=0.020..0.358 rows=35 loops=1)
  ->  Index Scan using data_point_index on data  (cost=0.00..1135.11 rows=1 width=297) (actual time=6369.327..7973.201 rows=33 loops=35)
        Index Cond: (point && polys.polygon)
        Filter: ((datasetid = ANY ('{1,3}'::integer[])) AND (readingdatetime >= '2009-01-01 00:00:00'::timestamp without time zone) AND (readingdatetime <= '2012-01-19 00:00:00'::timestamp without time zone))
Total runtime: 279065.540 ms

Best Answer

Effectively forcing the planner to do the thing you want might help. In this case, sub-setting the polygon table prior to executing the spatial join with the points table. You might be able to outwit the planner using "WITH" syntax:

WITH polys AS (
  SELECT * FROM area
  WHERE area.id in IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11)
)
SELECT * 
FROM polys JOIN data ON ST_Intersects(data.point, polys.polygon)
WHERE datasetid IN(3) 
AND (readingdatetime BETWEEN '1948-01-01' AND '2012-11-19');

The trouble with trying to play these games is that you are coding into your statement the assumption "my polygon list will always be more selective than my other query portions". Which might not be true for all parameterizations of your query, or for all applications of a particular query over a heterogeneously distributed dataset.

But it might work.

UPDATE: This goes even further down the dangerous road of assuming you know the selectivity of your clauses beforehand, this time we also take the attribute selection on the point table out and do it separately before the spatial join:

WITH polys AS (
  SELECT * FROM area
  WHERE area.id in IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,3,4,5,6,22,23,24,7,8,9,10,11)
),
WITH points AS (
  SELECT * FROM data
  WHERE datasetid IN(3) 
  AND (readingdatetime BETWEEN '1948-01-01' AND '2012-11-19')
)
SELECT * 
FROM polys JOIN points ON ST_Intersects(points, polys.polygon);
Related Question