Mercurial > gemma

changeset 2398:8481e6266691

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix corner case in area generation from stretch In case orthogonal lines of arbitrary length at the ends of the axis between two given distance marks cross the axis (e.g. because it describes a tight turn), the result of ST_Buffer with endcap=flat is not well defined. Test data have been amended to include such a case as well as appropriate data to test correctness of a fix. The fix includes construction of the afore mentioned orthogonal lines, splitting the given area by these lines and constructing the end result from those parts of the splitted area that intersect with the clipped axis. Due to numerical inaccuracy, the parts might overlap slightly and eventually cross the clipped axis where they should only touch. Therefore, a small buffer was introduced before testing intersection and intersecting parts are dissolved using ST_Union (which is necessary to revert splits inside the resulting polygon anyhow).
author Tom Gottfried <tom@intevation.de>
date Wed, 27 Feb 2019 18:00:35 +0100
parents b95234702ee9
children fd248ede0251
files schema/isrs_functions.sql schema/isrs_tests.sql schema/tap_tests_data.sql
diffstat 3 files changed, 49 insertions(+), 22 deletions(-) [+]
line wrap: on
line diff
--- a/schema/isrs_functions.sql	Wed Feb 27 16:26:15 2019 +0100
+++ b/schema/isrs_functions.sql	Wed Feb 27 18:00:35 2019 +0100
@@ -82,25 +82,29 @@
                     ST_Dump(ST_Transform(area, z)) AS a_dmp
                 WHERE ST_Intersects(a_dmp.geom, axis_substring.line)
             ),
+        rotated_ends AS (
+            SELECT ST_Collect(ST_Scale(
+                    ST_Translate(e,
+                        (ST_X(p1) - ST_X(p2)) / 2,
+                        (ST_Y(p1) - ST_Y(p2)) / 2),
+                    ST_Point(d, d), p1)) AS blade
+                FROM axis_substring, area_subset,
+                    LATERAL (SELECT i, ST_PointN(line, i) AS p1
+                        FROM (VALUES (1), (-1)) AS idx (i)) AS ep,
+                    ST_Rotate(ST_PointN(line, i*2), pi()/2, p1) AS ep2 (p2),
+                    ST_Makeline(p1, p2) AS e (e),
+                    LATERAL (SELECT (ST_MaxDistance(p1, area) / ST_Length(e))
+                            * 2) AS d (d)),
         range_area AS (
-            -- Create a buffer around the clipped axis, as large as it could
-            -- potentially be intersecting with the area polygon that
-            -- intersects with the clipped axis. Get the intersection of that
-            -- buffer with the area polygon, which can potentially
-            -- be a multipolygon.
-            SELECT (ST_Dump(ST_Intersection(
-                    ST_Buffer(
-                        axis_substring.line,
-                        ST_MaxDistance(
-                            axis_substring.line,
-                            area_subset.area),
-                        'endcap=flat'),
-                    area_subset.area))).geom
-                FROM axis_substring, area_subset)
+            -- Split area by orthogonal lines at the ends of the clipped axis
+            SELECT (ST_Dump(ST_CollectionExtract(
+                    ST_Split(area, blade), 3))).geom
+                FROM area_subset, rotated_ends)
         -- From the polygons returned by the last CTE, select only those
         -- around the clipped axis
-        SELECT ST_Collect(ST_Transform(range_area.geom, ST_SRID(area)))
+        SELECT ST_Multi(ST_Union(ST_Transform(range_area.geom, ST_SRID(area))))
             FROM axis_substring, range_area
-            WHERE ST_Intersects(range_area.geom, axis_substring.line)
+            WHERE ST_Intersects(ST_Buffer(range_area.geom, -0.0001),
+                axis_substring.line)
     $$
     LANGUAGE sql;
--- a/schema/isrs_tests.sql	Wed Feb 27 16:26:15 2019 +0100
+++ b/schema/isrs_tests.sql	Wed Feb 27 18:00:35 2019 +0100
@@ -42,7 +42,6 @@
     ) IS NULL,
     'ISRSrange_area returns NULL, if given area does not intersect with axis');
 
-\set test_area 'POLYGON((-1 1, 2 1, 2 -1, -1 -1, -1 1))'
 SELECT ok(
     ST_DWithin(
         (SELECT geom FROM waterway.distance_marks_virtual
@@ -50,8 +49,8 @@
         ST_Boundary(ISRSrange_area(isrsrange(
                 ('AT', 'XXX', '00001', '00000', 0)::isrs,
                 ('AT', 'XXX', '00001', '00000', 1)::isrs),
-            ST_SetSRID(:'test_area'::geometry,
-                4326)))::geography,
+            (SELECT ST_Collect(CAST(area AS geometry))
+                FROM waterway.waterway_area))),
         1)
     AND
     ST_DWithin(
@@ -60,11 +59,12 @@
         ST_Boundary(ISRSrange_area(isrsrange(
                 ('AT', 'XXX', '00001', '00000', 0)::isrs,
                 ('AT', 'XXX', '00001', '00000', 1)::isrs),
-            ST_SetSRID(:'test_area'::geometry,
-                4326)))::geography,
+            (SELECT ST_Collect(CAST(area AS geometry))
+                FROM waterway.waterway_area))),
         1),
     'Resulting polygon almost ST_Touches points corresponding to stretch');
 
+\set test_area 'POLYGON((-1 1, 2 1, 2 -1, -1 -1, -1 1))'
 SELECT ok(
     2 = ST_NumGeometries(
         ISRSrange_area(
--- a/schema/tap_tests_data.sql	Wed Feb 27 16:26:15 2019 +0100
+++ b/schema/tap_tests_data.sql	Wed Feb 27 18:00:35 2019 +0100
@@ -80,7 +80,12 @@
 
 INSERT INTO waterway.waterway_axis (wtwaxs, objnam) VALUES (
     ST_SetSRID(ST_CurveToLine(
-        'CIRCULARSTRING(0 0, 0.5 0.5, 1 0, 1.5 -0.2, 2 0)'::geometry),
+        'CIRCULARSTRING(0 0, 0.5 0.5, 0.6 0.4)'),
+        4326),
+    'testriver'
+), (
+    ST_SetSRID(ST_CurveToLine(
+        'CIRCULARSTRING(0.6 0.4, 1 0, 2 0)'),
         4326),
     'testriver'
 ), (
@@ -88,6 +93,24 @@
     'testriver'
 );
 
+-- Simulate waterway area as non-intersecting buffers around axis
+WITH RECURSIVE
+buffer AS (
+    SELECT id, ST_Buffer(wtwaxs, 10000, 'endcap=flat')::geometry AS buf
+        FROM waterway.waterway_axis),
+cleaned AS (
+    (SELECT ARRAY[id] AS ids, buf AS cbuf, buf AS others
+        FROM buffer ORDER BY id FETCH FIRST ROW ONLY)
+    UNION
+    (SELECT ids || id,
+            ST_Difference(buf, others),
+            ST_Union(buf, others)
+        FROM cleaned, buffer
+        WHERE id <> ALL(ids)
+        ORDER BY id ASC, ids DESC
+        FETCH FIRST ROW ONLY))
+INSERT INTO waterway.waterway_area (area) SELECT cbuf FROM cleaned;
+
 INSERT INTO users.templates (template_name, country, template_data)
     VALUES ('AT', 'AT', '\x'), ('RO', 'RO', '\x');