-- 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>

-- Clip an area 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 and the table waterway.waterway_axis to retrieve
-- perpendicular direction at these geo-locations.
-- Distance marks are assumed to be near the axis and 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
-- relevant part of the axis is used for clipping.
CREATE OR REPLACE FUNCTION ISRSrange_area(
    stretch isrsrange,
    area geometry
) RETURNS geometry
AS $$
    WITH
        -- Get coordinates of location codes
        points_geog AS (
            SELECT geom FROM waterway.distance_marks_virtual
                WHERE location_code = lower(stretch)
                    OR location_code = upper(stretch)),
        utm_zone AS (
            -- Find best matchting UTM zone
            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(
                    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
        points AS (
            SELECT ST_ClosestPoint(
                    wtwaxs,
                    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(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.
            -- 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 leads to unexpected results of the buffer with
            -- endcap=flat in the CTE below.
            SELECT ST_SimplifyPreserveTopology(ST_LineSubstring(
                        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
            -- problems with invalid/self-intersecting multipolygons
            SELECT ST_Union(a_dmp.geom) AS area
                FROM axis_substring, utm_zone,
                    ST_Dump(ST_Transform(area, z)) AS a_dmp
                WHERE ST_Intersects(a_dmp.geom, axis_substring.line)
            ),
        range_area AS (
            -- Create a buffer around the clipped axis, as large as it could
            -- potentially be intersecting with the area polygon that
            -- intersects with the clipped axis. Get the intersection of that
            -- buffer with the area polygon, which can potentially
            -- be a multipolygon.
            SELECT (ST_Dump(ST_Intersection(
                    ST_Buffer(
                        axis_substring.line,
                        ST_MaxDistance(
                            axis_substring.line,
                            area_subset.area),
                        'endcap=flat'),
                    area_subset.area))).geom
                FROM axis_substring, area_subset)
        -- From the polygons returned by the last CTE, select only those
        -- around the clipped axis
        SELECT ST_Collect(ST_Transform(range_area.geom, ST_SRID(area)))
            FROM axis_substring, range_area
            WHERE ST_Intersects(range_area.geom, axis_substring.line)
    $$
    LANGUAGE sql;