-- 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).
CREATE OR REPLACE FUNCTION ISRSrange_area(
    stretch isrsrange,
    area geometry
) RETURNS geometry
AS $$
    WITH
        -- Get coordinates of location codes
        from_geog AS (
            SELECT geom FROM waterway.distance_marks_virtual
                WHERE location_code = lower(stretch)),
        to_geog AS (
            SELECT geom FROM waterway.distance_marks_virtual
                WHERE location_code = upper(stretch)),
        utm_zone AS (
            -- Find best matchting UTM zone
            SELECT best_utm(ST_Collect(
                    from_geog.geom::geometry,
                    to_geog.geom::geometry)) AS z
                FROM from_geog, to_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
        from_point AS (
            SELECT ST_ClosestPoint(
                    wtwaxs,
                    ST_Transform(from_geog.geom::geometry, z)) AS geom
                FROM axis, from_geog, utm_zone),
        to_point AS (
            SELECT ST_ClosestPoint(
                    wtwaxs,
                    ST_Transform(to_geog.geom::geometry, z)) AS geom
                FROM axis, to_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(from_point.geom, to_point.geom) AS pts
                        FROM from_point, to_point) AS points
                WHERE ST_Covers(lines.line, points.pts)),
        axis_substring AS (
            -- Use linear referencing to clip axis between distance marks
            SELECT ST_LineSubstring(
                    axis_segment.line,
                    ST_LineLocatePoint(axis_segment.line, from_point.geom),
                    ST_LineLocatePoint(axis_segment.line, to_point.geom)
                ) AS line
            FROM axis_segment, from_point, to_point),
        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,
                            ST_Transform(area, z)),
                        'endcap=flat'),
                    ST_Transform(area, z)))).geom
                FROM axis_substring, utm_zone)
        -- From the polygons returned by the last CTE, select only the one
        -- around the clipped axis
        SELECT ST_Transform(range_area.geom, ST_SRID(area))
            FROM axis_substring, range_area
            WHERE ST_Intersects(range_area.geom, axis_substring.line)
    $$
    LANGUAGE sql;