Mercurial > gemma
view schema/isrs_functions.sql @ 2455:54c9fe587fe6
Subdivide SQL function to prepare for improved error handling
The context of an error (e.g. the function in which it occured)
can be inferred by the database client. Not doing all in one
statement will render the context more meaningful.
author | Tom Gottfried <tom@intevation.de> |
---|---|
date | Fri, 01 Mar 2019 18:38:02 +0100 |
parents | 7677a2850a2d |
children | 79f4a20e31c2 |
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; -- 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. CREATE OR REPLACE FUNCTION ISRSrange_axis(stretch isrsrange) RETURNS geometry AS $$ WITH RECURSIVE -- Get coordinates of location codes points_geog AS ( SELECT geom FROM waterway.distance_marks_virtual WHERE location_code = lower(stretch) OR location_code = upper(stretch)), utm_zone AS ( -- Find best matchting UTM zone SELECT best_utm(stretch) AS z), axis AS ( SELECT id, 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 points AS ( SELECT ST_ClosestPoint( wtwaxs, ST_Transform(points_geog.geom::geometry, z)) AS geom FROM points_geog, utm_zone, ( SELECT ST_Collect(wtwaxs) AS wtwaxs FROM axis) AS ax), 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_Intersects( ST_Buffer(axis.wtwaxs, 0.0001), points.geom) 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( ST_Boundary(refgeom), ST_Boundary(geom)) UNION -- Linestring to be added SELECT geom))) FROM axis_snapped AS axis_snapped (refids, refgeom), axis AS axis (id, geom), (SELECT ST_Collect(points.geom) AS pts FROM points) AS points WHERE id <> ALL(refids) AND NOT ST_Covers(ST_Buffer(refgeom, 0.0001), points.pts) ORDER BY ST_Distance(ST_Boundary(refgeom), ST_Boundary(geom)) FETCH FIRST ROW ONLY)), axis_segment AS ( -- Fetch end result from snapping SELECT wtwaxs AS line FROM axis_snapped WHERE array_length(ids, 1) = ( SELECT max(array_length(ids, 1)) FROM axis_snapped)) -- 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 leads to unexpected results of the buffer with -- endcap=flat in the CTE below. SELECT ST_SimplifyPreserveTopology(ST_LineSubstring( axis_segment.line, min(fractions.f), max(fractions.f)), 0.0001) AS line FROM axis_segment, LATERAL ( SELECT ST_LineLocatePoint(axis_segment.line, points.geom) AS f FROM points) AS fractions GROUP BY axis_segment.line $$ LANGUAGE sql STABLE PARALLEL SAFE; -- Clip an area to a stretch given by a pair of ISRS location codes. -- Uses ISRSrange_axis() to retrieve the respective clipped axis used to find -- perpendicular direction at geo-locations of ISRS codes. -- 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 -- relevant part of the axis is used for clipping. CREATE OR REPLACE FUNCTION ISRSrange_area( stretch isrsrange, area geometry ) RETURNS geometry AS $$ WITH axis_substring AS ( SELECT ISRSrange_axis(stretch) AS line), utm_zone AS ( SELECT best_utm(stretch) AS z), area_subset AS ( -- 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) AS area FROM axis_substring, utm_zone, LATERAL ( SELECT ST_MakeValid(ST_Transform(geom, z)) AS geom FROM ST_Dump(area)) AS a_dmp WHERE ST_Intersects(a_dmp.geom, axis_substring.line) ), 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 axis_substring, area_subset, LATERAL (SELECT i, ST_PointN(line, i) AS p1 FROM (VALUES (1), (-1)) AS idx (i)) AS ep, ST_Rotate(ST_PointN(line, i*2), pi()/2, p1) AS ep2 (p2), ST_Makeline(p1, p2) AS e (e), LATERAL (SELECT (ST_MaxDistance(p1, area) / 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, blade), 3))).geom FROM area_subset, 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))) FROM axis_substring, range_area WHERE ST_Intersects(ST_Buffer(range_area.geom, -0.0001), axis_substring.line) $$ LANGUAGE sql;