<div dir="ltr">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).<div><br></div><div>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).</div><div> </div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Apr 14, 2025 at 10:04 AM Martin Davis <<a href="mailto:mtnclimb@gmail.com" target="_blank">mtnclimb@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">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.<div><br></div><div>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.</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Mon, Apr 14, 2025 at 9:01 AM Sandro Santilli <<a href="mailto:strk@kbt.io" target="_blank">strk@kbt.io</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
I'm evaluating the use of the GEOS "grid-based" overlay operations<br>
available since version 3.9.0 as a way to reduce PostGIS Topology<br>
states in which vertices of incoming lines end up being closer than<br>
tolerance to segments of existing lines.<br>
<br>
Example of such problematic states:<br>
<br>
- <a href="https://trac.osgeo.org/postgis/ticket/5862" rel="noreferrer" target="_blank">https://trac.osgeo.org/postgis/ticket/5862</a><br>
- <a href="https://trac.osgeo.org/postgis/ticket/5786" rel="noreferrer" target="_blank">https://trac.osgeo.org/postgis/ticket/5786</a><br>
<br>
I thought using the "precise" overlay CAPI functions could reduce<br>
this problem by always finding intersections when facets are within<br>
the precision grid, so I tried the inputs of PostGIS ticket #5862:<br>
<br>
=# select ST_AsText(e1) e1, ST_AsText(e2) e2 from t5862_inputs;<br>
-[ RECORD 1 ]----------------------------------------------------------------------------------------------------------------------------------------------------<br>
e1 | LINESTRING(22.780107846871616 70.70515928614921,22.779899976871615 70.7046262461492)<br>
e2 | LINESTRING(22.79217056687162 70.70247684614921,22.779969266871618 70.70480392614921,22.780038556871617 70.7049816061492,22.796764346871615 70.7044482361492)<br>
<br>
Here you see how the internal points of the second line (e2) are at a<br>
distance which is smaller than 1e-14:<br>
<br>
=# select n, ST_Distance(e1, ST_PointN(e2, n)) from t5862_inputs, generate_series(1,4) n;<br>
n | st_distance<br>
---+------------------------<br>
1 | 0.012357374241807065<br>
2 | 4.855711461806118e-16<br>
3 | 2.8243441995579915e-15<br>
4 | 0.016671670112874255<br>
(4 rows)<br>
<br>
Asking GEOS 3.14.0dev-CAPI-1.20.0 to compute the intersection between<br>
the two lines with a precision of 1e-13 correctly returns those 2<br>
internal vertices:<br>
<br>
=# select ST_AsText( ST_Intersection(e1, e2, 1e-13) ) from t5862_inputs;<br>
st_astext<br>
---------------------------------------------------------------------------------<br>
LINESTRING(22.7800385568716 70.7049816061492,22.7799692668716 70.7048039261492)<br>
(1 row)<br>
<br>
But when a precision grid of 1e-14 is used, only one of those two<br>
internal points are returned (vertex 2, the closest). The other<br>
vertex, which was ~2.8e-15 distant, is not included in the<br>
intersection:<br>
<br>
=# select ST_AsText( ST_Intersection(e1, e2, 1e-14) ) from t5862_inputs;<br>
st_astext<br>
--------------------------------------------<br>
POINT(22.77996926687162 70.70480392614921)<br>
(1 row)<br>
<br>
How can this be explained ?<br>
<br>
--strk;<br>
<br>
Libre GIS consultant/developer 🎺<br>
<a href="https://strk.kbt.io/services.html" rel="noreferrer" target="_blank">https://strk.kbt.io/services.html</a><br>
</blockquote></div>
</blockquote></div>