Miércoles, 08 Junio 2016 17:16

Extrayendo datos del DERA con Python

Los Datos Espaciales de Referencia de Andalucía para escalas intermedias -DERA- "es un repertorio de bases cartográficas de diferente naturaleza geométrica (puntos, líneas, polígonos, imágenes raster) referidas al territorio andaluz". Más información en esta página del Instituto de Estadística y Cartografía de Andalucía.

Los datos pueden descargarse en un archivo comprimido zip por cada uno de los bloques. Dentro de cada zip la información se encuentra accesibles por capas en formato shapefile (.shp), en sistema de referencia geodésico ETRS89 y proyectadas en UTM huso 30.

Filtrado de datos con SIG

La extensión espacial de las capas es el ámbito geográfico de la comunidad autónoma de Andalucía. Esto significa que por ejemplo en la capa sv01_sanidad_centro_salud.shp se encuentran todos los centros de salud existentes en Andalucía y sus diferentes tipologías. Si queremos trabajar con una selección de datos, por ejemplo por provincia, deberemos:

  • Descarga el zip correspondiente desde la web del IECA
  • Descomprimirlo el zip
  • Abrir un Sistema de Información Geográfica (ArcGIS, QGIS, gvSIG...)
  • Cargar la capa (ej. sv01_sanidad_centro_salud.shp)

Normalmente las capa del DERA incluyen información del municipio y provincia, por lo que podemos hacer una consulta por el atributo que deseemos, seleccionar los datos y crear una nueva capa a partir de la selección.

Una vez que tenemos los datos en nuestro SIG, se nos pueden plantear algunas preguntas:

  • ¿Qué ocurre si la capa no dispone de un atributo como municipio o provincia que permita filtrarla? Esto puede suceder para capas de tipo lineal o poligonal cuya extensión supere un límite administrativo.
  • ¿Cómo extraer información de ámbitos geográficos definidos por nosotros? Este puede ser el caso de un barrio, o la delimitación de una mancomunidad. En este caso debemos contar con la capa de delimitación y a continuación seleccionar los elementos geométricos (ej. centros de salud) mediante una consulta espacial (superposición, dentro de, toca...).
  • ¿Y si quiero obtener todos los datos de un bloque temático, como por ejemplo "Servicios", para un barrio? En este caso, tendremos que repetir la operación anterior para cada una de las capas.

Trabajando con Python y GDAL

Algunos de los problemas que he planteado antes pueden solucionarse con las herramientas de creación de modelos de procesado integradas en los SIG. En el siguiente enlace a MappingGIS hay una sencilla guía de cómo hacerlo en QGIS.

Como sigo con el aprendizaje de Python que comencé con dxf2gmlcatastro he decidido crear un pequeño código que haga lo siguiente:

  • Descomprimir un archivo zip de capas Shape, en este caso del DERA.
  • Localizar una archivo poligonal que va a ser usado para "recortar" las capas Shape del zip
  • Crear una carpeta donde almacenar las nuevas capas recortadas.
  • Generar nuevas capas recortadas en una carpeta concreta y con un sufijo (_clip).
  • Borrar la carpeta donde se ha descomprimido el zip.

Los requisitos para ejecutar el código son

  • Tener instalado Python.
  • Tener instalada la librería GDAL/OGR.

Un ejemplo

$ python clipShapesZip.py 'G15_Patrimonio.zip' 'clip_area.shp' 'clipFolder'

Como nota importante comentar que si el geoproceso se aplica sobre una lineal o poligonal, la función hará su función y "recortará" las geometrías a partir del la capa indicada.

Este artículo continua la entrada anterior en la se expusieron algunos de los usos de los comandos gdalinfo y gdaltranslate de la librería GDAL. En el siguiente texto nos vamos a centrar en el comando gdalwarp. Este comando se usa para hacer reproyecciones del Sistema de Referencia de Coordenadas de una imagen georeferenciada.

Reproyección de una imagen georeferenada con gdalwarp

El uso básico es sencillo: partimos de la base de que la imagen ya posee un SRC definido, tras el comando indicaremos el SRC de salida usando el código EPSG (ej. -t_srs "EPSG:4326"), el fichero a reproyectar y el nombre la imagen proyectada.

$ gdalwarp  -t_srs "EPSG:4326"  -of GTiff img.tif img4326.tif

Las opciones de gdalwarp son muchos más variadas establecer un valor sin datos (- dstnodata) o usar una archivo vectorial para recortar la imagen de salida (-cutline). Todas las opciones pueden consultarse en la siguiente dirección http://www.gdal.org/gdalwarp.html

De ED50 a ETRS89

En España, a partir del 1 de enero de 2015  y según  el Real Decreto 1071/2007 , de 27 de julio, por el que se regula el sistema geodésico de referencia oficial en España, "...toda la cartografía y bases de datos de información geográfica y cartográfica producida o actualizada por las Administraciones Públicas deberá compilarse y publicarse conforme a lo que se dispone en este real decreto” o lo que es lo mismos debe estar encontrarse en el ETRS89. Más información en esta página del IGN.

Para facilitar la transformación del datum ED50 a ETRS89, el Instituto Geográfico Nacional ha generado dos rejillas (Península y Baleares) en formato NTV2. Estas rejillas pueden ser descargadas desde la web del IGN.  Las rejillas pueden añadidas  usadas en la mayoría de los SIG pero también podemos utilizarla directamente usando GDAL con el comando gdalwarp. En la web de Mappingis podréis encontrar un artículo de su uso en QGIS.

Si queremos trasformar con gdalwarp una imagen en ED50 UTM30N (EPGS:23030) a ETRS89 (EPGS:25830), lo primero será descargar la correspondiente rejilla a nuestro ordenador. Tras el comando debemos añadir los SRC de entrada (-s_srs) y salida (-t_srs), pero en esta ocasión añadiendo los parámetros de la librería PROJ.4 e indicando la ruta de la rejilla del IGN tras la opción nadgrids.

Suponiendo que hemos almacenado la rejilla de la península en nuestro disco duro C: en la carpeta "rejillas", el comando quedaría así:

$ gdalwarp -s_srs "+init=epsg:23030 +nadgrids=C/:rejilla/PENR2009.gsb +wktext" -t_srs "+init=epsg:25830 +nadgrids=null +wktext" -of GTiff img23030.tif img25830.tif

Un poco de programación: procesando de múltiples imágenes con gdalwarp

Como he comentado al principio, todo este textazo tiene en primer lugar una función didáctica para mi. El segundo aspecto por el que merece utilizar los comando de estas librerías es el hecho de poder crear nuestros propios fragmentos de código o scripts que nos ayuden a mejorar nuestra vida..al menos desde el punto de vista informático. Creo que en esta búsqueda de la felicidad, podría ser interesante hacer un un pequeño script que por ejemplo convirtiera a formato GeoTiff y mejorara el tamaño de, ¿porqué no?, 200 ficheros.

Creamos un script utilizando el shell de Linux con la extensión sh (ej. ed50etrs89.sh). El script está pensado para ser utilizado desde el terminal de Linux. Si alguno se lo trabaja para Windows o Mac, y quiere que lo ponga en la entrada avisad por correo.

El código hace lo siguiente:

  • Localiza las imágenes con la extensión JPG que se encuentran en el mismo directorio del archivo  ed50etrs89.sh
  • Reproyecta  y comprime las imágenes a ETRS89 UTM 30N usando la rejilla del IGN localizada en en mismo directorio del script.
  • Salva las imágenes con el sufijo “_25830”  y en formato GeoTiff
  • Usa gdalinfo para generar una archivo de texto con los metadatos de las imágenes reproyectadas.

 

Si es necesario le asignamos, los correspondientes permisos de ejecución con chmod

$ chmod +x ed50etrs89.sh

Y a continuación, lo ejecutamos

$ ./ed50etrs89.sh

GDAL (http://www.gdal.org/) es una biblioteca de software para la lectura y escritura de formatos de datos geoespaciales publicada bajo la MIT License por la Fundación OSGeo. Con esta librería se pueden realizar multitud de operaciones de transformación y procesamiento sobre gran variedad de datos raster y vectoriales.

En la siguiente serie de entradas vamos a tratar algunos los comandos destinados a obtener información de un raster (gdalinfo), convertir archivos a otros formatos (gdal_translade) y cambiar el sistema de referencia de coordenadas  de imágenes referencias  (gdalwarp).

Para saber más sobre esta librería se pueden consultar los siguientes artículos.

¿Terminal o SIG?

La librería GDAL es utilizada por gran número de paquetes geomáticos (https://trac.osgeo.org/gdal/wiki/SoftwareUsingGdal) como por ejemplo QGIS, gvSIG o ESRI ArcGIS 9.2+. Desde los menús de cualquiera de estos clientes SIG podemos acceder a las funciones de GDAL utilizando los formularios diseñados para ello. Si podemos utilizar sin complicaciones las funciones desde un Sistema de Información Geográfica ¿qué sentido tiene utilizar los comandos?. Para mi, la primera razón es el simple hecho de conocer, aprender y saber a manejar distintas herramientas vinculadas con la información geográfica. El segundo motivo, y no por ello el menos interesante, funcional y divertido es el intentar dar soluciones más efectivas mediante programación a operaciones que se repitan o en las que estén incluidas multitud de ficheros.

Instalación de GDAL

Comenzamos con la instalación de la librería. En mi, en un equipo con Ubuntu he añadido repositorio de fuentes, actualizado e instalado a librería usando los siguientes comandos.

$ sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
$ sudo apt-get update
$ sudo apt-get install gdal-bin

Una vez instalado podremos hacer una verificación mediante este comando.

$ gdalinfo --version

Para equipos con Windows lo más recomendable es utilizar el instalador OSGeo4W, aunque también está la alternativa de FWTools.

Como no me quiero meter en follones con la instalación en Mac por no haberla probado, aquí un dejo un enlace que puede ayudar.

Usando gdalinfo  para obtener información sobre un archivo

Podemos utilizar el comando gdalinfo para obtener gran cantidad de información de un determinado fichero de imagen. Si queremos podemos también comprobar los tipos de formatos soportados por GDAL con el comando gdalinfo –formats.

Si nos encontramos en el directorio donde se encuentra el archivo, desde Linux sólo debemos abrir nuestro terminal, escribir el comando seguido del nombre y extensión del fichero. En el caso de estar en otro directorio debemos añadir la ruta y el nombre más extensión del raster.

$ gdalinfo mapa.tif

o

$ gdalinfo utm.tif /home/user/mapa.tif

Entre otros datos obtendremos: tipo del fichero, tamaño, SRC, tamaño del píxel, metadatos, coordenadas de las esquinas y centro o la información de las bandas.

Si volvemos a escribir el comando, pero esta vez seguido del símbolo mayor que (>) con un nombre y la extensión txt (ej. utm_metados.txt), se genera un fichero con los datos. Esta opción es bastante interesante si queremos incorporar un fichero adjunto de metadatos a nuestra imagen.

$ gdalinfo mapa.tif > mapa_metados.txt

Cambios de formato con gdal_translate

La utilidad básica de gdal_translate es la de hacer una copia de un fichero existente en otro formato de los reconocidos por la librería.

Por ejemplo, al georreferenciar una imagen que sólo requiera escalado y traslación, lo normal es que se generar un archivo  de tipo world file con los datos del píxel y la coordenada de la esquina. Estos ficheros se sombran de igual manera que la imagen y tienen una extensión acabada en la letra w. (pnw, .jpgw/.jpegw o .wld). Si utilizamos este sistema de multiarchivo la información geográfica no es intrínseca al archivo, por lo que no podremos saber el Sistema de Referencia de coordenadas en que se encuentra la imagen. Para solventar este escollo es preferible trabajar con archivos en formato GeoTiff.

Para convertir el archivo a formato GeoTiff usamos gdal_translate y añadimos el fichero de origen y el nombre del nuevo fichero con su extensión.

$ gdal_translate -of GTiff mapa.jpg mapa.tif


Una vez realiza la conversión es interesante comprobar la diferencia de datos vinculados al fichero usando de nuevo el comando gdalinfo.

Lo normal es que utilicemos la librería GDAL con imágenes georreferenciados, pero puede ser también útil para pasar de formato otros archivos sin coordenadas previas. Por ejemplo, podemos convertir un plano en PDF a un archivo TIF, establecer la resolución de la imagen de salida o incluso obtener una zona concreta del PDF usando  -projwin.

$ gdal_translate in.pdf out300.tif  --config GDAL_PDF_DPI 300 -projwin 792 144 9450 5600.

Después utilizando un SIG o la misma librería podríamos asignarles las coordenadas para su georreferenciación.

Opciones de compresión

Podemos ir más allá usando los distintas opciones del comando. y por ejemplo definir las opciones de configuración para aplicar compresión a fichero y reducir su tamaño. Usando el comando y la opción “-co” se podrá asignar al nuevo fichero un algoritmo de compresión que puede reducir el peso del nuevo.

$ gdal_translate -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co TILED=YES  mapa.jpg mapa.tif

No voy a entrar en las diferentes algoritmos de compresión de imágenes que podemos utilizar pero sí os dejo unos enlaces que son de gran interés. Sólo indicar que existen opciones que nos permiten reducir casi a un 50% el peso de un fichero sin perder calidad en la imagen.

Sobre mí

SIGdeletras es Patricio Soriano y Patricio Soriano es SIGdeletras. Trabajo el campo las Tecnologías de la Información Geográfica y especialmente su aplicación en el ámbito del Administración Pública y el Patrimonio Cultural...  ¿Quieres saber más sobre mí?

 

Buscar