Mercurial > gemma
view schema/isrs_functions.sql @ 2425:4c16f5ad1905
Make area generation more robust
Valid geometries can become invalid geometries during transformation
e.g. due to very small 'minimum clearance' (see PostGIS docs) and
numerical inaccuracy of the transformation. Thus, try to repair
polygons of the input area after transformation and calculate
the union of the result polygons before transformation to avoid
'TopologyException' in the union operations.
The simulated waterway area of the changed axis in the test data
is affected here, but real waterway area data with complex shore
lines could likely be affected, too.
author | Tom Gottfried <tom@intevation.de> |
---|---|
date | Thu, 28 Feb 2019 16:36:49 +0100 |
parents | 8481e6266691 |
children | 48495bd3081d |
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, 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;