CODEWITHSUNDEEP
postgresqlpostgispolygonspatialintersection
TheRealJimShady
TheRealJimShady

Reputation: 917

Why does st_intersection return non-polygons?

I have two polygon layers. I want to run st_intersection on them, to give the result of the areas where they overlap as a new layer. The new layer should contain the attributes from both input layers. I found this image which seems to illustrate my desired end results.

enter image description here

My two input layers are both polygons:

SELECT      st_geometrytype(geom),
            COUNT(*)
FROM        a
GROUP BY    st_geometrytype(geom)
-- Result is 1368 st_polygons

SELECT      st_geometrytype(geom),
            COUNT(*)
FROM        b
GROUP BY    st_geometrytype(geom)
-- Result is 539548 st_polygons

The query I run is as below:

SELECT      a.*,
            b.*,
            st_intersection(a.geom, b.geom) as geom
FROM        a,b
WHERE       st_intersects(a.geom, b.geom)

However in the result I get not just polygons (which I expect), but lines, points, multipolygons and geometry collections. I guess because some of my input polygons share points but not true intersections perhaps?

Grateful for some advice please on how to deal with this, whether my query is correct, anything I can do to improve performance etc. Thanks.

Upvotes: 1

Views: 1507

Answers (1)

RafDouglas C. Tommasi
RafDouglas C. Tommasi

Reputation: 156

ST_intersect returns several geometry types, depending on the relative topology.

For example, running ST_intersect on two adjacent polygons returns the common part of the shared boundary.

While it ouptuts a single table (as you can verify in pgadmin, for example), in the Browser swatch of QGIS it will be shown as multiple tables of different geometry types (for example: POLYGON, MULTIPOLY, LINE, and POINT) but (somewhat confusingly) with the same name.

Visually, you can tell them apart observing the accompaining icons on the left:

enter image description here

You can however select which type of geometry you want, for example by adding a WHERE filter with ST_Dimension:

SELECT      a.*,
            b.*,
            st_intersection(a.geom, b.geom) as geom
FROM        a,b
WHERE       st_intersects(a.geom, b.geom) 
            AND ST_Dimension(st_intersects(a.geom, b.geom)) = 2;

or, for performance sake, re-write it in a fashion similar to:

SELECT clipped.* 
FROM (
         SELECT a.id, b."fieldName",
             (ST_Dump(ST_Intersection(a.geom, b.geom))).geom AS geom
         FROM "public"."table_A_name" AS a INNER JOIN "public"."table_B_name" AS b   
         ON ST_Intersects(a.geom, b.geom)
     ) AS clipped
WHERE ST_Dimension("clipped"."geom") = 2;

The latter solution creates an anonymous temporary table, which allows ST_Intersection to run only once.

You might have noticed thath the trick is in ST_Dimension("clipped"."geom") = 2.

ST_Dimensions which filters the outputs from ST_Intersection so as to keep only polygons (which have a topological dimension of 2).

Upvotes: 2

Related Questions