Intersection of touching geometries returns a polygon

Question

I'm using Postgis 2.3. I have two touching geometries that I got from an overlay (ST_Difference) operation:

id |                                              st_astext    
 3 | MULTIPOLYGON(((0 -2,0 4,3 4,6 1,0 -2)),((6 6,8 8,8 4,6 6))) 
 4 | POLYGON((2 6,7 9,7 7,6 6,7 5,7 0,4 0,6 1,4.33333333333333 2.66666666666667,5 4,2 6))

Yet when I check their intersection, I get not only linestrings, but also a polygon.

select st_AsText(ST_Intersection(a.the_geom, b.the_geom)) from renovation a, renovation b where a.id=4 and b.id=3;

GEOMETRYCOLLECTION(LINESTRING(7 7,6 6),LINESTRING(6 6,7 5),LINESTRING(4 0,6 1),POLYGON((6 1,4.33333333333333 2.66666666666667,4.33333333333333 2.66666666666667,6 1)))

This polygon is very weird; it's got two equal vertices, yet it is valid, and it has a non-zero area.

How does this happen, and how can I prevent that from happening or filter this unnecessary polygon out?

Answer 1

If you use the WKT presentations of the geometries the intersection returns only multilinestring

select ST_AsText(
 ST_Intersection(
 ST_GeomFromText(
 'MULTIPOLYGON(((0 -2,0 4,3 4,6 1,0 -2)),((6 6,8 8,8 4,6 6)))'),
 ST_GeomFromText(
 'POLYGON((2 6,7 9,7 7,6 6,7 5,7 0,4 0,6 1,4.33333333333333 2.66666666666667,5 4,2 6))')
 )
 );

"MULTILINESTRING((6 1,4 0),(6 6,7 7),(7 5,6 6))"

I could not re-produce your result even when geometries were stored into tables

create table intersection (id int, geom geometry);
create table intersection2 (id int, geom geometry);

insert into intersection values (1,ST_GeomFromText('MULTIPOLYGON(((0 -2,0 4,3 4,6 1,0 -2)),((6 6,8 8,8 4,6 6)))'));
insert into intersection2 values (1,ST_GeomFromText('POLYGON((2 6,7 9,7 7,6 6,7 5,7 0,4 0,6 1,4.33333333333333 2.66666666666667,5 4,2 6))'));

select ST_AsText(
ST_Intersection(a.geom, b.geom))
from intersection a, intersection2 b;

"MULTILINESTRING((6 1,4 0),(6 6,7 7),(7 5,6 6))"

The problem is probably in the inaccuracy of floating point computing. For example decimal number 2.66666666666667 cannot be presented accurately as a double precision number. The converter at https://www.exploringbinary.com/floating-point-converter/ gives 2.666666666666670071350608850480057299137115478515625as a result. It is quite possible that the vertices returned by ST_Difference are not exactly the same and therefore ST_Intersection can find also minimal polygons. Another issue is that the polygon member in the geometry collection is invalid.

A common approach to accept some tolerance in the location of vertices is to read the geometries through ST_SnapToGrid http://postgis.net/docs/manual-2.0/ST_SnapToGrid.html.