Mercurial > gemma
changeset 1629:9d51f022b8ee
Introduce SQL function to clip an area to a stretch
author | Tom Gottfried <tom@intevation.de> |
---|---|
date | Wed, 19 Dec 2018 17:59:31 +0100 |
parents | 8ad51ad5a9ee |
children | f6218f11062a |
files | schema/install-db.sh schema/isrs_functions.sql schema/isrs_tests.sql schema/run_tests.sh schema/tap_tests_data.sql |
diffstat | 5 files changed, 134 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/schema/install-db.sh Wed Dec 19 16:48:00 2018 +0100 +++ b/schema/install-db.sh Wed Dec 19 17:59:31 2018 +0100 @@ -125,6 +125,7 @@ -f "$BASEDIR/geonames.sql" \ -f "$BASEDIR/manage_users.sql" \ -f "$BASEDIR/auth.sql" \ + -f "$BASEDIR/isrs_functions.sql" \ -f "$BASEDIR/default_sysconfig.sql"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/schema/isrs_functions.sql Wed Dec 19 17:59:31 2018 +0100 @@ -0,0 +1,82 @@ +-- Clip an area to a stretch given by a pair of ISRS location codes. +-- Uses the table waterway.distance_marks_virtual to map ISRS location codes +-- to their geo-location and the table waterway.waterway_axis to retrieve +-- perpendicular direction at these geo-locations. +-- Distance marks are assumed to be near the axis and the area passed as +-- argument is assumed to intersect with the axis +-- (use e.g. waterway area or fairway dimensions). +CREATE OR REPLACE FUNCTION ISRSrange_area( + stretch isrsrange, + area geometry +) RETURNS geometry +AS $$ + WITH + -- Get coordinates of location codes + from_geog AS ( + SELECT geom FROM waterway.distance_marks_virtual + WHERE location_code = lower(stretch)), + to_geog AS ( + SELECT geom FROM waterway.distance_marks_virtual + WHERE location_code = upper(stretch)), + utm_zone AS ( + -- Find best matchting UTM zone + SELECT best_utm(ST_Collect( + from_geog.geom::geometry, + to_geog.geom::geometry)) AS z + FROM from_geog, to_geog), + axis AS ( + -- Transform and sew together contiguous axis chunks + SELECT ST_LineMerge(ST_Collect(ST_Transform( + wtwaxs::geometry, z))) AS wtwaxs + FROM waterway.waterway_axis, utm_zone), + -- In order to guarantee the following ST_Covers to work, + -- snap distance mark coordinates to axis + from_point AS ( + SELECT ST_ClosestPoint( + wtwaxs, + ST_Transform(from_geog.geom::geometry, z)) AS geom + FROM axis, from_geog, utm_zone), + to_point AS ( + SELECT ST_ClosestPoint( + wtwaxs, + ST_Transform(to_geog.geom::geometry, z)) AS geom + FROM axis, to_geog, utm_zone), + axis_segment AS ( + -- select the contiguous axis on which distance marks lie + SELECT line + FROM ( + SELECT (ST_Dump(wtwaxs)).geom AS line + FROM axis) AS lines, + (SELECT ST_Collect(from_point.geom, to_point.geom) AS pts + FROM from_point, to_point) AS points + WHERE ST_Covers(lines.line, points.pts)), + axis_substring AS ( + -- Use linear referencing to clip axis between distance marks + SELECT ST_LineSubstring( + axis_segment.line, + ST_LineLocatePoint(axis_segment.line, from_point.geom), + ST_LineLocatePoint(axis_segment.line, to_point.geom) + ) AS line + FROM axis_segment, from_point, to_point), + range_area AS ( + -- Create a buffer around the clipped axis, as large as it could + -- potentially be intersecting with the area polygon that + -- intersects with the clipped axis. Get the intersection of that + -- buffer with the area polygon, which can potentially + -- be a multipolygon. + SELECT (ST_Dump(ST_Intersection( + ST_Buffer( + axis_substring.line, + ST_MaxDistance( + axis_substring.line, + ST_Transform(area, z)), + 'endcap=flat'), + ST_Transform(area, z)))).geom + FROM axis_substring, utm_zone) + -- From the polygons returned by the last CTE, select only the one + -- around the clipped axis + SELECT ST_Transform(range_area.geom, ST_SRID(area)) + FROM axis_substring, range_area + WHERE ST_Intersects(range_area.geom, axis_substring.line) + $$ + LANGUAGE sql;
--- a/schema/isrs_tests.sql Wed Dec 19 16:48:00 2018 +0100 +++ b/schema/isrs_tests.sql Wed Dec 19 17:59:31 2018 +0100 @@ -28,3 +28,33 @@ $$, 22023, NULL, 'ISRS text input needs to have correct length'); + +SELECT ok( + ISRSrange_area(isrsrange( + ('AT', 'XXX', '00001', '00000', 0)::isrs, + ('AT', 'XXX', '00001', '00000', 1)::isrs), + ST_SetSRID('POLYGON((0 1, 0 2, 1 2, 1 1, 0 1))'::geometry, 4326) + ) IS NULL, + 'ISRSrange_area returns NULL, if given area does not intersect with axis'); + +SELECT ok( + ST_DWithin( + (SELECT geom FROM waterway.distance_marks_virtual + WHERE location_code = ('AT', 'XXX', '00001', '00000', 0)::isrs), + ST_Boundary(ISRSrange_area(isrsrange( + ('AT', 'XXX', '00001', '00000', 0)::isrs, + ('AT', 'XXX', '00001', '00000', 1)::isrs), + ST_SetSRID('POLYGON((-1 1, 2 1, 2 -1, -1 -1, -1 1))'::geometry, + 4326)))::geography, + 1) + AND + ST_DWithin( + (SELECT geom FROM waterway.distance_marks_virtual + WHERE location_code = ('AT', 'XXX', '00001', '00000', 1)::isrs), + ST_Boundary(ISRSrange_area(isrsrange( + ('AT', 'XXX', '00001', '00000', 0)::isrs, + ('AT', 'XXX', '00001', '00000', 1)::isrs), + ST_SetSRID('POLYGON((-1 1, 2 1, 2 -1, -1 -1, -1 1))'::geometry, + 4326)))::geography, + 1), + 'Resulting polygon almost ST_Touches points corresponding to stretch');
--- a/schema/run_tests.sh Wed Dec 19 16:48:00 2018 +0100 +++ b/schema/run_tests.sh Wed Dec 19 17:59:31 2018 +0100 @@ -28,7 +28,7 @@ -c 'SET client_min_messages TO WARNING' \ -c "DROP ROLE IF EXISTS $TEST_ROLES" \ -f tap_tests_data.sql \ - -c 'SELECT plan(45)' \ + -c 'SELECT plan(47)' \ -f isrs_tests.sql \ -f auth_tests.sql \ -f manage_users_tests.sql \
--- a/schema/tap_tests_data.sql Wed Dec 19 16:48:00 2018 +0100 +++ b/schema/tap_tests_data.sql Wed Dec 19 17:59:31 2018 +0100 @@ -64,6 +64,26 @@ 1, 'depth', 'testorganization', true ); +INSERT INTO waterway.distance_marks_virtual VALUES ( + ('AT', 'XXX', '00001', '00000', 0)::isrs, + ST_SetSRID('POINT(0 0)'::geometry, 4326), + 'someENC' +), ( + ('AT', 'XXX', '00001', '00000', 1)::isrs, + ST_SetSRID('POINT(1 0)'::geometry, 4326), + 'someENC' +); + +INSERT INTO waterway.waterway_axis (wtwaxs, objnam) VALUES ( + ST_SetSRID(ST_CurveToLine( + 'CIRCULARSTRING(0 0, 0.5 0.5, 1 0, 1.5 -0.2, 2 0)'::geometry), + 4326), + 'testriver' +), ( + ST_SetSRID('LINESTRING(0.5 0.5, 1 1)'::geometry, 4326), + 'testriver' +); + INSERT INTO users.templates (template_name, template_data) VALUES ('AT', '\x'), ('RO', '\x'); INSERT INTO users.user_templates