Grid based intersection, what to expect ?

Mon Apr 14 21:13:01 PDT 2025

One more option:  another approach used by OverlayNG to improve robustness
is Snapped Noding.   This is different to Snap-Rounding.  Instead of
rounding all coordinates to a precision grid, Snapping Noding uses a
distance tolerance and snaps vertices to lines or vertices within the
tolerance.  This improves robustness, but doesn't change vertices unless
required. (Note it's still not *fully* robust - for that, full
snap-rounding is required).

This is implemented in GEOS, but is not expose in the C API (and hence not
in PostGIS).  It could be exposed if required (with some coding, of course).

On Mon, Apr 14, 2025 at 10:12 AM Martin Davis <mtnclimb at gmail.com> wrote:

> Also, the input values have not been rounded to the requested precision.
> In general this should be done to use the overlay precision ops correctly.
> (The data will be rounded internally, but the results may not match what
> might be expected from the original inputs).
>
> The goal of the precision-based overlay ops is to support operations on a
> dataset which is stored using a given precision.  The input data
> should either be provided in the required precision, or reduced to it (e.g.
> using GEOSGeom_setPrecision).  All operations on the geometry should be
> done in a way which maintains that precision.  (One gap in the current API
> is that the spatial predicates do not support a precision value.  I'm
> hoping to provide that at some point).
>
>
> On Mon, Apr 14, 2025 at 10:04 AM Martin Davis <mtnclimb at gmail.com> wrote:
>
>> 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/6d18bb92/attachment.htm>