Mercurial > gemma
comparison 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 |
comparison
equal
deleted
inserted
replaced
2397:b95234702ee9 | 2398:8481e6266691 |
---|---|
80 SELECT ST_Union(a_dmp.geom) AS area | 80 SELECT ST_Union(a_dmp.geom) AS area |
81 FROM axis_substring, utm_zone, | 81 FROM axis_substring, utm_zone, |
82 ST_Dump(ST_Transform(area, z)) AS a_dmp | 82 ST_Dump(ST_Transform(area, z)) AS a_dmp |
83 WHERE ST_Intersects(a_dmp.geom, axis_substring.line) | 83 WHERE ST_Intersects(a_dmp.geom, axis_substring.line) |
84 ), | 84 ), |
85 rotated_ends AS ( | |
86 SELECT ST_Collect(ST_Scale( | |
87 ST_Translate(e, | |
88 (ST_X(p1) - ST_X(p2)) / 2, | |
89 (ST_Y(p1) - ST_Y(p2)) / 2), | |
90 ST_Point(d, d), p1)) AS blade | |
91 FROM axis_substring, area_subset, | |
92 LATERAL (SELECT i, ST_PointN(line, i) AS p1 | |
93 FROM (VALUES (1), (-1)) AS idx (i)) AS ep, | |
94 ST_Rotate(ST_PointN(line, i*2), pi()/2, p1) AS ep2 (p2), | |
95 ST_Makeline(p1, p2) AS e (e), | |
96 LATERAL (SELECT (ST_MaxDistance(p1, area) / ST_Length(e)) | |
97 * 2) AS d (d)), | |
85 range_area AS ( | 98 range_area AS ( |
86 -- Create a buffer around the clipped axis, as large as it could | 99 -- Split area by orthogonal lines at the ends of the clipped axis |
87 -- potentially be intersecting with the area polygon that | 100 SELECT (ST_Dump(ST_CollectionExtract( |
88 -- intersects with the clipped axis. Get the intersection of that | 101 ST_Split(area, blade), 3))).geom |
89 -- buffer with the area polygon, which can potentially | 102 FROM area_subset, rotated_ends) |
90 -- be a multipolygon. | |
91 SELECT (ST_Dump(ST_Intersection( | |
92 ST_Buffer( | |
93 axis_substring.line, | |
94 ST_MaxDistance( | |
95 axis_substring.line, | |
96 area_subset.area), | |
97 'endcap=flat'), | |
98 area_subset.area))).geom | |
99 FROM axis_substring, area_subset) | |
100 -- From the polygons returned by the last CTE, select only those | 103 -- From the polygons returned by the last CTE, select only those |
101 -- around the clipped axis | 104 -- around the clipped axis |
102 SELECT ST_Collect(ST_Transform(range_area.geom, ST_SRID(area))) | 105 SELECT ST_Multi(ST_Union(ST_Transform(range_area.geom, ST_SRID(area)))) |
103 FROM axis_substring, range_area | 106 FROM axis_substring, range_area |
104 WHERE ST_Intersects(range_area.geom, axis_substring.line) | 107 WHERE ST_Intersects(ST_Buffer(range_area.geom, -0.0001), |
108 axis_substring.line) | |
105 $$ | 109 $$ |
106 LANGUAGE sql; | 110 LANGUAGE sql; |