2

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?

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.

  • Thanks a lot! I actually tried to understand why you cannot reproduce the problem. Eventually I came to the conclusion that if I run the following query, everything works: UPDATE renovation SET the_geom=ST_GeomFromText(ST_AsText(the_geom)); Does it even make sense? – Nitzan Weissman Nov 28 '18 at 12:43
  • 1
    The AsText - GeomFromText conversion does not round-trip always because of the similar reason that all native binary geometries can't be expressed accurately as WKT. This time the slightly changed coordinates work for you but sometimes the effect of the change is the opposite. I would recommend ST_SnapToGrid with for example 0.001 m tolerance. The ST_AsText roundtrip actually does more or less the same with an accuracy of 15 decimals by default postgis.net/docs/ST_AsText.html. – user30184 Nov 28 '18 at 13:17

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.