[postgis-users] Need a method for "noding" a street network

[Stephen Woodbridge]

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?

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

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,
   -Steve

> -- second temp table with locus (index of intersection point on the line)
> -- to avoid updating the previous table
> -- we keep only intersection points occuring onto the line, not at one
> of its ends
> drop table if exists inter_loc;
> create temp table inter_loc as (
> select l1id, l2id, st_line_locate_point(l.geom, i.geom) as locus
> from intergeom i left join bdaways l on (l.id <http://l.id> = i.l1id)
> where st_line_locate_point(l.geom, i.geom) <> 0 and
> st_line_locate_point(l.geom, i.geom) <> 1
> );
>
> -- index on l1id
> create index inter_loc_id_idx on inter_loc(l1id);
>
> -- Then computes the intersection on the lines subset, which is much
> smaller than full set
> -- as there are very few intersection points
> drop table if exists res;
> create table res as
> with cut_locations as (
> select l1id as lid, locus
> from inter_loc
> -- then generates start and end locus for each line that have to be cut
> buy a location point
> UNION ALL
> select i.l1id  as lid, 0 as locus
> from inter_loc i left join bdaways b on (i.l1id = b.id <http://b.id>)
> UNION ALL
> select i.l1id  as lid, 1 as locus
> from inter_loc i left join bdaways b on (i.l1id = b.id <http://b.id>)
> order by lid, locus
> ),
> -- we generate a row_number index column for each input line
> -- to be able to self-join the table to cut a line between two
> consecutive locations
> loc_with_idx as (
> select lid, locus, row_number() over (partition by lid order by locus)
> as idx
> from cut_locations
> )
> -- finally, each original line is cut with consecutive locations using
> linear referencing functions
> select l.id <http://l.id>, 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 (lid) join bdaways l
> on (l.id <http://l.id> = loc1.lid)
> where loc2.idx = loc1.idx+1
> -- keeps only linestring geometries
> and geometryType(st_line_substring(l.geom, loc1.locus, loc2.locus)) =
> 'LINESTRING';
>
> The total time is 7727 ms and it generates 1865 new segments.
>
> I will see if some filtering clauses used here can be ported efficiently
> in the big query.
>
> Nicolas
>
>
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
>