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 <mtncl...@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 <mtncl...@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 <s...@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 >>> >>
