利用openstreetmap数据构建road network(road->road segments)思路

 

 

说明：

　　（原创）在做地图匹配工作时，会遇到这样的问题：

　　有时候一条road很长很长，它与很多路都有交点。那么我们构建一个合理的road network时，就需要将每一条road分成一个个road segmen(路段)。

参考：

　　(2)workshop.pgrouting.org/chapters/com2pgrouting.html

　　(1)workshops.boundlessgeo.com/tutorial-routing/

步骤：

第一阶段：

　　（1）中用到了pgrouting，那么我们首先分析pgrouting.sql源码。其中pgr_nodeNetwork发挥了很大作用。

　　　　　　

　　左图中只有1和2两条road, 右图中road1被分解为3和4两个road.

　　这正是pgr_nodeNetwork的功用，它将每条road从与其他road交点的地方分解成多个segment。

　　那么它是如何做到的呢？

　　

    执行该句之后，会把有交点的留下来，放在intergeom中。
--    -- First creates temp table with intersection points
    p_ret = 'create temp table intergeom on commit drop as (
        select l1.' || quote_ident(n_pkey) || ' as l1id, 
               l2.' || quote_ident(n_pkey) || ' as l2id, 
	       l1.' || quote_ident(n_geom) || ' as line,
	       pgr_startpoint(l2.' || quote_ident(n_geom) || ') as source,
	       pgr_endpoint(l2.' || quote_ident(n_geom) || ') as target,
               st_intersection(l1.' || quote_ident(n_geom) || ', l2.' || quote_ident(n_geom) || ') as geom 
        from ' || pgr_quote_ident(intab) || ' l1 
             join ' || pgr_quote_ident(intab) || ' l2 
             on (st_dwithin(l1.' || quote_ident(n_geom) || ', l2.' || quote_ident(n_geom) || ', ' || tolerance || '))'||
        'where l1.' || quote_ident(n_pkey) || ' <> l2.' || quote_ident(n_pkey)||' and 
	st_equals(Pgr_startpoint(l1.' || quote_ident(n_geom) || '),pgr_startpoint(l2.' || quote_ident(n_geom) || '))=false and 
	st_equals(Pgr_startpoint(l1.' || quote_ident(n_geom) || '),pgr_endpoint(l2.' || quote_ident(n_geom) || '))=false and 
	st_equals(Pgr_endpoint(l1.' || quote_ident(n_geom) || '),pgr_startpoint(l2.' || quote_ident(n_geom) || '))=false and 
	st_equals(Pgr_endpoint(l1.' || quote_ident(n_geom) || '),pgr_endpoint(l2.' || quote_ident(n_geom) || '))=false  )';
    raise debug '%',p_ret;	
    EXECUTE p_ret;	

　　上面的代码会创建一个临时表intergeom，其中的每一行都是 一条road 与 另一条road相交的case。具体包含了两条roadid, 其中一条road的linegeom,以及交点geom.  其中这个交点geom是st_intersection生成的。

执行后，会找到每条被中间分开的line的位置，是浮点表示的。
--HAD TO CHANGE THIS QUERY
    p_ret= 'create temp table inter_loc on commit drop as ( select * from (
        (select l1id, l2id, ' || vst_line_locate_point || '(line,source) as locus from intergeom)
         union
        (select l1id, l2id, ' || vst_line_locate_point || '(line,target) as locus from intergeom)) as foo
        where locus<>0 and locus<>1)';
    raise debug  '%',p_ret;	
    EXECUTE p_ret;	

    -- index on l1id
    create index inter_loc_id_idx on inter_loc(l1id);

 　　这一段代码中用到了st_linelocatepoint函数，例如road1和road2交于pointA, 那么st_linelocatepoint(road1,pointA)将返回一个(0~1)float值，表示pointA在road1上的位置。

　　

　　在这一段代码中，使用st_linesubstring来生成road segment。

　　该函数生成road segment的方式的问题在于使用st_intersection生成交点geom.

　　下面是我对st_intersection的测试

　　

　　可以看到line(0 2 0, 4 2 0)和line(2 0 1, 2 4 1)并不在一个面上，但st_intersection给出二者有交点。很奇怪。如果这个在有立交桥的情况下，下面的路与立交桥上的路将会有交点，由此会产生错误。

　　对这个问题，我也做了实际数据的测试，

　　1）在osm上下载了shp和osm两种格式的数据，直接用pgadmin可视工具导入beijing_line.shp得到bjlines数据表，利用osmosis导入osm格式得到多张表其中

          包括Nodes。

　　2）利用qgis工具导入bjlines表，从中我找到了gid=14228的line，它是京沪线rail，z-order为0；有找到了gid=13868的line，它是一条highway=residential

　　　  的路。这两者从图形上看他们是有交点的，但从真实内容看，二者不可能有交点（z-order不同，一个是公路一个是铁路）

　　3）利用st_intersection去查询他们，返回来一个交点，去nodes表中查询这个交点是没有的。说明st_intersection来分隔道路是有问题的。

　　

 

第二阶段

 　　由于存在第一阶段最后的顾虑，我看了osm2pgrouting的源码，其更加符合openstreetmap一书中所说的，如果两条road没有share node，那么两条road就没有交点。osm2pgrouting遵循这个原则，来进行分解road。

　　其原理是：每一条road都会引用不小于2个node，由此，对每一条road都去检索其引用的node的被引用数，如果该数大于1，就说明这个node被多个road引用了，此时从这个node就将原始road分解开。

　　（2）中关于osm2pgrouting源代码的mapconfig也是值得一看的，它和（1）中的关于过滤非road的内容有些重叠，但不失为一个好点。