40. Construcción de Geometrías¶

La capa nyc_subway_stations nos ha proporcionado muchos ejemplos interesantes, pero hay algo que llama la atención:

Aunque es una base de datos de todas las estaciones, ¡no permite visualizar fácilmente las rutas! En este capítulo utilizaremos funciones avanzadas de PostgreSQL y PostGIS para construir una nueva capa de rutas lineales a partir de la capa de puntos de las estaciones de metro.

Dos cuestiones dificultan nuestra tarea:

La columna routes de nyc_subway_stations tiene varios identificadores de ruta en cada fila, por lo que una estación que podría aparecer en varias rutas sólo aparece una vez en la tabla.
En relación con el problema anterior, no hay información sobre el orden de las rutas en la tabla de estaciones, por lo que, aunque es posible encontrar todas las estaciones de una ruta determinada, no es posible utilizar los atributos para determinar el orden en que los trenes pasan por las estaciones.

El segundo problema es el más difícil: dado un conjunto desordenado de puntos de una ruta, ¿cómo los ordenamos para que coincidan con la ruta real?.

Estas son las paradas del tren «Q»:

SELECT s.gid, s.geom
FROM nyc_subway_stations s
WHERE (strpos(s.routes, 'Q') <> 0);

En esta imagen, las paradas están etiquetadas con su clave primaria gid.

Si empezamos en una de las estaciones finales, la siguiente estación de la línea parece ser siempre la más cercana. Podemos repetir el proceso cada vez siempre que excluyamos de nuestra búsqueda todas las estaciones encontradas anteriormente.

Hay dos formas de ejecutar una rutina iterativa de este tipo en una base de datos:

Utilizando un lenguaje procedimental, como PL/PgSQL.
Uso recursivo de expresiones comunes.

Las expresiones comunes de tabla (CTE) tienen la virtud de no requerir una definición de función para ejecutarse. Aquí está la CTE para calcular la línea de ruta del tren “Q”, empezando por la parada más al norte (donde gid es 304).

WITH RECURSIVE next_stop(geom, idlist) AS (
    (SELECT
      geom,
      ARRAY[gid] AS idlist
    FROM nyc_subway_stations
    WHERE gid = 304)
    UNION ALL
    (SELECT
      s.geom,
      array_append(n.idlist, s.gid) AS idlist
    FROM nyc_subway_stations s, next_stop n
    WHERE strpos(s.routes, 'Q') != 0
    AND NOT n.idlist @> ARRAY[s.gid]
    ORDER BY ST_Distance(n.geom, s.geom) ASC
    LIMIT 1)
)
SELECT geom, idlist FROM next_stop;

El CTE consta de dos partes unidas entre sí:

La primera parte establece un punto de inicio para la expresión. Obtenemos la geometría inicial e inicializamos la matriz de identificadores visitados, utilizando el registro de «gid» 304 (el final de la línea).
La segunda parte itera hasta que no encuentra más registros. En cada iteración toma el valor devuelto de la iteración anterior a través de la autorreferencia a «next_stop». Buscamos en cada parada de la línea Q (strpos(s.rutas,”Q”)) que no hayamos añadido ya a nuestra lista de visitadas (NOT n.idlist @> ARRAY[s.gid]) y las ordenamos por su distancia al punto anterior, tomando sólo la primera (la más cercana).

Más allá de la propia CTE recursiva, hay una serie de características avanzadas de PostgreSQL que se utilizan aquí:

¡Estamos usando ARRAY! PostgreSQL soporta arreglos de cualquier tipo. En este caso tenemos un arreglo de enteros, pero también podríamos construir un arreglo de geometrías, o cualquier otro tipo de PostgreSQL.
Estamos utilizando array_append para construir nuestro vector de identificadores de estación visitados.
Estamos utilizando el operador de matrices @> («array contains») para averiguar cuál de las estaciones de tren Q ya hemos visitado. El operador @> requiere valores ARRAY en ambos lados, por lo que tenemos que convertir los números «gid» individuales en arrays con un solo elemento utilizando la sintaxis ARRAY[].

Al ejecutar la consulta, se obtiene cada geometría en el orden en que se encuentra (que es el orden de ruta), así como la lista de identificadores ya visitados. Envolver las geometrías en la función agregada ST_MakeLine de PostGIS convierte el conjunto de geometrías en una única línea, construida en el orden proporcionado.

WITH RECURSIVE next_stop(geom, idlist) AS (
    (SELECT
      geom,
      ARRAY[gid] AS idlist
    FROM nyc_subway_stations
    WHERE gid = 304)
    UNION ALL
    (SELECT
      s.geom,
      array_append(n.idlist, s.gid) AS idlist
    FROM nyc_subway_stations s, next_stop n
    WHERE strpos(s.routes, 'Q') != 0
    AND NOT n.idlist @> ARRAY[s.gid]
    ORDER BY ST_Distance(n.geom, s.geom) ASC
    LIMIT 1)
)
SELECT ST_MakeLine(geom) AS geom FROM next_stop;

La cual luce como:

Éxito!

Excepto por dos problemas:

Aquí sólo estamos calculando una ruta de metro, queremos calcular todas las rutas.
Nuestra consulta requiere un conocimiento a priori, el identificador inicial de la estación que sirve de semilla para el algoritmo de búsqueda que construye la ruta.

Abordemos primero el problema más difícil: averiguar cuál es la primera estación de una ruta sin echar un vistazo manual al conjunto de estaciones que la componen.

Nuestras paradas del tren «Q» pueden servir de punto de partida. ¿Qué caracteriza a las estaciones finales de la ruta?

Una respuesta es «son las estaciones más al norte (septentrionales) y más al sur ( meridionales)». Sin embargo, imagine que el tren «Q» circulara de este a oeste. ¿Se cumpliría la condición?

Una caracterización menos direccional de las estaciones finales es «son las estaciones más alejadas del centro de la ruta». Con esta caracterización no importa si la ruta discurre en dirección norte/sur o este/oeste, sólo que discurra más o menos en una dirección, sobre todo en los extremos.

Puesto que no existe una heurística al 100% para averiguar los puntos finales, probemos esta segunda regla.

Nota

Un fallo obvio de la regla «lo más lejos posible del centro» es una línea circular, como la Circle Line de Londres (Reino Unido). Afortunadamente, Nueva York no tiene líneas de este tipo!

To work out the end stations of every route, we first have to work out what routes there are! We find the distinct routes.

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
)
SELECT * FROM routes;

Note the use of two advanced PostgreSQL ARRAY functions:

string_to_array takes in a string and splits it into an array using a separator character. PostgreSQL supports arrays of any type, so it’s possible to build arrays of strings, as in this case, but also of geometries and geographies as we’ll see later in this example.
unnest takes in an array and builds a new row for each entry in the array. The effect is to take a «horizontal» array embedded in a single row and turn it into a «vertical» array with a row for each value.

The result is a list of all the unique subway route identifiers.

 route
-------
 1
 2
 3
 4
 5
 6
 7
 A
 B
 C
 D
 E
 F
 G
 J
 L
 M
 N
 Q
 R
 S
 V
 W
 Z
(24 rows)

We can build on this result by joining it back to the nyc_subway_stations table to create a new table that has, for each route, a row for every station on that route.

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
)
SELECT * FROM stops;

 gid |                      geom                      | route
-----+----------------------------------------------------+-------
   2 | 010100002026690000CBE327F938CD21415EDBE1572D315141 | 1
   3 | 010100002026690000C676635D10CD2141A0ECDB6975305141 | 1
  20 | 010100002026690000AE59A3F82C132241D835BA14D1435141 | 1
  22 | 0101000020266900003495A303D615224116DA56527D445141 | 1
                            ...etc...

Now we can find the center point by collecting all the stations for each route into a single multi-point, and calculating the centroid of that multi-point.

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
),
centers AS (
  SELECT ST_Centroid(ST_Collect(geom)) AS geom, route
  FROM stops
  GROUP BY route
)
SELECT * FROM centers;

The center point of the collection of “Q” train stops looks like this:

So the northern most stop, the end point, appears to also be the stop furthest from the center. Let’s calculate the furthest point for every route.

WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
),
centers AS (
  SELECT ST_Centroid(ST_Collect(geom)) AS geom, route
  FROM stops
  GROUP BY route
),
stops_distance AS (
  SELECT s.*, ST_Distance(s.geom, c.geom) AS distance
  FROM stops s JOIN centers c
  ON (s.route = c.route)
  ORDER BY route, distance DESC
),
first_stops AS (
  SELECT DISTINCT ON (route) stops_distance.*
  FROM stops_distance
)
SELECT * FROM first_stops;

We’ve added two sub-queries this time:

stops_distance joins the centers points back to the stations table and calculates the distance between the stations and center for each route. The result is ordered such that the records come out in batches for each route, with the furthest station as the first record of the batch.
first_stops filters the stops_distance output by only taking the first record for each distinct group. Because of the way we ordered stops_distance the first record is the furthest record, which means it’s the station we want to use as our starting seed to build each subway route.

Now we know every route, and we know (approximately) what station each route starts from: we’re ready to generate the route lines!

But first, we need to turn our recursive CTE expression into a function we can call with parameters.

CREATE OR REPLACE function walk_subway(integer, text) returns geometry AS
$$
WITH RECURSIVE next_stop(geom, idlist) AS (
    (SELECT
      geom AS geom,
      ARRAY[gid] AS idlist
    FROM nyc_subway_stations
    WHERE gid = $1)
    UNION ALL
    (SELECT
      s.geom AS geom,
      array_append(n.idlist, s.gid) AS idlist
    FROM nyc_subway_stations s, next_stop n
    WHERE strpos(s.routes, $2) != 0
    AND NOT n.idlist @> ARRAY[s.gid]
    ORDER BY ST_Distance(n.geom, s.geom) ASC
    LIMIT 1)
)
SELECT ST_MakeLine(geom) AS geom
FROM next_stop;
$$
language 'sql';

And now we are ready to go!

CREATE TABLE nyc_subway_lines AS
-- Distinct route identifiers!
WITH routes AS (
  SELECT DISTINCT unnest(string_to_array(routes,',')) AS route
  FROM nyc_subway_stations ORDER BY route
),
-- Joined back to stops! Every route has all its stops!
stops AS (
  SELECT s.gid, s.geom, r.route
  FROM routes r
  JOIN nyc_subway_stations s
  ON (strpos(s.routes, r.route) <> 0)
),
-- Collects stops by routes and calculate centroid!
centers AS (
  SELECT ST_Centroid(ST_Collect(geom)) AS geom, route
  FROM stops
  GROUP BY route
),
-- Calculate stop/center distance for each stop in each route.
stops_distance AS (
  SELECT s.*, ST_Distance(s.geom, c.geom) AS distance
  FROM stops s JOIN centers c
  ON (s.route = c.route)
  ORDER BY route, distance DESC
),
-- Filter out just the furthest stop/center pairs.
first_stops AS (
  SELECT DISTINCT ON (route) stops_distance.*
  FROM stops_distance
)
-- Pass the route/stop information into the linear route generation function!
SELECT
  ascii(route) AS gid, -- QGIS likes numeric primary keys
  route,
  walk_subway(gid, route) AS geom
FROM first_stops;

-- Do some housekeeping too
ALTER TABLE nyc_subway_lines ADD PRIMARY KEY (gid);

Here’s what our final table looks like visualized in QGIS:

As usual, there are some problems with our simple understanding of the data:

there are actually two “S” (short distance «shuttle») trains, one in Manhattan and one in the Rockaways, and we join them together because they are both called “S”;
the “4” train (and a few others) splits at the end of one line into two terminuses, so the «follow one line» assumption breaks and the result has a funny hook on the end.

Hopefully this example has provided a taste of some of the complex data manipulations that are possible combining the advanced features of PostgreSQL and PostGIS.

Introduction to PostGIS

40. Construcción de Geometrías¶

40.1. Ver también¶

Introduction to PostGIS

40. Construcción de Geometrías¶

40.1. Ver también¶

Navegación