Entrades classificades amb: distància

Implementación de las funcionalidades QGIS3 para realizar el cálculo local en el módulo CTE

En este post se mostrará el proceso de implementación de las nuevas funcionalidades que trae consigo la versión 3 de QGIS, las cuales permiten realizar operaciones complejas ejecutando simplemente un algoritmo. En este caso, se tomará como ejemplo el módulo CTE. Este módulo permite obtener los trayectos o caminos entre una dirección de Mataró y las N entidades más próximas.

Utilización manual de los algoritmos

En primer lugar, se debe habilitar el paquete de “Processing” en el caso de que no esté activado. Para ello se deberá ir a “Complementos” -> “Administrar e instalar complementos…” (Fig. 1).

Fig 1. Habilitar paquete “Processing”.

Tras haberlo habilitado aparecerá una nueva pestaña llamada “Procesos”, en la cual se deberá seleccionar “Caja de herramientas” para que se muestren en pantalla todos los algoritmos disponibles. Dicho elemento aparecerá en una barra lateral dentro de QGIS.

A continuación se debe probar si el algoritmo que se utilizará genera el resultado deseado. En el caso del CTE, el algoritmo que potencialmente puede interesarnos es el denominado “Ruta más corta (punto a capa)”, el cual se encuentra en la sección de análisis de redes (Fig. 2).

Fig 2. Caja de herramientas de procesos.

Dicho algoritmo necesita la capa vectorial que representa la red (en nuestro caso SegmentsXarxaCarrers), un punto de inicio y una capa vectorial con los puntos finales. Para realizar las pruebas se utilizarán (de forma totalmente arbitraria) un punto cualquiera del mapa de Mataró como punto de inicio y los IES como puntos finales. Por tanto, se deberán cargar en QGIS tanto SegmentsXarxaCarrers como los IES (Fig. 3).

Fig 3. Importación de los IES y SegmentsXarxaCarrers.

Abrimos el algoritmo, seleccionamos nuestra red, punto de inicio, puntos finales y el resto de los valores los dejaremos por defecto (Fig. 4). Acto seguido, ejecutamos.

Fig 4. Algoritmo “Ruta más corta (punto a capa)”.

Una vez haya finalizado la ejecución el algoritmo, se añadirá automáticamente el resultado a la leyenda (Fig. 5). Desde el punto arbitrario que se ha escogido en este caso particular el resultado sería el siguiente, en el cual se puede apreciar que produce el resultado deseado (buscar la ruta más corta a cada entidad).

Fig 5. Resultado del algoritmo (rutas más cortas).

Inclusión de los algoritmos en el código

Para poder utilizar los algoritmos proporcionados por QGIS e incluirlos en el módulo, será necesario realizar un “import processing” en el código para poder ejecutar la función “processing.run()”, la cual se encargará de ejecutar el algoritmo deseado. Dentro de esta función se le deberá indicar el ID del algoritmo que se vaya a utilizar y los parámetros que se le pasarán.

Para saber cuál es el ID del algoritmo, simplemente con poner el cursor sobre él nos aparecerá su nombre (Fig. 6).

Fig 6. ID del algoritmo.

Por otra parte, para saber qué tipos admite cada uno de los parámetros (y más información sobre el funcionamiento general del algoritmo) podemos usar la consola de QGIS e introducir el comando “processing.algorithmHelp(“native:shortestpathpointtolayer”)”.

En este módulo se ha querido preservar el cálculo que se realiza en el servidor mediante consultas SQL. Para poder mantener ambas opciones se ha optado simplemente por añadir un checkbox el cual propiciará que el módulo se ejecute en local si se ha marcado y, contrariamente, se ejecutará en el servidor cuando esté desmarcado (Fig. 7).

Fig 7. Checkbox de cálculo local en el módulo CTE.

Por tanto, en el código del módulo se ha escrito la siguiente función, la cual se encarga de ejecutar el algoritmo.

def calculo_Local(self,network_lyr,start_point,end_lyr):
    parameters={'INPUT':network_lyr,
                'STRATEGY':0,
                'DIRECTION_FIELD:'',
                'VALUE_FORWARD:'',
                'VALUE_BACKWARD:'',
                'VALUE_BOTH:'',
                'DEFAULT_DIRECTION':2,
                'SPEED_FIELD':'',
                'DEFAULT_SPEED':1,
                'TOLERANCE':0,
                'START_POINT':start_point,
                'END_POINTS':end_layer,
                'OUTPUT':'memory:'}
    return.processing.run('native:shortestpathpointtolayer',parameters)

A esta función se le pasan los siguientes parámetros:

  • network_lyr: es la capa vectorial que representa la red.
  • start_point: es el punto de inicio elegido por el usuario, es decir, la dirección que haya introducido.
  • end_lyr: es la capa vectorial que representa los puntos finales, es decir, las entidades que haya elegido el usuario.

El resto de parámetros son los que anteriormente se usaban por defecto y, teniendo en cuenta el uso que tiene este algoritmo en el módulo, nos interesa que sean esos valores fijos.

En el módulo CTE lo que se busca es mostrarle al usuario los N caminos más cortos, siendo N un número definido por el propio usuario. El algoritmo genera las rutas más cortas de todas las entidades por lo que, tras haber utilizado el algoritmo, se deberá realizar un filtrado de los caminos que tengan un menor coste según el número de entidades establecido por el usuario.

A partir de este punto, se representará el resultado utilizando el mismo código empleado cuando se realiza el cálculo en el servidor.

Cálculo por distancia/tiempo

El algoritmo utilizado, “Ruta más corta (punto a capa)”, puede funcionar definiendo el coste por distancia o por tiempo. Para realizar el cálculo utilizando la distancia lo que hará es comparar la longitud de los segmentos, mientras que por tiempo tomará el valor que le indiquemos de un campo en concreto o bien asignará una velocidad predeterminada.

En el caso del módulo CTE no representa ningún problema realizar el cálculo por distancia. Sin embargo, si se quisiera hacer el cálculo por tiempo, en la base de datos está definido un coste distinto dependiendo del sentido debido a que no se tarda el mismo tiempo en subir una pendiente que en bajarla. Desafortunadamente, este algoritmo está pensado para asignar una única velocidad para cada segmento sin tener en cuenta el sentido. Por tanto, todo lo que se ha descrito previamente sólo sirve para el caso de la distancia ya que en la actualidad este algoritmo no soporta esta dualidad de velocidades.

Diferencias rendimiento (ventajas/inconvenientes)

Las diferencias de rendimiento entre el cálculo en el servidor y la nueva opción de cálculo local varían en función de lo concurrido que esté el servidor en ese momento y del ordenador del cliente.

A continuación se muestra el tiempo que ha tardado el módulo en calcular cada uno de los destinos en servidor y local en unas pruebas de testeo.

Fig 8. Resultados del testeo. Eje X: puntos finales. Eje Y: tiempo de ejecución.

Tal y como se puede apreciar en el gráfico (Fig. 8), en el caso del cálculo en el servidor el tiempo de proceso aumenta de forma proporcional a la cantidad de entidades de destino que tiene que calcular. Por otra parte, el tiempo que tarda el cálculo local es sumamente sólido, siendo prácticamente una constante. Se tarda el mismo tiempo en calcular 4 destinos (CasalsAvisOficial) que 144 (ParadesBus).

Debido a que el cálculo local presenta una duración más estable, que de media tarda 5 segundos (mientras que el cálculo en el servidor tarda de media 7,9 segundos) y que libera carga en el servidor, esta nueva forma de realizar los cálculos representa una mejora en el módulo CTE y, por norma general, al usuario le interesará ejecutar este módulo marcando el cálculo local.

 

Cálculo en el servidor de los caminos más cortos en el módulo CTE

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.

Manual d’ús de l’usuari: Cerca de trajectes a entitats

Aquest document explica el funcionament del mòdul de ‘Cerca de Trajectes a Entitats’ per a QGIS. Aquest mòdul permet obtenir els trajectes o camins entre una adreça qualsevol de Mataró entrada per l’usuari i les N entitats més properes (també escollides per l’usuari). El cost pot estar en funció de la distància o el temps. Els desplaçaments obtinguts es mostren de forma ordenada de menys a més cost en forma de taula i de capa. Per a poder utilitzar-lo, el primer que s’ha de fer és executar el programa QGIS i un cop inicialitzat aquest, cal pitjar la icona següent CTE2o anar a Complementos -> CCU -> Cerca de Trajectes a Entitats i s’obrirà una finestra com la que podem veure a continuació a la figura 1.

CTE1

Figura 1: Imatge de l’estat inicial del mòdul CTE

Aquest mòdul es pot dividir bàsicament en tres parts diferenciades: origen (1) , destí (2), connexió i botons (3).

1. Connexió i altres botons

La primera part a mostrar inclou els botons de funcionament del mòdul (d’esquerra a dreta):

  • Connexió: pestanya desplegable on es mostren totes les connexions.
  • Botó ‘Inici’: aquest botó inicia el procés de consulta del mòdul.
  • Botó ‘Sortir’: aquest botó tanca el mòdul.

CTE5

Figura 2: Connexió i botons

2. Origen

En aquesta part del mòdul l’usuari troba les eines per escollir el punt d’inici, més concretament l’adreça. Seguidament es descriuran les seves parts, començant per la part superior:

CTE3

Figura 3: Origen

El primer element que trobem és una camp de text per cerca el nom del carrer a on es troba l’origen. Inicialment, en el gran requadre inferior, s’hi llisten tots els noms de carrers i a mesura que es va introduint un nom, la llista es va actualitzant automàticament amb els noms que s’assemblen al nom que introdueix l’usuari fins a trobar el desitjat. Una vegada, s’ha trobat el carrer, cal fer doble clic en el nom perquè quedi ben especificat en el camp de text.

Posteriorment, s’ha d’escriure el número del portal que desitgem, juntament amb la lletra. En el cas de la lletra no és necessari escollir-ne una.

Per acabar, a la dreta del camp de text per introduir el nom hi ha un botó amb una creu (X) que permet esborrar el contingut del camp de text.

3. Destí

En aquesta part del mòdul, l’usuari pot escollir la capa de destí i les diferents característiques dels camins: cost, número de camins i aproximació en cas de no existir l’adreça desitjada.

Començant per la part superior, el primer que trobem és una pestanya desplegable on es mostren les diferents capes de punts que poden ser destí. Just a sota hi trobem els modes de cerca. El primer que trobem són dos botons, de color verd i vermell respectivament.

El botó “Més proper” escull l’adreça amb el número de portal més proper a la adreça que s’ha seleccionat, independentment de si és parell o senar, i sempre que l’adreça escollida per l’usuari no existeixi.

En canvi, el botó “Més proper Parell/Senar” escull el portal amb el número senar o parell en funció de quin dels dos tipus sigui el número que s’hagi escollit prèviament, i sobretot que no existeixi.

A la dreta d’aquests dos botons, hi ha una spin box que permet escollir el nombre de camins desitjats, en un rang de 1 a 10 camins. En el cas que el nombre de camins sigui superior al nombre d’entitats d’una capa, el número de camins serà tants com entitats tingui la capa escollida.

Després trobem una pestanya desplegable on es pot triar el cost dels camins: distància o temps. En el cas d’escollir la opció de temps, hi ha un check box que permet incloure en el cost el cost de nusos.

A sota hi ha un gran requadre que s’utilitza per mostrar els resultats de la cerca que s’ha realitzat. La taula que es mostra té dos columnes: la primera mostra el nom de l’entitat de destí i la segona mostra el cost: si és distància en metres i si és en temps en minuts.

Entre la opció del cost i el requadre de resultats hi ha un espai buit que un cop s’ha fet la consulta hi apareix l’adreça de l’origen.

CTE4

Figura 4: Secció de destí

Primers passos en el ‘rutatge’ al QGIS

En aquest post mostraré els passos que he hagut de seguir per tal de crear el primer dels mòduls que utilitza funcions de rutatge en el QGIS. En aquest cas, volem traslladar el mòdul que calcula els 3 camins a les escoles bressol i llars d’infants més pròximes a cada portal de la ciutat de Mataró. En primer lloc presento les capes amb les quals es treballa. Llavors segueixo amb la creació del mòdul i la seva utilització.

Capes de treball

Actualment amb el GeoMedia treballem amb bases de dades de Microsoft Access majoritàriament. Aquestes incorporen tot tipus de dades, ja sigui taules amb informació sobre el padró o el cadastre, que no estan georeferenciades, o capes de punts, línies o polígons com poden ser els mapes de carrers, illes, cruïlles i serveis públics de la ciutat. Les dades en aquest format ens permeten una gran versatilitat d’ús i facilitat tan en l’ús com en la creació, modificació i eliminació. També són fàcilment exportables a qualsevol altre tipus de format per seguir treballant. Ocasionalment també treballem amb arxius SHAPE.

El QGIS permet gestionar formats raster i vectorials a través de les biblioteques GDAL i OGR, així com altres bases de dades. Una biblioteca GDAL i OGR és un conjunt de programes que estan formats per comandes, cada un amb moltes possibilitats d’ús.

Les capes més comunes amb les quals hem treballat són el Vector Layer, on utilitzem arxius SHAPE (.shp), i Delimited Text Layer on utilitzem arxius CSV (.csv).

Ambdos tipus de dades poden ser modificats fàcilment: tan crear, inserir, modificar o eliminar objectes i/o camps, la qual cosa els fa idonis treballar amb aquests tipus d’arxius ja que si volem fer proves per poder desenvolupar les consultes o els mòduls, ens són de gran utilitat.

Procés de creació

Per crear la interfície gràfica s’utilitza un creador de models que porta incorporats una sèrie de funcionalitats. Haurem d’anar a la caixa d’eines de processat i allà “crear model nou”.
fuc
I seguidament se’ns obra la següent finestra:
nou
Les dues imatges següents corresponen a la columna de l’esquerra de l’anterior imatge hi tenim tots els tipus de paràmetres amb els quals podem treballar, i si canviem a la pestanya de algoritmes, trobem totes les funcions amb les quals podrem treballar sobre les dades. Cal veure el funcionament d’aquestes per tal de poder-les utilitzar correctament.

fuc

dades

 

 

 

 

 

Seguidament, vam afegir el tipus de dades i buscar les funcionalitats adequades per construir la interfície gràfica.

dades2
Una vegada introduïdes les capes, va ser necessari buscar una funcionalitat que convertís el tipus de geometria dels objectes de les capes de punts ja que van sorgir problemes per què obteníem una capa buida. En el nostre cas vam utilitzar el mòdul “Convert Geometry type”.
convert
Seguidament vam començar la recerca d’entre les funcionalitats alguna que calculi una ruta entre dos punts. No existeix una funcionalitat que faci tal funció, així que vam seguir la cerca. Vam seguir buscant entre els plugins que permet instal·lar el QGIS. Vam trobar-ne un que si que fa aquesta funció però no ens era d’utilitat ja que no hi ha cap manera d’automatitzar el procés de càlcul de les distàncies. Això doncs, vam recórrer a l’últim recurs: crear un script en Python que realitzi la funció.

Vam trobar un script que calculava la ruta entre dos punts sobre un graf de carrers. Aquest script utilitza una funció interna de la API de QGIS anomenada “dijkstra” que et retorna un camí. El següent script és el que vam trobar a Internet:

##ruta=name
##ruta=name
##points=vector
##network=vector
##output=output vector

#Algorithm body
#==================================
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *
from processing.tools.vector import VectorWriter

point_layer = processing.getObject(points)
network_layer = processing.getObject(network)
writer = VectorWriter(output, None, [QgsField("order", QVariant.Int)],
network_layer.dataProvider().geometryType(), network_layer.crs())

# prepare graph
vl = network_layer
director = QgsLineVectorLayerDirector(vl,-1,'','','',3)
properter = QgsDistanceArcProperter()
director.addProperter( properter )
crs = vl.crs()
builder = QgsGraphBuilder( crs )

# prepare points
features = processing.features(point_layer)
point_count = point_layer.featureCount()
points = []
for f in features:
  points.append(f.geometry().asPoint())
tiedPoints = director.makeGraph( builder, points )
graph = builder.graph()
route_vertices = []

for i in range(0,point_count-1):
    progress.setPercentage(int(100 * i/ point_count))
    
    from_point = tiedPoints[i]
    to_point = tiedPoints[i+1]
    from_id = graph.findVertex(from_point)
    to_id = graph.findVertex(to_point)

    (tree,cost) = QgsGraphAnalyzer.dijkstra(graph,from_id,0)
    if tree[to_id] == -1:
        continue # ignore this point pair
    else:
        #collect all the vertices between the points
        route_points = []
        curPos = to_id 
        while (curPos != from_id):
            route_points.append(graph.vertex(
graph.arc(tree[curPos]).inVertex()).point())
           curPos = graph.arc( tree[ curPos ] ).outVertex()
        route_points.append(from_point)
    # add a feature
    fet = QgsFeature()
    fet.setGeometry(QgsGeometry.fromPolyline(route_points))
    fet.setAttributes([i])
    writer.addFeature(fet)
del writer

Un cop testejat i debuguejat, vam començar amb la modificació de l’script per tal d’obtenir el resultat desitjat. Aquest va ser un procés complex ja que s’havia de canviar el programa per dins. Després de forces entrebancs en el procés vam aconseguir un resultat força aproximat al que volíem. Van caldre forces hores per acabar de depurar el codi ja que teníem petits errors que costaven de detectar.

##dintreilla=vector
##EscolesBressol=vector
##network=vector
##output=output vector
#Algorithm body
#==================================
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import time
from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *
from processing.tools.vector import VectorWriter

start_time = time.time()
network_layer = processing.getObject(network)

inputPoint = processing.getObject(dintreilla)
features = processing.features(inputPoint)

inputPoint2 = processing.getObject(EscolesBressol)
features2 = processing.features(inputPoint2)

di= 0
eb= 0
id = -1
fields = []
fields.append (QgsField("ID", QVariant.Int))
fields.append (QgsField("Length", QVariant.Int))
fields.append(QgsField("From_Node",QVariant.String))
fields.append(QgsField("To_Node",QVariant.String))

writer = VectorWriter(output, None, 
fields, network_layer.dataProvider().geometryType(), network_layer.crs())

#Per buscar els 3 millors de cada punt

for fea1 in features:
    di=di+1
    #xx = fea1.geometry().asPoint().x()
    #yy = fea1.geometry().asPoint().y()
    #pStart = QgsPoint(xx, yy)
    from_node = fea1.attributes()
    inici = from_node[0]
    print inici
    fea2=None
    features2 = processing.features(inputPoint2)
    eb = -1
    vec = []
    for fea2 in features2:
        eb=eb+1
        id  = id + 1
        nom = fea2.attributes()
        desti = nom[2]
       #---------------------------------------------------------------
        vl = network_layer
        director = QgsLineVectorLayerDirector(vl,-1,'Cost','Cost_inver','',3)
        properter = QgsDistanceArcProperter()
        director.addProperter(properter)
        crs = vl.crs()
        builder = QgsGraphBuilder(crs ,True,0.001)

        # prepare points
        points = []
        points.append(fea1.geometry().asPoint())
        points.append(fea2.geometry().asPoint())

        tiedPoints = director.makeGraph( builder, points )
        graph = builder.graph()

        route_vertices = []
        for i in range(0,2-1):
            from_point = tiedPoints[i]
            to_point = tiedPoints[i+1]

            from_id = graph.findVertex(from_point)
            to_id = graph.findVertex(to_point)

            (tree,cost) = QgsGraphAnalyzer.dijkstra(graph,from_id,0)
            if tree[to_id] == -1:
                continue # ignore this point pair
            else:
                # collect all the vertices between the points
                route_points = []
                curPos = to_id 
                while (curPos != from_id):
                    route_points.append(graph.vertex(
graph.arc(tree[curPos]).inVertex()).point())
                    curPos = graph.arc(tree[curPos]).outVertex()

                route_points.append(from_point)

            # add a feature
            geom=QgsGeometry.fromPolyline(route_points)
            fet = QgsFeature()
            fet.setGeometry(QgsGeometry.fromPolyline(route_points))
            fet.setAttributes([id, geom.length(), inici, desti])
            vec.append(fet)
           
    if (len(vec) > 0):
        vec.sort(key=lambda vec: vec[1])
        for i in range (0,3):
            writer.addFeature(vec[i])
del writer
print("--- %s seconds ---" % (time.time() - start_time))

Finalment vam aconseguir un resultat que s’ajustava a les nostres necessitats, i així completar el procés de creació d’un mòdul amb el QGIS. Vam posar l’script en un mòdul per a python, vam posar-li totes les connexions necessàries.
python
I aquest en va ser el resultat final de tot el procés:
final

Procés d’utilització

Aquí es descriu el procés d’utilització del mòdul per trobar els 3 Camins més pròxims a una Escola Bressol i una d’infants. Aquest comença amb la preparació de les dades a la llegenda o panell de capes del QGIS. En el cas que ens ocupa necessitarem 4 capes SHAPE, 3 de punts i una de segments. Les 3 capes de punts són les Escoles Bressol, les Llars d’Infants  i els dintreilles, que són tots els portals de la ciutat.  I la capa de segments són el conjunts de carrers de la ciutat, és a dir, per on hem de trobar el camí.

Primer de tot, hem de tenir les capes en el tipus desitjat: SHAPE. En el cas que no estiguin ja en aquest format, cal transformar-les per tal de poder-hi treballar. Un cop les tinguem, les guardem per tal de poder-les agafar i emprar.

Tal com s’indica a la fotografia, s’afegeix cada capa via Capa -> Añadir capa -> Añadir capa vectorial:
afe3

afe

Seleccionem explorar i amb l’ajuda de la finestra, busquem els arxius SHAPE que volem posar.

afe2

Repetim l’acció 3 vegades més fins a aconseguir les 4 capes desitjades, tal i com es veu a la foto.

Un cop fet, anem al panell de la dreta de la pantalla on hi ha la “caja de herramientas de procesado” i busquem a l’apartat de Modelos -> CCU, un model anomenat “3EB més pròximes”.
afe4

Una vegada trobat, cal executar-lo. S’obrirà una pestanya amb el següent diàleg:

afe5

Posem a cada pestanya la capa que ens demani: a la primera hi posem la xarxa de carrers sobre la qual volem treballar, a la segona hi posem els Dintreilles o portals de la ciutat(assegurar-se de que sigui la versió “dintreilla_trajectes”) i finalment les escoles bressol o les llars d’infants. Haurem de repetir el procés per cada una de les capes: una per les EB i una altre per les LI.

Un cop tot estigui a punt, només cal executar el procés i esperar a obtenir el resultat. El temps d’espera pot variar segons l’ordinador on s’estigui executant aquest.

Quan finalitzi el procés cal guardar el resultat en un fitxer SHAPE, ja que està en un fitxer temporal i per tant, es perdria en el moment en què tanquem el programa.

Prova del mòdul

En aquest apartat es mostra un exemple de mostra del mòdul. No utilitzaré la capa del ‘dintreilla’ ja que es massa gran i per veure el funcionament, amb una simple capa en tenim prou.

En aquest cas utilitzaré una capa amb un punt (de color vermell) que serà el punt d’origen i una altre capa (de color verd) amb els possibles destins més propers. D’aquesta manera es pot veure ben clar el funcionament.
ex1

Posem en marxa el mòdul i el resultat que obtenim és el següent:

ex2

A la imatge es pot veure el graf de carrers de la ciutat de Mataró. S’hi pot veure 3 camins ressaltats de color verd clar. Tots tres tenen com a origen el punt vermell.

WFS: Interactuació amb mapes.

Descripció.

La necessitat de transferir cartografia per Internet ha sigut, i encara és, un problema difícil. La cartografia digital conté un gran volum de dades d’informació, i enviar aquestes quantitats de megabytes per Internet és una tasca pesada, i lenta.

És per això que es va començar a treballar amb WFS, sigles de “Web Feature Service”, que és un servei estàndard que ofereix una interfície que permet sol·licituds de comunicació, permetent interactuar amb mapes WMS (Web Map Service).

A través d’una URL, es pot accedir a les dades cartogràfiques que s’hagin publicat a la taula WFS, i fer consultes específiques, més endavant es mostra un exemple.

Per realitzar aquestes operacions s’utilitza el llenguatge GML, que deriva de XML,  és l’estàndard a través del qual es transmeten les comandes WFS.

Un cop s’ha fet la publicació de la taula WFS i es té la URL que apunta a ella, es poden veure les dades en alguns SIG (Sistemes d’Informació Geogràfica) que permetin el tractament d’aquest tipus de dades, en aquest cas al Geomedia.

En aquest cas, s’ha decidit a utilitzar la taula de trajectes i proximitat a les tres escoles bressol més properes.

Com fer la publicació?

Primerament, des del servidor on està instal·lat el GeoMediaWebMap Professional, s’obre el Server Configuration Utility i es crea un nou servei amb el botó “Add”.

Creació del servei.

El següent pas és seleccionar l’opció Manipulate Feature Web Service.

Es deixen els valors per defecte i es va prement “Next”, s’haurà d’introduir un nom pel WFS.

Nom del servei.

Quan es demani el tipus de base de dades, introduir Microsoft Access. La resta es deixa per defecte, i es finalitza, el servei s’ha creat.

En aquest punt, s’obre el Geomedia i es carrega el mapa que es vol publicar, s’obre el “GeomediaWebMap Publisher Administrator”.

Configurant el Geomedia.

Es selecciona el servei que s’acaba de crear i apareixerà una barra amb botons. Polsar sobre el cinqué botó “PublishandPopulatetheGeoWorkspace”, i seleccionar la primera opció “PublishtheGeoWorkspace contents to theMetaData”. Es selecciona “Sí”.

Botó “settings”.

Es selecciona el segon botó, “settings”. Aquí s’ha de seleccionar l’arxiu .csf que conté el sistema de coordenades. Després es prem l’últim botó per tancar el menú.

En aquest punt només queda inicialitzar el servei, per fer-ho, es torna al Publisher Server Configuration Utility, es selecciona el servei que s’ha creat i es prem “Initialize”.

Servei inicialitzat.

Provar el servei des del navegador introduint l següent línia de comanda:

http://geoportalccu.tecnocampus.cat/ProvaWFS/request.aspx?version=1.1.0&service=wfs&request=getcapabilities

Codi GML creat.

Interpretació gràfica del codi GML.

Veure així el codi no és interessant, el més interessant és veure-ho transformat en un mapa, al Geomedia per exemple, a continuació es veu com es pot configurar el geomedia per tal de veure la informació del servei WFS.

Primerament s’ha de crear una nova connexió, de tipus WFS, en aquest cas s’ha creat de WFS només lectura, per tal de que les dades no es puguin modificar.

Nova connexió.

S’ha d’introduir la URL abans esmentada, i s’accepta.

Un cop creada la connexió, ja es poden mostrar les seves entrades des de la llegenda.

Ventana de mapa del Geomedia amb les dades.

Polígons de Voronoi

En aquest article es descriurà el funcionament dels polígons de Voronoi, una eina molt important a l’hora d’estudiar àrees d’influència.

Els polígons de Voronoi es basen en la distància euclidiana, i són molt apropiats quan les dades són qualitatives. Es tracta de fer una partició del pla, a partir d’uns punts que anomenarem punts generadors.

Aquesta partició del pla en regions té la peculiaritat de que des de qualsevol punt de dins d’una regió determinada, la distància al punt generador corresponent és sempre menor que la distància a qualsevol altre punt generador extern. Per tant, les fronteres de les regions són equidistants de dos o mes punts generadors.

Inicialment, aquest polígons van ser creats per l’anàlisi de dades meteorològiques, però avui en dia s’utilitzen també per determinar zones d’influència, que és el que s’explicarà en aquest article.

Els polígons de Voronoi serveixen per dividir un espai en un número determinat de regions. S’especifiquen un conjunt de punts (punts generadors) i quan es fa el diagrama, aquests queden dividits pels polígons, un punt en cada regió. Les regions s’anomenen cel·les o polígons de Voronoi.

Cada polígon correspon a l’àrea d’influència, per dir-ho d’alguna manera, del punt  que conté.

Primerament, crec que és interessant posar un exemple per entendre millor per a què serveixen els polígons:

Suposem el cas de que s s’està estudiant els centres d’atenció primària (CAP) del terme municipal de Mataró, i es vol construir un altre i no se sap on. Gràcies a les zones d’influència creades mitjançant els polígons de Voronoi, es podrà situar més o menys el nou CAP.

El programa que s’ha utilitzat, a més del Geomedia, ha sigut el Global Mapper.

Global Mapper és una potent aplicació que combina eines de tractament de dades espacials amb un accés a gran varietat de formats d’arxius. És molt útil com a complement del Geomedia.

Amb el Geomedia s’han exportat els caps com una única entitat, i el contorn de Mataró s’ha aconseguit agafant el perfil de la unió de totes les illes que formen Mataró.

Terme municipal amb els CAPS

Tots els arxius exportats són de tipus ShapeFile, per tal de que Global Mapper els reconegui i es pugui treballar amb ells.

Un cop es té el terme municipal de Mataró amb els CAPS, és hora de passar a l’acció, el procés és simple, i el resultat és molt satisfactori.

Primerament es selecciona tot i es crea el diagrama de Voronoi des del menú d’anàlisi.

Apareixerà una finestra en la que s’haurà d’indicar que es volen allargar els límits uns 4000 metres, això permetrà que, en cas de que el diagrama no arribi a tocar el perímetre de Mataró, aquest s’allargui fins a tocar-lo.

Allargar límits

Un cop allargats, només cal dir que el límit fins on s’allarguen els polígons és l’àrea contenidora, és a dir, el terme municipal de Mataró. Això es configura a partir del botó “Bounds” i seleccionant l’última opció. Per últim s’accepta per crear els polígons.

Limitar polígons

La imatge de Mataró amb els seus CAPS quedaria dividida pels polígons de Voronoi, hi hauria una cel·la per cada CAP.

Aquests polígons resultants, permeten veure l’àrea d’influència de cada CAP.

Polígons creats

Observant la imatge, es pot veure les divisions que corresponen a cada CAP, hi ha una concentració més elevada al centre urbà, degut a que la població és notablement més elevada.

El CAP de dalt a la dreta, el de Rocafonda-Palau, té molta zona d’influència, això és perquè la zona Nord de Mataró no està tant urbanitzada com la resta.

Segons la població, i el número de places de cada CAP, fent ús dels diagrames de Voronoi, es podria situar un futur CAP.