Mercurial > gemma
changeset 5009:e8b2dc771f9e
Store axis as MultiLinestring
MultiLinestrings could already be imported but we stored them as
multiple Linestrings with identical attributes and even stored
Linestrings with self-intersections as multiple single Linestrings
with identical attributes. Avoid both by storing as MultiLinestring.
In passing, removed unnecessary processing steps in the INSERT
statemet for the sys_admin case and ensured that attempts to
convert to valid simple features are made after transformation,
which might lead to invalid features.
Since isrsrange_axis() relies on single Linestrings for linear
referencing, add an extra ST_Dump().
| author | Tom Gottfried <tom@intevation.de> |
|---|---|
| date | Wed, 11 Mar 2020 17:11:23 +0100 |
| parents | 0b97f5301a17 |
| children | 1d1be6bd5304 |
| files | pkg/imports/wx.go schema/gemma.sql schema/gemma_tests.sql schema/isrs_functions.sql schema/tap_tests_data.sql schema/updates/1423/01.axis_as_multilinestring.sql schema/updates/1423/02.adapt_func.sql schema/version.sql |
| diffstat | 8 files changed, 157 insertions(+), 49 deletions(-) [+] |
line wrap: on
line diff
--- a/pkg/imports/wx.go Wed Mar 11 15:54:37 2020 +0100 +++ b/pkg/imports/wx.go Wed Mar 11 17:11:23 2020 +0100 @@ -100,16 +100,18 @@ SELECT users.current_user_area_utm() AS a ) INSERT INTO waterway.waterway_axis (wtwaxs, objnam, nobjnam) -SELECT dmp.geom, $3, $4 +SELECT + ST_Multi(ST_Node(ST_CollectionExtract(ST_Transform(new_ax, 4326), 2))), + $3, $4 FROM ST_GeomFromWKB($1, $2::integer) AS new_line (new_line), - ST_Dump(ST_Transform(ST_CollectionExtract( + LATERAL (SELECT CASE WHEN pg_has_role('sys_admin', 'MEMBER') - THEN ST_Node(ST_Transform(new_line, - best_utm(ST_Transform(new_line, 4326)))) + THEN new_line ELSE ST_Intersection((SELECT a FROM resp), ST_Node(ST_Transform(new_line, (SELECT ST_SRID(a) FROM resp)))) - END, - 2), 4326)) AS dmp + END) AS new_ax (new_ax) + -- Do nothing if intersection is empty: + WHERE NOT ST_IsEmpty(new_ax) RETURNING id ` ) @@ -222,47 +224,34 @@ nobjnam = sql.NullString{String: *props.NObjNnm, Valid: true} } + var ls multiLineSlice switch feature.Geometry.Type { case "LineString": var l lineSlice if err := json.Unmarshal(*feature.Geometry.Coordinates, &l); err != nil { return err } - if err := storeLinestring( - ctx, - savepoint, - feedback, - l, - epsg, - props, - nobjnam, - &outside, - &features, - insertStmt); err != nil { - return err - } + ls = append(ls, l) case "MultiLineString": - var ls []lineSlice if err := json.Unmarshal(*feature.Geometry.Coordinates, &ls); err != nil { return err } - for _, l := range ls { - if err := storeLinestring( - ctx, - savepoint, - feedback, - l, - epsg, - props, - nobjnam, - &outside, - &features, - insertStmt); err != nil { - return err - } - } default: unsupported[feature.Geometry.Type]++ + continue + } + if err := storeLinestring( + ctx, + savepoint, + feedback, + ls, + epsg, + props, + nobjnam, + &outside, + &features, + insertStmt); err != nil { + return err } } return nil @@ -302,7 +291,7 @@ ctx context.Context, savepoint func(func() error) error, feedback Feedback, - l lineSlice, + l multiLineSlice, epsg int, props waterwayAxisProperties, nobjnam sql.NullString,
--- a/schema/gemma.sql Wed Mar 11 15:54:37 2020 +0100 +++ b/schema/gemma.sql Wed Mar 11 17:11:23 2020 +0100 @@ -620,7 +620,7 @@ CREATE TABLE waterway_axis ( id int PRIMARY KEY GENERATED BY DEFAULT AS IDENTITY, - wtwaxs geography(LINESTRING, 4326) NOT NULL + wtwaxs geography(MULTILINESTRING, 4326) NOT NULL CHECK(ST_IsSimple(CAST(wtwaxs AS geometry))), -- TODO: Do we need to check data set quality (DRC 2.1.6)? objnam varchar NOT NULL,
--- a/schema/gemma_tests.sql Wed Mar 11 15:54:37 2020 +0100 +++ b/schema/gemma_tests.sql Wed Mar 11 17:11:23 2020 +0100 @@ -26,8 +26,8 @@ SELECT throws_ok($$ INSERT INTO waterway.waterway_axis (wtwaxs, objnam) VALUES - (ST_GeogFromText('LINESTRING(0 0, 1 1)'), 'test'), - (ST_GeogFromText('LINESTRING(0 0, 1 1)'), 'test') + (ST_GeogFromText('MULTILINESTRING((0 0, 1 1))'), 'test'), + (ST_GeogFromText('MULTILINESTRING((0 0, 1 1))'), 'test') $$, 23505, NULL, 'No duplicate geometries can be inserted into waterway_axis');
--- a/schema/isrs_functions.sql Wed Mar 11 15:54:37 2020 +0100 +++ b/schema/isrs_functions.sql Wed Mar 11 17:11:23 2020 +0100 @@ -61,9 +61,11 @@ z = best_utm(stretch); CREATE TEMP TABLE axis AS - SELECT id, wtwaxs, ST_Boundary(wtwaxs) AS bdr - FROM (SELECT id, ST_Transform(wtwaxs::geometry, z) AS wtwaxs - FROM waterway.waterway_axis) AS axs; + 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)); CREATE INDEX axs_bdr ON axis USING GiST (bdr); ANALYZE axis;
--- a/schema/tap_tests_data.sql Wed Mar 11 15:54:37 2020 +0100 +++ b/schema/tap_tests_data.sql Wed Mar 11 17:11:23 2020 +0100 @@ -124,22 +124,23 @@ ); INSERT INTO waterway.waterway_axis (wtwaxs, objnam) VALUES ( - ST_SetSRID(ST_CurveToLine( - 'CIRCULARSTRING(0 0, 0.5 0.5, 0.6 0.4)'), + ST_SetSRID( + ST_Multi(ST_CurveToLine('CIRCULARSTRING(0 0, 0.5 0.5, 0.6 0.4)')), 4326), 'testriver' ), ( - ST_SetSRID(ST_CurveToLine('CIRCULARSTRING(0.6 0.4, 1 0, 1.5 -0.00001)'), + ST_SetSRID( + ST_Multi(ST_CurveToLine('CIRCULARSTRING(0.6 0.4, 1 0, 1.5 -0.00001)')), 4326), 'testriver' ), ( - ST_SetSRID('LINESTRING(0.5 0.5, 1 1)'::geometry, 4326), + ST_SetSRID('MULTILINESTRING((0.5 0.5, 1 1))'::geometry, 4326), 'testriver' ), ( - ST_SetSRID('LINESTRING(1.5 0, 1.55001 0)'::geometry, 4326), + ST_SetSRID('MULTILINESTRING((1.5 0, 1.55001 0))'::geometry, 4326), 'testriver' ), ( - ST_SetSRID('LINESTRING(1.55 0, 2 0)'::geometry, 4326), + ST_SetSRID('MULTILINESTRING((1.55 0, 2 0))'::geometry, 4326), 'testriver' );
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/schema/updates/1423/01.axis_as_multilinestring.sql Wed Mar 11 17:11:23 2020 +0100 @@ -0,0 +1,11 @@ +-- Cannot alter type of a column used in a trigger definition +DROP TRIGGER waterway_axis_wtwaxs_unique ON waterway.waterway_axis; + +ALTER TABLE waterway.waterway_axis + ALTER wtwaxs TYPE geography(MULTILINESTRING, 4326) + USING ST_Multi(CAST(wtwaxs AS geometry)); + +-- Re-create trigger +CREATE CONSTRAINT TRIGGER waterway_axis_wtwaxs_unique + AFTER INSERT OR UPDATE OF wtwaxs ON waterway.waterway_axis + FOR EACH ROW EXECUTE FUNCTION prevent_st_equals('wtwaxs');
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/schema/updates/1423/02.adapt_func.sql Wed Mar 11 17:11:23 2020 +0100 @@ -0,0 +1,105 @@ +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)); + 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;