Bases de Datos Geoespaciales
En esta práctica vamos a hacer análisis de rutas utilizando la extensión pgrouting para PostGis. Lo primero que tenemos que hacer es instalar dicha extensión:
pgrouting-pg92-binaries-2.0.0w64gcc48.zip y copia los archivos en
la ruta de instalación de postgres: C:\Program Files\PostgreSQL\9.2Para probar que la extensión quedó instalada correctamente, desde una consola de sql (conectada a alguna de las bases de datos que hemos utilizado), de pgAdmin escribe:
CREATE EXTENSION pgrouting;
Si la consulta no regresa ningún error, la extensión quedó instalada correctamente.
Lo primero que tenemos que hacer es, como siempre, subir nuestro
shape a la base de datos. Los datos que vamos a usar están en el
archivo red_utm.shp. Súbelos en una table que se llame calles
##Parte II: Creación de la topología y pesos
Una vez que tenemos las calles que nos interesan, vamos a usar pgrouting para crear la topología de red sobre las calles. Lo primero que necesitamos es agregar dos campos para almacenar los nodos de orígen y destino de cada segmento:
alter table calles add column source integer;
alter table calles add column target integer;
Ahora, vamos a llamar a la función select pgr_createTopology('lines', tolerancia, 'geom', 'id');, para crear los nodos y asignar los
identificadores correspondientes. Los argumentos de la función son los
siguientes:
En nuestro caso:
select pgr_createTopology('calles', 0.0001, 'geom', 'gid');
Como pueden ver, esta función crea la tabla calles_vertices_pgr,
idealmente esta tabla contiene todos los nodos de la red, examínenla
en Qgis.
Ahora, vamos a asignar algunos pesos a las calles, para eso podemos usar (igual que en la práctica de redes de Análisis Espacial), la categoría vial y estimar una velocidad promedio de recorrido a partir de eso. Primero agregamos las columnas que nos faltan:
alter table calles add column speed float;
alter table calles add column cost float;
Ahora vamos a popular los valores de speed utilizando la columna catvial (esto es un ejemplo, fíjate a qué tipo de calle corresponde cada categoría y pon un valor razonable):
update calles set speed =
case
when catvial = 'CUARTO ORDEN'then 10
when catvial = 'TERCER ORDEN' then 20
when catvial = 'SEGUNDO ORDEN' then 30
when catvial = 'PRIMER ORDEN' then 50
else null
end
Ahora podemos calcular la columna costo usando el tiempo de viaje:
update red_calles_utm set cost = ((st_length(geom)/1000)/speed)*(60)
Finalmente, para terminar esta parte del ejercicio, vamos a calcular una ruta usando dos algoritmos diferentes, primero vamos a usar [Dijkstra] (http://docs.pgrouting.org/2.0/en/src/dijkstra/doc/index.html#pgr-dijkstra):
select c.gid, c.geom from calles c,
(SELECT seq, id1 AS node,
id2 AS edge, cost FROM pgr_dijkstra(
' SELECT gid AS id,
source::integer,
target::integer,
cost::double precision AS cost
FROM calles', 300032, 241417, true, true)) as ruta
where c.gid = ruta.edge
Ahora vamos a utilizar A*, para este algoritmo (heurístico), necesitamos agregar cuatro nuevas columnas y popularlas:
ALTER TABLE calles
ADD COLUMN x1 double precision,
ADD COLUMN y1 double precision,
ADD COLUMN x2 double precision,
ADD COLUMN y2 double precision;
UPDATE calles SET
x1 = ST_X(ST_startPoint(ST_GeometryN(geom,1))),
y1 = ST_Y(ST_startPoint(ST_GeometryN(geom,1))),
x2 = st_x(st_endpoint(ST_GeometryN(geom,1))),
y2 = st_y(st_endpoint(ST_GeometryN(geom,1)));
Ahora sí, podemos utilizar el algoritmo A*:
select c.gid, c.geom from calles c,
(SELECT seq, id1 AS node, id2 AS edge, cost
FROM pgr_astar(
' SELECT
gid AS id,
source::integer,
target::integer,
cost::double precision AS cost,
x1,
y1,
x2,
y2
FROM calles', 162867, 163952, true, false) ) as ruta
where c.gid =ruta.edge
Jueguen un rato con los nodos de inicio y fin, con lo algoritmos de rutas, investiguen y, finalmente, intenten contestar las siguientes preguntas:
En esta parte de la práctica, vamos a utilizar una red de calles extraida de OpenStreetMap, estos son datos contribuidos por usuarios, una especie de wikipedia para cartografía digital. Una ventaja de OSM es que desde un principio fue pensado como una fuente de datos para calcular rutas, de modo que su estructura permite construir una red topológica de manera natural (por suspuesto, tiene la desventaja de ser VGI, lo que nos puede hacer dudar de su precisión, validez, etc.).
El proceso para importar la red de OSM a postgres es demasiado largo como para hacerlo en el taller, entonces trabajaremos a partir de un respaldo de una base preparada con anticipación. De cualquier modo, si te interesa saber cómo utilizar los datos de OSM en pgrouting, el proceso involucra dos etapas:
Para importar los datos de esta práctica necesitas crear una nueva
base de datos, digamos, red_osm. No es necesario que le agregues las
extensiones de PostGis y pgrouting, el respaldo ya las incluye (claro,
sólo si están ya instaladas en el servidor). Una vez que hayas creado
la base de datos puedes, desde pgAdmin, dar botón derecho y
seleccionar la opción "Restaurar", navega hasta el archivo
osm_mex.backup y selecciónalo. Listo! tenemos una base de datos
lista para trabajar.
La base de datos que acabamos de crear tiene las siguientes tablas:
Schema | Name | Type | Owner
--------+--------------------------+----------+-------
public | classes | table | user
public | nodes | table | user
public | relation_ways | table | user
public | relations | table | user
public | types | table | user
public | way_tag | table | user
public | ways | table | user
public | ways_car | table | user
public | ways_vertices_pgr | table | user
las tablas ways y ways_vertices_pgr son las que contienen los segmentos y los nodos respectivamente. Dentro de la tabla ways vas a encontrar las columnas source y target cuyo significado ya debes de conocer bien. Además puedes notar que hay una columna llamada r_cost, esta representa el costo de recorrer la calle en sentido contrario (sí, esta red tiene los sentidos de la calle bien hechos!). La columna to_cost, que viene vacía sirve para que nosotros almacenemos un costo por defecto para la red.
Ahora bien, la red que importamos desde OSM contiene segmentos que no
corresponden a calles (ríos, canales, etc.) o bien segmentos por donde
no pueden circular automóviles, entonces, para el primer ejercicio
vamos a utilizar la tabla ways_car que es un extracto de ways que
contiene sólo los segmentos que corresponden a calles por donde pueden
circular automóviles.
Para ir agarrando familiaridad con la red, calculemos una ruta usando los algoritmos que ya conocemos:
select c.gid, c.the_geom from ways c,
(SELECT * FROM pgr_dijkstra(
' SELECT
gid AS id,
source::integer,
target::integer,
cost_s::double precision AS cost,
reverse_cost_s::double precision AS reverse_cost
FROM ways',
36198, 2064, directed := false)) as ruta
where c.gid = ruta.edge
Antes de usar el siguiente algoritmo, vamos a ver lo que nos regresa Dijkstra (antes de unirlo con la geometría):
SELECT * FROM pgr_dijkstra(
' SELECT gid AS id,
source::integer,
target::integer,
cost_s::double precision AS cost,
reverse_cost_s::double precision AS reverse_cost
FROM ways',
36198, 2064, directed := false)
seq | path_seq | node | edge | cost | agg_cost
-----+----------+--------+--------+-------------------+------------------
1 | 1 | 36198 | 55424 | 4.54431199114601 | 0
2 | 2 | 57113 | 5516 | 0.569442072315959 | 4.54431199114601
3 | 3 | 20607 | 30515 | 31.9701056592301 | 5.11375406346197
4 | 4 | 62719 | 55423 | 2.23134388147577 | 37.0838597226921
5 | 5 | 30394 | 5515 | 2.70381549531985 | 39.3152036041679
6 | 6 | 26561 | 5514 | 3.87482926327961 | 42.0190190994877
7 | 7 | 12693 | 5513 | 13.6339632399199 | 45.8938483627673
8 | 8 | 34693 | 5512 | 0.62428992858025 | 59.5278116026872
9 | 9 | 1248 | 5511 | 0.569880735703789 | 60.1521015312674
10 | 10 | 10730 | 5510 | 16.2132749069428 | 60.7219822669712
11 | 11 | 52478 | 5509 | 66.1974711848956 | 76.9352571739141
12 | 12 | 64210 | 30514 | 10.1381957558963 | 143.13272835881
13 | 13 | 15936 | 80335 | 2.7172145482802 | 153.270924114706
14 | 14 | 17178 | 80334 | 13.3013325757714 | 155.988138662986
15 | 15 | 79450 | 55422 | 11.0730604379678 | 169.289471238758
16 | 16 | 32573 | 5508 | 17.3805523793542 | 180.362531676725
17 | 17 | 16538 | 30513 | 7.73287188110276 | 197.74308405608
18 | 18 | 5355 | 52159 | 0.935288256531237 | 205.475955937182
19 | 19 | 8199 | 240078 | 5.00239753364591 | 206.411244193714
20 | 20 | 185946 | 234915 | 7.24769785967377 | 211.413641727359
21 | 21 | 175987 | 237612 | 6.27383758313728 | 218.661339587033
Como pueden ver nos regresa los ids de los segmentos y nodos que vamos
atravesando, así como los costos y costos acumumulados. Lo que es importante
ver es la columna path_seq, que nos permite ordenar los segmentos en caso de
que el algoritmo no los regrese ordenados.
select c.gid, c.the_geom from ways c,
(SELECT seq, id1 AS node, id2 AS edge, cost
FROM pgr_astar(' SELECT
gid AS id,
source::integer,
target::integer,
cost::double precision AS cost,
reverse_cost::double precision AS reverse_cost,
x1,
y1,
x2,
y2
FROM ways',
36198, 2064, false, true)) as ruta
where c.gid = ruta.edge
Como puedes ver, hay dos diferencias con lo que hicimos el ejercicio anterior:
Estamos llamando a la función pgr_dijkstra con el argumento
directed := false, para indicarle que queremos usar la versión no-dirigida
de la gráfica, en lugar de pasar los argumentos false,false, como en el
ejercicio anterior. En principio estas dos formas son equivalentes,
pero al llamarlo como en este ejercicio nos aseguramos de usar la última
versión de la función
(aquí puedes ver la documentación).
Estamos usando la columna reverse_cost para indicarle al
algoritmo cuál es el costo de recorrer el segmento en sentido
opuesto. En caso de que la vía sea de un sólo sentido, el costo en
reversa es negativo para evitar que el algoritmo seleccione ese arco.
Ahora, vamos a utilizar como costo el tiempo de recorrido asumiendo una velocidad constante para cada tipo de via. Por ejemplo, utilicemos la velocidad máxima para cada segmento como base para calcular el tiempo de recorrido:
select (st_length(the_geom::geography)/1000)/maxspeed_forward as tiempo
from ways_car limit 100
Nota: cuando hacemos
st_length(the_geom::geography) estamos calculando la distancia del
segmento sobre el esferoide.
Con la consulta anterior tenemos el tiempo de recorrido (en horas) para cada segmento, ahora, esto lo podemos usar directamente como costo en el algoritmo de ruta:
select c.gid, c.the_geom from ways c,
(SELECT seq, id1 AS node, id2 AS edge, cost
FROM pgr_astar(
'SELECT
gid AS id,
source::integer,
target::integer,
(st_length(the_geom::geography)/1000)/maxspeed_forward::double precision AS cost,
reverse_cost::double precision AS reverse_cost,
x1,
y1,
x2,
y2
FROM ways',
36198, 2064, true, true)) as ruta
where c.gid = ruta.edge
Como pueden ver, la ruta en este caso es igual con ambos costos. Compliquemos las cosas un poco, supongamos que estamos en hora pico y que las velocidades se ven modificadas de la siguiente forma:
Primero vamos a calcular la nueva velocidad máxima para cada tipo de segmento:
select class_id,
case
when class_id in(101,102,103) then maxspeed_forward/8
when class_id in(106,107,108) then maxspeed_forward/4
else maxspeed_forward/2
end
from ways
Para simplificar las consultas siguientes, vamos a poner estos valores en una nueva columna:
alter table ways add column velocidad_pico float;
update ways set velocidad_pico =
case
when class_id in(101,102,103) then maxspeed_forward/8
when class_id in(106,107,108) then maxspeed_forward/4
else maxspeed_forward/2
end;
Ahora sí, vamos a calcular la ruta usando las nuevas velocidades (lo único que necesitamos cambiar es la velocidad que vamos a usar):
select c.gid, c.the_geom from ways c,
(SELECT seq, id1 AS node, id2 AS edge, cost
FROM pgr_astar(
'SELECT
gid AS id,
source::integer,
target::integer,
(st_length(the_geom::geography)/1000)/velocidad_pico::double precision AS cost,
reverse_cost::double precision AS reverse_cost,
x1,
y1,
x2,
y2
FROM ways',
36198, 2064, true, true)) as ruta
where c.gid = ruta.edge
Comparen las dos rutas y los tiempos de traslado en cada caso.
Ahora vamos a usar pgrouting para resolver el problema de encontrar el camino óptimo para un repartidor que tiene que visitar varias localizaciones en su ruta y regresar al lugar de origen. Matemáticamente, el problema consiste en encontrar un ciclo hamiltoniano mínimo en una gráfica dirigida y con pesos.
El primer paso es definir cuales son los lugares por donde debe pasar el agente, para esto vamos a seleccionar un conjunto de nodos de la red que servirán como los puntos de reparto y un nodo que será la base del repartidor.
Para resolver el problema vamos a utilizar el algoritmo pgr_tsp de pgrouting. Este algoritmo es bastante más complejo y funciona un poco diferente que los que hemos usado, en lugar de regrasarte la ruta entre todos los puntos, regresa el orden en el que estos deben ser visitados de acuerdo a su distancia euclidiana o a una matriz de distancia que nosotros definamos. Idealmente, la matriz de distancia la podríamos construir tomando todas las distancias entre nuestros nodos de interes, calculadas usando Dijkstra, por ejemplo. Sin embargo, para simplificar el problema, vamos a utilizar la versión más simple del algoritmo. De la documentación podemos ver que lo que necesitamos para correr el algoritmo es:
Entonces, nuestra consulta queda de la siguiente manera:
SELECT seq, id1, id2, round(cost::numeric, 2) AS cost
FROM pgr_tsp(
'select id::int, st_x(the_geom) as x, st_y(the_geom) as y FROM
ways_vertices_pgr where id in (36104,2099,26248,25170)',
36104, 25170)
Al ejecutar la consulta, lo que nos regresa es una tabla con la secuencia en la que tenemos que recorrer los nodos y el costo estimado en cada segmento del recorrido. Ahora visualicemos la secuencia:
select p.id, orden.seq, p.the_geom
from (SELECT seq, id1, id2,round(cost::numeric, 2) AS cost
FROM pgr_tsp(
'select id::int,
st_x(the_geom) as x,
st_y(the_geom) as y
FROM ways_vertices_pgr where
id in (36104,2099,26248,25170)', 36104, 25170)) as orden
join ways_vertices_pgr p
on p.id = orden.id2
Como pueden ver, el algoritmo no nos regresa la ruta que debemos seguir, sin embargo es posible obtenerla usando alguno de los algoritmos de ruteo que conocemos. El ejercicio es obtener la ruta completa del circuito y dibujarla en QGIS.