[postgis-users] Closing polylines

[Nicolas Ribot]

Hi George,

Based on the dataset you've uploaded, here are some steps taken to merge
the contours.
Process is using the same principle as the one described in my previous
emails, with some changes to cope with the fact that contours are opened at
map border and that several distinct contours have to be rebuilt for a
given elevation.

1°) Defines the map frame: extent enclosing the working data.

create table frame as (
select st_exteriorRing(st_setSRID(st_extent(geom)::geometry, 2193)) as geom
from contour_100
);

2°) Table of individual segments from contour lines: its dumps
MULTILINESTRINGS into LINESTRINGS

drop table if exists seg;
create table seg as (
select elevint as elev, (st_dump(geom)).path[1] as path,
(st_dump(geom)).geom as geom
from contour_100
);

-- adds unique gid column to seg
alter table seg add column gid serial primary key;
-- spatial index on geom and index on evel
create index seg_geom_gist on seg using gist(geom);
create index seg_elev_idx on seg (elev);

3°) Identify lines from seg table that can be removed based on their
topology: These lines are kept in tables for future use

  • closed lines should be removed (represent a closed contour):

drop table if exists seg_closed;
create table seg_closed as (
select * from seg where st_isclosed(geom)
);

delete from seg where st_isclosed(geom);

• lines touching the map extent at start AND end points: These lines should
be closed with map extent boundary, or left opened for later processing,
but represent, for our dataset, connected contours:

drop table if exists seg_border;
create table seg_border as (
select s.*
from seg s, frame f
where st_dwithin(st_startpoint(s.geom), f.geom, 0.001)
and st_dwithin(st_endpoint(s.geom), f.geom, 0.001)
);

delete from seg s using seg_border sb
where s.gid = sb.gid;

4°) Creates the table giving, for each segment, the 2 closest segments,
that are closer to 25 meters.
This threshold value is choosen to avoid segments to far away from each
other, for instance if a segment is cut by the frame.
This figure should be adapted according to the dataset, maybe based on
contour equidistance.
Segments touching the map frame will have only one closest segment.

drop table if exists tmp;
create table tmp as (
with closest as (
select s1.elev as elev, s1.gid, s1.path as path, s2.gid as closest,
s2.geom,
st_distance(st_collect(st_startpoint(s1.geom), st_endpoint(s1.geom)),
st_collect(st_startpoint(s2.geom), st_endpoint(s2.geom))) as dist,
row_number() over (
partition by s1.gid
order by s1.gid,
st_distance(st_collect(st_startpoint(s1.geom), st_endpoint(s1.geom)),
st_collect(st_startpoint(s2.geom), st_endpoint(s2.geom)))) as r
from seg s1, seg s2
where s1.elev = s2.elev
and s1.gid <> s2.gid
) select distinct * from closest
where r < 3 and dist < 25
order by elev, gid, r
);

5°) A first iteration to merge all lines touching map frame.
It will be our "seed" to start iteration: each segment touching the frame
will be completed to its neighbour, and so forth:

drop table if exists merged_contour1;
create table merged_contour1 as (
with recursive tab as (
-- begins with first segment of each contour: makes a line with it and its
closest segment
-- storing already processed segments as a stop conditition
select s.gid, s.elev, t.closest, array[s.gid, t.closest] as gids,
-- handle segments orientation to guarantee that closest segments will be
oriented
-- such that calling st_makeline(g1, g2) will connect the shortest distance
st_makeShortestLine(s.geom, t.geom) as geom, 1 as rank
from seg s, tmp t, frame f
where s.gid = t.gid
and s.elev = t.elev
-- takes lines touching the map frame at one end point
and (st_dwithin(st_startpoint(s.geom), f.geom, 0.001)
     or st_dwithin(st_endpoint(s.geom), f.geom, 0.001))

UNION ALL

-- takes closest segment from current one
select tab.gid, tab.elev, tmp.closest, gids || tmp.closest,
st_makeShortestLine(tab.geom, tmp.geom) as geom, rank+1
from tmp, tab
where tmp.elev = tab.elev
and tmp.gid = tab.closest
and not (tmp.closest = any(gids))
) select distinct on (gid) gid, elev, rank, gids, geom from (
select * from tab
order by gid, rank desc
) as foo
);

6°) These segments are then removed from original table, to keep only
segments entirely contained in the map frame, to process them later.
Merged_contour1.gids is the array of all gid used to form the contour. Use
it to remove segments already merged:

delete from seg s using merged_contour1 m
where s.gid = any(gids);

0 rows left: all segments from input table are treated.
If there are still rows in seg (for instance, contours completely inside
the map frame), one should rebuild a new association table based on these
rows, and relaunch an iteration to merge these segments.
A distinct clause on geometries/elevation at the end of iteration could
help to remove duplicate contours.

7°) The final query to group together all contour sharing the same
elevation, from our working tables merged_contour1,  seg_closed and
seg_border:

create table all_contour as (
with all_c as (
select elev, geom from merged_contour1
UNION
SELECT elev, geom from seg_border
UNION
select elev, geom from seg_closed
) select elev, st_union(geom)
from all_c
group by elev
);

This all_contour table contains one MULTILINESTRING for each elevation.
Each individual LINESTRING in the collection is either closed, or connected
and opened at the map frame.
It's up to you to either leave contours opened, or to close them with the
map frame, according to your needs.

The attached image shows the final result, with the attributes of contours.

Nicolas
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20120518/785e92d7/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screen Shot 2012-05-18 at 3.18.07 PM.png
Type: image/png
Size: 118757 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20120518/785e92d7/attachment.png>