Is it possible for you to share the dataset or a subset of it, to give it a try ?
Nicolas On 8 May 2013 22:04, Nicolas Ribot <[email protected]> wrote: > > > > On 8 May 2013 21:50, Stephen Woodbridge <[email protected]> wrote: > >> I tried to change my first with to just query and existing table, like: >> >> with lines as ( >> select id as gid, * from bdaways >> /* >> >> select 1 as gid, 'MULTILINESTRING((0 1,2 1))'::geometry as geom >> union all >> select 2 as gid, 'MULTILINESTRING((1 0,1 2))'::geometry as geom >> union all >> select 3 as gid, 'LINESTRING(1 1.5,2 2)'::geometry as geom >> */ >> ), >> >> but when I run it I get: >> >> ERROR: line_locate_point: 2st arg isnt a point >> >> ********** Error ********** >> >> ERROR: line_locate_point: 2st arg isnt a point >> SQL state: XX000 >> >> This looks like st_intersection(l1.geom, l2.geom) is not returning a >> point. > > > Yes it may return an empty geometryCollection. It could be filtered in the > cut_locations cte by adding a where clause geometryType = 'POINT'. > > >> >> On 5/8/2013 3:32 PM, Nicolas Ribot wrote: >> >>> Hi, >>> >>> Glad it helps. >>> >>> CTE are supported since 8.4 according to the doc: >>> http://www.postgresql.org/**docs/8.4/static/queries-with.**html<http://www.postgresql.org/docs/8.4/static/queries-with.html> >>> >> >> ok, cool, will give that a try also. >> >> >> Otherwise, the CTE's can be replaced by subqueries. >>> The stored procedure may also help by creating temp tables instead of >>> chaining subqueries, though I don't know if it will run faster. >>> >> >> One thing I noticed is that CTE's can not be indexed, so I might get >> better performance is I create tables and index them based on the needs of >> successive queries in a stored procedure. I'll need to play with this a bit >> to figure out what works best. >> >> >> Using topology should be pretty simple in your case: build a new >>> topology based on the lines table, then query the topology.edge table, >>> keeping initial line gid. It may be worth trying it. >>> >> >> I need to find the how to for working with topologies. I have seem some >> in the past and just need to give it a try now that I have a real use case >> for it. >> >> -Steve >> >> Nicolas >>> >>> >>> On 8 May 2013 21:01, Stephen Woodbridge <[email protected] >>> <mailto:woodbri@swoodbridge.**com <[email protected]>>> wrote: >>> >>> Hi Nicolas, >>> >>> Wow! thank you for an excellent example. This is very help. Since I >>> want this to work on pg 8.4+, I'll convert this into a stored >>> procedure since I can't use CTE subqueries. >>> >>> Now I have some work cut out to do on this. :) >>> >>> Thanks again, >>> -Steve >>> >>> >>> On 5/8/2013 2:39 PM, Nicolas Ribot wrote: >>> >>> Hi Stephen, >>> >>> Building a topology would definitively help in this situation, >>> though it >>> may take some time on very large dataset I guess. >>> If you plan to use some topological functions on the dataset in >>> addition >>> with pgRouting functions, it may be worth the effort. >>> >>> Concerning st_union and its magic "segmentize" feature, would it >>> be >>> possible to divide the initial set of lines into smaller areas >>> and >>> process these subsets to avoid filling up the memory ? >>> >>> Looking at this subject recently (cutting lines by points, cf. >>> http://trac.osgeo.org/postgis/**__wiki/__** >>> UsersWikiSplitPolygonWithPoint**__s<http://trac.osgeo.org/postgis/__wiki/__UsersWikiSplitPolygonWithPoint__s> >>> <http://trac.osgeo.org/**postgis/wiki/** >>> UsersWikiSplitPolygonWithPoint**s<http://trac.osgeo.org/postgis/wiki/UsersWikiSplitPolygonWithPoints> >>> >) >>> >>> I >>> found that linear referencing functions can help in such a case. >>> >>> The principle is to get the location index of intersection >>> points for >>> each line, and then to cut this line by its locations, using >>> st_line_substring. >>> It appears to be very efficient, using st_dwithin to trigger >>> spatial >>> index, then joining on the lines primary keys, which should be >>> fast. >>> >>> In your usecase, intersection nodes between lines have to be >>> identified >>> before their locations can be computed. >>> >>> Concerning the tolerance, I'm pretty sure snapping the input >>> dataset to >>> a grid would help to run a precise st_intersection between lines. >>> >>> Based on the linestring sample data, here is the query using >>> linear >>> referencing. It uses CTE subqueries to identify each step: >>> >>> with lines as ( >>> select 1 as gid, 'MULTILINESTRING((0 1,2 >>> 1))'::geometry as geom >>> union all >>> select 2 as gid, 'MULTILINESTRING((1 0,1 >>> 2))'::geometry as geom >>> union all >>> select 3 as gid, 'LINESTRING(1 1.5,2 2)'::geometry as >>> geom >>> ), >>> -- multilinestrings are dumped into simple objects >>> -- if multilinestrings have several parts, one should generate a >>> unique >>> id based >>> -- on their gid and path into the collection. >>> dumped_lines as ( >>> select gid, (st_dump(l.geom)).geom >>> from lines l >>> ), >>> -- This query computes the locations, for each input line, of the >>> intersection points with other lines. >>> -- this will be used to cut lines based on these locations. >>> -- to be able to cut lines from their beginning to their end, we >>> generate the 0 and 1 location index >>> cut_locations as ( >>> select l1.gid as lgid, st_line_locate_point(l1.geom, >>> st_intersection(l1.geom, l2.geom)) as locus >>> from dumped_lines l1 join dumped_lines l2 on (st_dwithin(l1.geom, >>> l2.geom, 0.01)) >>> where l1.gid <> l2.gid >>> -- then generates start and end locus for each line, to be able >>> to cut them >>> UNION ALL >>> select l.gid as lgid, 0 as locus >>> from dumped_lines l >>> UNION ALL >>> select l.gid as lgid, 1 as locus >>> from dumped_lines l >>> order by lgid, locus >>> ), >>> -- This query generates a row_number index column for each input >>> line >>> and intersection point. >>> -- This index will be used to self-join the table to cut a line >>> between >>> two consecutive locations >>> -- (idx, idx+1) pairs. >>> -- window function is used to generate the index inside each >>> line partition >>> loc_with_idx as ( >>> select lgid, locus, row_number() over (partition by lgid order >>> by locus) >>> as idx >>> from cut_locations >>> ) >>> -- finally, each original line is cut with consecutive locations >>> using >>> linear referencing function. >>> -- a filtering is done to eliminate points produced when lines >>> connect >>> at their ends >>> select l.gid, loc1.idx as sub_id, st_line_substring(l.geom, >>> loc1.locus, >>> loc2.locus) as geom , >>> st_geometryType(st_line___**substring(l.geom, loc1.locus, >>> >>> loc2.locus)) as type >>> from loc_with_idx loc1 join loc_with_idx loc2 using (lgid) join >>> dumped_lines l on (l.gid = loc1.lgid) >>> where loc2.idx = loc1.idx+1 >>> -- filter out point geometries occuring if intersection point is >>> at >>> line's start or end point. >>> -- there must be a faster way to filter out theses geometries. >>> and st_geometryType(st_line___**substring(l.geom, loc1.locus, >>> >>> loc2.locus)) >>> <> 'ST_Point'; >>> >>> >>> A new unique ID key can be computed based on line gid and subgid >>> generated by the query. >>> Initial line attributes can be moved to the new segments using >>> the line >>> gid key. >>> >>> Nicolas >>> >>> >>> On 8 May 2013 16:27, Stephen Woodbridge <[email protected] >>> <mailto:woodbri@swoodbridge.**com <[email protected]>> >>> <mailto:woodbri@swoodbridge.__**com >>> >>> <mailto:woodbri@swoodbridge.**com <[email protected]>>>> >>> wrote: >>> >>> Hi all, >>> >>> This question comes up reasonably often on the pgRouting >>> list and >>> has been posted he on occasion under titles like "How to >>> break >>> streets at intersections?" >>> >>> It seems to me that this would be a good function to create >>> in >>> either postgis or pgrouting. >>> >>> THE PROBLEM: >>> >>> I have a table of 10's of thousands of street segments to >>> 10's of >>> millions of street segments. These street segments are >>> LINSTRING or >>> MULTILINESTRING geometries with some arbitrary number of >>> attribute >>> columns. The geometries may cross one another and are not >>> noded >>> correctly for use with pgRouting. >>> >>> THE RESULTS: >>> >>> We want to process the table and create a new table with >>> the same >>> structure (see comment about primary key below), and in the >>> new >>> table all the geometries are broken at intersections and >>> all the new >>> pieces of the original segment that have been broken have >>> the >>> original attributes propagated to them. So if the original >>> segment >>> has column foo='abc' and was split into 3 new segments, >>> each of the >>> three new segments would also have foo='abc'. The exception >>> to this >>> might be that the new table needs a new primary column as >>> the old >>> primary key will now be duplicated for the multiple parts. >>> >>> POTENTIAL SOLUTIONS: >>> >>> 1. I think one way to do this would be to create a topology >>> and load >>> the table into it, then extra a new table from the topology. >>> Although I'm not sure of the specifics for doing this or the >>> efficency of doing it this way. >>> >>> 2. Another way seems to be using a query like: >>> >>> select (st_dump(bar.the_geom)).* from ( >>> select st_union(foo.the_geom) as the_geom from mytable >>> foo >>> ) as bar; >>> >>> And then taking each of the dump.geom objects and using >>> st_contains >>> to find which original segment it belonged to so we can >>> move the >>> attributes to the new segment. This method also loose any >>> association to the original record and forces the use of >>> st_contains >>> to re-associate the new segments to the original segments. >>> >>> My concern with this is that the st_union has to load the >>> whole >>> table which may be 10's of millions of street segments and >>> this will >>> likely be a memory problem. Also running the st_contains() >>> does not >>> seems to me to be optimal. >>> >>> 3. Is there a good recipe for doing this somewhere that I >>> have not >>> found? or other better approaches to this problem? >>> >>> What would be the best way to add tolerance to the problem? >>> using >>> snap to grid? >>> >>> Thoughts on how to do this efficiently? >>> >>> Since I'm working on the pgRouting 2.0 release I thought >>> this might >>> be a nice function to add to that if we can come up with a >>> generic >>> way to do this. >>> >>> Thanks, >>> -Steve >>> >>> >>> -- Example to demonstrate st_union above >>> select st_astext((st_dump(bar.the____**_geom)).geom) from ( >>> >>> >>> select st_union(foo.the_geom) as the_geom from ( >>> select 'MULTILINESTRING((0 1,2 1))'::geometry as >>> the_geom >>> union all >>> select 'MULTILINESTRING((1 0,1 2))'::geometry as >>> the_geom >>> union all >>> select 'LINESTRING(1 1.5,2 2)'::geometry as >>> the_geom >>> ) as foo >>> ) as bar; >>> >>> "LINESTRING(1 1.5,2 2)" >>> "LINESTRING(1 0,1 1)" >>> "LINESTRING(1 1,1 1.5)" >>> "LINESTRING(1 1.5,1 2)" >>> "LINESTRING(0 1,1 1)" >>> "LINESTRING(1 1,2 1)" >>> ______________________________**_____________________ >>> >>> postgis-users mailing list >>> [email protected] >>> >>> <mailto:postgis-users@lists.**osgeo.org<[email protected]> >>> > >>> <mailto:postgis-users@lists.__**osgeo.org <http://osgeo.org> >>> >>> <mailto:postgis-users@lists.**osgeo.org<[email protected]> >>> >> >>> http://lists.osgeo.org/cgi-___**_bin/mailman/listinfo/postgis-** >>> ____users<http://lists.osgeo.org/cgi-____bin/mailman/listinfo/postgis-____users> >>> <http://lists.osgeo.org/cgi-__**bin/mailman/listinfo/postgis-_** >>> _users<http://lists.osgeo.org/cgi-__bin/mailman/listinfo/postgis-__users> >>> > >>> >>> <http://lists.osgeo.org/cgi-__**bin/mailman/listinfo/postgis-_** >>> _users<http://lists.osgeo.org/cgi-__bin/mailman/listinfo/postgis-__users> >>> <http://lists.osgeo.org/cgi-**bin/mailman/listinfo/postgis-** >>> users <http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users>>> >>> >>> >>> >>> >>> >>> >>> ______________________________**___________________ >>> postgis-users mailing list >>> [email protected] <mailto:postgis-users@lists.** >>> osgeo.org <[email protected]>> >>> http://lists.osgeo.org/cgi-__**bin/mailman/listinfo/postgis-_** >>> _users<http://lists.osgeo.org/cgi-__bin/mailman/listinfo/postgis-__users>< >>> http://lists.osgeo.org/cgi-**bin/mailman/listinfo/postgis-**users<http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> >>> > >>> >>> >>> ______________________________**___________________ >>> postgis-users mailing list >>> [email protected] <mailto:postgis-users@lists.** >>> osgeo.org <[email protected]>> >>> http://lists.osgeo.org/cgi-__**bin/mailman/listinfo/postgis-_** >>> _users<http://lists.osgeo.org/cgi-__bin/mailman/listinfo/postgis-__users> >>> >>> <http://lists.osgeo.org/cgi-**bin/mailman/listinfo/postgis-**users<http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> >>> > >>> >>> >>> >>> >>> ______________________________**_________________ >>> postgis-users mailing list >>> [email protected] >>> http://lists.osgeo.org/cgi-**bin/mailman/listinfo/postgis-**users<http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> >>> >>> >> ______________________________**_________________ >> postgis-users mailing list >> [email protected] >> http://lists.osgeo.org/cgi-**bin/mailman/listinfo/postgis-**users<http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users> >> > >
_______________________________________________ postgis-users mailing list [email protected] http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
