Mercurial > gemma
changeset 5513:68358e4603c8
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 | 6738655809eb |
| children | 147a2585e136 |
| files | schema/isrs_functions.sql schema/updates/1466/01.fix_ISRSrange_axis.sql schema/updates/1466/02.fix_bn_areas.sql schema/version.sql |
| diffstat | 4 files changed, 152 insertions(+), 2 deletions(-) [+] |
line wrap: on
line diff
--- a/schema/isrs_functions.sql Fri Oct 15 15:14:33 2021 +0200 +++ b/schema/isrs_functions.sql Tue Oct 19 13:12:39 2021 +0200 @@ -65,7 +65,8 @@ geom AS wtwaxs, ST_Boundary(geom) AS bdr FROM waterway.waterway_axis, - ST_Dump(ST_Transform(wtwaxs::geometry, z)); + ST_Dump(ST_Transform(wtwaxs::geometry, z)) + WHERE validity @> current_timestamp; CREATE INDEX axs_bdr ON axis USING GiST (bdr); ANALYZE axis;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/schema/updates/1466/01.fix_ISRSrange_axis.sql Tue Oct 19 13:12:39 2021 +0200 @@ -0,0 +1,106 @@ +CREATE OR REPLACE FUNCTION ISRSrange_axis( + stretch isrsrange, + tolerance float + -- in m, up to which linestrings will be connected at their boundary +) RETURNS geometry +AS $$ +DECLARE z int; +DECLARE result_geom geometry; +BEGIN + -- Find best matchting UTM zone + z = best_utm(stretch); + + CREATE TEMP TABLE axis AS + SELECT row_number() OVER () AS id, + geom AS wtwaxs, + ST_Boundary(geom) AS bdr + FROM waterway.waterway_axis, + ST_Dump(ST_Transform(wtwaxs::geometry, z)) + WHERE validity @> current_timestamp; + CREATE INDEX axs_bdr ON axis USING GiST (bdr); + ANALYZE axis; + + WITH RECURSIVE + -- In order to guarantee the following ST_Covers to work, + -- snap distance mark coordinates to axis + points0 AS ( + SELECT ST_ClosestPoint( + wtwaxs, + ST_Transform(geom, z)) AS geom + FROM ST_Dump(ISRSrange_points(stretch)), ( + SELECT ST_Collect(wtwaxs) AS wtwaxs + FROM axis) AS ax), + -- Ensure two distinct points on axis have been found + points AS ( + SELECT geom + FROM points0 + WHERE 2 = (SELECT count(DISTINCT geom) FROM points0)), + axis_snapped AS ( + -- Iteratively connect non-contiguous axis chunks + -- to find the contiguous axis on which given distance marks lie + (SELECT ARRAY[id] AS ids, wtwaxs + FROM axis, points + WHERE ST_Intersects( + ST_Buffer(axis.wtwaxs, 0.0001), points.geom) + FETCH FIRST ROW ONLY) + UNION + -- Connect endpoint of next linestring with closest + -- endpoint of merged linestring until a contiguous + -- linestring connecting both distance marks is build up + (SELECT refids || id, + ST_LineMerge(ST_Collect(ARRAY( + -- Linestring build up so far + SELECT refgeom + UNION + -- Fill eventual gap + SELECT ST_ShortestLine( + ST_Boundary(refgeom), bdr) + UNION + -- Linestring to be added + SELECT geom))) + FROM axis_snapped AS axis_snapped (refids, refgeom), + axis AS axis (id, geom, bdr), + (SELECT ST_Collect(points.geom) AS pts + FROM points) AS points + WHERE id <> ALL(refids) + AND ST_DWithin( + ST_Boundary(refgeom), bdr, tolerance) + AND NOT ST_Covers(ST_Buffer(refgeom, 0.0001), points.pts) + ORDER BY ST_Boundary(refgeom) <-> bdr + FETCH FIRST ROW ONLY)), + axis_segment AS ( + -- Fetch end result from snapping + SELECT wtwaxs AS line + FROM axis_snapped, + (SELECT ST_Collect(points.geom) AS pts + FROM points) AS points + -- Return end result only if both distance marks were connected + WHERE ST_Covers(ST_Buffer(wtwaxs, 0.0001), points.pts)) + -- Use linear referencing to clip axis between distance marks. + -- Simplification is used to work-around the problem, that + -- ST_LineSubstring might generate very small line segments at an + -- end of the resulting linestring, that significantly differ from + -- the direction of the input linestring due to finite precision + -- of the calculation. The generated small segment of the + -- resulting line would lead e.g. to unexpected results in an area + -- generated by ISRSrange_area(). + SELECT ST_SimplifyPreserveTopology(ST_LineSubstring( + axis_segment.line, min(fractions.f), max(fractions.f)), + 0.0001) AS line + INTO STRICT result_geom + FROM axis_segment, LATERAL ( + SELECT ST_LineLocatePoint(axis_segment.line, points.geom) AS f + FROM points) AS fractions + GROUP BY axis_segment.line; + + -- Drop temporary table to avoid side effects on PostgreSQL's MVCC, + -- because otherwise subsequent invocations of the function will not see + -- changes on the underlying waterway.waterway_axis that might have + -- occured. + DROP TABLE axis; + + RETURN result_geom; +END; + $$ + LANGUAGE plpgsql + PARALLEL RESTRICTED;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/schema/updates/1466/02.fix_bn_areas.sql Tue Oct 19 13:12:39 2021 +0200 @@ -0,0 +1,43 @@ +DO $$ +DECLARE + bns CURSOR FOR SELECT * FROM waterway.bottlenecks FOR UPDATE; + call_stack text; +BEGIN + FOR cur_bn IN bns LOOP + BEGIN + UPDATE waterway.bottlenecks bn SET area = ISRSrange_area( + ISRSrange_axis(stretch, + COALESCE(( + -- Guess tolerance from the last successful + -- bottleneck import owned by a waterway_admin of + -- the country matching the bottleneck_id + SELECT DISTINCT ON (usr.country) + CAST(substring(msg FROM '((\d*\.)?\d+)$') + AS float) + FROM import.import_logs log + JOIN import.imports imp + ON log.import_id = imp.id + JOIN users.list_users usr USING (username) + WHERE starts_with(log.msg, + 'Tolerance used to snap waterway axis:') + AND imp.kind = 'bn' + AND imp.state IN('accepted', 'reviewed') + AND usr.rolname = 'waterway_admin' + AND usr.country = substring( + bn.bottleneck_id FROM 1 FOR 2) + ORDER BY usr.country, imp.changed DESC), + -- Use default tolerance if originally used cannot + -- be determined + 5)), + (SELECT ST_Collect(CAST(area AS geometry)) + FROM waterway.waterway_area)) + WHERE CURRENT OF bns; + EXCEPTION + WHEN no_data_found THEN + GET STACKED DIAGNOSTICS call_stack = PG_EXCEPTION_CONTEXT; + RAISE WARNING '% (%): %, CONTEXT: %', + cur_bn.bottleneck_id, cur_bn.validity, + SQLERRM, call_stack; + END; + END LOOP; +END $$;