Tuesday, 6 March 2018

geoprocessing - Select features that do NOT intersect in PostGIS


This seems to me like such a simple question (and it probably is) yet I cannot seem to find an example that gives me the answer. Using PostGIS, I just want to select points that fall outside of polygons. Ultimately this is the inverse of the ST_Intersects, as far as I can see it.


Example: I have a taxlot layer and a address point layer. I assume I should be using the ST_Intersects, but how do I tell it to do the reverse selection? I thought maybe adding a NOT statement in front of the code below, but that did not work.


CREATE table t_intersect AS
SELECT
hp.gid,
hp.st_address,
hp.city,
hp.st_num,

hp.the_geom
FROM
public.parcel as par,
public.housepoints as hp
WHERE
ST_Intersects(hp.the_geom,par.the_geom);

Answer



The reason it doesn't work with "not intersects" is that you only compare geometries in pairs; there will be the same problem with disjoint. Every housepoint will disjoint some parcels even if it intersects one parcel.


underdark's suggestion doesn't have that problem. There is also another trick that probably will make more effective use of indexes:


CREATE TABLE t_intersect AS

SELECT
hp.gid,
hp.st_address,
hp.city,
hp.st_num,
hp.the_geom
FROM
public.housepoints AS hp LEFT JOIN
public.parcel AS par ON
ST_Intersects(hp.the_geom,par.the_geom)

WHERE par.gid IS NULL;

The idea is to join them with st_intersects and get the rows where parcel id is not present.


The indexes needed here are a spatial index and an index on gid in parcels (assuming that id in parcels table is called gid too).


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...