--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/schema/isrs_functions.sql	Wed Dec 19 17:59:31 2018 +0100
@@ -0,0 +1,82 @@
+-- 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;