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.

Como he venido comentando en las dos últimas entradas, en el repositorio de GitHub de SIGdeletras podéis encontrar el script dxf2gmlcatastro permite transformar un archivo de parcela catastral CAD con extensión DXF al formato GML establecido por Catastro para conseguir su validación gráfica en la Sede Electrónica del Catastro.

Tras verlo con algunos compañeros, presentarlo en las última reunión de Geoinquietos Córdoba y recibir varios correos de personas interesadas en usarlo, la conclusión a la que he llegado es la de que por muy útil que sea una herramienta si la gente no sabe utilizarla ya se está perdiendo el sentido básico que motiva su creación.

En esta línea, y poco después de publicar el script, Óscar Martínez de MásqueSIG escribió una estupenda entrada en la que explicaba cómo usar dxf2gmlcatastro en gvSIG. Siguiendo las reflexiones de Óscar, podemos decir que integrar código en un SIG como gvSIG, o en nuestro caso QGIS, permite:

  • Al poder integrar el código en un SIG que trabaja con Python y GDAL, el uso inicial de la librería es más fácil.
  • No estamos limitados a la instalación en un sistema operativo concreto ya tanto gvSIG como QGIS son SIG multiplataforma.
  • Podemos mejorarlo usando las funcionalidades que nos ofrece el propio SIG (cuadros de dialogo, interfaz visual...).
  • Se podría integrar en la barra de herramientas o incluso convertirlo en una extensión.

En conclusión y como dice Óscar …”que pueda llegar a más gente, que al final es para lo que lo hacemos”.

Definir la variable PYTHONPATH

La variable PYTHONPATH es utilizada en QGIS para acceder a los módulos de Python. Esta variable ya se define durante la instalación apuntando a una carpeta denominada Python dentro del directorio de instalación del programa, o en la carpeta .qgis2 en Linux. A pesar de ello, vamos a añadir una ruta nueva más accesible donde vamos a guardar nuestro módulo.

Los pasos a seguir son los siguientes:

  • Abrir QGIS.
  • Ir al menú Configuración>Opciones.
  • En la pestaña Sistema, buscamos el apartado Entorno.
  • Activamos la opción “Usar variables personalizadas…”
  • Pinchamos en el botón “Añadir” y definimos la variable con las siguientes opciones:
    • Aplicar: “Poner a continuación”
    • Variable: PYTHONPATH
    • Valor: Carpeta donde vamos a guardar nuestros archivos Python (ej. C:\scriptsqgis)
  • Pinchamos en Aceptar y reiniciamos QGIS.

 

Usar dxf2gmlcatastro en QGIS.

Lo primero que debemos hacer es descarga el código, o hacer un git clone,  de dxf2gmlcatastro desde el repositorio de GitHub y copiar los archivos en la carpeta definida en PYTHONPATH (ej. C:\scriptsQGIS).

Lo primero es comprobar que QGIS accede a la nueva ruta que hemos añadido a PYTHONPATH. Una vez ejecutado de nuevo QGIS,  vamos a abrir la consola de Python integrada (Complementos > Consola de Python) y escribimos import dxf2gmlcatastro y después pulsamos "Intro". Si hemos seguido correctamente los pasos anteriores no nos debería devolver ningún error.

A continuación vamos a ejecutar la función crea_gml de dxf2gmlcatastro que convertirá nuestro archivo DXF al GML de Catastro. Dentro de la función tendremos de añadir los argumentos que nos indiquen:

  • Dónde se encuentra el archivo DXF (ej. 'C:\carpeta\archivoparcela.dxf')
  • El lugar y el nombre del GML(ej. 'C:\carpeta\gmlcatastro.gml')
  • El código EPSG del Sistema de Referencia de Coordenadas del DXF. (Los SRC admitidos son 25828, 25829, 25830 y 25831)
dxf2gmlcatastro.crea_gml('C:\carpeta\archivoparcela.dxf', 'C:\carpeta\gmlcatastro.gml', '25830')

Si todo ha ido correcto podemos cargar el nuevo GML en QGIS y comprobar que se ha generado correctamente y que incluye los atributos según el modelo de Inspire.

Crear un script con PyQGIS personalizado

Una vez que disponemos de dxf2gmlcatastro podemos crear otros fragmentos de código para integrarlo en QGIS. Aquí os dejo un ejemplo que carga el GML en QGIS.

Para crear este archivo hemos usado el editor de código de QGIS y una vez terminado lo hemos salvado en nuestro equipo con el nombre catastroQGIS.py. Sólo deberemos cargar el archivo en el terminal de Python de QGIS y cambiar la ruta y nombre de los archivos para utilizarlo cada vez que queramos.


El archivo catastroQGIS.py se encuentra en la carpeta ejemplo del repositorio GitHub

 

Dentro del programa de la 11ª reunión de Geoinquietos Córdoba voy a realizar de mini taller donde explicar cómo se usa el script Pyhton dxf2gmlcatastro que estoy generando para convertir un archivo CAD DXF al formato GML definido por tras la Resolución conjunta Catastro-Registro.

Desde hace unos meses he empezado a aprender Python y su aplicación temas geo. Este ejercicio de programación me está sirviendo bastante para ir adquiriendo los conocimientos básicos de este lenguaje. Además de las lecturas y pruebas, la ayuda de Marcos Ortega de Indavelopers está siendo crucial.

Un rápido sondeo ha puesto de manifiesto que la mayoría de las personas que asistente a las reuniones de geoinquietos usan Windows. Aunque todo la programación está hecha en Linux (Ubuntu), he montado una máquina virtual con Virtual Box para ver los pasos a dar para poder ejecutar el script en este sistema operativo. Existe una manera más rápida de configurar GDAL para Python que es instalando OSGeo4W y usar la shell incorporada. Pero como siempre es un buen momento para aprender he decido tomar el camino más largo y aprovechar para escribir esta entrada.

Instalación de Python

Para instalar Python es necesario entrar en el apartado de descargas de la página de Pyhton) y bajarse la versión del lenguaje que se quiera instalar según el sistema operativo. Para esta guía hemos utilizado Python 3.4.4 - 2015-12-21 Windows x86 MSI installer.

Una vez descargada, ejecutamos el archivo e instalamos según la configuración por defecto. El único aspecto al que en principio hay que atender es a marcar, si no lo está, la casilla que añade python.exe a nuestro PATH.

Terminada la instalación, accedemos a Sistema>Configuración avanzada del sistema y pinchamos en el botón Variables del entorno. Desde aquí, podemos comprobar que se ha añadido la ruta a la carpeta de Python, en nuestro caso "C:\Python34", para la variable PATH. Editamos esta variable y añadimos también las rutas C:\Python34\Lib\site-packages\ y C:\Python34\Scripts.

Instalación de GDAL

Para poder usar la biblioteca geoespacial GDAL con Python es necesario tener instalados los binarios con anterioridad. Esta parte de la entrada recoge la guía en inglés del siguiente enlace.

Lo primero es saber la versión del compilador que estamos utilizando. Para ello buscamos en nuestros programas IDLE (Python GUI) dentro de la carpeta de Python y lo ejecutamos. En la información de inicial podemos apreciar la versión del compilador (MSC v.1600). En este mismo texto obtendremos la arquitectura de nuestra máquina (32 o 64 bits).

Recopilados estos datos, entramos en la web www.gisinternals.com para acceder a los archivos de GDAL según nuestra configuración (ej. release-1600-gdal-1-11-3-mapserver-6-4-2). Una vez en la página de descarga, localizaremos el enlace del core de GDAL (Generic installer for the GDAL core components), descargamos e instalamos. Para nuestro equipo el archivo descargado ha sido gdal-111-1600-core.msi

Tras la instalación, tendremos que volver a editar de nuevo las variables del entorno y añadiremos a PATH, precedida de punto y coma, la dirección de la carpeta de GDAL de nuestro equipo (ej. ;C:\Program Files (x86)\GDAL).

A continuación debemos añadir dos variables nuevas:

  • GDAL_DATA que apunte a la carpeta gdal-data (C:\Program Files (x86)\GDAL\gdal-data)
  • GDAL_DRIVER_PATH que apunte a la carpeta gdalplugins (C:\Program Files (x86)\GDAL\gdalplugins)

Para comprobar que hemos realizado correctamente la instalación accederemos al Símbolo del sistema como administrador y escribiremos gdalinfo --version. Si todo es correcto obtrendremos la versión de GDAL instalada

Instalación de Python bindings

Para poder utilizar GDAL con Python vamos a instalar los bindings que se encuentran en la misma web que hemos descargado los archivos binarios de GDAL. Un binding es una adaptación de una biblioteca para ser usada en un lenguaje de programación distinto de aquel en el que ha sido escrita.

Como partimos de una instalación de Python 3, descagaremos e instalamos el archivo GDAL-1.11.3.win32-py3.3.msi

Para terminar, vamos a acceder de nuevo la GUI de Python IDLE y escribiremos el siguiente código importando GDAL en Python

import gdal

Si no nos devuelve ningún mensaje de error, hemos terminado nuestra instalación.

Si queréis empezar trastear con Python y GDAL os recomiendo seguir los ejemplos de Python GDAL/OGR Cookbook.. También puede puede revisar el curso "Geoprocessing with Python using Open Source GIS" de la Utah State University.

Estas últimas semanas conocidos y amigos que se dedican a la arquitectura, la tasación inmobiliaria me han hecho consultas sobre la generación del archivo GML que debe añadirse a la descripción grafíca que describe informáticamente las parcelas catastrales.

Laboralmente no he tenido que realizar trabajos vinculados con Catastro, por para enterarme un poco de que iba el tema tuve que leer algunos documentos técnicos colgados en la web de la Dirección General de Catastro, blog especializados y algún que otro foro.

Según he podido, leer cuando la certificación catastral descriptiva y gráfica no coincide con la realidad es necesario elaborar una representación gráfica alternativa a la que se debe añadir un XML con contenido geográfico que no es otro que el famoso GML.

Este mismo archivo es también accesible, como adjunto, en las certificaciones catastrales descriptivas y gráficas, desde la consulta interactiva de un bien inmueble y desde servicios web WFS.

Mi primera respuesta a fue que se usara un SIG de escritorio para pasar el DFX al GML pero parece que no es tan sencilla el esquema no está disponible aún en la mayoría de los SIG (parece que en la última versión de gvSIG 2.3 está disponible). El formato de parcela catastral debe cumplir el estándar INSPIRE cadastral parcel definido en INSPIRE Data Specification on Cadastral Parcels - Guidelines version 3.0.1. L

¿Cómo generar el GML?

La Dirección General del Catastro ha colgado de su web una serie de guías técnicas en las que se explica entre otros cómo generar un GML de Parcela Catastral. En el mismo apartado hay un par de ejemplos con explicaciones comentadas de archivos GML validados para parcela catastral y para edificio.

En la guía se definen los siguientes pasos:

  • Paso 1: Descarga de un fichero DXF de la sede electrónica del Catastro conteniendo la cartografía de la zona en la que se desea intervenir.
  • Paso 2: Edición del fichero con AUTOCAD (...u otro programa CAD o incluso SIG) y generación de un nuevo fichero DXF.
  • Paso 3: Modificación manual del fichero DXF generado y obtención de coordenadas.
  • Paso 4: Generación del fichero GML para adaptarlo al formato de parcela catastral.
  • Paso 5: Validación del fichero GML en la Sede Electrónica del Catastro.

Todo este proceso es bastante "artesano" y lo primero que se me ocurrió es que algunos de estos pasos (3 y 4) se podría automatizar con un poco de programación con el fin de ganar tiempo y evitar errores que se puedan generar por la edición manual. Tras algunas lectura y documentación pensé en generar un pequeño código en Python usando la librería GDAL y compartirlo en GitHub para que pueda ser utilizado y espero que mejorado.

¿Cómo funciona dxf2gmlcatastro.py?

Descargar dxf2gmlcatastro

En primer lugar debemos descargar/descomprimir o hacer un git clone de los archivos python disponibles en el repositorio dxf2gmlcatastro en GitHub.

Instalar Python y la libraría GDAL.

En Ubuntu Python viene instalado por defecto. De todas formas si queremos comprobarlo y ver la versión instalada sólo hay que ejecutar desde terminal el siguiente comando

$ python --version
Python 2.7.6

La librería GDAL es la que se encargará de todas las operaciones de acceso y lectura del archivo DXF. Para usarla con Python he instalado python-gdal.

$ sudo apt-get install python-gdal

Todo el trabajo se ha realizado en Ubuntu. Para Windows, tras instalar Python, podría utilizarse el administrador de paquetes Python Pip.

Ejecutar dxf2gmlcatastro.py (v.2*)

* Agradecer al inestimable ayuda y PR de Marcos Ortega de Indavelopers

Desde la consola ejecutamos el comando dxf2gmlcatastro.py indicando a continuación la ubicación del archivo DXF, el nombre del nuevo archivo GML y el código EPSG del Sistema de Referencia de Coordendas del archivo DXF.

$ pyhon dxf2gmlcatastro.py archivocad.dxf archivogmlcatasro.gml 25830

Toda la información junto a archivos de ejemplo puede consultarse en el repositorio GitHub.

2do

  • Permitir elegir el SRC del GML. (v.2)
  • Investigar qué es el "Identificativo local de la parcela" y se si debería solicitar al ejecutar el script.
  • Poder elegir el archivo DXF a transformar. (v.2)
  • Poder elegir el archivo GML a crear.(v.2)
  • Generar un GML de varias parcelas catastrales.
  • Crear un script para edificio.
  • Probarlo en otros Sistemas Operativos. (MacOS)
  • Probarlo en QGIS

Referencias

Durante las pasadas fiestas navideñas, tuve tenido la oportunidad de realizar el curso "Introducción a Scripting en gvSIG 2.1" organizado por la Asociación gvSIG mediante la Plataforma de Capacitación a Distancia gvSIG-Training.

El curso tenía como objetivo “el de dar a conocer un poco más sobre el potencial de la programación geoespacial, la posibilidad de crear nuevas herramientas, nuevos geoprocesos o análisis de datos, que nos aumentarán la potencia de gvSIG adaptándose a nuestras necesidades. También la automatización de tareas, que nos podrían generar un ahorro de tiempo y de trabajo considerable.”

El lenguaje de programación manejado ha sido Python, utilizado mediante las herramientas scripting incluidas en la última versión gvSIG la 2.1, en mi caso, corriendo sobre Xubuntu.

Mi experiencia, partiendo que mis conocimientos sobre Python eran nulos, ha sido bastante satisfactoria. Tanto ha sido así, que aprovechando la opción que ofrecía la plataforma, me animé completar el curso para obtener el Certificado de Aprovechamiento, correspondiente a 30 créditos del programa de Certificación de gvSIG.

Además de lo aprendido, uno de los puntos fuertes del curso ha sido el buen hacer del tutor Óscar Martínez (@masquesig), responsable de la web Másquesig en la que se puede encontrar bastante información de calidad sobre scripting en gvSIG y QGIS. Óscar está también disfrutando de la beca Google Summer of Code 2014 desarrollando el proyecto sobre “Acceso a los geoprocesos desde el módulo de scripting en gvSIG 2.x″ que permite la ejecución de los algoritmos de Sextante y gvSIG desde la consola de Jython o desde el Scripting Composer, ubicados ambos dentro del módulo de Scripting.

Además de la página web de Óscar Martínez, se puede consultar esta entrada resumen de su participación en las las 10as Jornadas Internacionales gvSIG el pasado Diciembre en Valencia con el Taller de Scripting en las 10as Jornadas Internacionales de gvSIG.

Como parte de los trabajos para conseguir el certificado, además del abono de las tasas, se incluye el desarrollo de un proyecto final, en que se apliquen los conocimientos aprendidos durante el curso. Una vez revisado el código por parte del tutor, los proyectos presentados serán incorporados en el apartado Scripts de la web gvSIG Outreach quedando así disponibles para toda la comunidad gvSIG.

Mi proyecto, y gracias a las orientaciones del profesor y a las aportaciones de los compañeros del curso, consistió en la creación de un script que generara un archivo en formato txt, con información de la capa (nombre, geometría, SRC, envolvente) para cada una de las capas (shape) cargadas en la vista activa y lo guarda en un directorio seleccionado. De forma opcional y utilizando la librería gvpy desarrollada por Óscar Martínez hace una copia de las capas en la vista utilizando el comando export en el mismo directorio.

Seguro que profundizando un poco más en el material entregado en el curso podré mejorarlo. He apuntado algunas ideas que podían mejorarlo y que podrían convertirlo en script de más utilidad para trabajos que requieran hacer copias de grandes conjuntos de datos incorporando una información o metadatos básicos. Aquí dejo algunas:

  • Comprobar que la librería gvpy está instalada y en el caso contrario devolver un mensaje de aviso. Esta parte ya está solucionada incluyendo estas líneas de código

    try:
        import gvpy
    except:
        commonsdialog.msgbox("ERROR: cannot find gvpy module. Please visit https://github.com/oscar9/gvpy for more information", "Error",commonsdialog.FORBIDEN)
        sys.exit()
    
  • Si la carpeta elegida tiene datos, permitir seleccionar otra o incluso borrar los existentes. Para ello tengo que buscar ejemplos utilizando la librería "os" de Python

  • La generación del txt y la exportación se hacen sobre todas las capas cargadas en la vista activa. Quizás sería interesante crear un dialog en el listado de capas para poder realizar una selección de las capas con las que trabajaría el script.

  • Una opción interesante sería la de incluir la opción de crear un archivo comprimido con todos los ficheros creados.

  • Se podría incluir la opción de transformar otros Sistema de Referencia de Coordenadas o incluso a otros formatos (ej. GeoJSON)

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