Hi George,<br><br>Based on the dataset you've uploaded, here are some steps taken to merge the contours.<div>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.</div>
<div><br>1°) Defines the map frame: extent enclosing the working data.</div><div><div><br></div><div><font face="'courier new', monospace">create table frame as (</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select st_exteriorRing(st_setSRID(st_extent(geom)::geometry, 2193)) as geom </font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>from contour_100</font></div><div><font face="'courier new', monospace">);</font></div><div><br></div>
<div>2°) Table of individual segments from contour lines: its dumps MULTILINESTRINGS into LINESTRINGS</div><div><br></div><div><font face="'courier new', monospace">drop table if exists seg;</font></div><div><font face="'courier new', monospace">create table seg as (</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select elevint as elev, (st_dump(geom)).path[1] as path, (st_dump(geom)).geom as geom</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>from contour_100</font></div>
<div><font face="'courier new', monospace">);</font></div><div><font face="'courier new', monospace"><br></font></div><div><font face="'courier new', monospace">-- adds unique gid column to seg</font></div>
<div><font face="'courier new', monospace">alter table seg add column gid serial primary key; </font></div><div><font face="'courier new', monospace">-- spatial index on geom and index on evel</font></div>
<div><font face="'courier new', monospace">create index seg_geom_gist on seg using gist(geom);</font></div><div><font face="'courier new', monospace">create index seg_elev_idx on seg (elev);</font></div><div>
<br></div><div>3°) Identify lines from seg table that can be removed based on their topology: These lines are kept in tables for future use</div><div><br></div><div> • closed lines should be removed (represent a closed contour): </div>
<div><br></div><div><font face="'courier new', monospace">drop table if exists seg_closed;</font></div><div><font face="'courier new', monospace">create table seg_closed as (</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select * from seg where st_isclosed(geom)</font></div>
<div><font face="'courier new', monospace">);</font></div><div><font face="'courier new', monospace"><br></font></div><div><font face="'courier new', monospace">delete from seg where st_isclosed(geom);</font></div>
<div><br></div><div>• 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:</div>
<div><br></div><div><font face="'courier new', monospace">drop table if exists seg_border;</font></div><div><font face="'courier new', monospace">create table seg_border as (</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select s.*</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>from seg s, frame f </font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>where st_dwithin(st_startpoint(s.geom), f.geom, 0.001)</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>and st_dwithin(st_endpoint(s.geom), f.geom, 0.001)</font></div><div><font face="'courier new', monospace">);</font></div>
<div><font face="'courier new', monospace"><br></font></div><div><span style="font-family:'courier new',monospace">delete from seg s using seg_border sb</span></div><div><font face="'courier new', monospace">where s.gid = sb.gid;</font></div>
<div><br></div><div>4°) Creates the table giving, for each segment, the 2 closest segments, that are closer to 25 meters.</div><div>This threshold value is choosen to avoid segments to far away from each other, for instance if a segment is cut by the frame.</div>
<div>This figure should be adapted according to the dataset, maybe based on contour equidistance.</div><div>Segments touching the map frame will have only one closest segment. </div><div><br></div><div><font face="'courier new', monospace">drop table if exists tmp;</font></div>
<div><font face="'courier new', monospace">create table tmp as (</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>with closest as (</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select s1.elev as elev, s1.gid, s1.path as path, s2.gid as closest,</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>s2.geom,</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>st_distance(st_collect(st_startpoint(s1.geom), st_endpoint(s1.geom)), st_collect(st_startpoint(s2.geom), st_endpoint(s2.geom))) as dist,</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>row_number() over (</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>partition by s1.gid </font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>order by s1.gid,</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>st_distance(st_collect(st_startpoint(s1.geom), st_endpoint(s1.geom)), st_collect(st_startpoint(s2.geom), st_endpoint(s2.geom)))) as r</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>from seg s1, seg s2</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>where s1.elev = s2.elev</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>and s1.gid <> s2.gid</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>) select distinct * from closest </font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>where r < 3 and dist < 25</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>order by elev, gid, r</font></div>
<div><font face="'courier new', monospace">);</font></div><div><br></div><div>5°) A first iteration to merge all lines touching map frame. </div><div>It will be our "seed" to start iteration: each segment touching the frame will be completed to its neighbour, and so forth:</div>
<div><br></div><div><font face="'courier new', monospace">drop table if exists merged_contour1;</font></div><div><font face="'courier new', monospace">create table merged_contour1 as (</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>with recursive tab as (</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>-- begins with first segment of each contour: makes a line with it and its closest segment</font></div><div>
<font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>-- storing already processed segments as a stop conditition</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select s.gid, s.elev, t.closest, array[s.gid, t.closest] as gids,</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>-- handle segments orientation to guarantee that closest segments will be oriented</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>-- such that calling st_makeline(g1, g2) will connect the shortest distance</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>st_makeShortestLine(s.geom, t.geom) as geom, 1 as rank</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>from seg s, tmp t, frame f</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>where s.gid = t.gid</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>and s.elev = t.elev</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>-- takes lines touching the map frame at one end point</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>and (st_dwithin(st_startpoint(s.geom), f.geom, 0.001)</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span> or st_dwithin(st_endpoint(s.geom), f.geom, 0.001))</font></div><div><font face="'courier new', monospace"><br>
</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>UNION ALL</font></div><div><font face="'courier new', monospace"><br></font></div><div>
<font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>-- takes closest segment from current one </font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select tab.gid, tab.elev, tmp.closest, gids || tmp.closest, </font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>st_makeShortestLine(tab.geom, tmp.geom) as geom, rank+1</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>from tmp, tab</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>where tmp.elev = tab.elev</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>and tmp.gid = tab.closest</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>and not (tmp.closest = any(gids)) </font></div><div><span class="Apple-tab-span" style="font-family:'courier new',monospace;white-space:pre"> </span><span style="font-family:'courier new',monospace">) select distinct on (gid) gid, elev, rank, gids, geom from (</span></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select * from tab </font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>order by gid, rank desc</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>) as foo</font></div><div><font face="'courier new', monospace">);</font></div><div><br></div><div>6°) These segments are then removed from original table, to keep only segments entirely contained in the map frame, to process them later.</div>
<div>Merged_contour1.gids is the array of all gid used to form the contour. Use it to remove segments already merged:</div><div><br></div><div><font face="'courier new', monospace">delete from seg s using merged_contour1 m</font></div>
<div><font face="'courier new', monospace">where s.gid = any(gids);</font></div><div><br></div><div>0 rows left: all segments from input table are treated.</div><div>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. </div>
<div>A distinct clause on geometries/elevation at the end of iteration could help to remove duplicate contours.</div><div><br></div><div>7°) The final query to group together all contour sharing the same elevation, from our working tables merged_contour1, seg_closed and seg_border:</div>
<div><br></div><div><font face="'courier new', monospace">create table all_contour as (</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>with all_c as (</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select elev, geom from merged_contour1</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>UNION</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>SELECT elev, geom from seg_border</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>UNION</font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>select elev, geom from seg_closed</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>) select elev, st_union(geom) </font></div>
<div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>from all_c</font></div><div><font face="'courier new', monospace"><span class="Apple-tab-span" style="white-space:pre"> </span>group by elev</font></div>
<div><font face="'courier new', monospace">);</font></div><div><br></div></div><div>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.</div>
<div>It's up to you to either leave contours opened, or to close them with the map frame, according to your needs.</div><div><br></div><div>The attached image shows the final result, with the attributes of contours.</div>
<div><br></div><div>Nicolas</div>