<div dir="ltr">Thanks Sandro, I will try to add these snapping steps in the process.<div><br></div><div>The following code:</div><div><br></div><div><div>select dropTopology('topo1');</div><div>SELECT topology.CreateTopology('topo1',4326, 0.000001);</div>


<div>select st_createTopoGeo('topo1', st_collect(geom))</div><div>from ways;</div><div><br></div><div>Still runs after  1790 sec.</div><div><br></div><div style>Nicolas</div><div class="gmail_extra"><br><br><div class="gmail_quote">

On 9 May 2013 18:02, Sandro Santilli <span dir="ltr"><<a href="mailto:strk@keybit.net" target="_blank">strk@keybit.net</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Note that ST_Intersect may sligthly move the inputs while attempting<br>



to catch robustness issues, so it's recommended to snap the final<br>
set of intersection points to a grid and remove duplicates, then<br>
finally snap each input segment to the so-found points to ensure<br>
the are found on the lines they initially were computed from.<br>
<br>
Out of curiosity: how much does ST_CreateTopoGeo takes when<br>
passed the input as a collection ?<br>
<br>
--strk;<br>
<div><div><br>
On Thu, May 09, 2013 at 05:52:54PM +0200, Nicolas Ribot wrote:<br>
> On 9 May 2013 15:43, Stephen Woodbridge <<a href="mailto:woodbri@swoodbridge.com" target="_blank">woodbri@swoodbridge.com</a>> wrote:<br>
><br>
> > On 5/9/2013 7:20 AM, Nicolas Ribot wrote:<br>
> ><br>
> >> (I spammed this thread a bit with image attachment....)<br>
> >><br>
> >> Hi Steve,<br>
> >><br>
> >> We the given dataset, my approach is indeed slow compared to st_union<br>
> >> approach (though precision for the st_dwithin clause must be adapted to<br>
> >> the current dataset. I took the following precision: 0.000001)<br>
> >><br>
> >> The st_union method generates 18322 segments in 7318 ms, though the<br>
> >> final association between original  lines and new segment is not done<br>
> >> here.<br>
> >><br>
> >> With the query I gave, the st_dwithin part takes 11.7 sec on a recent<br>
> >> laptop machine (1.8 Ghz Intel Core I7, 1024 mb of ram for shared_buffer,<br>
> >> 512 for work_mem)...<br>
> >><br>
> >> The complete query returns 17292 segments in 17956 ms.<br>
> >><br>
> >> As the lines are almost already noded, it generates a lot of<br>
> >> intersection points coincident with one line ends.<br>
> >><br>
> >> As you noted, intermediate temp tables may help here:<br>
> >><br>
> >> I decomposed the query into intermediate steps and the performance is<br>
> >> about the same as with st_union :<br>
> >><br>
> >> -- First creates temp table with intersection points<br>
> >> drop table  if exists intergeom;<br>
> >> create temp table intergeom as<br>
> >> select <a href="http://l1.id" target="_blank">l1.id</a> as l1id, <a href="http://l2.id" target="_blank">l2.id</a> as l2id,<br>
> >> st_intersection(l1.geom, l2.geom) as geom<br>
> >> from bdaways l1 join bdaways l2 on (st_dwithin(l1.geom, l2.geom,<br>
> >> 0.000001))<br>
> >> where <a href="http://l1.id" target="_blank">l1.id</a>  <> <a href="http://l2.id" target="_blank">l2.id</a> ;<br>
> >><br>
> >> -- keeps only true intersection points<br>
> >> -- must handle the case where lines intersects at a linestring...<br>
> >><br>
> ><br>
> > Would it make sense to take all the geometryType(geom) <> 'LINESTRING' and<br>
> > just add the end points to the intergeom table. I think this would then add<br>
> > break points at the ends of overlapping segments insure that they get<br>
> > divided at those points and it would also add extraneous points at the<br>
> > other end, but these would get filter out later. So add these two lines<br>
> > here?<br>
> ><br>
><br>
> You meant, geometryType(geom) = 'LINESTRING ?<br>
> Yes it would make sense. I will try.<br>
><br>
><br>
> > insert intergeom (l1id, l2id, geom)<br>
> >        values (l1id, l2id, st_startpoint(geom))<br>
> >        where geometryType(geom) <> 'LINESTRING';<br>
> ><br>
> > insert intergeom (l1id, l2id, geom)<br>
> >        values (l1id, l2id, st_endpoint(geom))<br>
> >        where geometryType(geom) <> 'LINESTRING';<br>
> ><br>
> ><br>
> ><br>
> >  delete from intergeom where geometryType(geom) <> 'POINT';<br>
> >><br>
> ><br>
> > I think some OSM data from a small central american country might make a<br>
> > good test case. I have not tried any of these:<br>
> ><br>
</div></div>> > <a href="http://downloads.cloudmade." target="_blank">http://downloads.cloudmade.</a>**com/americas/central_america<<a href="http://downloads.cloudmade.com/americas/central_america" target="_blank">http://downloads.cloudmade.com/americas/central_america</a>><br>



<div><div>> ><br>
> > but I notice that if you click into a country like Belize or Costa Rica,<br>
> > there is a shapefile for the country near the bottom of the list or you can<br>
> > pick just a city for a smaller data set.<br>
> ><br>
> > I'll start packaging this into a stored procedure. This is an awesome job.<br>
> ><br>
><br>
> Thanks.<br>
> I tested the costa rica highways dataset and found the following results:<br>
><br>
> The query runs in about 28 seconds to process 29187 ways that are not at<br>
> all connected between them.<br>
> It generates 48415 segments.<br>
><br>
> As an extra steps, all original ways that did not need to be cut by<br>
> intersection points have to be added to the final results: these lines<br>
>  are already clean but were dismissed during the cutting process:<br>
><br>
> insert into res (gid, sub_id, geom, type)<br>
> select ways.gid, 1 as sub_id, ways.geom, geometryType(ways.geom)<br>
>  from ways<br>
> where not exists (<br>
> select res.gid from res where res.gid = ways.gid<br>
>  );<br>
><br>
> The result looks good, as all ways seems to be segmented (cf.<br>
> <a href="http://sd-38177.dedibox.fr/img1.png" target="_blank">http://sd-38177.dedibox.fr/img1.png</a> and <a href="http://sd-38177.dedibox.fr/img2.png" target="_blank">http://sd-38177.dedibox.fr/img2.png</a>:<br>
> black: initial lines, blue: new segments. Lines ends are shown as circles).<br>
><br>
> A problem remains with non-simple or dirty lines as they self-intersect or<br>
> have several ends and so are not noded properly.<br>
> Though it should be simple to process them by identifying the problem. here<br>
> is an example of such a line: <a href="http://sd-38177.dedibox.fr/img3.png" target="_blank">http://sd-38177.dedibox.fr/img3.png</a><br>
><br>
> I go looking at the procedure ;)<br>
><br>
>  Nicolas<br>
</div></div><div><div>_______________________________________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org" target="_blank">postgis-users@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users" target="_blank">http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users</a><br>
</div></div></blockquote></div><br></div></div></div>