# HG changeset patch
# User Tom Gottfried <tom@intevation.de>
# Date 1634840430 -7200
# Node ID 3cfbc5769e8ba7276c4e7bf7a0222e910b06d06d
# Parent  05db984d3db1d025051a77ee18af80c545d159a8
Make invalid output of ISRSrange_area() less likely

Since ST_Transform() might produce an invalid polygon, add an extra
ST_MakeValid().

diff -r 05db984d3db1 -r 3cfbc5769e8b schema/isrs_functions.sql
--- a/schema/isrs_functions.sql	Thu Oct 21 19:50:39 2021 +0200
+++ b/schema/isrs_functions.sql	Thu Oct 21 20:20:30 2021 +0200
@@ -198,7 +198,8 @@
                 FROM rotated_ends)
         -- From the polygons returned by the last CTE, select only those
         -- around the clipped axis
-        SELECT ST_Multi(ST_Transform(ST_Union(range_area.geom), ST_SRID(area)))
+        SELECT ST_Multi(ST_MakeValid(ST_Transform(
+                ST_Union(range_area.geom), ST_SRID(area))))
             INTO result_geom
             FROM range_area
             WHERE ST_Intersects(ST_Buffer(range_area.geom, -0.0001), axis);
diff -r 05db984d3db1 -r 3cfbc5769e8b schema/updates/1468/01.makevalid_isrsrange_area.sql
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/schema/updates/1468/01.makevalid_isrsrange_area.sql	Thu Oct 21 20:20:30 2021 +0200
@@ -0,0 +1,51 @@
+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_MakeValid(ST_Transform(
+                ST_Union(range_area.geom), ST_SRID(area))))
+            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;
diff -r 05db984d3db1 -r 3cfbc5769e8b schema/version.sql
--- a/schema/version.sql	Thu Oct 21 19:50:39 2021 +0200
+++ b/schema/version.sql	Thu Oct 21 20:20:30 2021 +0200
@@ -1,1 +1,1 @@
-INSERT INTO gemma_schema_version(version) VALUES (1467);
+INSERT INTO gemma_schema_version(version) VALUES (1468);