Ejemplos de SQL geoespacial (I)

En un post anterior, hicimos una breve introducción al SQL geoespacial, mencionando el Dimensionally Extended nine-Intersection Model (DE-9IM), su implementación mediante predicados espaciales y, finalmente, se explicó el concepto de índice geoespacial. En este post y el siguiente, veremos algunos ejemplos específicos de consultas SQL geoespaciales.

Usaremos para estos dos posts unos datos de prueba, correspondientes a Colombia, que pueden ser descargados de aquí. Tras descomprimir el fichero, veremos que contiene cuatro archivos CSV.

Cargaremos estos datos en CartoDB, lo cuál nos ahorrará instalar ningún tipo de software en nuestra máquina, para realizar los ejercicios. Dado que CartoDB se basa en PostgreSQL junto con la extensión espacial PostGIS, podemos usar la API SQL de PostGIS

Carga de datos en CartoDB

Para cargar los ficheros de datos, simplemente los arrastramos al Dashboard, y CartoDB se encargará de crear la tabla en PostgreSQL automáticamente, y llevarnos a la pantalla donde veremos sus datos. En la siguiente imagen, se puede ver una captura de la tabla barrios_de_bogota, que será ampliamente utilizada en los ejemplos

Si pinchamos en Map View, en la parte superior, podemos ver ya directamente nuestros datos sobre un mapa

Adicionalmente, ya se han creado índices espaciales sobre la columna geométrica, así que ya estaríamos listos para comenzar a ejecutar consultas con estos datos.

Repetiremos la operación con el resto de ficheros, y ya estamos listos para trabajar.

Uso de predicados geoespaciales simples

Veremos a continuación el uso práctico de los predicados geoespaciales mencionados en el post anterior, aplicados a los datos específicos que hemos cargado en CartoDB

Intersects

Vamos a realizar una prueba del predicado ST_Intersects en la tabla barrios_de_bogota, formada por polígonos. Para ello, primero etiquetaremos los polígonos con su id, y así podremos referenciarlos en las consultas.

Etiquetar polígonos en CartoDB es muy sencillo. Basta con seguir estos pasos

Desde la vista de mapa de la tabla, elegimos el icono de wizards en la barra lateral

En el primer wizard (simple), seleccionamos gid en el desplegable Label Text. Elegimos un tamaño de letra suficientemente grande como para que se vea bien (en nuestro ejemplo, 18)

Con esto, ya veremos nuestros polígonos etiquetados con su atributo gid.

Como podemos ver en la captura, los polígonos que representan los barrios de Bogotá intersectan unos con otros. Concretamente, vemos que el polígono con gid 16 intersecta con los polígonos con gid 8 y 15

Comprobemos este hecho con una consulta SQL geoespacial, que ejecutaremos en la consola SQL de CartoDB

WITH poligono_16 as (select the_geom from barrios_de_bogota where gid = 16)

SELECT gid
FROM barrios_de_bogota
WHERE ST_Intersects(the_geom, poligono_16.the_geom)
AND gid != 16

El resultado, como esperábamos, consiste en los ids 8 y 15.

Hay dos detalles a destacar sobre esta consulta:

Por definición, un polígono intersecta consigo mismo. Por esa razón, se añade la condición gid != 16, excluyendo al propio polígono de la comparación

Mediante el uso de la cláusula WITH, evitamos ejecutar una consulta varias veces, precalculándola antes de ejecutar la consulta principal, y utilizando su resultado en la misma.  La consulta alternativa a ésta habría sido la siguiente

SELECT gid
FROM barrios_de_bogota
WHERE ST_Intersects(the_geom, (select the_geom from barrios_de_bogota where gid = 16))
AND gid != 16

Nótese que esta segunda consulta ejecuta la instrucción SELECT más interna una vez por cada fila de la tabla barrios_de_bogota. Algo que no es necesario. En este caso, al ser la tabla bastante pequeña, realmente no notamos gran diferencia en cuanto a tiempos de procesamiento, pero en el caso de grandes tablas, con millones de registros, optimizaciones como ésta pueden significar la diferencia entre llegar a obtener o no un resultado en un tiempo razonable.

Contains

Veamos a continuación un ejemplo del predicado Contains, para responder a la pregunta ¿En que barrio se encuentra el Museo del 20 de Julio?

Ejecutaremos la consulta, una vez más, desde la consola SQL de la tabla barrios_de_bogota

SELECT b.name from barrios_de_bogota b, points p
WHERE ST_Contains(b.the_geom, p.the_geom) and
p.name = ‘Museo del 20 de Julio’

El resultado de la consulta es San Cristobal, tal y como esperábamos.

En esta consulta, hemos usado dos tablas. La ya conocida barrios_de_bogota, y la tabla points, que contiene POIs (puntos de interés) de Colombia, entre los que se encuentra el mencionado Museo del 20 de Julio

Distance

Veremos en este apartado uno de los ejemplos más comunes: cálculo de distancias. Con una pequeña variante. Aplicaremos el predicado ST_DWithin, implementado por PostGIS, que nos sirve para calcular elementos geométricos a una distancia máxima de uno dado.

Supongamos que estamos interesados en saber qué puntos de interés se encuentran a una distancia máxima de 2km de la oficina de turismo bogotanismo.com, en la ciudad de Bogotá. Para ello, ejecutaremos esta consulta SQL

WITH oficina_turismo as (SELECT ST_Transform(the_geom, 21818) as the_geom

               FROM points

               WHERE name=’Bogotanisimo.com’)

SELECT name

FROM points a, oficina_turismo b

WHERE

       name is not null and

       name != ” and

       name != ‘Bogotanisimo.com’ and

       ST_DWithin(

            ST_Transform(a.the_geom, 21818),

            b.the_geom,

            2000

          )

El resultado de la consulta es la panadería Los Hornitos

En esta consulta también hemos utilizado la cláusula WITH, para evitar repetir una consulta. Además, hemos usado el predicado ST_DWithin, que nos permite calcular los puntos a una distancia máxima del punto que nos interesa. Por último, destacamos el uso de la función ST_Transform, de PostGIS. Esta función sirve para transformar geometrías de un sistema de coordenadas a otro. Es necesario porque las geometrías almacenadas en nuestras tablas utilizan un sistema de coordenadas geográficas, en el que la unidad de medida son los grados, en lugar de los metros. Si queremos medir distancias en metros, debemos usar un sistema de coordenadas que utilice metros, como es el caso del sistema de coordenadas con SRID 21818, usado en Colombia.

Crosses

Para finalizar, veremos el uso del predicado ST_Crosses de PostGIS, que nos permitirá obtener de manera sencilla los nombres de todos los barrios de Bogotá por los que pasa el Río Bogotá. La consulta la ejecutaremos desde la consola SQL de la tabla barrios_de_bogota, y es la siguiente

SELECT b.name
FROM barrios_de_bogota b JOIN waterways w
ON ST_Crosses(b.geom, w.geom)
WHERE w.name = ‘Rio Bogotá’

El resultado consta de varias filas, con los nombres de los barrios deseados

Terminamos, por ahora, con estos ejercicios de SQL geoespacial. En el próximo post, veremos algunos sencillos ejemplos de análisis espacial usando los predicados espaciales que nos proporciona PostGIS.