Mercurial > gemma
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 |
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; |