Saturday, 12 September 2015

pgrouting - Is it possible to make a routable graph from Polygons?


I would like to know if it is possible to make a graph with pgrouting 1.05 on polygons. What I mean is that I have to make a graph on a database without roads but with polygon geometry. I would like to make a graph on this geometry.


Thanks for your help and I'm sorry for my english



Answer



Okay, I looked further into the idea of Steve above. I'll try to demontrate his idea in a QGIS/PostGIS/pgrouting environment. You will get results such as this:


enter image description here


First, let's assume you have a geodata table with your shelves/obstacles looking like this (I made them up for this purpose):


enter image description here


Make sure your shelf data is in a projected coordinate system with meters as units (I used a Mercator CRS(SRID=3785). If it's in lat/lng you have to reproject your data.


Then, open up your database connection and import your shelves geodata (either use a QGIS plugin like SPIT or the command line tool shp2pgsql). Now, create a point grid. Use the function from here: How to create regular point grid inside a polygon in Postgis? :



CREATE OR REPLACE FUNCTION makegrid(geometry, integer)
RETURNS geometry AS
'SELECT ST_Collect(ST_SetSRID(ST_POINT(x,y),ST_SRID($1))) FROM
generate_series(floor(st_xmin($1))::int, ceiling(st_xmax($1)-st_xmin($1))::int, $2) as x
,generate_series(floor(st_ymin($1))::int, ceiling(st_ymax($1)-st_ymin($1))::int,$2) as y
where st_intersects($1,ST_SetSRID(ST_POINT(x,y),ST_SRID($1)))'
LANGUAGE sql;

Following SQL gives you a table containing a point grid with spacing of one meter between each point:


CREATE SEQUENCE seq;

CREATE TABLE grid AS
SELECT nextval('seq') AS gid, (ST_Dump(makegrid(ST_Transform(ST_SetSRID(ST_Extent(the_geom),3785),3785),1))).geom AS the_geom FROM shelves;

This grid only covers the extent of your shelves file:


enter image description here


Now remove every grid point which is intersected by/contained by a shelf (you could also do this directly in the grid query above) and create an updated table:


CREATE TABLE open_space_grid AS 
SELECT
grid.*
FROM

grid left join
shelves ON
ST_Intersects(grid.the_geom,shelves.the_geom)
WHERE shelves.gid is null;

this should give you this:


enter image description here


Now you have to connect the points to form a network. This query connects each point pair with a distance of one meter, forming a grid network.


CREATE TABLE edges AS
SELECT nextval('seq') AS gid, the_geom

FROM (
select ST_SetSRID(ST_MakeLine(a.the_geom, b.the_geom),3785) AS the_geom
FROM open_space_grid AS a,open_space_grid AS b) AS tmp WHERE ST_Length(ST_Transform(ST_SetSRID(tmp.the_geom,3785),3785)) = 1;

Then throw out all duplicates:


DELETE FROM edges 
WHERE gid IN (
SELECT e1.gid
FROM edges e1, edges e2
WHERE st_equals(e1.the_geom, e2.the_geom)

AND e1.gid < e2.gid
);

Now we are here:


enter image description here


Now create the routable network (Everything here taken from http://underdark.wordpress.com/2011/02/07/a-beginners-guide-to-pgrouting/)


CREATE OR REPLACE VIEW net_ext AS
SELECT *, ST_StartPoint(ST_SetSRID(the_geom,3785)) as startpoint, ST_EndPoint(ST_SetSRID(the_geom,3785)) as endpoint
FROM edges;


CREATE TABLE node AS
SELECT row_number() OVER (ORDER BY foo.p)::integer AS id,
foo.p AS the_geom
FROM (
SELECT DISTINCT net_ext.startpoint AS p FROM net_ext
UNION
SELECT DISTINCT net_ext.endpoint AS p FROM net_ext
) foo
GROUP BY foo.p;


CREATE TABLE network AS
SELECT a.*, b.id as start_id, c.id as end_id, ST_Length(ST_MakeLine(a.startpoint, a.endpoint)) as length
FROM net_ext AS a
JOIN node AS b ON a.startpoint = b.the_geom
JOIN node AS c ON a.endpoint = c.the_geom;

Now try it out (e.g. execute the query below in QGIS via pgQuery). Input your start and end node id's (in my case, those are 1 and 590):


SELECT * 
FROM network
JOIN

(SELECT * FROM shortest_path('
SELECT gid AS id,
start_id::int4 AS source,
end_id::int4 AS target,
length AS cost
FROM network',
1,
590,
false,
false)) AS route

ON
network.gid = route.edge_id;

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...