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
[email protected]
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
