En este post se explicará cómo se realiza el cálculo en el servidor en el módulo CTE, el cual permite obtener los trayectos o caminos entre una dirección de Mataró y las N entidades más próximas.
Este cálculo se divide en 2 fases: cálculo de los caminos más cortos y la selección de los segmentos finales. Esto es debido a que la función pgr_withPointsKSP (la cual se utiliza para calcular los caminos más cortos) tiene asociado el problema de que devuelve los caminos del grafo de nodo a nodo, no teniendo en cuenta que en la mayor parte de los casos tanto el punto de inicio como los finales estarán situados en un punto intermedio de un segmento (una arista). Para obtener más información sobre este tema se puede consultar el post “Cobertura mitjançant el graf de trams de carrers (GTC)” de Josep López Xarbau.
Cálculo de los caminos más cortos
En esta primera fase se utiliza la función de PGRouting, pgr_withPointsKSP. Esta función busca los N caminos más cortos utilizando el algoritmo de Yen.
Todo el código que se presentará a continuación pertenece al módulo CTE, el cual está programado en Python y se comunica con la base de datos de Postgres.
En primer lugar lo que nos interesa realizar es unir todos los puntos involucrados (tanto los de destino como el de origen) en una misma tabla para poder prepararlos para la función pgr_withPointsKSP.
De este modo, primero se borra y se crea una nueva tabla para guardar dichos puntos.
drop = 'DROP TABLE IF EXISTS NecessaryPoints_'+Fitxer+';'
try:
cur.execute(drop)
conn.commit()
except:
print ("DROP TABLE ERROR 1")
create = 'CREATE TABLE NecessaryPoints_'+Fitxer+' (\n'
create += "\tpid serial primary key,\n"
create += "\tthe_geom geometry,\n"
create += "\tentitatID int8,\n"
create += "\tedge_id BIGINT,\n"
create += "\tfraction FLOAT,\n"
create += "\tnewPoint geometry);"
try:
cur.execute(create)
conn.commit()
except:
print ("CREATE TABLE NecessaryPoints ERROR")
Acto seguido, se añaden los puntos a la tabla.
insert = 'INSERT INTO NecessaryPoints_'+Fitxer
insert += ' (entitatID,the_geom) (SELECT 0, ST_Centroid("geom") the_geom from "dintreilla"'
insert += ' WHERE "Carrer_Num_Bis" = \''+CNB+'\');\n'
insert += 'INSERT INTO NecessaryPoints_'+Fitxer
insert += ' (entitatID, the_geom) (SELECT "id", ST_Centroid("geom") the_geom FROM"'
insert += self.dlg.comboCapaDesti.currentText() + '" ORDER BY "id");'
try:
cur.execute(insert)
conn.commit()
except:
print ("Insert Points NecessaryPoints ERROR")
Después se añade el id del tramo más próximo a cada punto, los puntos proyectados sobre el grafo y la fracción de segmento en la que se encuentran.
update = 'UPDATE NecessaryPoints_'+Fitxer
update += ' SET "edge_id"=tram_proper."tram_id"'
update += ' FROM (SELECT distinct on(Poi."pid") Poi."pid" AS Punt_id,Sg."id" AS Tram_id,'
update += ' ST_Distance(Sg."the_geom",Poi."the_geom") AS dist '
update += 'FROM "Xarxa_Prova" as Sg,NecessaryPoints_'+Fitxer+' AS Poi '
update += 'ORDER BY Poi."pid",ST_Distance(Sg."the_geom",Poi."the_geom"),Sg."id") tram_proper'
update += ' WHERE NecessaryPoints_'+Fitxer+'."pid"=tram_proper."punt_id";\n'
update += 'UPDATE NecessaryPoints_'+Fitxer
update += ' SET fraction = ST_LineLocatePoint(e.the_geom, NecessaryPoints_'+Fitxer+'.the_geom),'
update += 'newPoint = ST_LineInterpolatePoint(e."the_geom",'
update += ' ST_LineLocatePoint(e."the_geom", NecessaryPoints_'+Fitxer+'."the_geom"))'
update += ' FROM "Xarxa_Prova" AS e WHERE NecessaryPoints_'+Fitxer+'."edge_id" = e."id";\n'
try:
cur.execute(update)
conn.commit()
except:
print ("Update Points NecessaryPoints ERROR")
Ahora ya está todo preparado para poder realizar el cálculo. Se hace una consulta para poder generar una sentencia SQL que haga la búsqueda de todos los caminos más cortos a todos los puntos necesarios y después se añaden a una tabla llamada “Resultat”.
select = 'select * from NecessaryPoints_'+Fitxer+' order by pid'
cur.execute(select)
vec = cur.fetchall()
create = 'create local temp table "Resultat" as SELECT * FROM (\n'
for x in range (0,len(vec)):
if x < len(vec) and x >= 2:
create += 'UNION\n'
if x != 0:
if vec[x][4] == 1.0 or vec[x][4] == 0.0:
create +='select '+ str(x) +' AS routeID,'+ str(vec[x][2])
create +=' AS entitatID, * FROM pgr_withPointsKSP'
create +='(\'SELECT id, source, target, cost, reverse_cost '
create +='FROM "Xarxa_Prova" ORDER BY id\','
create +='\'SELECT pid, edge_id, fraction FROM NecessaryPoints_'
create +=Fitxer+'\',-1,' + str(vec[x][2])+',1)\n'
else:
create += 'select '+ str(x) +' AS routeID,'+ str(vec[x][2])
create +=' AS entitatID, * FROM pgr_withPointsKSP'
create +='(\'SELECT id, source, target, cost, reverse_cost '
create +='FROM "Xarxa_Prova" ORDER BY id\','
create +='\'SELECT pid, edge_id, fraction FROM NecessaryPoints_'
create +=Fitxer+'\',-1,-' + str(vec[x][0]) +',1)\n'
create += ')QW ORDER BY routeID, seq;'
drop = 'DROP TABLE IF EXISTS "Resultat";'
try:
cur.execute(drop)
conn.commit()
except:
print ("DROP TABLE ERROR 2")
try:
cur.execute(create)
conn.commit()
except:
print ("CREATE TABLE Resultat global ERROR")
A continuación se deberá solucionar el problema que se ha comentado en la introducción del post, es decir, se seleccionarán los segmentos que son inicio y final para añadirlos al resultado final.
Selección de los segmentos finales
Lo primero que se debe hacer es borrar y crear la tabla “Segments finals”, en la cual figurarán todos los caminos posibles que son principio y/o final.
drop = "DROP TABLE IF EXISTS \"SegmentsFinals\";"
try:
cur.execute(drop)
conn.commit()
except:
print ("DROP TABLE ERROR 1")
create = "CREATE local temp TABLE \"SegmentsFinals\" (\n"
create += "\trouteid int8,\n"
create += "\tedge int8,\n"
create += "\t\"edgeAnt\" int8,\n"
create += "\tfraction FLOAT,\n"
create += "\t\"ordreTram\" int8,\n"
create += "\t\"cutEdge\" geometry);"
try:
cur.execute(create)
conn.commit()
except:
print ("CREATE TABLE SegmentsFinals ERROR")
Después se realiza una consulta que determinará qué segmentos son inicio y final.
select = 'SELECT routeid, node, edge FROM "Resultat" ORDER BY routeid, path_seq;'
try:
cur.execute(select)
vec = cur.fetchall()
conn.commit()
except:
print ("SELECT Resultat ERROR")
insert = ''
for x in range (len(vec)):
if vec[x][1] < 0:
if vec[x][1] != -1:
insert +='INSERT INTO "SegmentsFinals" (routeid, edge, "edgeAnt", "ordreTram") '
insert +='VALUES (' + str(vec[x][0]) + ', ' + str(vec[x-1][2]) + ', '
insert +=str(vec[x-2][2]) + ', ' + str(2) +');\n'
else:
insert +='INSERT INTO "SegmentsFinals" (routeid, edge, "edgeAnt", "ordreTram") '
insert +=VALUES (' + str(vec[x][0]) + ', ' + str(vec[x][2]) + ', '
insert +=str(vec[x+1][2]) + ', ' + str(1) + ');\n'
try:
cur.execute(insert)
conn.commit()
except:
print ("INSERT TABLE SegmentsFinals ERROR")
Se realiza un UPDATE para poder añadir la fracción de segmento en la cual se encuentra el punto.
select = 'SELECT routeid, edge, "ordreTram" FROM "SegmentsFinals" ORDER BY routeid, "ordreTram";'
try:
cur.execute(select)
vec = cur.fetchall()
conn.commit()
except:
print ("SELECT SegmentsFinals ERROR")
update = ''
for x in range(len(vec)):
ruta = vec[x][0]
edge = vec[x][1]
ordre = vec[x][2]
if ordre == 1:
update +='UPDATE "SegmentsFinals" s SET fraction = n.fraction FROM NecessaryPoints_'
update +=Fitxer+' n where n.edge_id = '+str(edge)+' AND s.edge ='+str(edge)
update +=' AND s."ordreTram" = 1 AND s.routeid = '+str(ruta)+' AND n.entitatid = 0;\n'
else:
update +='UPDATE "SegmentsFinals" s SET fraction = n.fraction FROM NecessaryPoints_'
update +=Fitxer+' n where n.edge_id = '+str(edge)+' AND s.edge ='+str(edge)
update +=' AND s."ordreTram" = 2 and s.routeid = '+str(ruta)
update +=' AND n.pid = '+str(ruta+1)+';\n'
try:
cur.execute(update)
conn.commit()
except:
print ("UPDATE TABLE SegmentsFinals ERROR")
A continuación se realiza una consulta para escoger y añadir el trozo de tramo que corresponde a cada inicio y final. Posteriormente se hace un UPDATE del campo de geometría de la tabla “SegmentsFinals” con los tramos ya recortados.
select = 'SELECT * FROM "SegmentsFinals" ORDER BY routeid;'
try:
cur.execute(select)
vec = cur.fetchall()
conn.commit()
except:
print ("SELECT SegmentsFinals ERROR")
updateSegment = ''
for x in range(len(vec)):
ordre = vec[x][4]
fraction = vec[x][3]
edgeAnt = vec[x][2]
edge = vec[x][1]
selectTouch ='SELECT ST_Touches((SELECT ST_Line_Substring("Xarxa_Prova"."the_geom",0,'
selectTouch +=str(fraction)+') AS geom FROM "Xarxa_Prova" WHERE"id"='+str(edge)+'),'
selectTouch +='(SELECT the_geom as geom FROM "Xarxa_Prova" WHERE "id"='+str(edgeAnt)+'));'
try:
cur.execute(selectTouch)
resposta = cur.fetchall()
conn.commit()
except:
print ("SELECT TOUCH ERROR")
if edgeAnt != -1:
if resposta[0][0]:
updateSegment +='UPDATE"SegmentsFinals" sf SET "cutEdge" = '
updateSegment +='ST_Line_Substring(s."the_geom",0,'+str(fraction)+') '
updateSegment +='FROM "Xarxa_Prova" s '
updateSegmnet +='WHERE sf."edge"='+str(edge)+' AND s."id"='+str(edge)
updateSegment +=' AND sf."routeid" = '+str(vec[x][0])+';\n'
else:
updateSegment +='UPDATE "SegmentsFinals" sf SET "cutEdge" = '
updateSegment +='ST_Line_Substring(s."the_geom",'+str(fraction)+',1) '
updateSegment +='FROM "Xarxa_Prova" s '
updateSegment +='WHERE sf."edge"='+str(edge)+' and s."id"='+str(edge)
updateSegment +=' and sf."routeid" = '+str(vec[x][0])+';\n'
else:
if ordre == 1:
fractForward = vec[x+1][3]
else:
fractForward = vec[x-1][3]
if fraction >= fractForward:
updateSegment +='UPDATE "SegmentsFinals" sf SET "cutEdge" = '
updateSegment +='ST_Line_Substring(s."the_geom",'+str(fractForward)+','
updateSegment +=str(fraction)+') FROM "Xarxa_Prova" s '
updateSegment +='WHERE sf."ordreTram" = '+ str(ordre)+' and sf."edge"='
updateSegment +=str(edge)+' and s."id"='+str(edge)+' and sf."routeid" = '
updateSegment +=str(vec[x][0])+';\n'
else:
updateSegment +='UPDATE "SegmentsFinals" sf SET "cutEdge" = '
updateSegment +='ST_Line_Substring(s."the_geom",'+str(fraction)+','
updateSegment +=str(fractForward)+') FROM "Xarxa_Prova" s '
updateSegment +='WHERE sf."ordreTram" = '+ str(ordre)+' and sf."edge"='
updateSegment +=str(edge)+' and s."id"='+str(edge)+' and sf."routeid" = '
updateSegment +=str(vec[x][0])+';\n'
try:
cur.execute(updateSegment)
conn.commit()
except:
print ("UPDATE TABLE SegmentsFinals Geometries ERROR")
Se añade y se actualiza el campo de geometría en la tabla “Resultat”.
alter = 'ALTER TABLE "Resultat" ADD COLUMN newEdge geometry;\n'
alter += 'UPDATE"Resultat" r SET newedge = s.the_geom FROM "Xarxa_Prova" s WHERE s.id = r.edge;'
try:
cur.execute(alter)
conn.commit()
except:
print ("ALTER and UPDATE TABLE Resultat Geometries ERROR")
Acto seguido se actualizan los tramos recortados en la tabla “Resultat”.
update = 'UPDATE "Resultat" r SET newedge = s."cutEdge" FROM "SegmentsFinals" s '
update += 'WHERE s."routeid" = r.routeid AND s.edge = r.edge;'
try:
cur.execute(update)
conn.commit()
except:
print ("ALTER and UPDATE TABLE Resultat Geometries ERROR")
Se seleccionan los N caminos más próximos a la dirección indicada en función del límite introducido por el usuario.
limit = self.getLimit()
select ='SELECT e."'+ nomCamp[0][0] +'" AS NomEntitat, r.agg_cost AS Cost, r.entitatID '
select +='FROM "Resultat" r JOIN "' + self.dlg.comboCapaDesti.currentText()
select += '" e ON r.entitatID = e.id WHERE r.edge = -1 ORDER BY 2 ASC limit ' + str(limit) + ';'
try:
cur.execute(select)
vec = cur.fetchall()
conn.commit()
except:
print ("SELECT resultats ERROR")
Finalmente se borrará y se creará una tabla para obtener todos los tramos para cada camino óptimo escogido y, al mismo tiempo, se añade la información obtenida en el SELECT anterior.
createTrams = 'DROP TABLE IF EXISTS "TramsNous_'+Fitxer+'";\n'
createTrams += 'CREATE TABLE "TramsNous_'+Fitxer+'" AS SELECT * FROM (\n'
rowCount = self.dlg.taulaResultat.rowCount()
self.dlg.taulaResultat.setColumnCount(2)
self.dlg.taulaResultat.setHorizontalHeaderLabels(['Entitat', lbl_Cost])
if self.dlg.comboCost.currentText() == 'Distancia':
rnd = 0
else:
rnd = 1
for x in range (rowCount,len(vec)):
self.dlg.taulaResultat.insertRow(x)
self.dlg.taulaResultat.setItem(x, 0, QTableWidgetItem(str(vec[x][0])))
self.dlg.taulaResultat.setItem(x, 1, QTableWidgetItem(str(round(vec[x][1],rnd))))
if x < len(vec) and x >= 1:
createTrams += 'UNION\n'
createTrams +='SELECT entitatid, \'' + str(vec[x][0].replace("'","''"))
createTrams +='\' AS "NomEntitatDesti" ,'+str(round(vec[x][1]))
createTrams +=' AS agg_cost, ST_Union(newedge) AS the_geom from "Resultat" '
createTrams +='WHERE entitatid = '+str(vec[x][2])+' GROUP BY entitatid\n'
createTrams += ")total ORDER BY agg_cost ASC;"
QApplication.processEvents()
try:
cur.execute(createTrams)
conn.commit()
except:
print ("create trams ERROR")
A partir de este punto, lo único que quedará por hacer es presentar en pantalla el resultado.
Conclusión
Como conclusión final del post, cabe destacar que este método permite mejorar y refinar el resultado que se obtiene de una de las funciones de la librería pgrouting. De este modo, se ha podido conseguir que al calcular la distancia más corta de uno a varios puntos se obtenga el punto de segmento concreto en el cual se encuentra el inicio o el final, en vez de quedarse simplemente con el nodo más cercano a dichos puntos.
