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

[Stephen Woodbridge]

Nicolas,

Here is the code in a stored procedure (attached). Interestingly it runs 
in 28sec in the stored procedure and 55 sec when I run your SQL. But I 
think there is a problem with the logic in that both are only returning 
2130 row res table.

-Steve



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 <http://l1.id> as l1id, l2.id <http://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 <http://l1.id> <> l2.id <http://l2.id>;
>
> -- keeps only true intersection points
> -- must handle the case where lines intersects at a linestring...
> delete from intergeom where geometryType(geom) <> 'POINT';
>
> -- 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
>

-------------- next part --------------
?
CREATE OR REPLACE FUNCTION pgr_node_table(intab text, n_pkey text, n_geom text,  outtab text, tol double precision)
    RETURNS text AS
$BODY$
DECLARE

BEGIN

    -- First creates temp table with intersection points
    drop table  if exists intergeom;
    EXECUTE 'create temp table intergeom as 
        select l1.' || quote_ident(n_pkey) || ' as l1id, 
               l2.' || quote_ident(n_pkey) || ' as l2id, 
               st_intersection(l1.' || quote_ident(n_geom) || ', l2.' || quote_ident(n_geom) || ') as geom 
        from ' || quote_ident(intab) || ' l1 
             join ' || quote_ident(intab) || ' l2 
             on (st_dwithin(l1.' || quote_ident(n_geom) || ', l2.' || quote_ident(n_geom) || ', ' || tol || '))
        where l1.' || quote_ident(n_pkey) || ' <> l2.' || quote_ident(n_pkey);

    -- must handle the case where lines intersects at a linestring...

    EXECUTE 'insert into intergeom (l1id, l2id, ' || quote_ident(n_geom) || ')
        select l1id, l2id, st_startpoint(' || quote_ident(n_geom) || ')
        from intergeom where geometryType(' || quote_ident(n_geom) || ') <> ''LINESTRING'' ';

    EXECUTE 'insert into intergeom (l1id, l2id, ' || quote_ident(n_geom) || ')
        select l1id, l2id, st_endpoint(' || quote_ident(n_geom) || ')
        from intergeom where geometryType(' || quote_ident(n_geom) || ') <> ''LINESTRING'' ';

    -- keeps only true intersection points
    EXECUTE 'delete from intergeom where geometryType(' || quote_ident(n_geom) || ') <> ''POINT'' ';

    -- 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;
    EXECUTE 'create temp table inter_loc as (
        select l1id, l2id, st_line_locate_point(l.' || quote_ident(n_geom) || ', i.' || quote_ident(n_geom) || ') as locus
        from intergeom i left join ' || quote_ident(intab) || ' l on (l.' || quote_ident(n_pkey) || ' = i.l1id)
        where st_line_locate_point(l.' || quote_ident(n_geom) || ', i.' || quote_ident(n_geom) || ') <> 0 
              and st_line_locate_point(l.' || quote_ident(n_geom) || ', i.' || quote_ident(n_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
   EXECUTE 'drop table if exists ' || quote_ident(outtab);
   EXECUTE 'create table ' || quote_ident(outtab) || ' 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 ' || quote_ident(intab) || ' b on (i.l1id = b.' || quote_ident(n_pkey) || ')
           UNION ALL
           select i.l1id  as lid, 1 as locus
           from inter_loc i left join ' || quote_ident(intab) || ' b on (i.l1id = b.' || quote_ident(n_pkey) || ')
           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.' || quote_ident(n_pkey) || ', loc1.idx as sub_id, st_line_substring(l.geom, loc1.locus, loc2.locus) as ' || quote_ident(n_geom) || ' 
       from loc_with_idx loc1 join loc_with_idx loc2 using (lid) join ' || quote_ident(intab) || ' l on (l.' || quote_ident(n_pkey) || ' = loc1.lid)
       where loc2.idx = loc1.idx+1
           -- keeps only linestring geometries
           and geometryType(st_line_substring(l.' || quote_ident(n_geom) || ', loc1.locus, loc2.locus)) = ''LINESTRING'' ';

    RETURN 'OK';
END;
$BODY$
    LANGUAGE 'plpgsql' VOLATILE STRICT COST 100;

select pgr_node_table('bdaways', 'id', 'geom',  'noded', 0.000001);
-- Total query runtime: 28263 ms. only created table of 2130