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
signature.asc
Description: PGP signature
