A lo largo de la última sesión del Master GIS Web y a la par que desplegábamos algunos pesados procesos para el tratamiento masivo de información, nos hemos planteado y hemos resuelto, de la mano de uno de nuestros alumnos , Alberto Callejas, un interesante problema de geoprocesos que nos ha permitido valorar el rendimiento de dos de las herramientas más usadas para realizar geoprocesos.

El resultado del trabajo de Alberto, nos ha servido además para incluirlo en el pull de procesos que usamos en varios proyectos en los que actualmente estamos involucrados. Es decir, nuestro alumno Alberto ha resuelto un problema REAL y PRÁCTICO de una forma brillante.

El problema

Debido a la alta demanda de petróleo y gas se va a proceder a realizar una serie de extracciones usando el tan polémico método de “fracking” en ciertos lugares previamente determinados cuyos terrenos van a ser destinados al efecto.

Por parte del organismo estatal competente (ente ficticio) se quiere estudiar los posibles daños a estructuras y usos del suelo que se puede provocar. Para ello se solicita un estudio de impacto a un técnico especializado en análisis espaciales en entornos SIG.

Datos de partida

  • Puntos –> Capa de pozos
  • Polígonos –> Capa de usos del suelo

Resultado buscado

Se quiere obtener los polígonos resultantes de la intersección del buffer de los pozos con la capa de usos del suelo. Además es útil saber el porcentaje del área del polígono resultante de la intersección respecto al área total del polígono de ese uso del suelo.

Alternativas:

1 – GIS de escritorio ArcGis o gvSIG: Buffer pozos + Intersect con Usos + Cálculo de porcentajes

El proceso lleva más de un día y se puede ejecutar con una condiciones de memoria virtual disponible muy específicas. Inviable si la base de datos crece. Que crecer

2 – Uso de funciones ST dentro de la base de datos PostGIS

El proceso lleva 15 minutos con un consumo mucho menor de recursos. Óptimo. Más abajo os detallamos el proceso del Análisis con funciones ST de PostGIS

 

Ejemplo gráfico

imagen1_PostGIS

El técnico aplica un buffer de 400 metros de distancia a cada pozo (área de influencia) e interseca a los distintos usos del suelo contenidos en cada buffer.

imagen2_PostGIS

imagen3_PostGIS

Buffer de 400m al pozo con identificador 270393187.

 imagen4_PostGIS

 

Ejemplo de intersección realizada a un pozo con identificador único igual a 270393187.

imagen5_PostGIS

 

Usos del suelo dentro del área de influencia del pozo.

 

Geoproceso realizado con funciones ST de PostGIS

El software utilizado en este caso es Postgresql y, además debido a la naturaleza espacial de la consulta se ha creado la extensión PostGIS en la base de datos para poder utilizar las geometrías de las tablas para realizar una consulta espacial.

PostgreSQL es un sistema de gestión de bases de datos objeto-relacional, distribuido bajo licencia BSD y con su código fuente disponible libremente. Es el sistema de gestión de bases de datos de código abierto más potente del mercado y en sus últimas versiones no tiene nada que envidiarle a otras bases de datos comerciales. PostgreSQL utiliza un modelo cliente/servidor y usa multiprocesos en vez de multihilos para garantizar la estabilidad del sistema. Un fallo en uno de los procesos no afectará el resto y el sistema continuará funcionando (fuente: http://www.postgresql.org.es/sobre_postgresql).

Pasos que vamos a seguir:

Paso 1 – Cargar las capas en la base de datos

Paso 2 – Escribir los campos que queremos que resulten del proceso

Paso 3 – Intersección de buffer y usos del suelo

Paso 4 – Union con los índices espaciales

Paso 5 –  Cálculo de superficies y porcentajes

Paso 1 – Cargar las capas en la base de datos

El primer paso es cargar los shapefiles en la base de datos PostGIS. Para simplificar, hemos cargado en la base de datos la capa de pozos con un buffer hecho. Pero podríamos haber lanzado esta operación también en base de datos.

  • sscc_30n: Tabla con usos del suelo. 2377 registros. Tiene un campo campo “geometry” (en el que se estipula el tipo de geometría) de tipo “MultiPolygon” y con un SRID = 32630.
  • buf_pozos_30n: Tabla con los pozos. con un campo “geometry” de tipo “MultiPolygon” y el mismo SRID (1)  que la tabla “sscc_30n”. Esta tabla es el resultado de aplicar un buffer de distancia igual a 400 metros sobre cada pozo con coordenadas X,Y

imagen6_PostGIS

Imágenes de tablas pozos_buffer y usos del suelo

Paso 2 – Escribir los campos que queremos que resulten del proceso

Como resultado esperado buscamos una tabla con los campos gid (integer), id_local (bigint), id_poligono (bigint), porcentaje (double precision) y geometry (multipolygon).

Paso 3 Intersección de buffer y usos del suelo

ST_Intersection: Para intersectar

ST_Astext: Para quedarnos con datos interpretables a primera vista

STX_Extract: Para devolver sólo polígonos en la intersección

i) ST_Intersection + ST_Astext: Primero tenemos que saber qué parcelas “intersecan” con el buffer de la tabla buf_pozos_30n. Para ello usamos la función ST_Intersection(geom1, geom2). Esta función devuelve las geometrías de la tabla 1 que intersectan con las geometrías de la tabla 2. Utiliza como argumentos las geometrías de las dos tablas.

Esta función devuelve no sólo las geometrías de tipo superficial que intersectan sino también las que no (valores nulos mostrados como GEOMETRY COLLECTION EMPTY).Con la función **ST_Astext recuperamos las geometrías como texto, lo cual facilita su interpretación, ya que de otro modo se mostraría codificado.

imagen7_PostGIS

ii) STX_Extract Para obtener sólo las geometrías de tipo superficial necesitamos usar la función STX_Extract. La función STX_Extract extrae de la geometría pasada como primer argumento las geometrías que tengan la dimensión específicada en ndimensions (0-puntuales, 1-lineales, 2- superficiales). En este caso usamos la dimension 2 ya que se trata de superficies. Si la geometría no contiene ninguna geometría de la dimension ndimension entonces devuelve null.

 

Para eliminar los valores nulos y sólo obtener los registros con geometría de ndimension igual a 2 debemos usar una concatenación externa (JOIN) entre las dos tablas. En este caso Inner Join.

El INNER JOIN permite emparejar filas de distintas tablas de forma más eficiente que con el producto cartesiano cuando una de las columnas de emparejamiento está indexada. Ya que en vez de hacer el producto cartesiano completo y luego seleccionar la filas que cumplen la condición de emparejamiento, para cada fila de una de las tablas busca directamente en la otra tabla las filas que cumplen la condición, con lo cual se emparejan sólo las filas que luego aparecen en el resultado.

Esto nos permite obtener los registros de cada tabla que tienen un campo en común, evitando así la obtención de valores nulos.

En este caso el campo común es el campo geometría (MultiPolygon en ambas tablas) por lo tanto hemos de indexar los campos geometry de las tablas sscc_30n y buf_pozos_30n.

Ahora con los índices espaciales creados y usando Inner Join:

imagen8_PostGIS

En nuestro caso partimos en un producto cruzado entre la tabla pozos y uso del suelo, o lo que es lo mismo 140.018 * 2.377 = 332.822.786 parejas de registros.

La razón por la que usa el Inner Join + ST_Intersects(b.geom, s.geom) en la consulta es porque es la única forma de “llamar” a los índices espaciales anteriormente creados. Estas líneas indican al programa que limite la consulta a los registros de las dos tablas cuyas geometrías intersecan, dejando el resto de registros fuera de la consulta. Esto agiliza muchísimo la consulta y más aún cuando se trabaja con un número elevado de registros.

Los registros devueltos son las geometrías de ndimension = 2 resultado de la intersección entre las geometrías de las dos tablas. No devuelve valores nulos ya que usa como argumentos sólo las geometrías de ndimension 2 presentes en las dos tablas.

Ahora tenemos que calcular el área de dichas intersecciones y dividirla entre el área de la parcela original para hallar el porcentaje: ** quitamos ST_Astext de la sentencia.

Campo porcentaje:

imagen9_PostGIS

Vamos a por el campo geometría:

imagen10_PostGIS