Rutas en PostGIS con la nueva versión de Pgrouting (parte 4/4): Nodificando una red

routing1=# select 
    st_length (st_transform (linea, 900913)) as mercator, 
    st_length (st_transform (linea, 32630)) as utmwgs84huso30,                
    st_length (linea::geography) as elipsoidalreal  
    from (
       select st_geomfromtext 
      ('LINESTRING (-0.4013 39.4549,-0.3414 39.4805)', 4326) as linea    
    ) as tabla;
    mercator     | utmwgs84huso30  |  elipsoidalreal  
------------------+-----------------+------------------
 7621.67507015997 | 5887.6082569941 | 5886.25654225072

 

La distancia mercator es 1764 metros mayor que la utmwgs84, por el contrario esta última se aproxima a la real en menos de 1 metro. Como vamos a tener que nodificar la red y calcular las intersecciones necesitamos trabajar en un sistema proyectado para que PostGIS realiza correctamente los cálculos, tenemos dos opciones:

  • Reproyectar las geometrías de la capa a una proyección local adecuada con capacidad métrica como UTM (siempre que los datos OSM queden contenidos dentro de un único huso, en caso contrario habría que seleccionar otro tipo de proyección como una Lambert).
  • Reproyectar las geometrías a una proyección de ámbito mundial como una mercator (EPSG: 3857) con lo cual no hay que preocuparse de la extensión de la cartografía (puede incluso considerar varios continentes) y utilizar el comando ST_Length_spheroid para el cálculo de costes de los tramos.

En nuestro ejemplo vamos a seguir la primera opción al tener la cartografía una extensión que cabe dentro de un único huso UTM (y la latitud no es superior a 80 grados). Vamos a cargar el fichero ejes.sql que contiene la cartografía del núcleo de Valencia pero esta vez no tiene intersecciones y está en 4326 como ya hemos comentado

routing1=# \d ejes                                        
 Columna |           Tipo            |                Modificadores
---------+---------------------------+------------------------------
 gid     | integer                   | not null valor por omisión
                                    nextval('ejes_gid_seq'::regclass)
 name    | text                      |
 oneway  | integer                   |
 geom    | geometry(LineString,4326) |
Índices:
    "ejes_pkey" PRIMARY KEY, btree (gid)

 

Asignamos el CRS 32630 a todas las geometrías de la capa ‘ejes’ y actualizamos la entrada en geometry_colums (proceso válido para cualquier versión de PostGIS):

routing1=# select updategeometrysrid ('ejes','geom', 32630);
--Reproyección de las geometrías al CRS 32630.
routing1=# update ejes 
      set geom = st_transform (st_setsrid (geom, 4326), 32630);

 

Ahora calculamos las intersecciones de la capa ejes con el comando pgr_nodenetwork. Se creará una nueva tabla ejes_noded con los nuevos tramos segmentados. Además se añadirán dos campos source y target vacíos a la nueva tabla.

routing1=# select pgr_nodenetwork ('ejes', 0.01, 'gid', 'geom');
NOTICE:   Total New segments: 46841
NOTICE:   New Table: public.ejes_noded