# HG changeset patch
# User Tom Gottfried <tom@intevation.de>
# Date 1554883327 -7200
# Node ID 524f1382558a7aa8a4e5734a34c4acfe763f1407
# Parent  7ee9a3ef90d4b39bac694848f9e8adfa471fbc27
Make snapping of many axis segments faster

Since a CTE is materialized anyhow, make the one for the UTM projected
axis a temporary table, extended by the precalculated geometry
boundaries and a spatial index on these, which then can be used in
geometry comparisons within the recursive CTE. This makes snapping
several hundreds of segments faster by more than an order of magnitude
at the expense of longer runtime for snapping only very few segments,
because of the additional cost for index creation. This could be
further optimized by first checking if there is a direct match without
any snapping, but since the function is currently only used in import
processes, a runtime of a few seconds seems acceptable.

diff -r 7ee9a3ef90d4 -r 524f1382558a schema/isrs_functions.sql
--- a/schema/isrs_functions.sql	Tue Apr 09 18:48:16 2019 +0200
+++ b/schema/isrs_functions.sql	Wed Apr 10 10:02:07 2019 +0200
@@ -54,22 +54,27 @@
     -- in m, up to which linestrings will be connected at their boundary
 ) RETURNS geometry
 AS $$
+DECLARE z int;
 DECLARE result_geom geometry;
 BEGIN
+    -- Find best matchting UTM zone
+    z = best_utm(stretch);
+
+    CREATE TEMP TABLE axis AS
+        SELECT id, wtwaxs, ST_Boundary(wtwaxs) AS bdr
+            FROM (SELECT id, ST_Transform(wtwaxs::geometry, z) AS wtwaxs
+                FROM waterway.waterway_axis) AS axs;
+    CREATE INDEX axs_bdr ON axis USING GiST (bdr);
+    ANALYZE axis;
+
     WITH RECURSIVE
-        utm_zone AS (
-            -- Find best matchting UTM zone
-            SELECT best_utm(stretch) AS z),
-        axis AS (
-            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
         points0 AS (
             SELECT ST_ClosestPoint(
                     wtwaxs,
                     ST_Transform(geom, z)) AS geom
-                FROM ST_Dump(ISRSrange_points(stretch)), utm_zone, (
+                FROM ST_Dump(ISRSrange_points(stretch)), (
                     SELECT ST_Collect(wtwaxs) AS wtwaxs
                         FROM axis) AS ax),
         -- Ensure two distinct points on axis have been found
@@ -96,19 +101,19 @@
                         UNION
                         -- Fill eventual gap
                         SELECT ST_ShortestLine(
-                                ST_Boundary(refgeom), ST_Boundary(geom))
+                                ST_Boundary(refgeom), bdr)
                         UNION
                         -- Linestring to be added
                         SELECT geom)))
                 FROM axis_snapped AS axis_snapped (refids, refgeom),
-                    axis AS axis (id, geom),
+                    axis AS axis (id, geom, bdr),
                     (SELECT ST_Collect(points.geom) AS pts
                         FROM points) AS points
                 WHERE id <> ALL(refids)
                     AND ST_DWithin(
-                        ST_Boundary(refgeom), ST_Boundary(geom), tolerance)
+                        ST_Boundary(refgeom), bdr, tolerance)
                     AND NOT ST_Covers(ST_Buffer(refgeom, 0.0001), points.pts)
-                ORDER BY ST_Distance(ST_Boundary(refgeom), ST_Boundary(geom))
+                ORDER BY ST_Boundary(refgeom) <-> bdr
                 FETCH FIRST ROW ONLY)),
         axis_segment AS (
             -- Fetch end result from snapping
@@ -135,11 +140,17 @@
                 FROM points) AS fractions
         GROUP BY axis_segment.line;
 
+    -- Drop temporary table to avoid side effects on PostgreSQL's MVCC,
+    -- because otherwise subsequent invocations of the function will not see
+    -- changes on the underlying waterway.waterway_axis that might have
+    -- occured.
+    DROP TABLE axis;
+
     RETURN result_geom;
 END;
     $$
     LANGUAGE plpgsql
-    STABLE PARALLEL SAFE;
+    PARALLEL RESTRICTED;
 
 -- Clip an area to a stretch given by a geometry representing an axis (e.g.
 -- the output of ISRSrange_axis()).