Checking if a polygon contains multiple point geometries using PostGIS?

Question

This is not a a duplicate of the "Count Points in Polygon with Postgis", since in the ST_Contains(geometry geomA, geometry geomB), this case has two geometries (start_geom and end_geom) in geomB, while the other question is about only one geometry.

Input polygon table us_county:

Column  |            Type             |                        Modifiers                        
----------+-----------------------------+--------------------------------
 gid      | integer                     | not null default nextval('us_county_gid_seq'::regclass)
 statefp  | character varying(2)        | | 
 geoid    | character varying(5)        | 
 name     | character varying(100)      | 
 geom     | geometry(MultiPolygon,4269) |

Input point table trips:

Column    |         Type         |                      Modifiers                      
--------------+----------------------+-----------------------------------
 id           | integer              | 
 gid          | integer              | not null default nextval('trips_gid_seq'::regclass)
 start_geom  | geometry(Point,4326) | 
 end_geom    | geometry(Point,4326) | 

I want to use to check if a polygon in us_county could contain both start_geom and end_geom of a trip points.

I've tried the following queries:

1. this query returns an empty result, which is not right

select us_county.name, count(trips.gid) 

from trips, us_county

where st_contains(us_county.geom, trips.start_geom) and 
      st_contains(us_county.geom, trips.end_geom)

group by us_county.name;

2. a second thought is to use ST_Collect() function:

select us_county.name, count(trips.gid) 

from trips, us_county 

where st_contains(ST_Transform(us_county.geom, 4326),  
ST_Collect(trips.start_geom, trips.end_geom))

Group by us_county.name;

But this turns out to be very slow, for 3k polygon and 50k points, the query run over 40min.

Answer 1

You can combine the points with ST_Collect:

SELECT us_county.name, count(trips.gid) 
FROM trips, us_county
WHERE ST_Contains(us_county.geom, 
                  ST_Collect(trips.start_geom, trips.end_geom))
GROUP BY us_county.name;

Answer 2

It looks like your main issue is the different projection systems used. You should be able to re-project on the fly to complete the query using ST_Transform:

select us_county.name, count(trips.gid) 
from trips, us_county
where st_contains(ST_Transform(us_county.geom,4326), trips.start_geom) and 
      st_contains(ST_Transform(us_county.geom,4326), trips.end_geom)
group by us_county.name;

Answer 3

One way that you could get this to work and speed up the query would be to "region tag" the geometries of the points table first. This would look something like:

ALTER TABLE trips
ADD COLUMN start_county_id char(5) -- assuming this is a fips code
ADD COLUMN end_county_id char(5);

UPDATE trips
SET start_county_id = a.geoid
FROM us_county AS a
WHERE ST_Intersects(trips.start_geom, a.geom);

UPDATE trips
SET end_county_id = a.geoid
FROM us_county AS a
WHERE ST_Intersects(trips.end_geom, a.geom);

The updates will be expensive, but this is often acceptable since it only has to occur once. Then you all you have to do is match the start/end county id's to the county fips code. You can put an index on the start_county_id and end_county_id columns, and the query should be much faster. Something like:

SELECT a.geoid, c.cnt
FROM us_counties AS a,
     (SELECT b.start_county_id, b.end_county_id, COUNT(b.id) as cnt
     FROM trips as b
     WHERE b.start_count_id = b.end_county_id
     GROUP BY b.start_county_id) as c
WHERE a.geoid = b.start_county_id

Spatial queries by nature are expensive. I have a points table with 20 million + rows and a polygon table of 75,000, and this technique cuts queries (very similar to yours) down to just a few seconds.