Grid based intersection, what to expect ?

[Martin Davis]

I suspect that the anomalies you're seeing are because the precision grid
size is so small.  A precision of 1e-14 combined with input values which
are on the order of 1e2 means that you are asking for 16 digits of decimal
precision, which is at or over the precision that can be represented using
double-precision FP.

The precision-based overlay ops are only able to support "reasonable"
precision for a given data magnitude.  I would say 14 digits of total
precision is the very most that can be evaluated, and using 12 or even 10
is safer.  Note that a precision of 10 digits allows representing earth
coordinates to about 1 mm, which should be plenty for real-world cases.

On Mon, Apr 14, 2025 at 9:01 AM Sandro Santilli <strk at kbt.io> wrote:

>
> 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 --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/geos-devel/attachments/20250414/3048b4ed/attachment.htm>