Mercurial > gemma
view schema/isrs_functions.sql @ 2419:f374ca26a337 staging_consolidation
removed depency on old component
author | Thomas Junk <thomas.junk@intevation.de> |
---|---|
date | Thu, 28 Feb 2019 14:03:07 +0100 |
parents | 8481e6266691 |
children | 4c16f5ad1905 |
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> -- Clip an area 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 and the table waterway.waterway_axis to retrieve -- perpendicular direction at these geo-locations. -- Distance marks are assumed to be near the axis and 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 -- 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(ST_Collect(geom::geometry)) AS z FROM points_geog), axis AS ( -- Transform and sew together contiguous axis chunks SELECT ST_LineMerge(ST_Collect(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 axis, points_geog, utm_zone), axis_segment AS ( -- select the contiguous axis on which distance marks lie SELECT line FROM ( SELECT (ST_Dump(wtwaxs)).geom AS line FROM axis) AS lines, (SELECT ST_Collect(points.geom) AS pts FROM points) AS points WHERE ST_Covers(ST_Buffer(lines.line, 0.0001), points.pts)), axis_substring AS ( -- 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), 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, ST_Dump(ST_Transform(area, z)) 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_Union(ST_Transform(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;