8

When clipping the polygons of one table with polygons in another, ST_Intersection can return a set of results that can be handled with ST_Dump. The returned multiple geometries are not necessarily ST_Polygon but also ST_LineString (probably also a point). So when running a query

INSERT INTO c (geom)
  (SELECT (ST_Dump(ST_Intersection(a.geom,b.geom))).geom
  FROM a INNER JOIN b ON ST_Intersects(a.geom, b.geom));

when trying to populate the table "c" with clipped polygons, it fails with ERROR: Geometry type (LineString) does not match column type (Polygon)

I did another nested SELECT statement so only the polygon geometries came through, like:

INSERT INTO c (geom)
  (SELECT geom FROM
    (SELECT (ST_Dump(ST_Intersection(a.geom,b.geom))).geom
    FROM a INNER JOIN b ON ST_Intersects(a.geom, b.geom))) AS cl
    WHERE ST_GeometryType(cl.geom)='ST_Polygon');

As this is a bit cumbersome I am wondering if there is a more elegant solution for dropping invalid geometries?

  • BTW, I put <!-- language: lang-sql --> before the code blocks but there is still no highligting. Any hints for a noob? – Robert Špendl Mar 10 '14 at 23:14
  • To format code simply paste it in, select it and use the Code Sample button {} above the Question editing window. – PolyGeo Mar 10 '14 at 23:30
  • Hm, I did that, the code is marked as "code", but the keywords are still not colored. – Robert Špendl Mar 11 '14 at 0:01
  • I don't normally try to prettify SQL code here but I just found a Meta posting that looks useful at meta.stackexchange.com/questions/90521/… - have you left a "blank line between the comment and the code"? – PolyGeo Mar 11 '14 at 0:32
  • An intersection of two geometries may yield any type of geometry, so I think your current solution is the most direct. – Mike T Mar 11 '14 at 4:18
8

This might be a good spot to use an SQL-language function. Here's a quick one that should work for this situation:

CREATE OR REPLACE FUNCTION PolygonalIntersection(a geometry, b geometry)
RETURNS geometry AS $$
SELECT ST_Collect(geom)
FROM 
(SELECT (ST_Dump(ST_Intersection(a, b))).geom 
UNION ALL
-- union in an empty polygon so we get an 
-- empty geometry instead of NULL if there
-- is are no polygons in the intersection
SELECT ST_GeomFromText('POLYGON EMPTY')) SQ
WHERE ST_GeometryType(geom) = 'ST_Polygon';
$$ LANGUAGE SQL;

This will retain the polygonal components of an intersection, but throw away everything else. It always returns a MultiPolygon, even if you have one or no components.

WITH 
      square   as (SELECT ST_GeomFromText('POLYGON ((0 0, 0  1,  1  1,  1  0, 0 0))') AS geom),
biggersquare   as (SELECT ST_GeomFromText('POLYGON ((0 0, 0 10, 10 10, 10  0, 0 0))') AS geom),
adjacentsquare as (SELECT ST_GeomFromText('POLYGON ((0 0, 1  0,  1 -1, -1 -1, 0 0))') AS geom)   

SELECT ST_AsText(PolygonalIntersection(square.geom, biggersquare.geom))
  FROM square, biggersquare;
--"MULTIPOLYGON(((0 0,0 1,1 1,1 0,0 0)))"

SELECT ST_AsText(PolygonalIntersection(square.geom, adjacentsquare.geom))
  FROM square, adjacentsquare;
--"MULTIPOLYGON(EMPTY)"
1

Very good answer from @dbaston. However, returning an empty geometry instead of null might cause issues because returned empty geometry has no srid. St_Intersection might return MultiPolygon as well. This updated function was really useful for me:

CREATE OR REPLACE FUNCTION PolygonalIntersection(a geometry, b geometry)
RETURNS geometry AS $$
SELECT ST_Collect(geom)
FROM 
(SELECT (ST_Dump(ST_Intersection(a, b))).geom 
) SQ
WHERE ST_GeometryType(geom) = 'ST_Polygon' OR ST_GeometryType(geom) = 'ST_MultiPolygon';
$$ LANGUAGE SQL;

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.