Grid based intersection, what to expect ?

[Sandro Santilli]

I'm evaluating the use of the GEOS "grid-based" overlay operations
available since version 3.9.0 as a way to reduce PostGIS Topology
states in which vertices of incoming lines end up being closer than
tolerance to segments of existing lines.

Example of such problematic states:

  - https://trac.osgeo.org/postgis/ticket/5862
  - https://trac.osgeo.org/postgis/ticket/5786

I thought using the "precise" overlay CAPI functions could reduce
this problem by always finding intersections when facets are within
the precision grid, so I tried the inputs of PostGIS ticket #5862:

  =# select ST_AsText(e1) e1, ST_AsText(e2) e2 from t5862_inputs;
  -[ RECORD 1 ]----------------------------------------------------------------------------------------------------------------------------------------------------
  e1 | LINESTRING(22.780107846871616 70.70515928614921,22.779899976871615 70.7046262461492)
  e2 | LINESTRING(22.79217056687162 70.70247684614921,22.779969266871618 70.70480392614921,22.780038556871617 70.7049816061492,22.796764346871615 70.7044482361492)

Here you see how the internal points of the second line (e2) are at a
distance which is smaller than 1e-14:

  =# select n, ST_Distance(e1, ST_PointN(e2, n)) from t5862_inputs, generate_series(1,4) n;
   n |      st_distance
  ---+------------------------
   1 |   0.012357374241807065
   2 |  4.855711461806118e-16
   3 | 2.8243441995579915e-15
   4 |   0.016671670112874255
  (4 rows)

Asking GEOS 3.14.0dev-CAPI-1.20.0 to compute the intersection between
the two lines with a precision of 1e-13 correctly returns those 2
internal vertices:

  =# select ST_AsText( ST_Intersection(e1, e2, 1e-13) ) from t5862_inputs;
                                      st_astext
  ---------------------------------------------------------------------------------
   LINESTRING(22.7800385568716 70.7049816061492,22.7799692668716 70.7048039261492)
  (1 row)

But when a precision grid of 1e-14 is used, only one of those two
internal points are returned (vertex 2, the closest). The other
vertex, which was ~2.8e-15 distant, is not included in the
intersection:

  =# select ST_AsText( ST_Intersection(e1, e2, 1e-14) ) from t5862_inputs;
                   st_astext
  --------------------------------------------
   POINT(22.77996926687162 70.70480392614921)
  (1 row)

How can this be explained ?

--strk;

  Libre GIS consultant/developer 🎺
  https://strk.kbt.io/services.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 659 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/geos-devel/attachments/20250414/d4668848/attachment.sig>