SELECT * FROM incendio1;

Podemos ver en la parte inferior de PgAdmin 4 que nos aparece la tabla que hemos seleccionado. Si nos fijamos bien, aparece una columna con el nombre de “geom” con los datos geoespaciales de nuestro polígono. A la derecha de la tabla, verás un icono azul con un ojo en el centro. Clica en él para ver las geometrías de la tabla. ¿Alucinante, no?

Perfecto, unas líneas atrás hemos comentado que deberíamos agrupar todas las tablas en una. Para ello, podemos utilizar Union all. El código sería el siguiente:

SELECT * 
INTO   tabla_incendios 
FROM   (SELECT * FROM   incendio1 
	UNION ALL 
	SELECT * FROM   incendio2
	UNION ALL 
	SELECT * FROM   incendio3
	UNION ALL  
	SELECT * FROM   incendio4
	UNION ALL  
	SELECT * FROM   incendio5
	UNION ALL  
	SELECT * FROM   incendio6) incendio

El uso de la CTE en SQL para el cálculo de la media entre geometrías

Ahora que ya hemos agrupado todos los datos en una misma tabla, podemos crear un polígono personalizado a partir de la intersección de sus bordes con los demás polígonos que hemos importado. Es decir, tal y como podemos ver en la imagen de abajo, todos los polígonos están superpuestos. Las intersecciones entre todos los bordes de todos los polígonos crean polígonos más pequeños. Cada una de estas geometrías tiene n polígonos por encima o debajo de ella. Es decir, con los que interseca. Nosotros queremos obtener aquel polígono que se ha creado a partir de todos estos pequeños polígonos donde n >= 3 (con 3 o más intersecciones), y así obtener una media personalizada.

Para realizar esta operación, debemos hacer más de un cálculo, lo que significa que tendremos que operar con resultados que no existen en el sistema de manera inherente. Una solución sería trabajar a partir de vistas, pero el problema es que estas pasan a ser permanentes en el sistema y eso no sería eficiente. Una buena aproximación es el uso de la Common Table Expression (CTE). La CTE nos permite crear apartados (como subquerys), que solo existirán durante la ejecución de esta consulta. Con el siguiente código podemos obtener el resultado esperado.

CREATE TABLE media_incendios AS 
WITH bordes AS (
    SELECT (ST_Dump(
     		ST_UnaryUnion(
     		ST_Collect(
     		ST_ExteriorRing(i.geom))))).geom AS geom
    FROM tabla_incendios i
),
 partes AS (
    SELECT (ST_Dump(ST_Polygonize(geom))).geom FROM bordes
),
partes_contadas AS (
    SELECT 
	 partes.geom AS geom, 
	 COUNT(*)
    FROM partes
    JOIN tabla_incendios i
    ON ST_Intersects(i.geom, ST_PointOnSurface(partes.geom))
    GROUP BY partes.geom
),
resultado AS (
SELECT ST_Union(geom) AS geometria
FROM partes_contadas
WHERE count > 3
) SELECT 
	geometria, 
	round((ST_Area(geometria)/10000)::numeric, 3) AS area 
  FROM resultado

En la primera línea de código especificamos el nombre de la tabla que vamos a crear y que va a ser el resultado de todo lo que se encuentra por debajo del código. Con WITH empezamos nuestra CTE.  A continuación, usamos:

• ST_ExteriorRing para devolver en Linestring los bordes de los polígonos. 
• ST_Collect para recoger todas las geometrías.
• ST_UnaryUnion para partir todas las líneas de las geometrías en sus respectivos puntos de unión. 
• ST_Dump para separar cada línea que ha sido partida en una sola geometría. 

El resultado de “bordes” es una tabla con todas las líneas de los perímetros de los polígonos que hemos partido. Si observamos las geometrías que se han creado, podemos ver lo siguiente:

En esta parte del código, con ST_Polygonize, hacemos una poligonización de las partes de los polígonos que son delimitados por todas las líneas que hemos creado anteriormente. 

La poligonización consiste en crear todos los polígonos que sean posibles utilizando todas las combinaciones entre las líneas obtenidas anteriormente.

Por otra parte, el uso de la función ST_Dump implica que todas las intersecciones entre esos perímetros se descompongan en líneas individuales. Algo imprescindible para la poligonización.

El resultado espacial sería el siguiente:

Ahora que ya tenemos todos los polígonos creados, podemos contar para cada uno de ellos el número de geometrías que tiene superpuestas (o con las que interseca). En primer lugar, para realizar este cálculo utilizaremos:

• ST_PointOnSurface para obtener una geometría de tipo punto dentro de cada polígono. Esta geometría, la utilizamos para intersecarla con los polígonos. 
• ST_Intersects donde especificamos como parámetros las geometrías de “tabla_incendios” y de “partes”. 

A nivel descriptivo, obtenemos un punto de cada polígono (ST_PointOnSurface) y lo intersecamos con todos los polígonos. La función COUNT(*) junto con la cláusula GROUP BY “partes.geom” nos proporciona el número de intersecciones de cada polígono.

El resultado es una tabla con todas las geometrías de la tabla “partes” con la cantidad de solapamientos con otras geometrías contados.

Por último, ya podemos juntar todas las partes en una sola dependiendo del valor que le especifiquemos a n y obtener el polígono deseado. Para ello, utilizaremos:

• ST_Union para juntar todas las geometrías que tienen un solapamiento más grande de tres. 

Una vez realizados todos los cálculos para obtener la media, es interesante calcular la extensión de nuestra zona de estudio. Para ello, en la última parte del código hemos utilizado:

• ST_Area para calcular el área que dividiremos por 10.000 para así obtener el resultado en km².
• ROUND para redondear el resultado de ST_Area en tres decimales. 

El resultado es un polígono que se ha creado a partir de calcular la media personalizada sobre el número de solapamientos de todas las geometrías con SQL. Al realizarlo con SQL podemos observar que el tiempo de espera de ejecución del código es muy bajo. Por lo tanto, nuestra aproximación al problema de la media ha sido solucionado de una manera eficiente y rápida. A continuación, podemos ver el resultado de “media_incendios”.

Cargar el resultado en QGIS desde PgAdmin 4

El resultado que hemos obtenido en el paso anterior lo podemos cargar a QGIS fácilmente. Para ello, nos vamos a nuestro proyecto de QGIS -> Layer -> Add Layer -> Add PostGis Layers donde veremos esta ventana:

Le damos al botón Connect y tendremos todas las capas de nuestra base de datos en nuestro proyecto de QGIS. Si lo deseamos, podemos ver todas las capas con las que hemos trabajado y marcar el borde de la capa final de nuestra operación con un color más intenso. En nuestro caso, el polígono resultante de la media es el que tiene el borde rojo más intenso.

Si te ha gustado el post, te invitamos a que no pierdas la motivación y sigas tu trayectoria en el mundo de las bases de datos. Sigue con nuestras píldoras de información en nuestro blog, y no olvides echarle un vistazo a nuestro Máster profesional en SIG, con el cual podrás especializarte en analista o programador SIG, unas de las posiciones más demandadas del mercado.