changeset 2431:48495bd3081d

Construct stretch area between distance marks also from non-contiguous axis Before, both distance marks were required to be next to the same linestring of the waterway axis.
author Tom Gottfried <tom@intevation.de>
date Thu, 28 Feb 2019 19:44:38 +0100
parents 2a93a8649751
children 851c904416f6
files schema/isrs_functions.sql schema/isrs_tests.sql schema/run_tests.sh schema/tap_tests_data.sql
diffstat 4 files changed, 60 insertions(+), 13 deletions(-) [+]
line wrap: on
line diff
--- a/schema/isrs_functions.sql	Thu Feb 28 17:30:26 2019 +0100
+++ b/schema/isrs_functions.sql	Thu Feb 28 19:44:38 2019 +0100
@@ -26,7 +26,7 @@
     area geometry
 ) RETURNS geometry
 AS $$
-    WITH
+    WITH RECURSIVE
         -- Get coordinates of location codes
         points_geog AS (
             SELECT geom FROM waterway.distance_marks_virtual
@@ -37,9 +37,7 @@
             SELECT best_utm(ST_Collect(geom::geometry)) AS z
                 FROM points_geog),
         axis AS (
-            -- Transform and sew together contiguous axis chunks
-            SELECT ST_LineMerge(ST_Collect(ST_Transform(
-                    wtwaxs::geometry, z))) AS wtwaxs
+            SELECT id, ST_Transform(wtwaxs::geometry, z) AS wtwaxs
                 FROM waterway.waterway_axis, utm_zone),
         -- In order to guarantee the following ST_Covers to work,
         -- snap distance mark coordinates to axis
@@ -47,16 +45,49 @@
             SELECT ST_ClosestPoint(
                     wtwaxs,
                     ST_Transform(points_geog.geom::geometry, z)) AS geom
-                FROM axis, points_geog, utm_zone),
-        axis_segment AS (
-            -- select the contiguous axis on which distance marks lie
-            SELECT line
-                FROM (
-                    SELECT (ST_Dump(wtwaxs)).geom AS line
-                        FROM axis) AS lines,
+                FROM points_geog, utm_zone, (
+                    SELECT ST_Collect(wtwaxs) AS wtwaxs
+                        FROM axis) AS ax),
+        axis_snapped AS (
+            -- Iteratively connect non-contiguous axis chunks
+            -- to find the contiguous axis on which given distance marks lie
+            (SELECT ARRAY[id] AS ids, wtwaxs
+                FROM axis, points
+                WHERE ST_Intersects(
+                    ST_Buffer(axis.wtwaxs, 0.0001), points.geom)
+                FETCH FIRST ROW ONLY)
+            UNION
+            -- Connect endpoint of next linestring with closest
+            -- endpoint of merged linestring until a contiguous
+            -- linestring connecting both distance marks is build up
+            (SELECT refids || id,
+                    ST_LineMerge(ST_Collect(ARRAY(
+                        -- Linestring build up so far
+                        SELECT refgeom
+                        UNION
+                        -- Fill eventual gap
+                        SELECT ST_MakeLine(
+                            ST_ClosestPoint(
+                                ST_Boundary(refgeom), ST_Boundary(geom)),
+                            ST_ClosestPoint(
+                                ST_Boundary(geom), ST_Boundary(refgeom)))
+                        UNION
+                        -- Linestring to be added
+                        SELECT geom)))
+                FROM axis_snapped AS axis_snapped (refids, refgeom),
+                    axis AS axis (id, geom),
                     (SELECT ST_Collect(points.geom) AS pts
                         FROM points) AS points
-                WHERE ST_Covers(ST_Buffer(lines.line, 0.0001), points.pts)),
+                WHERE id <> ALL(refids)
+                    AND NOT ST_Covers(ST_Buffer(refgeom, 0.0001), points.pts)
+                ORDER BY ST_Distance(ST_Boundary(refgeom), ST_Boundary(geom))
+                FETCH FIRST ROW ONLY)),
+        axis_segment AS (
+            -- Fetch end result from snapping
+            SELECT wtwaxs AS line
+                FROM axis_snapped
+                WHERE array_length(ids, 1) = (
+                    SELECT max(array_length(ids, 1)) FROM axis_snapped)),
         axis_substring AS (
             -- Use linear referencing to clip axis between distance marks.
             -- Simplification is used to work-around the problem, that
--- a/schema/isrs_tests.sql	Thu Feb 28 17:30:26 2019 +0100
+++ b/schema/isrs_tests.sql	Thu Feb 28 19:44:38 2019 +0100
@@ -100,3 +100,12 @@
                             0)),
                 4326))),
     'Self-intersecting multipolygon leads to one polygon in result');
+
+SELECT ok(
+    ISRSrange_area(
+        isrsrange(
+            ('AT', 'XXX', '00001', '00000', 0)::isrs,
+            ('AT', 'XXX', '00001', '00000', 2)::isrs),
+        (SELECT ST_Collect(CAST(area AS geometry))
+            FROM waterway.waterway_area)) IS NOT NULL,
+    'Area generated from non-matching distance mark and non-contiguous axis');
--- a/schema/run_tests.sh	Thu Feb 28 17:30:26 2019 +0100
+++ b/schema/run_tests.sh	Thu Feb 28 19:44:38 2019 +0100
@@ -28,7 +28,7 @@
     -c 'SET client_min_messages TO WARNING' \
     -c "DROP ROLE IF EXISTS $TEST_ROLES" \
     -f tap_tests_data.sql \
-    -c 'SELECT plan(57)' \
+    -c 'SELECT plan(58)' \
     -f isrs_tests.sql \
     -f auth_tests.sql \
     -f manage_users_tests.sql \
--- a/schema/tap_tests_data.sql	Thu Feb 28 17:30:26 2019 +0100
+++ b/schema/tap_tests_data.sql	Thu Feb 28 19:44:38 2019 +0100
@@ -76,6 +76,10 @@
     ('AT', 'XXX', '00001', '00000', 1)::isrs,
     ST_SetSRID('POINT(1 0)'::geometry, 4326),
     'someENC'
+), (
+    ('AT', 'XXX', '00001', '00000', 2)::isrs,
+    ST_SetSRID('POINT(1.6 0)'::geometry, 4326),
+    'someENC'
 );
 
 INSERT INTO waterway.waterway_axis (wtwaxs, objnam) VALUES (
@@ -89,6 +93,9 @@
 ), (
     ST_SetSRID('LINESTRING(0.5 0.5, 1 1)'::geometry, 4326),
     'testriver'
+), (
+    ST_SetSRID('LINESTRING(1.5 0.1, 2 0)'::geometry, 4326),
+    'testriver'
 );
 
 -- Simulate waterway area as non-intersecting buffers around axis