Mercurial > gemma
view schema/updates/1469/01.fix_isrsrange_area.sql @ 5559:ce9a9a1bf92f
Make invalid output of ISRSrange_area() less likely, next try
Since ST_MakeValid() might return a collection of lower-to-equal dimension
geometries, distill only the polygons from it. This should prevent
respective errors when trying to save the result to a column of type
MultiPolygon.
author | Tom Gottfried <tom@intevation.de> |
---|---|
date | Thu, 02 Dec 2021 12:37:33 +0100 |
parents | |
children |
line wrap: on
line source
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_CollectionExtract(ST_MakeValid(ST_Transform( ST_Union(range_area.geom), ST_SRID(area))), 3)) 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;