ST_intersection returns erroneous polygons

Question

I am trying to perform a clip by attribute in PostGIS between two shapefile containing 3000 polygons. Both shapefiles are free of errors (checked with topology checker), both shapefile are using EPSG:3857.

The following query works well for 99.9% of polygons, but for some it creates incomplete polygons or misses some when loaded in QGIS. The polygons with errors do have the same name used as a clipping attribute. I really don't understand why some polygons fail to be generated compared to others.

create table poly_20m_checked (gid serial PRIMARY KEY, polyID text, geometry geometry);
    insert into poly_20m_checked(polyID, geometry)
    select
    a.name as polyID,
    ST_Intersection(a.geom,b.geom) as geometry
    FROM donut20m a, symm20m b
    WHERE a.name=b.name
    Group by polyID, geometry

donut20m: symm20m: result (one polygon is missing):

Sample files that contains 4 polygons, among which one that returns an error (moved to make location anonymous): https://drive.google.com/open?id=0B7mqg7Gq1NGVVVVsbTZqMlVhZnc

The reason I am doing it in PostGIS is that I don't think there's a way to clip polygons in 2 shapefiles by using attribute in QGIS. What could be causing errors? I am very open to suggestion to get around the problem. I thought about spliting polygons by ID and then clip each pair of shapefiles but haven't really managed to code it.

EDIT: download for sample polygons. EDIT2: Tried to simplify geometries by 1m, which makes thing a bit better, but still bugy

Answer 1

Unfortunately, the ST_Intersection functions occasionally gives a GeometryCollection where you wouldn't expect it. This is the case with your data. When you inspect the output with ST_AsText you will see that the second donut from the top outputs both a polygon and a linestring. This linestring is bogus, since it has only 2 identical coordinates, but due to some rounding process, postgis (or likely GEOS) sees this as a seperate entity. This especially happens with polygons that are supposed to be precisely overlapping.

The solution is simple, always use a ST_Dump on the output of ST_Intersection. Example:

select
    a.name as polyID,
    (ST_Dump(ST_Intersection(a.geom,b.geom))).geom as geometry
    FROM tmp.sampledonut20m   a, tmp.samplesymm20m b
    WHERE a.name=b.name
    Group by polyID, geometry

Use a filter later on if you want to get rid of the linestrings. There is a nice example with a function from Dan Baston over here: PostGIS ST_Intersection of polygons can return lines