-- This is Free Software under GNU Affero General Public License v >= 3.0
-- without warranty, see README.md and license for details.

-- SPDX-License-Identifier: AGPL-3.0-or-later
-- License-Filename: LICENSES/AGPL-3.0.txt

-- Copyright (C) 2018, 2019 by via donau
--   – Österreichische Wasserstraßen-Gesellschaft mbH
-- Software engineering by Intevation GmbH

-- Author(s):
--  * Tom Gottfried <tom@intevation.de>
--  * Sascha Wilde <wilde@intevation.de>

CREATE OR REPLACE FUNCTION best_utm(stretch isrsrange) RETURNS integer
AS $$
    SELECT best_utm(ST_Collect(geom::geometry))
        FROM waterway.distance_marks_virtual
        WHERE location_code IN (lower(stretch), upper(stretch))
    $$
    LANGUAGE sql
    STABLE PARALLEL SAFE;

-- Return a multipoint with coordinates of stretch end points
CREATE OR REPLACE FUNCTION ISRSrange_points(
    stretch isrsrange
) RETURNS geometry
AS $$
DECLARE result_geom geometry;
BEGIN
    SELECT geom
        INTO STRICT result_geom
        FROM (
            SELECT ST_Collect(CAST(geom AS geometry)) AS geom
                FROM waterway.distance_marks_virtual
                WHERE location_code IN(lower(stretch), upper(stretch))) AS pts
        -- Ensure both have been found
        WHERE ST_NumGeometries(geom) = 2;

    RETURN result_geom;
END;
    $$
    LANGUAGE plpgsql
    STABLE PARALLEL SAFE;

-- Clip waterway axis 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.
-- Distance marks are assumed to be near the axis.
-- Returns the axis geometry transformed to the best matching UTM zone.
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_DWithin(axis.wtwaxs, points.geom, 0.0001)
                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(refbdr, bdr)
                        UNION
                        -- Linestring to be added
                        SELECT geom)))
                FROM axis_snapped AS axis_snapped (refids, refgeom),
                    axis AS axis (id, geom, bdr),
                    ST_Boundary(refgeom) AS refbdr (refbdr)
                WHERE id <> ALL(refids)
                    AND ST_DWithin(refbdr, bdr, tolerance)
                    -- Stop if refgeom goes through both distance marks
                    AND NOT 0.0001 >= ALL(SELECT refgeom <-> geom FROM points)
                ORDER BY refbdr <-> bdr
                FETCH FIRST ROW ONLY)),
        axis_segment AS (
            -- Fetch end result from snapping
            SELECT wtwaxs AS line
                FROM axis_snapped
                -- Return end result only if both distance marks were connected
                WHERE 0.0001 >= ALL(SELECT wtwaxs <-> geom FROM points))
        -- 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;

-- Clip an area to a stretch given by a geometry representing an axis (e.g.
-- the output of ISRSrange_axis()).
-- Clipping is done by cutting the area in perpendicular direction at
-- the ends of the given axis.
-- The area passed as argument is assumed to intersect with the axis
-- (use e.g. waterway area or fairway dimensions).
-- If a multipolygon is passed, the union of the polygons intersecting with the
-- axis is used for clipping.
CREATE OR REPLACE FUNCTION ISRSrange_area(
    axis geometry,
    area geometry
) RETURNS geometry
AS $$
DECLARE
    area_subset geometry;
    result_geom geometry;
BEGIN
    -- In case area is a multipolygon, process the union of those
    -- polygons, which intersect with the axis. The union is to avoid
    -- problems with invalid/self-intersecting multipolygons
    SELECT ST_Union(a_dmp.geom)
        INTO STRICT area_subset
        FROM (SELECT ST_MakeValid(ST_Transform(geom, ST_SRID(axis)))
                FROM ST_Dump(area)) AS a_dmp (geom)
        WHERE ST_Intersects(a_dmp.geom, axis)
        HAVING ST_Union(a_dmp.geom) IS NOT NULL;

    WITH
        rotated_ends AS (
            SELECT ST_Collect(ST_Scale(
                    ST_Translate(e,
                        (ST_X(p1) - ST_X(p2)) / 2,
                        (ST_Y(p1) - ST_Y(p2)) / 2),
                    ST_Point(d, d), p1)) AS blade
                FROM (SELECT i, ST_PointN(axis, i) AS p1
                        FROM (VALUES (1), (-1)) AS idx (i)) AS ep,
                    ST_Rotate(ST_PointN(axis, i*2), pi()/2, p1) AS ep2 (p2),
                    ST_Makeline(p1, p2) AS e (e),
                    LATERAL (
                        SELECT (ST_MaxDistance(p1, area_subset) / ST_Length(e))
                            * 2) AS d (d)),
        range_area AS (
            -- Split area by orthogonal lines at the ends of the clipped axis
            SELECT (ST_Dump(ST_CollectionExtract(
                    ST_Split(area_subset, blade), 3))).geom
                FROM rotated_ends)
        -- From the polygons returned by the last CTE, select only those
        -- around the clipped axis
        SELECT ST_Multi(ST_CollectionExtract(ST_MakeValid(ST_Transform(
                ST_Union(range_area.geom), ST_SRID(area))), 3))
            INTO result_geom
            FROM range_area
            WHERE ST_Intersects(ST_Buffer(range_area.geom, -0.0001), axis);

    RETURN result_geom;
END;
    $$
    LANGUAGE plpgsql
    STABLE PARALLEL SAFE;


-- Check if a given string looks like an ISRS code
CREATE OR REPLACE FUNCTION is_ISRSstring(str text) RETURNS boolean
AS $$
BEGIN
  str = upper(str);
  RETURN (SELECT str SIMILAR TO '[A-Z]{2}[A-Z0-9]{13}[0-9]{5}')
         AND is_country(substring(str from 1 for 2));
END;
    $$
    LANGUAGE plpgsql
    IMMUTABLE PARALLEL SAFE;