ST_Intersection error

[ezimanyi]

I obtained the following error when I apply twice the ST_Intersection function

select st_astext(st_intersection(geometry 'Linestring(1 1,3 3.5,3 3)', geometry 'Polygon((0 0,0 3,3 3,3 0,0 0))'));
-- GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(1 1,2.6 3))

select st_astext(st_intersection(geometry 'Linestring(1 1,3 3.5,3 3)', geometry 'GEOMETRYCOLLECTION(POINT(3 3),LINESTRING(1 1,2.6 3))'));
-- MULTIPOINT(1 1,3 3)

select postgis_full_version();
POSTGIS="3.0.3 6660953" [EXTENSION] PGSQL="130" GEOS="3.9.1dev-CAPI-1.14.1" PROJ="7.0.0" LIBXML="2.9.4" LIBJSON="0.12.1"

The result above MULTIPOINT(1 1,3 3) is an error. I have tested also with 2.5.5 and GEOS 3.8 and the result is the same.

[mdavis]

This is not a bug, it is actually expected behaviour. The reason is because the floating point representation used internally for numbers is not able to precisely represent certain decimal numbers (even apparently simple ones like 2.6). The usual term for this is "roundoff error".

Simplifying the example makes this more clear:

SELECT ST_AsTExt(ST_Intersection('LINESTRING(1 1,3 3.5)','LINESTRING(1 1,2.6 3)'));
-- POINT(1 1)

This is actually the correct result, since the point (2.6 3) cannot be represented exactly internally, and so it does not lie exactly on the first line. This can be seen from:

SELECT ST_Intersects('LINESTRING (1 1, 3 3.5)', 'POINT (2.6 3)');
-- f

In this case the intersection and intersects operations are consistent. They just don't match what you think they should be. But you are computing in exact decimal arithmetic!

One way to mitigate this issue is to use the new ST_Intersection function that takes a gridsize. Using a small gridsize will cause this operation to return the "expected" result, since it causes the point to snap to the line.

Incidentally, there are also situations where the operations themeselves are internally inconsistent, also due to roundoff error. This is a classic example:

WITH a AS (
  SELECT 'LINESTRING(0.2 0.2, 1.2 1.2)'::geometry AS l1,
         'LINESTRING(1.2 0.2, 0.2 1.2)'::geometry AS l2
)
SELECT ST_AsText(ST_Intersection(l1, l2)) AS intersection_l1_l2,
       ST_Intersects(ST_Intersection(l1, l2), l1) AS l1_intersects_intersection,
       ST_Intersects(ST_Intersection(l1, l2), l2) AS l2_intersects_intersection
FROM a;  
---
-[ RECORD 1 ]--------------+---------------
intersection_l1_l2         | POINT(0.7 0.7)
l1_intersects_intersection | t
l2_intersects_intersection | f