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

Improve this question
2
  • I think it will be hard to debug this without your sharing the actual geometries donut20m and symm20m. Commented Jun 5, 2017 at 14:25
  • Good idea, here is a sample that contains 4 polygons, among which one that returns an error: drive.google.com/open?id=0B7mqg7Gq1NGVVVVsbTZqMlVhZnc Commented Jun 6, 2017 at 2:41

1 Answer 1

Reset to default
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

Improve this answer
6
  • 1
    you can also use the st_collectionextract Commented Jun 7, 2017 at 0:32
  • 1
    and isnt group by a geometry not the right way to go? Commented Jun 7, 2017 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? Commented Jun 7, 2017 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' ? Commented Jun 7, 2017 at 8:07
  • Yes sorry, I do mean more than one polygon Commented Jun 8, 2017 at 3:11

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.