Note that ST_Intersect may sligthly move the inputs while attempting to catch robustness issues, so it's recommended to snap the final set of intersection points to a grid and remove duplicates, then finally snap each input segment to the so-found points to ensure the are found on the lines they initially were computed from.
Out of curiosity: how much does ST_CreateTopoGeo takes when passed the input as a collection ? --strk; On Thu, May 09, 2013 at 05:52:54PM +0200, Nicolas Ribot wrote: > On 9 May 2013 15:43, Stephen Woodbridge <[email protected]> wrote: > > > On 5/9/2013 7:20 AM, Nicolas Ribot wrote: > > > >> (I spammed this thread a bit with image attachment....) > >> > >> Hi Steve, > >> > >> We the given dataset, my approach is indeed slow compared to st_union > >> approach (though precision for the st_dwithin clause must be adapted to > >> the current dataset. I took the following precision: 0.000001) > >> > >> The st_union method generates 18322 segments in 7318 ms, though the > >> final association between original lines and new segment is not done > >> here. > >> > >> With the query I gave, the st_dwithin part takes 11.7 sec on a recent > >> laptop machine (1.8 Ghz Intel Core I7, 1024 mb of ram for shared_buffer, > >> 512 for work_mem)... > >> > >> The complete query returns 17292 segments in 17956 ms. > >> > >> As the lines are almost already noded, it generates a lot of > >> intersection points coincident with one line ends. > >> > >> As you noted, intermediate temp tables may help here: > >> > >> I decomposed the query into intermediate steps and the performance is > >> about the same as with st_union : > >> > >> -- First creates temp table with intersection points > >> drop table if exists intergeom; > >> create temp table intergeom as > >> select l1.id as l1id, l2.id as l2id, > >> st_intersection(l1.geom, l2.geom) as geom > >> from bdaways l1 join bdaways l2 on (st_dwithin(l1.geom, l2.geom, > >> 0.000001)) > >> where l1.id <> l2.id ; > >> > >> -- keeps only true intersection points > >> -- must handle the case where lines intersects at a linestring... > >> > > > > Would it make sense to take all the geometryType(geom) <> 'LINESTRING' and > > just add the end points to the intergeom table. I think this would then add > > break points at the ends of overlapping segments insure that they get > > divided at those points and it would also add extraneous points at the > > other end, but these would get filter out later. So add these two lines > > here? > > > > You meant, geometryType(geom) = 'LINESTRING ? > Yes it would make sense. I will try. > > > > insert intergeom (l1id, l2id, geom) > > values (l1id, l2id, st_startpoint(geom)) > > where geometryType(geom) <> 'LINESTRING'; > > > > insert intergeom (l1id, l2id, geom) > > values (l1id, l2id, st_endpoint(geom)) > > where geometryType(geom) <> 'LINESTRING'; > > > > > > > > delete from intergeom where geometryType(geom) <> 'POINT'; > >> > > > > I think some OSM data from a small central american country might make a > > good test case. I have not tried any of these: > > > > http://downloads.cloudmade.**com/americas/central_america<http://downloads.cloudmade.com/americas/central_america> > > > > but I notice that if you click into a country like Belize or Costa Rica, > > there is a shapefile for the country near the bottom of the list or you can > > pick just a city for a smaller data set. > > > > I'll start packaging this into a stored procedure. This is an awesome job. > > > > Thanks. > I tested the costa rica highways dataset and found the following results: > > The query runs in about 28 seconds to process 29187 ways that are not at > all connected between them. > It generates 48415 segments. > > As an extra steps, all original ways that did not need to be cut by > intersection points have to be added to the final results: these lines > are already clean but were dismissed during the cutting process: > > insert into res (gid, sub_id, geom, type) > select ways.gid, 1 as sub_id, ways.geom, geometryType(ways.geom) > from ways > where not exists ( > select res.gid from res where res.gid = ways.gid > ); > > The result looks good, as all ways seems to be segmented (cf. > http://sd-38177.dedibox.fr/img1.png and http://sd-38177.dedibox.fr/img2.png: > black: initial lines, blue: new segments. Lines ends are shown as circles). > > A problem remains with non-simple or dirty lines as they self-intersect or > have several ends and so are not noded properly. > Though it should be simple to process them by identifying the problem. here > is an example of such a line: http://sd-38177.dedibox.fr/img3.png > > I go looking at the procedure ;) > > Nicolas _______________________________________________ postgis-users mailing list [email protected] http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
