[GIS] Checking if a polygon contains multiple point geometries using PostGIS

postgis

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.

Best Answer

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;

Related Solutions

[GIS] Problem PostGIS Spatial Query ST_Contains ST_Disjoint ST_Intersects

You have a set of polygon, they cover an area of the map.

You have a set of points, some are inside the polygons and some are outside.

The st_intersect will compare all points against all polygons and only return points that are "on top" one or more polygons.

The st_disjoint will return a point if a point is not "on top" of a polygon. So for a set of points against a set of polygons, so unless the polygons overlap, each point is going to be "outside" of one or more polygons.

You can have a nested query of the intersect, and another inverted nested query of the intersect to get your t/f result

Adding a solution

There's a fiddle that goes with this here, since it was a PostgreSQL fiddle rather than PostGIS I've change the fiddle test to be a simple integer equality test.

So the first half of the sql is your intersect test:

    SELECT distinct points.id
    FROM points, multipolygon
    WHERE st_intersects(points.geom, multipolygon.geom)

I add a TRUE column to match your sample, by adding a wrapping SQL select with a hardcoded value of TRUE:

SELECT
   id,
   TRUE AS intersects_at_least_one_polygon
FROM 
(
    SELECT distinct points.id
    FROM points, multipolygon
    WHERE st_intersects(points.geom, multipolygon.geom)
) intersecting ;

To find the points that are NOT intersecting the polygons, well that is the total list of point id, minus the above points. So I select all the points, the outer select below, and then say "WHERE id NOT IN" and just repeat your SQL, plus I add a hardcode a FALSE:

SELECT
   points.id,
   FALSE AS intersects_at_least_one_polygon
FROM
   points WHERE id NOT IN
   (
        SELECT distinct points.id
        FROM points, multipolygon
        WHERE st_intersects(points.geom, multipolygon.geom)
   )
   order by id
   ;


SELECT
   id,
   TRUE AS intersects_at_least_one_polygon
FROM 
(
    SELECT distinct points.id
    FROM points, multipolygon
    WHERE st_intersects(points.geom, multipolygon.geom)
) intersecting
UNION
SELECT
   points.id,
   FALSE AS intersects_at_least_one_polygon
FROM
   points WHERE id NOT IN
   (
        SELECT distinct points.id
        FROM points, multipolygon
        WHERE st_intersects(points.geom, multipolygon.geom)
   )
   order by id
   ;

I'm certain that there are more elegant, efficient, solutions to this, but hope you're happy with a simple working version.

[GIS] Importing MultiPolygon into PostGIS

Your MultiPolygon is invalid. You can test it with ogrinfo:

ogrinfo -ro -dialect sqlite -sql "select IsValid(geometry) from
OGRGeoJSON" multipolygon.json
INFO: Open of `multipolygon.json'
      using driver `GeoJSON' successful.
GEOS warning: Self-intersection at or near point -115.907104 39.162109000000001

Layer name: SELECT
Geometry: None
Feature Count: 1
Layer SRS WKT:
(unknown)
isvalid(geometry): Integer (0.0)
OGRFeature(SELECT):0
  isvalid(geometry) (Integer) = 0

You can try to correct it in PostGIS with ST_MakeValid http://postgis.net/docs/ST_MakeValid.html. However, it is not simple

Step 1 Correct the faulty geometry:

create table valid as select st_makevalid(geometry) as geometry from multipolygon;

Step 2 Test:

select * from valid where ST_Intersects(ST_PointFromText('POINT( -116.024551 38.485773 )',4326), geometry);

Result: error

ERROR:  Relate Operation called with a LWGEOMCOLLECTION type.  This is unsupported.
HINT:  Change argument 2: 'GEOMETRYCOLLECTION(POLYGON((-115.000793 38.677307,-115.000793 38.499878,-115....'

Pity that ST_Intersects does not handle geometrycollections. So you must get rid of the collection. First part of the collection is the polygon so you can take just that.

Step 3 Keep just the polygon part

update valid set geometry=ST_GeometryN(geometry,1);

Step 4 Test again

select * from valid where ST_Intersects(ST_PointFromText('POINT( -116.024551 38.485773 )',4326), geometry);

Result: Success!

Best Answer

Related Solutions

[GIS] Problem PostGIS Spatial Query ST_Contains ST_Disjoint ST_Intersects

[GIS] Importing MultiPolygon into PostGIS

Related Question