Mercurial > gemma
view schema/updates/1466/01.fix_ISRSrange_axis.sql @ 5736:55892008ec96 default tip
Fixed a bunch of corner cases in WG import.
author | Sascha Wilde <wilde@sha-bang.de> |
---|---|
date | Wed, 29 May 2024 19:02:42 +0200 |
parents | 68358e4603c8 |
children |
line wrap: on
line source
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;