4

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: donut symm20m: symm result (one polygon is missing): result

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

3

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

  • 1
    you can also use the st_collectionextract – ziggy Jun 7 '17 at 0:32
  • 1
    and isnt group by a geometry not the right way to go? – ziggy Jun 7 '17 at 0:32
  • Excellent, thank you for your answer. it's working perfectly now on the larger dataset. It does make some polygon as multipolygon but a little dissolve by attribute in QGIS resolves that nicely. Is there a way to do that in PostGIS to save an import/export step? – Leo Jun 7 '17 at 6:52
  • I guess multipolygon is correct in some cases, since your donut is being cut into more than 1 piece. Or do you mean 'more than 1 polygon' instead of 'multipolygon' ? – tilt Jun 7 '17 at 8:07
  • Yes sorry, I do mean more than one polygon – Leo Jun 8 '17 at 3:11

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.