Mercurial > gemma
view schema/isrs_functions.sql @ 5520:05db984d3db1
Improve performance of bottleneck area calculation
Avoid buffer calculations by replacing them with simple distance comparisons
and calculate the boundary of the result geometry only once per iteration.
In some edge cases with very large numbers of iterations, this reduced
the runtime of a bottleneck import by a factor of more than twenty.
author | Tom Gottfried <tom@intevation.de> |
---|---|
date | Thu, 21 Oct 2021 19:50:39 +0200 |
parents | 68358e4603c8 |
children | 3cfbc5769e8b |
line wrap: on
line source
-- This is Free Software under GNU Affero General Public License v >= 3.0 -- without warranty, see README.md and license for details. -- SPDX-License-Identifier: AGPL-3.0-or-later -- License-Filename: LICENSES/AGPL-3.0.txt -- Copyright (C) 2018, 2019 by via donau -- – Österreichische Wasserstraßen-Gesellschaft mbH -- Software engineering by Intevation GmbH -- Author(s): -- * Tom Gottfried <tom@intevation.de> -- * Sascha Wilde <wilde@intevation.de> CREATE OR REPLACE FUNCTION best_utm(stretch isrsrange) RETURNS integer AS $$ SELECT best_utm(ST_Collect(geom::geometry)) FROM waterway.distance_marks_virtual WHERE location_code IN (lower(stretch), upper(stretch)) $$ LANGUAGE sql STABLE PARALLEL SAFE; -- Return a multipoint with coordinates of stretch end points CREATE OR REPLACE FUNCTION ISRSrange_points( stretch isrsrange ) RETURNS geometry AS $$ DECLARE result_geom geometry; BEGIN SELECT geom INTO STRICT result_geom FROM ( SELECT ST_Collect(CAST(geom AS geometry)) AS geom FROM waterway.distance_marks_virtual WHERE location_code IN(lower(stretch), upper(stretch))) AS pts -- Ensure both have been found WHERE ST_NumGeometries(geom) = 2; RETURN result_geom; END; $$ LANGUAGE plpgsql STABLE PARALLEL SAFE; -- Clip waterway axis 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. -- Distance marks are assumed to be near the axis. -- Returns the axis geometry transformed to the best matching UTM zone. CREATE OR REPLACE FUNCTION ISRSrange_axis( stretch isrsrange, tolerance float -- in m, up to which linestrings will be connected at their boundary ) RETURNS geometry AS $$ DECLARE z int; DECLARE result_geom geometry; BEGIN -- Find best matchting UTM zone z = best_utm(stretch); CREATE TEMP TABLE axis AS SELECT row_number() OVER () AS id, geom AS wtwaxs, ST_Boundary(geom) AS bdr FROM waterway.waterway_axis, ST_Dump(ST_Transform(wtwaxs::geometry, z)) WHERE validity @> current_timestamp; CREATE INDEX axs_bdr ON axis USING GiST (bdr); ANALYZE axis; WITH RECURSIVE -- In order to guarantee the following ST_Covers to work, -- snap distance mark coordinates to axis points0 AS ( SELECT ST_ClosestPoint( wtwaxs, ST_Transform(geom, z)) AS geom FROM ST_Dump(ISRSrange_points(stretch)), ( SELECT ST_Collect(wtwaxs) AS wtwaxs FROM axis) AS ax), -- Ensure two distinct points on axis have been found points AS ( SELECT geom FROM points0 WHERE 2 = (SELECT count(DISTINCT geom) FROM points0)), axis_snapped AS ( -- Iteratively connect non-contiguous axis chunks -- to find the contiguous axis on which given distance marks lie (SELECT ARRAY[id] AS ids, wtwaxs FROM axis, points WHERE ST_DWithin(axis.wtwaxs, points.geom, 0.0001) FETCH FIRST ROW ONLY) UNION -- Connect endpoint of next linestring with closest -- endpoint of merged linestring until a contiguous -- linestring connecting both distance marks is build up (SELECT refids || id, ST_LineMerge(ST_Collect(ARRAY( -- Linestring build up so far SELECT refgeom UNION -- Fill eventual gap SELECT ST_ShortestLine(refbdr, bdr) UNION -- Linestring to be added SELECT geom))) FROM axis_snapped AS axis_snapped (refids, refgeom), axis AS axis (id, geom, bdr), ST_Boundary(refgeom) AS refbdr (refbdr) WHERE id <> ALL(refids) AND ST_DWithin(refbdr, bdr, tolerance) -- Stop if refgeom goes through both distance marks AND NOT 0.0001 >= ALL(SELECT refgeom <-> geom FROM points) ORDER BY refbdr <-> bdr FETCH FIRST ROW ONLY)), axis_segment AS ( -- Fetch end result from snapping SELECT wtwaxs AS line FROM axis_snapped -- Return end result only if both distance marks were connected WHERE 0.0001 >= ALL(SELECT wtwaxs <-> geom FROM points)) -- Use linear referencing to clip axis between distance marks. -- Simplification is used to work-around the problem, that -- ST_LineSubstring might generate very small line segments at an -- end of the resulting linestring, that significantly differ from -- the direction of the input linestring due to finite precision -- of the calculation. The generated small segment of the -- resulting line would lead e.g. to unexpected results in an area -- generated by ISRSrange_area(). SELECT ST_SimplifyPreserveTopology(ST_LineSubstring( axis_segment.line, min(fractions.f), max(fractions.f)), 0.0001) AS line INTO STRICT result_geom FROM axis_segment, LATERAL ( SELECT ST_LineLocatePoint(axis_segment.line, points.geom) AS f FROM points) AS fractions GROUP BY axis_segment.line; -- Drop temporary table to avoid side effects on PostgreSQL's MVCC, -- because otherwise subsequent invocations of the function will not see -- changes on the underlying waterway.waterway_axis that might have -- occured. DROP TABLE axis; RETURN result_geom; END; $$ LANGUAGE plpgsql PARALLEL RESTRICTED; -- Clip an area to a stretch given by a geometry representing an axis (e.g. -- the output of ISRSrange_axis()). -- Clipping is done by cutting the area in perpendicular direction at -- the ends of the given axis. -- The area passed as argument is assumed to intersect with the axis -- (use e.g. waterway area or fairway dimensions). -- If a multipolygon is passed, the union of the polygons intersecting with the -- axis is used for clipping. 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_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; -- Check if a given string looks like an ISRS code CREATE OR REPLACE FUNCTION is_ISRSstring(str text) RETURNS boolean AS $$ BEGIN str = upper(str); RETURN (SELECT str SIMILAR TO '[A-Z]{2}[A-Z0-9]{13}[0-9]{5}') AND is_country(substring(str from 1 for 2)); END; $$ LANGUAGE plpgsql IMMUTABLE PARALLEL SAFE;