Mercurial > gemma
changeset 5521:3cfbc5769e8b
Make invalid output of ISRSrange_area() less likely
Since ST_Transform() might produce an invalid polygon, add an extra
ST_MakeValid().
| author | Tom Gottfried <tom@intevation.de> |
|---|---|
| date | Thu, 21 Oct 2021 20:20:30 +0200 |
| parents | 05db984d3db1 |
| children | 728b58946c34 |
| files | schema/isrs_functions.sql schema/updates/1468/01.makevalid_isrsrange_area.sql schema/version.sql |
| diffstat | 3 files changed, 54 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- 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);
--- /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;