Mercurial > gemma

comparison schema/updates/1466/01.fix_ISRSrange_axis.sql @ 5513:68358e4603c8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Use current axis only for calculating bottleneck areas This is a fixup of rev. cf25b23e3eec, which introduced historic data for the waterway axis but missed to take this into account in the calculation of bottleneck areas, leading to sometimes excessive runtime and bad results due to multiple (almost) equal axis geometries being considered as candidates in the bottleneck stretch. The database migration tries to recalculate all bottleneck areas, while some might fail that did not fail on import. A warning message is emitted for these and the area is left untouched.
author Tom Gottfried <tom@intevation.de>
date Tue, 19 Oct 2021 13:12:39 +0200
parents
children
comparison
equal deleted inserted replaced
5512:6738655809eb 5513:68358e4603c8
1 CREATE OR REPLACE FUNCTION ISRSrange_axis(
2 stretch isrsrange,
3 tolerance float
4 -- in m, up to which linestrings will be connected at their boundary
5 ) RETURNS geometry
6 AS $$
7 DECLARE z int;
8 DECLARE result_geom geometry;
9 BEGIN
10 -- Find best matchting UTM zone
11 z = best_utm(stretch);
12
13 CREATE TEMP TABLE axis AS
14 SELECT row_number() OVER () AS id,
15 geom AS wtwaxs,
16 ST_Boundary(geom) AS bdr
17 FROM waterway.waterway_axis,
18 ST_Dump(ST_Transform(wtwaxs::geometry, z))
19 WHERE validity @> current_timestamp;
20 CREATE INDEX axs_bdr ON axis USING GiST (bdr);
21 ANALYZE axis;
22
23 WITH RECURSIVE
24 -- In order to guarantee the following ST_Covers to work,
25 -- snap distance mark coordinates to axis
26 points0 AS (
27 SELECT ST_ClosestPoint(
28 wtwaxs,
29 ST_Transform(geom, z)) AS geom
30 FROM ST_Dump(ISRSrange_points(stretch)), (
31 SELECT ST_Collect(wtwaxs) AS wtwaxs
32 FROM axis) AS ax),
33 -- Ensure two distinct points on axis have been found
34 points AS (
35 SELECT geom
36 FROM points0
37 WHERE 2 = (SELECT count(DISTINCT geom) FROM points0)),
38 axis_snapped AS (
39 -- Iteratively connect non-contiguous axis chunks
40 -- to find the contiguous axis on which given distance marks lie
41 (SELECT ARRAY[id] AS ids, wtwaxs
42 FROM axis, points
43 WHERE ST_Intersects(
44 ST_Buffer(axis.wtwaxs, 0.0001), points.geom)
45 FETCH FIRST ROW ONLY)
46 UNION
47 -- Connect endpoint of next linestring with closest
48 -- endpoint of merged linestring until a contiguous
49 -- linestring connecting both distance marks is build up
50 (SELECT refids || id,
51 ST_LineMerge(ST_Collect(ARRAY(
52 -- Linestring build up so far
53 SELECT refgeom
54 UNION
55 -- Fill eventual gap
56 SELECT ST_ShortestLine(
57 ST_Boundary(refgeom), bdr)
58 UNION
59 -- Linestring to be added
60 SELECT geom)))
61 FROM axis_snapped AS axis_snapped (refids, refgeom),
62 axis AS axis (id, geom, bdr),
63 (SELECT ST_Collect(points.geom) AS pts
64 FROM points) AS points
65 WHERE id <> ALL(refids)
66 AND ST_DWithin(
67 ST_Boundary(refgeom), bdr, tolerance)
68 AND NOT ST_Covers(ST_Buffer(refgeom, 0.0001), points.pts)
69 ORDER BY ST_Boundary(refgeom) <-> bdr
70 FETCH FIRST ROW ONLY)),
71 axis_segment AS (
72 -- Fetch end result from snapping
73 SELECT wtwaxs AS line
74 FROM axis_snapped,
75 (SELECT ST_Collect(points.geom) AS pts
76 FROM points) AS points
77 -- Return end result only if both distance marks were connected
78 WHERE ST_Covers(ST_Buffer(wtwaxs, 0.0001), points.pts))
79 -- Use linear referencing to clip axis between distance marks.
80 -- Simplification is used to work-around the problem, that
81 -- ST_LineSubstring might generate very small line segments at an
82 -- end of the resulting linestring, that significantly differ from
83 -- the direction of the input linestring due to finite precision
84 -- of the calculation. The generated small segment of the
85 -- resulting line would lead e.g. to unexpected results in an area
86 -- generated by ISRSrange_area().
87 SELECT ST_SimplifyPreserveTopology(ST_LineSubstring(
88 axis_segment.line, min(fractions.f), max(fractions.f)),
89 0.0001) AS line
90 INTO STRICT result_geom
91 FROM axis_segment, LATERAL (
92 SELECT ST_LineLocatePoint(axis_segment.line, points.geom) AS f
93 FROM points) AS fractions
94 GROUP BY axis_segment.line;
95
96 -- Drop temporary table to avoid side effects on PostgreSQL's MVCC,
97 -- because otherwise subsequent invocations of the function will not see
98 -- changes on the underlying waterway.waterway_axis that might have
99 -- occured.
100 DROP TABLE axis;
101
102 RETURN result_geom;
103 END;
104 $$
105 LANGUAGE plpgsql
106 PARALLEL RESTRICTED;