Grid based intersection, what to expect ?

[Martin Davis]

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/8227ebf1/attachment-0001.htm>