CREATE OR REPLACE FUNCTION ISRSrange_axis(
    stretch isrsrange,
    tolerance float
    -- 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 row_number() OVER () AS id,
                geom AS wtwaxs,
                ST_Boundary(geom) AS bdr
            FROM waterway.waterway_axis,
                ST_Dump(ST_Transform(wtwaxs::geometry, z))
            WHERE validity @> current_timestamp;
    CREATE INDEX axs_bdr ON axis USING GiST (bdr);
    ANALYZE axis;

    WITH RECURSIVE
        -- 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)), (
                    SELECT ST_Collect(wtwaxs) AS wtwaxs
                        FROM axis) AS ax),
        -- Ensure two distinct points on axis have been found
        points AS (
            SELECT geom
                FROM points0
                WHERE 2 = (SELECT count(DISTINCT geom) FROM points0)),
        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_ShortestLine(
                                ST_Boundary(refgeom), bdr)
                        UNION
                        -- Linestring to be added
                        SELECT geom)))
                FROM axis_snapped AS axis_snapped (refids, refgeom),
                    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), bdr, tolerance)
                    AND NOT ST_Covers(ST_Buffer(refgeom, 0.0001), points.pts)
                ORDER BY ST_Boundary(refgeom) <-> bdr
                FETCH FIRST ROW ONLY)),
        axis_segment AS (
            -- Fetch end result from snapping
            SELECT wtwaxs AS line
                FROM axis_snapped,
                    (SELECT ST_Collect(points.geom) AS pts
                        FROM points) AS points
                -- Return end result only if both distance marks were connected
                WHERE ST_Covers(ST_Buffer(wtwaxs, 0.0001), points.pts))
        -- Use linear referencing to clip axis between distance marks.
        -- Simplification is used to work-around the problem, that
        -- ST_LineSubstring might generate very small line segments at an
        -- end of the resulting linestring, that significantly differ from
        -- the direction of the input linestring due to finite precision
        -- of the calculation. The generated small segment of the
        -- resulting line would lead e.g. to unexpected results in an area
        -- generated by ISRSrange_area().
        SELECT ST_SimplifyPreserveTopology(ST_LineSubstring(
                    axis_segment.line, min(fractions.f), max(fractions.f)),
                0.0001) AS line
        INTO STRICT result_geom
        FROM axis_segment, LATERAL (
            SELECT ST_LineLocatePoint(axis_segment.line, points.geom) AS f
                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
    PARALLEL RESTRICTED;