As user30184 noted, you have to use ST_Crosses() instead of ST_Intersects().
Additional, your expression works on one line compared to one polygon only: The black line touches Polygon 1 and doesn't cross Polygon 1. In more detail:
SELECT x FROM line, polygons WHERE ...
means: Build a case from every pair of a line and a polygon, so number of cases is n_lines x n_polygons. Do the WHERE clause for every case and keep cases where it is true. Result: blue line, P1: false, blue line, P2: true, black line, P1: true, black line, P2: false.
Maybe you want something like this:
SELECT lines.geom
FROM lines, polygons
WHERE ST_Touches(lines.geom, polygons.geom) AND
NOT EXISTS (SELECT 1 FROM polygons p2 WHERE ST_Crosses(lines.geom, p2.geom));
To SELECT
- those that are within others (courtyards):
SELECT a.*
FROM <buildings> AS a
WHERE EXISTS (
SELECT 1
FROM <buildings> AS b
WHERE a.<id> <> b.<id>
AND ST_Within(a.geom, b.geom)
);
- those that contain others (houses):
SELECT a.*
FROM <buildings> AS a
WHERE EXISTS (
SELECT 1
FROM <buildings> AS b
WHERE a.<id> <> b.<id>
AND ST_Contains(a.geom, b.geom)
);
Note that the EXISTS
construct is slightly more performant and more elegant than a JOIN
in cases where multiple courtyards exists for a single house.
To UPDATE
those that contain others with the footprint of those contained (subtract courtyards from houses; produce Polygons with inner rings):
UPDATE <houses> AS a
SET geom = (
SELECT ST_Difference(a.geom, ST_Union(b.geom))
FROM <houses> AS b
WHERE a.<id> <> b.<id>
AND ST_Contains(a.geom, b.geom)
)
WHERE EXISTS (
SELECT 1
FROM <buildings> AS b
WHERE a.<id> <> b.<id>
AND ST_Contains(a.geom, b.geom)
);
Oviously, this will alter data in your table - make sure that is what you want. Note that this fails in cases where the difference results in a multi-part Polygon - this would need to get handled either by changing the geom
column definition, or by deleting the originals and reinserting single-part Polygons from a similar query.
To DELETE
those that are contained by others (courtyards):
DELETE --test with SELECT *
FROM <houses> AS a
WHERE EXISTS (
SELECT 1
FROM <buildings> AS b
WHERE a.<id> <> b.<id>
AND ST_Within(a.geom, b.geom)
);
Obviously, this will erase data from your table - so make sure that is what you want before running this query; test a DELETE
by running a SELECT *
instead and make sure the resulting rows are what you intend to delete!
Finally, to CREATE
a new TABLE
with only those that contain others, with the footprint of those contained removed:
CREATE TABLE <houses_with_holes> AS (
SELECT a.<id>,
a.<column_1>,
...
a.<column_n>,
ST_SetSRID(ST_Difference(a.geom, ST_Union(b.geom)), <SRID>)::GEOMETRY(POLYGON, <SRID>) AS geom
FROM <houses> AS a
JOIN <houses> AS b
ON a.<id> <> b.<id> AND ST_Contains(a.geom, b.geom)
GROUP BY
1, 2, ..., n
);
This is probably the best way to test the outcome. This query can also get improved to handle cases where the difference results in multi-part Polygons.
Best Answer
Since you said you get a group of intersecting polygons for each polygon you're interested in, you may want to create what is referred to as a "polygon overlay".
This isn't exactly what Adam's solution is doing. To see the difference, take a look at this picture of an ABC intersection:
I believe Adam's solution will create an "AB" polygon that covers both the area of "AB!C" and "ABC", as well as an "AC" polygon that covers "AC!B" and "ABC", and a "BC" polygon that is "BC!A" and "ABC". So the "AB", "AC", and "BC" output polygons would all overlap the "ABC" area.
A polygon overlay produces non-overlapping polygons, so AB!C would be one polygon and ABC would be one polygon.
Creating a polygon overlay in PostGIS is actually pretty straightforward.
There are basically three steps.
Step 1 is extract the linework [Note that I'm using the exterior ring of the polygon, it does get a little more complicated if you want to correctly handle holes]:
Step 2 is to "node" the linework (produce a node at every intersection). Some libraries like JTS have "Noder" classes you can use to do this, but in PostGIS the ST_Union function does it for you:
Step 3 is to create all the possible non-overlapping polygons that can come from all those lines, done by the ST_Polygonize function:
You could save the output of each of those steps into a temp table, or you can combine them all into a single statement:
I'm using ST_Dump because the output of ST_Polygonize is a geometry collection, and it is (usually) more convenient to have a table where each row is one of the polygons that makes up the polygon overlay.