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;