changeset 1629:9d51f022b8ee

Introduce SQL function to clip an area to a stretch
author Tom Gottfried <tom@intevation.de>
date Wed, 19 Dec 2018 17:59:31 +0100
parents 8ad51ad5a9ee
children f6218f11062a
files schema/install-db.sh schema/isrs_functions.sql schema/isrs_tests.sql schema/run_tests.sh schema/tap_tests_data.sql
diffstat 5 files changed, 134 insertions(+), 1 deletions(-) [+]
line wrap: on
line diff
--- a/schema/install-db.sh	Wed Dec 19 16:48:00 2018 +0100
+++ b/schema/install-db.sh	Wed Dec 19 17:59:31 2018 +0100
@@ -125,6 +125,7 @@
        -f "$BASEDIR/geonames.sql" \
        -f "$BASEDIR/manage_users.sql" \
        -f "$BASEDIR/auth.sql" \
+       -f "$BASEDIR/isrs_functions.sql" \
        -f "$BASEDIR/default_sysconfig.sql"
 
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/schema/isrs_functions.sql	Wed Dec 19 17:59:31 2018 +0100
@@ -0,0 +1,82 @@
+-- Clip an area to a stretch given by a pair of ISRS location codes.
+-- Uses the table waterway.distance_marks_virtual to map ISRS location codes
+-- to their geo-location and the table waterway.waterway_axis to retrieve
+-- perpendicular direction at these geo-locations.
+-- Distance marks are assumed to be near the axis and the area passed as
+-- argument is assumed to intersect with the axis
+-- (use e.g. waterway area or fairway dimensions).
+CREATE OR REPLACE FUNCTION ISRSrange_area(
+    stretch isrsrange,
+    area geometry
+) RETURNS geometry
+AS $$
+    WITH
+        -- Get coordinates of location codes
+        from_geog AS (
+            SELECT geom FROM waterway.distance_marks_virtual
+                WHERE location_code = lower(stretch)),
+        to_geog AS (
+            SELECT geom FROM waterway.distance_marks_virtual
+                WHERE location_code = upper(stretch)),
+        utm_zone AS (
+            -- Find best matchting UTM zone
+            SELECT best_utm(ST_Collect(
+                    from_geog.geom::geometry,
+                    to_geog.geom::geometry)) AS z
+                FROM from_geog, to_geog),
+        axis AS (
+            -- Transform and sew together contiguous axis chunks
+            SELECT ST_LineMerge(ST_Collect(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
+        from_point AS (
+            SELECT ST_ClosestPoint(
+                    wtwaxs,
+                    ST_Transform(from_geog.geom::geometry, z)) AS geom
+                FROM axis, from_geog, utm_zone),
+        to_point AS (
+            SELECT ST_ClosestPoint(
+                    wtwaxs,
+                    ST_Transform(to_geog.geom::geometry, z)) AS geom
+                FROM axis, to_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,
+                    (SELECT ST_Collect(from_point.geom, to_point.geom) AS pts
+                        FROM from_point, to_point) AS points
+                WHERE ST_Covers(lines.line, points.pts)),
+        axis_substring AS (
+            -- Use linear referencing to clip axis between distance marks
+            SELECT ST_LineSubstring(
+                    axis_segment.line,
+                    ST_LineLocatePoint(axis_segment.line, from_point.geom),
+                    ST_LineLocatePoint(axis_segment.line, to_point.geom)
+                ) AS line
+            FROM axis_segment, from_point, to_point),
+        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,
+                            ST_Transform(area, z)),
+                        'endcap=flat'),
+                    ST_Transform(area, z)))).geom
+                FROM axis_substring, utm_zone)
+        -- From the polygons returned by the last CTE, select only the one
+        -- around the clipped axis
+        SELECT ST_Transform(range_area.geom, ST_SRID(area))
+            FROM axis_substring, range_area
+            WHERE ST_Intersects(range_area.geom, axis_substring.line)
+    $$
+    LANGUAGE sql;
--- a/schema/isrs_tests.sql	Wed Dec 19 16:48:00 2018 +0100
+++ b/schema/isrs_tests.sql	Wed Dec 19 17:59:31 2018 +0100
@@ -28,3 +28,33 @@
     $$,
     22023, NULL,
     'ISRS text input needs to have correct length');
+
+SELECT ok(
+    ISRSrange_area(isrsrange(
+            ('AT', 'XXX', '00001', '00000', 0)::isrs,
+            ('AT', 'XXX', '00001', '00000', 1)::isrs),
+        ST_SetSRID('POLYGON((0 1, 0 2, 1 2, 1 1, 0 1))'::geometry, 4326)
+    ) IS NULL,
+    'ISRSrange_area returns NULL, if given area does not intersect with axis');
+
+SELECT ok(
+    ST_DWithin(
+        (SELECT geom FROM waterway.distance_marks_virtual
+            WHERE location_code = ('AT', 'XXX', '00001', '00000', 0)::isrs),
+        ST_Boundary(ISRSrange_area(isrsrange(
+                ('AT', 'XXX', '00001', '00000', 0)::isrs,
+                ('AT', 'XXX', '00001', '00000', 1)::isrs),
+            ST_SetSRID('POLYGON((-1 1, 2 1, 2 -1, -1 -1, -1 1))'::geometry,
+                4326)))::geography,
+        1)
+    AND
+    ST_DWithin(
+        (SELECT geom FROM waterway.distance_marks_virtual
+            WHERE location_code = ('AT', 'XXX', '00001', '00000', 1)::isrs),
+        ST_Boundary(ISRSrange_area(isrsrange(
+                ('AT', 'XXX', '00001', '00000', 0)::isrs,
+                ('AT', 'XXX', '00001', '00000', 1)::isrs),
+            ST_SetSRID('POLYGON((-1 1, 2 1, 2 -1, -1 -1, -1 1))'::geometry,
+                4326)))::geography,
+        1),
+    'Resulting polygon almost ST_Touches points corresponding to stretch');
--- a/schema/run_tests.sh	Wed Dec 19 16:48:00 2018 +0100
+++ b/schema/run_tests.sh	Wed Dec 19 17:59:31 2018 +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(45)' \
+    -c 'SELECT plan(47)' \
     -f isrs_tests.sql \
     -f auth_tests.sql \
     -f manage_users_tests.sql \
--- a/schema/tap_tests_data.sql	Wed Dec 19 16:48:00 2018 +0100
+++ b/schema/tap_tests_data.sql	Wed Dec 19 17:59:31 2018 +0100
@@ -64,6 +64,26 @@
         1, 'depth', 'testorganization', true
     );
 
+INSERT INTO waterway.distance_marks_virtual VALUES (
+    ('AT', 'XXX', '00001', '00000', 0)::isrs,
+    ST_SetSRID('POINT(0 0)'::geometry, 4326),
+    'someENC'
+), (
+    ('AT', 'XXX', '00001', '00000', 1)::isrs,
+    ST_SetSRID('POINT(1 0)'::geometry, 4326),
+    'someENC'
+);
+
+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),
+        4326),
+    'testriver'
+), (
+    ST_SetSRID('LINESTRING(0.5 0.5, 1 1)'::geometry, 4326),
+    'testriver'
+);
+
 INSERT INTO users.templates (template_name, template_data)
     VALUES ('AT', '\x'), ('RO', '\x');
 INSERT INTO users.user_templates