Serious rounding issues

[postgispaul]

I have 3 polygons (in 4326 projection), and when I calculate st_difference on those, there appears to be great rounding-issues.

Query which results in erroneous polygon: select st_difference(t1.geom,st_collect(t2.geom)) from error as t1 join error t2 on true where t1.id=1 and t2.id in (600,1351) group by t1.geom

Screenshot is in attachment. Green filled polygon is the result of said query, black outlines are the contours of the polygons in the table.

Table should be filled with the contents of polygons.txt from attachment

Versions: "POSTGIS="2.2.0 r14208" GEOS="3.5.0-CAPI-1.9.0 r4090" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.0.1, released 2015/09/15 GDAL_DATA not found" LIBXML="2.7.8" LIBJSON="0.12" TOPOLOGY RASTER" "PostgreSQL 9.4.1, compiled by Visual C++ build 1800, 64-bit"

[postgispaul]

sources

[postgispaul]

Just some extra information: Geometry is in 4326 (WGS84), which does seem to have some drawbacks, but the results I'm seeing here are only when this particular polygon (1351) is added to the mix. I did have a lot of other polygons (about 3000 of them) in that area, but none of them produced those rounded/truncated results. (Was quite a struggle to pinpoint this issue on that specific polygon)

[dbaston]

I tried to reproduce this issue in JTS, which threw the familiar "This method does not support GeometryCollection arguments" exception. It seems that the result of ST_Collect on polygons 600 and 1351 is not a valid MultiPolygon (the two polygons overlap).

The comments in GEOS suggest that it should be blocking this operation, as JTS does, in geometry.h: ​https://github.com/libgeos/libgeos/blob/597a93df71c763205f5a021d35594b14db791a3f/include/geos/geom/Geometry.h#L626

For this particular case, the rounding problem can be resolved by using ST_Union instead of ST_Collect.

[robe]

dbaston,

Thanks so much for looking into this. I hadn't even noticed that he was using ST_Collect instead of ST_Union. Yah ST_Collect produces invalid polygons and all bets are off when you feed invalid geometries to GEOS.

[robe]

postgispaul,

Can you try changing your ST_Collect to ST_Union as dbaston suggested. I'm going to close this out. If you still have issue feel free to reopen.

[dbaston]

robe, is there a GEOS issue here? Based on the code comments, it seems like GEOS should be rejecting these inputs, but it's not.

[postgispaul]

Thanks for investigating!

I can vaguely remember that an older version of postgis (somewhere in the latest few of the 2.1 series) gave me an error about invalid polygons, but that was in one of my many experiments with st_collect/st_union/st_makevalid combo's.

Will clarify tomorrow when I have access to the database.

FWIW, if the conclusion of dbaston is correct, I think giving an error has my preference over giving incorrect results

[robe]

Yah me too unfortunately that may be tricky to do since the error would need to bubble up from GEOS and would not be isolated to this function so hard to tell where to putit as some functions do work okay with invalid polygons like ST_MakeValid

[postgispaul]

Just tested, and st_union works indeed correctly (on both Postgis 2.2.0 and 2.1.7)

[dbaston]

I was a bit unclear in my earlier comment. Here is some info on the output of ST_Collect:

SELECT ST_GeometryType(ST_Collect(geom)), ST_IsValid(ST_Collect(geom)) FROM bugreport WHERE id in (600, 1351);
    st_geometrytype    | st_isvalid 
-----------------------+------------
 ST_GeometryCollection | t
(1 row)

So there is no invalid geometry involved. What I meant to say what the output of ST_Collect, if it were a MultiPolygon (which it is not), would be invalid. I think the issue here is how ST_Difference should be handling GeometryCollections, not how it should be handling invalid geometries.

The PostGIS docs for ST_Difference say "Do not call with a GeometryCollection as an argument", without any further elaboration, and the GEOS headers for Geometry::difference say that it does not support GeometryCollections. So it seems that either:

• PostGIS should be calling errorIfGeometryCollection(), before passing this on to GEOS, or
• GEOS should be throwing an exception, which PostGIS should handle and report an error to the user. PostGIS is already handling errors from GEOS correctly in this case, so it seems that GEOS is not throwing an exception.