3

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.

  • Possible duplicate of Count Points in Polygon with Postgis – Midavalo Apr 19 '16 at 7:35
  • 1
    Is it only a typo in a question that you create fields start_geom and end_geom but query for pickup_geom and dropoff_geom? When you write it fails, how does it fail? Do you get an error from PostGIS or is the result something else than what you'd like to get? Edit this into your question. – user30184 Apr 19 '16 at 8:26
  • @user30184, thanks for correcting me on this. I have made edits. – enaJ Apr 19 '16 at 8:31
2

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;
  • I was thinking using ST_Collect, too. But the operation took a long time and actually never finished (I've tried run it for more than 40min and no signs for finish, so I killed the operation). – enaJ Apr 19 '16 at 18:15
  • I have 3k polygons and 50k points in total. Is that the problem of the big data or the bugs in SQL scripts? – enaJ Apr 19 '16 at 18:16
  • 1
    40 minutes sounds like a really long time for a join of 50k points against TIGER counties. Maybe post your query plan in the question? – dbaston Apr 19 '16 at 20:45
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;
  • That is true the projection systems were mixed, I've corrected that though it still doesn't work. – enaJ Apr 19 '16 at 18:04
  • What error do you get? Or do you just get a count of? – HeikkiVesanto Apr 19 '16 at 19:56
  • this query returns nothing. But I was expecting it return a list of non-zero values – enaJ Apr 19 '16 at 20:03
  • And if you add them to QGIS with the projections mentioned they overlap? – HeikkiVesanto Apr 19 '16 at 20:12
  • yes, the point and polygon are overlapping with each other – enaJ Apr 19 '16 at 20:45
2

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.

  • @enaJ No problem. I haven't tested the last query on a problem exactly like yours, so it may be a little off, but you get the idea. You can also index the start_county_id and end_county_id columns to further speed up the calculations. B-tree indices are the default in postgresql and probably more appropriate for this type of data. Once again, the indexing may take a long time but only has to be performed once. – haff Apr 20 '16 at 23:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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