annotate schema/isrs_functions.sql @ 2398:8481e6266691

Fix corner case in area generation from stretch In case orthogonal lines of arbitrary length at the ends of the axis between two given distance marks cross the axis (e.g. because it describes a tight turn), the result of ST_Buffer with endcap=flat is not well defined. Test data have been amended to include such a case as well as appropriate data to test correctness of a fix. The fix includes construction of the afore mentioned orthogonal lines, splitting the given area by these lines and constructing the end result from those parts of the splitted area that intersect with the clipped axis. Due to numerical inaccuracy, the parts might overlap slightly and eventually cross the clipped axis where they should only touch. Therefore, a small buffer was introduced before testing intersection and intersecting parts are dissolved using ST_Union (which is necessary to revert splits inside the resulting polygon anyhow).
author Tom Gottfried <tom@intevation.de>
date Wed, 27 Feb 2019 18:00:35 +0100
parents 4aabbf324e55
children 4c16f5ad1905
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1720
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
1 -- This is Free Software under GNU Affero General Public License v >= 3.0
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
2 -- without warranty, see README.md and license for details.
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
3
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
4 -- SPDX-License-Identifier: AGPL-3.0-or-later
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
5 -- License-Filename: LICENSES/AGPL-3.0.txt
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
6
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
7 -- Copyright (C) 2018, 2019 by via donau
1720
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
8 -- – Österreichische Wasserstraßen-Gesellschaft mbH
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
9 -- Software engineering by Intevation GmbH
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
10
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
11 -- Author(s):
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
12 -- * Tom Gottfried <tom@intevation.de>
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
13 -- * Sascha Wilde <wilde@intevation.de>
29a0fe218af9 Add missing license header
Tom Gottfried <tom@intevation.de>
parents: 1717
diff changeset
14
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
15 -- Clip an area to a stretch given by a pair of ISRS location codes.
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
16 -- Uses the table waterway.distance_marks_virtual to map ISRS location codes
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
17 -- to their geo-location and the table waterway.waterway_axis to retrieve
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
18 -- perpendicular direction at these geo-locations.
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
19 -- Distance marks are assumed to be near the axis and the area passed as
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
20 -- argument is assumed to intersect with the axis
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
21 -- (use e.g. waterway area or fairway dimensions).
2232
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
22 -- If a multipolygon is passed, the union of the polygons intersecting with the
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
23 -- relevant part of the axis is used for clipping.
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
24 CREATE OR REPLACE FUNCTION ISRSrange_area(
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
25 stretch isrsrange,
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
26 area geometry
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
27 ) RETURNS geometry
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
28 AS $$
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
29 WITH
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
30 -- Get coordinates of location codes
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
31 points_geog AS (
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
32 SELECT geom FROM waterway.distance_marks_virtual
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
33 WHERE location_code = lower(stretch)
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
34 OR location_code = upper(stretch)),
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
35 utm_zone AS (
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
36 -- Find best matchting UTM zone
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
37 SELECT best_utm(ST_Collect(geom::geometry)) AS z
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
38 FROM points_geog),
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
39 axis AS (
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
40 -- Transform and sew together contiguous axis chunks
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
41 SELECT ST_LineMerge(ST_Collect(ST_Transform(
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
42 wtwaxs::geometry, z))) AS wtwaxs
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
43 FROM waterway.waterway_axis, utm_zone),
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
44 -- In order to guarantee the following ST_Covers to work,
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
45 -- snap distance mark coordinates to axis
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
46 points AS (
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
47 SELECT ST_ClosestPoint(
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
48 wtwaxs,
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
49 ST_Transform(points_geog.geom::geometry, z)) AS geom
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
50 FROM axis, points_geog, utm_zone),
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
51 axis_segment AS (
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
52 -- select the contiguous axis on which distance marks lie
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
53 SELECT line
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
54 FROM (
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
55 SELECT (ST_Dump(wtwaxs)).geom AS line
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
56 FROM axis) AS lines,
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
57 (SELECT ST_Collect(points.geom) AS pts
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
58 FROM points) AS points
1723
50548a6df009 Fixed fix of fix.
Sascha L. Teichmann <sascha.teichmann@intevation.de>
parents: 1722
diff changeset
59 WHERE ST_Covers(ST_Buffer(lines.line, 0.0001), points.pts)),
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
60 axis_substring AS (
1811
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
61 -- Use linear referencing to clip axis between distance marks.
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
62 -- Simplification is used to work-around the problem, that
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
63 -- ST_LineSubstring might generate very small line segments at an
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
64 -- end of the resulting linestring, that significantly differ from
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
65 -- the direction of the input linestring due to finite precision
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
66 -- of the calculation. The generated small segment of the
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
67 -- resulting line leads to unexpected results of the buffer with
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
68 -- endcap=flat in the CTE below.
c21f72775f6b Fix numerical problems in stretch area generation
Tom Gottfried <tom@intevation.de>
parents: 1723
diff changeset
69 SELECT ST_SimplifyPreserveTopology(ST_LineSubstring(
2373
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
70 axis_segment.line, min(fractions.f), max(fractions.f)),
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
71 0.0001) AS line
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
72 FROM axis_segment, LATERAL (
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
73 SELECT ST_LineLocatePoint(axis_segment.line, points.geom) AS f
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
74 FROM points) AS fractions
4aabbf324e55 Fix cutting of axis segment between two distance marks
Tom Gottfried <tom@intevation.de>
parents: 2232
diff changeset
75 GROUP BY axis_segment.line),
2232
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
76 area_subset AS (
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
77 -- In case area is a multipolygon, process the union of those
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
78 -- polygons, which intersect with the axis. The union is to avoid
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
79 -- problems with invalid/self-intersecting multipolygons
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
80 SELECT ST_Union(a_dmp.geom) AS area
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
81 FROM axis_substring, utm_zone,
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
82 ST_Dump(ST_Transform(area, z)) AS a_dmp
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
83 WHERE ST_Intersects(a_dmp.geom, axis_substring.line)
7936b46a88d4 Handle invalid multipolygons in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 1983
diff changeset
84 ),
2398
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
85 rotated_ends AS (
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
86 SELECT ST_Collect(ST_Scale(
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
87 ST_Translate(e,
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
88 (ST_X(p1) - ST_X(p2)) / 2,
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
89 (ST_Y(p1) - ST_Y(p2)) / 2),
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
90 ST_Point(d, d), p1)) AS blade
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
91 FROM axis_substring, area_subset,
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
92 LATERAL (SELECT i, ST_PointN(line, i) AS p1
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
93 FROM (VALUES (1), (-1)) AS idx (i)) AS ep,
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
94 ST_Rotate(ST_PointN(line, i*2), pi()/2, p1) AS ep2 (p2),
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
95 ST_Makeline(p1, p2) AS e (e),
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
96 LATERAL (SELECT (ST_MaxDistance(p1, area) / ST_Length(e))
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
97 * 2) AS d (d)),
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
98 range_area AS (
2398
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
99 -- Split area by orthogonal lines at the ends of the clipped axis
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
100 SELECT (ST_Dump(ST_CollectionExtract(
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
101 ST_Split(area, blade), 3))).geom
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
102 FROM area_subset, rotated_ends)
1983
f9f1babe52ae Fix area generation from multipolygon input
Tom Gottfried <tom@intevation.de>
parents: 1811
diff changeset
103 -- From the polygons returned by the last CTE, select only those
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
104 -- around the clipped axis
2398
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
105 SELECT ST_Multi(ST_Union(ST_Transform(range_area.geom, ST_SRID(area))))
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
106 FROM axis_substring, range_area
2398
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
107 WHERE ST_Intersects(ST_Buffer(range_area.geom, -0.0001),
8481e6266691 Fix corner case in area generation from stretch
Tom Gottfried <tom@intevation.de>
parents: 2373
diff changeset
108 axis_substring.line)
1629
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
109 $$
9d51f022b8ee Introduce SQL function to clip an area to a stretch
Tom Gottfried <tom@intevation.de>
parents:
diff changeset
110 LANGUAGE sql;