--- 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()).