Importando información raster con el comando raster2pgsql de PostGIS en Debian

En un post anterior, se describió la manera de introducir información vectorial a una base de datos con componente espacial PostGIS mediante el comando shp2pgsql-gui. Esta información puede ser desplegada visualmente porque GIS’s de escritorio como QGIS, gvSIG o OpenJUMP tienen funcionalidad para permitir conectar directamente a la base de datos.

En el caso de raster, la posibilidad de incorporar esta información en modo visual para Debian desaparece con el comando raster2pgsql (no existe el equivalente raster2pgsql-gui) y el antiguo plugin de visualización de información raster en bases PostGIS ha quedado inservible al pasar a las versiones 2.x de QGIS. No obstante, es posible la visualización para aquellos que usen gvSIG.

Para utilizar el comando y verificar que los resultados son comparables, vamos a utilizar el mismo ejemplo del libro “PostGis Cookbook” que sugiere bajar las imágenes de tmin y tmax de http://www.worldclim.org/current con una resolución de 10 minutos de arco. Estas imágenes las desempaqueté en la carpeta worldclim de mi home de usuario y las renombré para que las numeración reflejase valores de dos dígitos (…01… en lugar de …1…, etc) para los 12 meses del año.

Para cargar el raster a la base de datos PostGis voy a usar, en consola de bash y ubicado en mi home de usuario (observen que hay una referencia relativa a tmax01.bil ubicado en worldclim), una adaptación del sugerido en el libro porque el que está allí no me funciona:

raster2pgsql -I -C -F -t 100x100 -s 4326 worldclim/tmax01.bil chp01.tmax01 | psql -d postgis_cookbook -W -f tmax01.sql

Lo que significa cada modificador del comando (-I, -C, etc) pueden averiguarlo con un raster2pgsql -?. No obstante, voy a señalar que -s hace referencia al SRID de la proyección del ráster (EPSG: 4326, long/lat WGS84) y -W hace posible la solicitud del password del propietario de la base de datos (a mi, particularmente, el sugerido -U nombre_usuario no me funciona en Debian).

Después de cargar exitosamente el ráster a la base de datos con la instrucción anterior, vamos a conectarnos a ésta y a explorar su contenido.

su
Contraseña: *********
# su - zeito
psql postgis_cookbook
Contraseña: *********
psql (9.4.2)
conexión SSL (protocolo: TLSv1.2, cifrado: ECDHE-RSA-AES256-GCM-SHA384, bits: 256, compresión: desactivado)
Digite «help» para obtener ayuda.

postgis_cookbook=#

Con el comando:

postgis_cookbook=# \d chp01.tmax01

podremos averiguar las columnas, tipo, modificadores, índices y restricciones que fueron creadas en la base de datos.

Como el ráster se teseló en bloques de 100 x 100 píxeles (indicado por la opción -t de raster2pgsql) entonces el comando siguiente devuelve el valor de 198:

postgis_cookbook=# SELECT count(*) FROM chp01.tmax01;
count
-------
198
(1 row)

Para entender de donde sale este valor empleamos, en consola de bash, el comando gdalinfo siguiente:

gdalinfo PG":host=localhost port=5432 dbname=postgis_cookbook user=zeito password=zeito1955 schema='chp01' table='tmax01' mode=2"

Driver: PostGISRaster/PostGIS Raster driver
Files: none associated
Size is 2160, 900
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.166666666666667,-0.166666666666667)
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000, -60.0000000) (180d 0' 0.00"W, 60d 0' 0.00"S)
Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
Lower Right ( 180.0000000, -60.0000000) (180d 0' 0.00"E, 60d 0' 0.00"S)
Center      (   0.0000000,  15.0000000) (  0d 0' 0.00"E, 15d 0' 0.00"N)
Band 1 Block=100x100 Type=Int16, ColorInterp=Gray
  NoData Value=-9999

El tamaño del ráster es de 900 filas por 2160 columnas (1944000 píxeles) por lo que dividido entre los 100 x 100 píxeles (10000) de los bloques con los que está almacenado en la base de datos nos da 194.4. El sistema los acomoda en 198 bloques para adaptarlo a una configuración rectangular. Más adelante veremos como los distribuye finalmente utilizando el comando ogr2ogr para producir una rejilla (temp_grid.shp) que permite visualizar la distribución de estos bloques.

ogr2ogr temp_grid.shp PG:"host=localhost port=5432 dbname='postgis_cookbook' user='zeito' password='su_password'" -sql "SELECT rid, filename, ST_Envelope(rast) as the_geom FROM chp01.tmax01"

La última columna del archivo temp_grid.shp es más angosta que las demás (por la diferencia entre 194.4 y 198); tal como se presenta a continuación:

grid_temp

Para finalizar, vamos primero a incluir, mediante shp2pgsql-gui, el shapefile countries.shp en la base de datos PostGis según los parámetros que se encuentran en la imagen siguiente:

shp2pgsql1

Esto es para “intersectar” el centroide de cada feature (países) en el shapefile con el ráster y seleccionar los 10 países cuyo valor de temperatura máxima sea la más baja en el mes de enero. El código y el resultado es el siguiente:

SELECT * FROM (
SELECT c.name, ST_Value(t.rast, ST_Centroid(c.the_geom))/10 as tmax_jan FROM chp01.tmax01 AS t
JOIN chp01.countries AS c
ON ST_Intersects(t.rast, ST_Centroid(c.the_geom))
) AS foo
ORDER BY tmax_jan LIMIT 10;

                  name                  | tmax_jan 
----------------------------------------+----------
 Greenland                              |    -29.8
 Russia                                 |    -27.9
 Canada                                 |    -27.4
 Svalbard                               |    -13.6
 Mongolia                               |    -11.7
 Kazakhstan                             |    -11.4
 Tajikistan                             |      -10
 Korea, Democratic People's Republic of |     -8.5
 Kyrgyzstan                             |     -7.9
 Finland                                |     -6.8
(10 filas)

Observen que los values en el ráster se dividen entre 10 para obtener los valores respectivos de temperatura centígrada (ver metadata del ráster).

Por otra parte, si en lugar de ordenar por tmax_jan ordenamos por nombre de país (name), entonces aparecen “huecos” en la información para algunos de ellos (los archipiélagos); tal como se presenta a continuación:

SELECT * FROM (
SELECT c.name, ST_Value(t.rast, ST_Centroid(c.the_geom))/10 as tmax_jan FROM chp01.tmax01 AS t
JOIN chp01.countries AS c
ON ST_Intersects(t.rast, ST_Centroid(c.the_geom))
) AS foo
ORDER BY name LIMIT 10;

        name         | tmax_jan 
---------------------+----------
 Afghanistan         |     -5.4
 Åland Islands       |     -0.8
 Albania             |      8.8
 Algeria             |     17.5
 American Samoa      |         
 Andorra             |      0.6
 Angola              |     26.7
 Anguilla            |     28.1
 Antigua and Barbuda |         
 Argentina           |     33.6
 Armenia             |     -5.2
 Aruba               |     29.1
 Australia           |     37.4
 Austria             |     -2.2
 Azerbaijan          |      6.5
 Bahamas             |         
 Bahrain             |     20.7
 Bangladesh          |     25.5
 Barbados            |     28.3
 Belarus             |     -3.9
(20 filas)

Esto es producto del hecho de tomar como criterio valido de intersección el centroide. Para los países insulares el centroide probablemente caerá en cualquiera de los océanos y en el raster esto corresponde a ‘no data’.

Esta entrada fue publicada en Postgres+postgis, SIG, Software Libre. Guarda el enlace permanente.

2 respuestas a Importando información raster con el comando raster2pgsql de PostGIS en Debian

  1. Oriel Romero dijo:

    Hola Jose , la pregunta es en realidad para el admin de django, quisiera saber si puedo cargar coordenadas del tipo gg:mm:ss,ssss y convertirlas a grados decimales. gracias.

    • No he usado django pero por lo que veo en su documentación lo que subyace es Python. Independientemente del administrador de base de datos con componente espacial que pienses usar (PostGIS?), el algoritmo de conversión es:

      gg+mm/60+ss.ssss/3600 = grados decimales

      e incluirlo en un script algo muy sencillo.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s