4

In PostGIS I have two tables:

Table that contains polygons:

table a.provinces

gid | geom       |   name    |
------------------------------
1   | 010604 ... | Champagne |

Table that contains points:

table a.stations

gid | geom     | province_short |
---------------------------------
1   | 0132 ... | ch             |

They are not connected with foreign key!

I would like to get:

  1. Only points that are inside of the polygon
  2. Only points that are outside of the polygon

SRID: 3857

Something like this: enter image description here

Can it be done just like this or I have to modify database?

  • 2
    Use ST_Intersects and ST_Disjoint. Read documents and tutorials. – user30184 Sep 12 '16 at 17:28
  • 1
    Welcome to gis.stackexchange! Please note that a good question on this site is expected to show some degree of research on your part, i.e. what you have tried and - if applicable - code so far. For more info, you can check our faq. – underdark Sep 12 '16 at 17:47
8
-- get all the points inside polygons
SELECT poi.*
FROM polygons pol
JOIN points poi ON (ST_Within(poi.geom, pol.geom));


-- get all the points not in polygon
-- if polygons overlap it might be better to go for NOT EXISTS clause
SELECT poi.*
FROM polygons pol
LEFT JOIN points poi ON (ST_Within(poi.geom, pol.geom))
WHERE pol.gid IS NULL;
  • We have trouble with this both in PostGIS and SQL Server... and I can't this sample to work either... – DPSSpatial Sep 12 '16 at 17:39
  • Thanks for fast response! I have a problem with part --get all the points not in polygon...NOTICE: gserialized_gist_joinsel: jointype 1 not supported. I guess it is because of "LEFT" part. – vayacondios2015 Sep 12 '16 at 21:44
  • It's just a notice, not an error. – Michal Zimmermann Sep 13 '16 at 5:26
  • I made views from your querries and then I made visualization of first one in QGIS, but second view itself wasnt shown in PostGIS layer at all and in that case table view was empty. I guessed all of this its happening because of the notice (not error, sorry about that), but maybe I am wrong. Any suggestion? – vayacondios2015 Sep 13 '16 at 14:21
  • 1
    Views not working in QGIS are not related to that notice in any way. It's just QGIS doesn't handle spatial views very well. – Michal Zimmermann Sep 14 '16 at 7:58
6

All points within the polygon:

SELECT pts.id 
FROM pts JOIN polys 
ON ST_Intersects(pts.geom, polys.geom)

Now the hard one. Pairwise disjointness is not helpful, and neither is it particularly easy to optimize searches for disjointness. So instead, run the intersects test again, but keep around the rows that fail the test, using a left join:

SELECT pts.id
FROM pts LEFT JOIN polys
ON ST_Intersects(pts.geom, polys.geom)
WHERE polys.id IS NULL

Yes, it looks insane, but it works, and it's basically the fastest way to get the answer.

1

This isn't an answer so much as it is an echo of the problem at hand - using ST_Disjoint to test if features intersect another feature doesn't seem to give the expected results:

Given the following query which tests which schools are in a particular boundary:

select distinct sch.* 
, ST_Disjoint(sch.geom, b.geom) as out
from dpsdata."Schools_Current"  as sch
, dpsdata."SchoolBoundaries_All" as b
where b.high_num = '451'

The results are not what would answer the question orginally answered above.

While features that 'do' intersect are given a 'f' value

enter image description here

...features that are given a 't' value are also found in the resultset that 'do not' intersect:

enter image description here

Is there something wrong with the query I'm using with ST_Disjoint?

  • ST_Disjoint compares all points to all polygons, and unless all the polygons are the same (and they're not, in your case), there'll always be points disjoint from some polygons. The result in the first picture is correct, you're asking PostGIS to give you points not lying in the green polygon - all of them lie in that polygon though, so you get false as a result. – Michal Zimmermann Sep 14 '16 at 8:02
  • @MichalZimmermann I've updated the image slightly - sill the same issue with the features getting 'T' in the 2nd image - they should not be features that intersect the polygon, isn't that right? – DPSSpatial Sep 14 '16 at 14:22
  • There is another discussion on this topic here: gis.stackexchange.com/questions/210184/… – DPSSpatial Sep 16 '16 at 15:54
0

Method 4

1) Create two worksheets for stations by running the query:

SELECT * INTO stations_0 FROM stations; SELECT * INTO stations_1 FROM stations;

2) Delete the stations that are included in the polygons:

DELETE FROM stations_0 as t1 WHERE EXISTS (SELECT 1 FROM provinces as t2 WHERE ST_Intersects(t2.geom, t1.geom));

and do not enter polygons:

DELETE FROM stations_1 as t1 WHERE NOT EXISTS (SELECT 1 FROM provinces as t2 WHERE ST_Intersects(t2.geom, t1.geom));

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.