PostGIS Raster¶

Introducción¶

Desde la versión 2.0 de PostGIS, es posible cargar y manipular datos de naturaleza ráster en una base de datos espacial, gracias a PostGIS Raster: un nuevo conjunto de tipos y funciones que dotan a PostGIS de la posibilidad de manipular datos ráster.

El objetivo de PostGIS Raster es la implementación de un tipo de datos RASTER lo más parecido posible al tipo GEOMETRY de PostGIS , y ofrecer un único conjunto de funciones SQL que operen de manera transparente tanto en coberturas vectoriales como en coberturas de tipo ráster.

Como ya hemos mencionado, PostGIS Raster es parte oficial de PostGIS 2.0, de manera que no es necesario instalar ningún software adicional. Cuando una base de datos PostgreSQL es activada con PostGIS, ya es al mismo tiempo también activada con PostGIS Raster.

Tipo de datos Raster¶

El aspecto más importante a considerar sobre el nuevo tipo de datos RASTER definido por PostGIS Raster es que tiene significado por si mismo. O de otra forma: una columna de tipo RASTER es una cobertura raster completa, con metadatos y posiblemente geolocalizada, pese a que pertenezca a una cobertura raster mayor.

Tradicionalmente, las bases de datos espaciales con soporte para ráster, han permitido cargar y teselar coberturas ráster para operar con ellas. Normalmente, se almacenaban los metadatos de la cobertura completa por un lado (geolocalización, tamaño de píxel, srid, extensión, etc) y las teselas por otro, como simples chunks de datos binarios adyacentes. La filosofía tras PostGIS Raster es diferente. También permite la carga y teselado de coberturas completas, pero cada tesela por separado, contiene sus propios metadatos y puede ser tratada como un objeto raster individual. Además de eso, una tabla de PostGIS Raster cargada con datos pertenecientes a una misma cobertura:

• Puede tener teselas de diferentes dimensiones (alto, ancho).
• Puede tener teselas no alineadas con respecto a la misma rejilla.
• Puede contener huecos o teselas que se solapan unas con otras.

Este enfoque hace a PostGIS Raster una herramienta muy poderosa, aunque también tiene algunos problemas inherentes (como el bajo rendimiento en aplicaciones de visualización de datos ráster en tiempo real). Para saber más sobre PostGIS Raster el mejor sitio donde acudir es http://trac.osgeo.org/postgis/wiki/WKTRaster

En los siguientes apartados, veremos algunas de las operaciones que se pueden realizar con PostGIS Raster

Procesando y cargando raster con GDAL VRT¶

Como ya hemos visto en el capítulo de importación y exportación, la manera de cargar datos ráster en PostGIS Raster es a través de raster2pgsql. En este capítulo, trabajaremos con datos sobre las temperaturas en Colombia durante el mes de Septiembre del año 2010.

Tenemos datos relativos a todas las zonas del planeta y varios meses disponibles en la url http://www.worldclim.org. Las zonas que nos interesan a nosotros son 2: zona 23 y zona 33. La zona 23 contiene los datos de la mayor parte de Colombia, y el resto está en la zona 33. Las url de descarga de datos de esas zonas son:

• http://www.worldclim.org/tiles.php?Zone=23
• http://www.worldclim.org/tiles.php?Zone=33

Los datos en si son ficheros GeoTIFF con una sola banda. Aquí vemos ambos fragmentos visualizados en QGIS:

. image:: _images/temperaturas_colombia.png

scale:	50 %

Podríamos cargar los ficheros GeoTIFF por separado en una misma tabla, mediante dos llamadas a raster2pgsql, pero lo que vamos a hacer es construir un raster virtual en formato VRT, y cargar ese ráster en PostGIS Raster con una única llamada a raster2pgsql. Lo haremos en 3 pasos:

1. Construir fichero VRT a partir de los ficheros TIFF con gdalbuildvrt:

   # gdalbuildvrt tmean9.vrt tmean9_*.tif
   

2. Crear el fichero SQL a partir del fichero VRT generado en el paso anterior:

   # raster2pgsql -I -C -F -t 36x36 -P -M -s 4326 tmean9.vrt > tmean9.sql
   

3. Cargar el fichero SQL en PostGIS Raster:

   # psql -d taller_semana_geomatica -f tmean9.sql -U postgres
   

Con eso ya tendríamos ambos ficheros cargados en PostGIS Raster

Un problema que tenemos ahora es que los ficheros GeoTIFF tienen muchos más datos que simplemente los datos de Colombia. Es por eso que, para los siguientes apartados, vamos a recortar una parte del fichero GeoTIFF, que comprenda solo Colombia, y cargaremos esa parte en PostGIS Raster, con el nombre tmean9_colombia.

El recorte se hará usando la herramienta clipper de QGIS, que no es más que una interfaz gráfica de usuario para llamar a gdal_translate pasándole las coordenadas de inicio y fin y la altura y anchura del rectángulo a obtener. A pesar de que esta operación se realizará en clase, no se profundizará en ella, por exceder de los límites del curso. Se realiza únicamente para justificar la existencia de la tabla utilizada en los ejemplos posteriores.

La imagen recortada queda así (se ha aplicado un pseudo-color a la capa para apreciar el contorno de Colombia dentro del fichero GeoTiff):

. image:: _images/temperaturas_recorte_colombia.png

scale:	50 %

Bastante más pequeña y manejable.

De todas formas, aun podemos afinar más esta operación de recorte. En un apartado posterior veremos como utilizar la geometría que define los límites de Colombia como molde para quedarnos únicamente con la porción del ráster contenida dentro de esos límites.

Obtención de metadatos y estadísticas de una capa PostGIS Raster¶

Mediante consultas SQL, es posible obtener metadatos y estadísticas de las capas ráster almacenadas.

#SELECT
        r_table_name,
        r_raster_column,
        srid,
        scale_x,
        scale_y,
        blocksize_x,
        blocksize_y,
        same_alignment,
        regular_blocking,
        num_bands,
        pixel_types,
        nodata_values,
        out_db,
        ST_AsText(extent) AS extent
FROM raster_columns WHERE r_table_name = 'tmean9_colombia';

# SELECT
        rid,
        (ST_BandMetadata(rast, 1)).*
FROM tmean9_colombia
WHERE rid = 1265;

#  rid | pixeltype | nodatavalue | isoutdb | path
 ------+-----------+-------------+---------+------
  1266 | 32BF      |             | f       |

# WITH stats AS (
        SELECT
                (ST_SummaryStats(rast, 1)).*
        FROM tmean9_colombia
        WHERE rid = 1266
)
SELECT
        count,
        sum,
        round(mean::numeric, 2) AS mean,
        round(stddev::numeric, 2) AS stddev,
        min,
        max
FROM stats;

# count |  sum   |  mean  |  stddev | min | max
 -------+--------+--------+---------+-----+-----
   1296 | 326501 | 251.93 |    7.21 | 223 | 263

MapAlgebra sobre capas PostGIS Raster¶

En el apartado anterior, vimos como los valores de temperaturas de la capa ráster estaban escalados por 100. Vamos a cambiar todos estos valores usando una expresión de MapAlgebra. Para ello, añadiremos una nueva banda con los valores cambiados:

# UPDATE tmean9_colombia SET
        rast = ST_AddBand(
                rast,
                ST_MapAlgebraExpr(rast, 1, '32BF', '[rast] / 100.', -9999),
                1
        );

En la llamada a MapAlgebra, hemos especificado que la banda de salida tendrá un tamaño de píxel de 32BF y un valor NODATA de -9999. Con la expresión [rast] / 100, convertimos cada valor de píxel a su valor previo al escalado.

Tras ejecutar esa consulta, el resultado es éste:

# ERROR:  new row for relation "tmean9_colombia" violates check constraint "enforce_out_db_rast"
DETAIL:  Failing row contains (1, 0100000200563C2A37C011813F18FD8BFEC51081BF00000000426E54C0000000..., tmean9_colombia.tif)

Como vemos, la consulta no ha funcionado. El problema es que, cuando cargamos esta capa ráster usando raster2pgsql, especificamos el flag -C. Este flag activa una serie de restricciones en nuestra tabla, para garantizar que todas las columnas de tipo RASTER tienen los mismos atributos (más información en http://postgis.net/docs/manual-2.0/RT_AddRasterConstraints.html).

El mensaje de error nos dice que hemos violado una de esas restricciones. Concretamente la restricción de out-db. A primera vista, puede parecer extraño, porque nosotros no estamos especificando que la nueva banda sea de tipo out-db. El problema es que esta restricción solo funciona con una banda, y si se intenta añadir una segunda banda a un ráster que ya tiene una, la restricción lo hace fallar.

La solución a nuestro problema pasa por:

1. Eliminar las restricciones de la tabla mediante DropRasterConstraints
2. Volver a ejecutar la consulta
3. Volver a activar las restricciones (OJO: Es una operación costosa en datos raster muy grandes)

Las consultas a ejecutar son las siguientes:

# SELECT DropRasterConstraints('tmean9_colombia', 'rast'::name);
# UPDATE tmean9_colombia SET rast = ST_AddBand(rast, ST_MapAlgebra(rast, 1, '32BF', '[rast] / 100.', -9999), 1);
# SELECT AddRasterConstraints('tmean9_colombia', 'rast'::name);

Y el resultado es:

# droprasterconstraints
-----------------------
t

# UPDATE 2950

# addrasterconstraints
----------------------
t

Ahora comprobaremos que una nueva banda ha sido añadida a nuestro ráster:

# SELECT
        (ST_Metadata(rast)).numbands
FROM tmean9_colombia
WHERE rid = 1266;

Devuelve:

# numbands
----------
2

¿Y cuáles son los detalles de esas dos bandas?:

# WITH stats AS (
        SELECT
                1 AS bandnum,
                (ST_SummaryStats(rast, 1)).*
        FROM tmean9_colombia
        WHERE rid = 1266
        UNION ALL
        SELECT
                2 AS bandnum,
                (ST_SummaryStats(rast, 2)).*
        FROM tmean9_colombia
        WHERE rid = 1266
)

SELECT
        bandnum,
        count,
        round(sum::numeric, 2) AS sum,
        round(mean::numeric, 2) AS mean,
        round(stddev::numeric, 2) AS stddev,
        round(min::numeric, 2) AS min,
        round(max::numeric, 2) AS max
FROM stats
ORDER BY bandnum;

El resultado es:

# bandnum | count |    sum    |  mean  | stddev |  min   |  max
 ---------+-------+-----------+--------+--------+--------+--------
        1 |  1296 | 326501.00 | 251.93 |   7.21 | 223.00 | 263.00
        2 |  1296 |   3265.01 |   2.52 |   0.07 |   2.23 |   2.63

Vemos que el valor en la banda 2 ha sido corregido, dividiendo los valores de temperaturas entre 100. Ahora las temperaturas tienen sentido como grados centígrados

Clip de datos ráster usando geometrías¶

Una de las grandes ventajas de poder tener datos de naturaleza ráster y vectorial cargados en PostGIS es que se puede operar con ellos mediante la utilización de la misma API SQL. En este ejemplo, veremos como recortar un raster usando una geometría como modelo.

Trabajaremos con los datos ráster de temperaturas, y con los datos vectoriales de Colombia. Como vemos en esta imagen (coloreada con pseudocolor en QGIS 2.0), el ráster ocupa bastante más extensión que Colombia:

. image:: _images/raster_with_vector.png

scale:	50 %

Lo que queremos es recortar la parte del ráster que queda dentro de los límites de Colombia. Y lo haremos únicamente con consultas SQL. Posteriormente, volcaremos ese ráster recortado a disco, en formato GeoTIFF.

La consulta que se queda solamente con la parte del ráster comprendida dentro de los límites de Colombia es:

# CREATE TABLE tmean9_colombia_clip AS
SELECT t.rid, t.rast, c.admin_name
FROM tmean9_colombia t JOIN co c ON ST_Intersects(t.rast, c.geom)

Con esa consulta hemos logrado crear una tabla con datos ráster únicamente comprendidos dentro de los límites de Colombia. Para visualizar esa tabla, tenemos dos opciones. Ambas requieren de GDAL 2.0

• Volcar el contenido de la tabla a disco, a formato GeoTIFF, mediante el uso de gdal_translate http://www.gdal.org/gdal_translate.html
• Instalar en QGIS el plugin de visualización de PostGIS Raster. El problema es que aun no se ha portado el plugin a la versión 2.0 de QGIS

Elegimos la primera opción, por no requerir la instalación de ningún software adicional. La orden que debemos ejecutar es:

# gdal_translate PG:"host=localhost port=5432 dbname=taller_semana_geomatica user=postgres password=postgres table=tmean9_colombia_clip mode=2" tmean9_colombia_clip.tif

Y el aspecto de este ráster recortado una vez colocado sobre el mapa y coloreado con pseudocolor en QGIS 2.0 es:

. image:: _images/postgis_raster_clipped.png

scale:	30 %

Combinando raster y geometrías para análisis espacial¶

Vamos a ver ahora cuáles fueron las temperaturas máximas, mínimas y medias de todos los barrios de Bogotá durante el mes de Septiembre. Para ello, usaremos nuevamente la API SQL de PostGIS y PostGIS Raster junto con las funciones de agregación de PostgreSQL.

La consulta a realizar es la siguiente:

# WITH stats AS (
        SELECT rast, (ST_SummaryStats(rast, 2)).*
        FROM tmean9_colombia_clip
)

SELECT
        b.name,
        ROUND(AVG(s.mean::numeric), 2) AS tmean,
        ROUND(AVG(s.min::numeric), 2) as tmin,
        ROUND(AVG(s.max::numeric), 2) as tmax
FROM stats s JOIN barrios_de_bogota b ON ST_Intersects(s.rast, b.geom)
GROUP BY b.name
ORDER BY b.name

El resultado es el siguiente:

#         name      | tmean | tmin | tmax
----------------+-------+------+------
 Antonio Nariño |  1.19 | 0.63 | 2.04
 Barrios Unidos |  1.25 | 0.79 | 1.72
 Bosa           |  1.39 | 0.66 | 2.23
 Chapinero      |  1.25 | 0.79 | 1.72
 Ciudad Bolívar |  1.19 | 0.63 | 2.04
 Ciudad Kennedy |  1.19 | 0.63 | 2.04
 Engativá       |  1.25 | 0.79 | 1.72
 Fontibón       |  1.25 | 0.79 | 1.72
 Los Mártires   |  1.19 | 0.63 | 2.04
 Puente Aranda  |  1.19 | 0.63 | 2.04
 Rafael Uribe   |  1.19 | 0.63 | 2.04
 San Cristóbal  |  1.19 | 0.63 | 2.04
 Santa Fé       |  1.19 | 0.63 | 2.04
 Suba           |  1.30 | 0.96 | 1.41
 Sumapáz        |  1.10 | 0.46 | 2.19
 Teusaquillo    |  1.19 | 0.63 | 2.04
 Tunjuelito     |  1.19 | 0.63 | 2.04
 Usaquén        |  1.25 | 0.79 | 1.72
 Usme           |  1.08 | 0.56 | 2.08