changeset 2373:4aabbf324e55

Fix cutting of axis segment between two distance marks ST_LineSubstring() requires it's second argument to be smaller than the third, which is now ensured.
author Tom Gottfried <tom@intevation.de>
date Thu, 21 Feb 2019 15:56:07 +0100
parents 20e4efa64320
children 6efd7ecd3a7d
files schema/isrs_functions.sql
diffstat 1 files changed, 17 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/schema/isrs_functions.sql	Thu Feb 21 15:51:56 2019 +0100
+++ b/schema/isrs_functions.sql	Thu Feb 21 15:56:07 2019 +0100
@@ -4,7 +4,7 @@
 -- SPDX-License-Identifier: AGPL-3.0-or-later
 -- License-Filename: LICENSES/AGPL-3.0.txt
 
--- Copyright (C) 2018 by via donau
+-- Copyright (C) 2018, 2019 by via donau
 --   – Österreichische Wasserstraßen-Gesellschaft mbH
 -- Software engineering by Intevation GmbH
 
@@ -28,18 +28,14 @@
 AS $$
     WITH
         -- Get coordinates of location codes
-        from_geog AS (
+        points_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)),
+                WHERE location_code = lower(stretch)
+                    OR 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),
+            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(
@@ -47,24 +43,19 @@
                 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 (
+        points 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),
+                    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,
-                    (SELECT ST_Collect(from_point.geom, to_point.geom) AS pts
-                        FROM from_point, to_point) AS points
+                    (SELECT ST_Collect(points.geom) AS pts
+                        FROM points) AS points
                 WHERE ST_Covers(ST_Buffer(lines.line, 0.0001), points.pts)),
         axis_substring AS (
             -- Use linear referencing to clip axis between distance marks.
@@ -76,11 +67,12 @@
             -- resulting line leads to unexpected results of the buffer with
             -- endcap=flat in the CTE below.
             SELECT ST_SimplifyPreserveTopology(ST_LineSubstring(
-                    axis_segment.line,
-                    ST_LineLocatePoint(axis_segment.line, from_point.geom),
-                    ST_LineLocatePoint(axis_segment.line, to_point.geom)
-                ), 0.0001) AS line
-            FROM axis_segment, from_point, to_point),
+                        axis_segment.line, min(fractions.f), max(fractions.f)),
+                    0.0001) AS line
+            FROM axis_segment, LATERAL (
+                SELECT ST_LineLocatePoint(axis_segment.line, points.geom) AS f
+                    FROM points) AS fractions
+            GROUP BY axis_segment.line),
         area_subset AS (
             -- In case area is a multipolygon, process the union of those
             -- polygons, which intersect with the axis. The union is to avoid