Desde el año 2005 se ha desarrollado el proyecto de investigación aplicada denominado: “Ramón Llull”, para la investigación del tráfico de información sobre activos inmuebles, patrocinado por el Colegio Notarial de Valencia y el Consejo General del Notariado, bajo la dirección de Antonio Jiménez Clar y Alfredo Ramón Morte, profesores de las áreas de Derecho Civil y Análisis Geográfico Regional, al frente de un grupo de investigadores del Laboratorio de Geomática del Instituto Interuniversitario de Geografía de la Universidad de Alicante.

Entre los objetivos conseguidos, destacaremos la creación de un conjunto de aplicaciones informáticas capaces de facilitar a los profesionales, técnicos, juristas y en especial, a los notarios, herramientas basadas en Sistemas de Información Geográfica para la gestión de información geográfica catastral, muy amigables pero con altas prestaciones en cuanto su exactitud geométrica y fiabilidad, acordes con los requerimientos de la propia gestión cartográfica de la Dirección General de Catastro, cuya colaboración ha sido decisiva, más allá de los diversos servicios telemáticos que ya presta la Oficina Virtual de Catastro.

Estos desarrollos han conseguido que profesionales sin perfil técnico sean capaces de comprobar la geometría de las parcelas catastrales implicadas en diversos tipos de operaciones de alteración básicas (agregación, segregación, división o deslinde), obra nueva y división horizontal, mediante el uso “transparente” de matrices de comparación topológica y el consiguiente diagnóstico automático sobre pertinencia geométrica / cartográfica de la situación a la que se quiere otorgar escritura pública. La comparación entre la finca resultante de la alteración propuesta por el interesado / propietario y la parcela catastral de referencia, incorpora en el proceso los criterios de calidad y exactitud exigidos por la Dirección General de Catastro para este tipo de operaciones, tomando siempre la cartografía catastral como base de cada operación y como destino final del resultado, para su fácil actualización.

De esta manera, se facilita la comunicación entre el interesado, el técnico proyectista de la alteración y el notario, sobre una plataforma informática de gestión de la información que para poder resumir sus prestaciones, destacaremos lo siguiente: automatiza los flujos de intercambio de datos “geo-gráficos” y alfanuméricos entre usuarios y formatos, los procesos de cotejo cartográfico y la generación automática de los componentes necesarios para construir la propuesta de alteración en el formato digital de intercambio de información de Catastro (FXCC), los formularios en formato de documentos tipo: 902 y 903, según el caso, así como las salidas gráficas a pantalla/papel de la Cartografía en proyección UTM y en ETRS89 o ED50 (según se va asumiendo el cambio de elipsoide en el SIGCA), incorporando el plano a la propia escritura, dando mayor seguridad jurídica a la alteración de la información sobre la propiedad y generando un documento muy útil en la mejora de la descripción de la finca resultante para su correspondiente asiento en el Registro de la Propiedad.

COMUNICACIÓN y SENCILLEZ en las aplicaciones cliente

Las aplicaciones cliente de escritorio se han desarrollado en Lenguaje C#, basadas en componentes de código abierto como framework de rendering SharpMap o la librería de análisis topológico / geométrico NTS Net Topology Suite. Mediante una solicitud HTTP a la Oficina Virtual de Catastro, se descarga en local la descripción alfanumérica y gráfica en formato FXCC (en el futuro, la necesidad del proceso electrónico hará que la propia Oficina Virtual de Catastro genere una “certificación electrónica” que supere esta dualidad de formas de información, mientras tanto, el formato de CAD va por un lado y el de texto por otro…). A partir de las parcelas demandadas por los usuarios, en el ámbito organizativo de un expediente de alteración catastral, la aplicación convierte el CAD en geometrías OGC, para su gestión topológica y su visualización mediante capas gráficas, evaluando de forma automática las relaciones espaciales entre las parcelas catastrales y la información técnica del proyecto de alteración, con el uso de predicados y matrices de cálculo topológico acordes con los criterios de calidad cartográfica dictados por el propio Catastro.

La clave del éxito de los programas del “Ramón Llull” reside en la capacidad de comunicación segura y efectiva entre personal con perfiles de formación jurídica y técnica, desde un entorno de trabajo intuitivo y en ello, no debemos quitar protagonismo a los servicios de mapa de la Infraestructura de Datos Espacial de España (IDEE). De hecho, la incorporación de los servicios WMS de la propia Cartografía de Catastro o del Plan Nacional de Ortofotogrametría Aérea (PNOA) cumplen el papel de fondo cartográfico de referencia básica sobre el que revisar la propuesta técnica de alteración de información catastral por parte de todos los implicados.

No obstante, los servicios de descarga de datos (como los límites de las parcelas catastrales implicadas en el proyecto), la conversión de las geometrías de CAD a formatos de datos GIS, la detección de errores topológicos conflictivos para la incorporación de la propuesta en el SIGCA de Catastro y la posibilidad de revisar el proyecto o expediente de alteración desde un entorno de trabajo gráfico con herramientas diseñadas específicamente para estos fines, convierten a la herramienta en un auténtico “banco de pruebas” previo al otorgamiento de su escritura pública y Declaración Catastral en la Oficina Virtual de Catastro. En este sentido, ha sido muy interesante el retorno de información de los primeros grupos de usuarios y el apoyo de los colegios profesionales de Arquitectos, Arquitectos Técnicos, Topógrafos o Geógrafos, así como el del grupo de trabajo sobre “Delimitación de la propiedad en la Comunidad Valenciana” (promovido por la Universidad Politécnica de Valencia) integrado por un equipo con una visión multidisciplinar (Derecho Civil, Ingeniería Cartográfica, Topógrafos, Geógrafos, Notarios, Registradores de la Propiedad…)

SEGURIDAD y ROBUSTEZ en las aplicaciones de servidor

La gestión de los expedientes de alteración, el acceso con reconocimiento de certificados y firmas electrónicas, así como los permisos de cada usuario, según su perfil profesional, quedan administrados desde un Servidor con PostgreSQL / PostGIS como backend, mientras se desarrolla el trabajo de los distintos agentes y hasta que se finaliza el expediente, si cabe, con el otorgamiento de una escritura y la remisión de una declaración catastral.

El conjunto de aplicaciones del “Ramón Llull” aseguran su robustez, competencia y seguridad a través de un entorno de trabajo con acceso certificado desde un cliente local, un servicio WEB de asistencia técnica – formación on-line – actualización de versiones y un Servidor de Datos dedicado a gestionar todos los expedientes de alteración catastral elaborados por las aplicaciones cliente, para dar la necesaria seguridad al trámite y a los procesos de trabajo de técnicos o juristas. Gracias a ello, el proyecto consigue dar respuesta a las necesidades planteadas desde la reciente entrada en vigor de la Ley de Economía Sostenible cuyo articulado supone un cambio para la Ley del Catastro Inmobiliario, concediendo un gran papel al uso de la referencia catastral como código clave de identificación del inmueble, promueve la remisión electrónica de la información a Catastro, la aceleración de los trámites en la resolución de alteraciones o el intercambio de información digital con registradores y notarios.

La NECESIDAD de ESTÁNDARES y CONOCIMIENTO ABIERTO

La vocación universitaria de difusión del conocimiento, asumida en las líneas de trabajo del propio Laboratorio de Geomática, ha determinado que las aplicaciones del proyecto constituyen el germen de un proyecto de software libre. Desde junio de 2012 se ha habilitado un lugar WEB para la descarga de versiones, documentación y asistencia técnica para los usuarios, dentro de una fase experimental de puesta en explotación. La remisión de los primeros expedientes resueltos está sirviendo para ensayar los protocolos de entrega en la Oficina Virtual, en colaboración con las Delegaciones Territoriales de Catastro, y para estudiar futuros métodos electrónicos de descarga y remisión de información más ágiles y eficaces (generación de esquemas XML o generación automática de nuevas referencias catastrales). Todo ello podrá ser posible gracias a la firma electrónica y la seguridad del servidor de expedientes, aunque por ahora, se sigue dependiendo mucho de formatos de CAD y de texto, así como del proceso administrativo analógico, por lo que la herramienta está preparada para generar la información electrónica en estos formatos digitales (FXCC) y para la impresión en PDF de los formularios tipo.
Todo ello sin perder de vista el requisito de que la base operativa sea open source y software libre, con un grupo de desarrollo e investigación detrás en continua evolución, actualizando e implementando herramientas que permitan hacer uso de otro tipo de servicios interoperables, como la capacidad de proceso de los servicios WPS que podrían suponer un futuro muy distinto, incluso la sustitución de la aplicación cliente por un simple navegador y el alojamiento directo de la información de un expediente concluido y de la escritura pública en la Oficina Virtual de Catastro, por ejemplo.
En la actualidad, también se está llevando a cabo una iniciativa programada de formación de usuarios, dentro de las tareas de divulgación del proyecto. La primera edición de estos cursos sobre las aplicaciones del proyecto Ramón Llull, celebrada en la Universidad de Alicante durante el pasado mes de marzo, tuvo una gran acogida y resultó de gran utilidad para el alumnado y para el propio equipo de desarrollo, en espera de las próximas ediciones que tendrán lugar en el Colegio Notarial de Valencia y en la Sede Universitaria de Benissa, los próximos meses de mayo, junio y julio.
Para más información sobre el proyecto Ramón Llull, manuales, materiales de los cursos ya celebrados, así como las próximas ediciones de cursos, se pueden consultar la web: http://www.arsmagna.es o contactar mediante correo electrónico: proyecto.ramonllull@ua.es

Geomondrian es un servidor OLAP con soporte espacial, lo que permite utilizar la componente espacial a la hora de realizar consultas en el almacén.

Este artículo describe los pasos necesarios que hay que realizar para instalar y probar un servidor Geomondrian  en un computador. El objetivo no es tener un sistema en producción, sino que se pueda instalar fácilmente el sistema y se pueda “jugar” un poco para ver la potencia que puede llegar a tener.

Instalación

Partimos de una máquina con una Ubuntu 11.10 de 64 bits. Instalamos el Tomcat 6 en Ubuntu. Para ello ejecutamos el siguiente comando:

sudo apt-get install tomcat6 tomcat6-admin

Editamos el fichero /etc/tomcat6/tomcat-users.xml y añadimos el usuario admin con contraseña admin:

 <role rolename="manager"/>
 <role rolename="admin"/>
 <user username="admin" password="admin" roles="manager,admin"/>

Reiniciamos el tomcat:

sudo /etc/init.d/tomcat6 restart

En este punto ya tenemos una instancia de Tomcat corriendo en nuestra máquina junto con la aplicación de administración que se va a utilizar para instalar el Geomondrian.

Descargamos el Geomondrian 1.0 desde su página web. Se trata de un archivo con extensión war, que puede ser fácilmente desplegado en Tomcat.

Abrimos un navegador y vamos a la página de administración de Tomcat (http://localhost:8080/manager/html). La primera vez que accedemos tenemos que escribir el nombre de usuario y la contraseña que hemos creado previamente.

Buscamos la opción WAR file to deploy y seleccionamos el archivo geomondrian.war que nos hemos descargado previamente. A continuación hacemos click en deploy.

Si todo ha ido bien aparecerá una nueva entrada en la lista de aplicaciones que tiene el Tomcat instaladas:

Hacemos click en el enlace o bien escribimos la url http://localhost:8080/geomondrian/ y ya podemos acceder a la ventana principal del Geomondrian.

Por defecto la aplicación incluye una base de datos embebida que contiene algunos ejemplos sin información espacial. En la web de Geomondrian hay un script de creación de una base de datos con información espacial pero se trata de una base de datos Postgis, por lo que tenemos que instalar previamente Postgres y Postgis.

Para ello ejecutamos desde una consola el siguiente comando:

 sudo apt-get -install postgresql-9.1 postgresql-contrib-9.1 postgis postgresql-9.1-postgis
 

Cambiamos el password del usuario postgres a postgres.

 sudo -u postgres psql
 postgres=# ALTER USER postgres WITH ENCRYPTED password 'postgres';
 

A continuación vamos a crear una base de datos plantilla y vamos a instalar Postgis en ella.

sudo su postgres
createdb postgistemplate
psql -d postgistemplate -f /usr/share/postgresql/9.1/contrib/postgis-1.5/postgis.sql
psql -d postgistemplate -f /usr/share/postgresql/9.1/contrib/postgis-1.5/spatial_ref_sys.sql

Una vez que ya tenemos la base de datos plantilla con todas las funciones de Postgis instaladas vamos a crear una base de datos a partir de esta plantilla que será usada por el Geomondrian.

 createdb -T postgistemplate simple_geofoodmart
 

Ahora vamos a instalar el script creación de la base de datos OLAP que hemos tenido que descargar previamente desde la web de Geomondrian. Una vez descargado, vamos al directorio donde lo hemos descargado y ejecutamos lo siguiente:

 psql -h localhost -U postgres -W -d simple_geofoodmart -f simple_geofoodmart.sql
 

A continuación tenemos que editar la cadena de conexión para que el Geomondrian pueda conectar con la base de datos espacial. Para ello abrimos el fichero /var/lib/tomcat6/webapps/geomondrian/WEB-INF/web.xml, buscamos el servlet GeoMDXQueryServlet y editamos la propiedad conectString:

 <servlet>
 <servlet-name>GeoMDXQueryServlet</servlet-name>
 <servlet-class>
 mondrian.web.servlet.MdxQueryServlet
 </servlet-class>
 <init-param>
 <param-name>connectString</param-name>
 <param-value>Provider=mondrian;Jdbc=jdbc:postgresql_postGIS://localhost/simple_geofoodmart?user=postgres&amp;password=postgres;Catalog=/WEB-INF/queries/simple_geofoodmart.xml;JdbcDrivers=org.postgis.DriverWrapper;</param-value>
 </init-param>
 </servlet>
 

Llegados a este punto ya tenemos un servidor SOLAP en funcionamiento que está conectado con una base de datos con soporte espacial.

Estructura de la base de datos

Antes de empezar a probar la base de datos, hay que entenderla. Existe poca información en la web sobre la estructura de la misma, pero se puede consultar el xml que utiliza el Geomondrian para explotarla que se encuentra en el fichero /var/lib/tomcat6/webapps/geomondrian/WEB-INF/queries/simple_geofoodmart.xml.

Se trata de una base de datos de ventas de supermercados que tiene un único cubo Ventas (Sales) para el que se han definido 5 dimensiones: Tienda (Store), Producto (Product), Promociones (Promotions), Clientes (Customers) y Tiempo (Time). El objetivo es estudiar las ventas en función de estas dimensiones:

El cubo ventas, tiene 3 medidas que son Ventas (Store Sales), Coste (Store Cost) y Unidades (Unit Sales) que son las variables que van a ser estudiadas.

La dimensión Tienda establece una jerarquía espacial de objetos (Tienda -> Ciudad ->Estado -> País) que tiene un atributo de tipo geométrico. Este atributo representa el polígono que describe el contorno del objeto y será estudiado en un apartado posterior.

Consulta de la base de datos

Desde la ventana principal del Geomondrian, hay que seleccionar la opción JPivot pivot table with simple_geofoodmart spatial cube para acceder a la consola de administración de  Geomondrian que contiene la base de datos espacial.

Lo primero que se muestra es una ventana en la que se han seleccionado dos dimensiones (Store y Product) y en la que se muestran las tres medidas que hay en el cubo.

Seleccionando el primer botón de la derecha se abre el navegador de OLAP que sirve para configurar el cubo.

En esta ventana se configuran las columnas y las filas que se quieren mostrar y los filtros que hay que aplicar en alguna de las dimensiones. Los iconos que hay delante de las dimensiones las mueven de un lugar a otro. Por ejemplo, si queremos realizar un análisis por cliente y queremos que aparezca como una columna en la tabla, simplemente seleccionaremos el icono de columna que aparece delante de Customers.

En la figura se puede observar como en Time aparece un  Year=1997, lo que indica que se están filtrando registros que no sean del año 1997. Se puede pulsar en en nombre de la dimensión y se abrirá una nueva ventana donde se pueden cambiar las condiciones del filtro.

Sobre la tabla que muestra los resultados se pueden realizar operaciones de drill down y drill up de forma que cada vez vemos más o menos detalle. La siguiente imagen se ha generado haciendo operaciones de drill down hasta que se han visualizado las ventas de bebida de la tienda 15 de Seattle.

En la imagen se puede ver como no hay datos de ejemplo para las tiendas de México y Canadá.

El lenguaje MDX

Otra opción interesante es la de MDX, que se puede configurar seleccionado el segundo botón de la derecha de la barra de menú. Las expresiones multidimensionales (MDX es el acrónimo de MultiDimensional eXpressions) es un lenguaje de consulta para bases de datos multidimensionales sobre cubos OLAP, se utiliza en Business Intelligence para generar reportes para la toma de decisiones basados en datos históricos, con la posibilidad de cambiar la estructura, o permitiendo rotar el cubo (wikipedia).

Conociendo la sintaxis de MDX, se pueden crear cubos a medida sin necesidad de utilizar el interfaz de usuario. Además, existe una equivalencia entre la consulta que hay en la ventana de MDX y el interfaz gráfico de forma que si se modifica uno, el otro se actualiza inmediatamente.

La siguiente imagen muestra un ejemplo de consulta MDX desde el interfaz web de la aplicación. El interfaz permite la edición manual de la consulta.

Gráficas

También se pueden añadir gráficas para una consulta determinada. Las gráficas son interactivas, de forma que a medida que el usuario navega haciendo operaciones de drill up y del drill down, la gráfica se refresca con las dimensiones seleccionadas.

La siguiente imagen muestra un ejemplo de ello. Se han ido haciendo operaciones de drill down y en la gráfica han ido apareciendo nuevas columnas dinámicamente.

Existen varios tipos de gráficas que el usuario puede configurar, como la gráfica de barras, la de tartas, etc.

El soporte espacial

Todo lo que se ha comentado hasta este apartado es tan válido para el Geomondrian como para el Mondrian ya que, hasta este momento, no se ha hecho ninguna mención al soporte espacial en los datos.

En el apartado anterior se comentó la importancia del lenguaje MDX y su mapeo con el interfaz gráfico de la aplicación. En realidad, la diferencia que aporta Geomondrian frente al Mondrian es que extiende el lenguaje MDX para que se puedan utilizar operadores espaciales. A este nuevo lenguaje lo llama GeoMDX.

Ya hemos visto que existen algunas tablas de la base de datos que tienen un atributo de tipo geometría. Ese tipo de datos no existe en Postgres, y es la extensión espacial Postgis la que lo añade. La siguiente imagen muestra una consulta sobre la tabla de ciudades realizada con el pgAdmin.

Se puede observar como existe una columna geometry que tiene un valor binario que no representa nada. Si se intenta utilizar un atributo de este tipo en un Mondrian lo que puede ocurrir es que o bien no procese el campo al no reconocer su tipo o bien lo muestre con su valor en binario que, para un observador humano, no dice nada.

A continuación vamos a intentar visualizar el mismo atributo desde el Geomondrian. Para ello vamos a seleccionar que nos muestre los atributos de las tablas de dimensiones y vamos a hacer drill down hasta llegar al nivel de ciudad. La siguiente figura muestra el resultado de San Diego y San Francisco.

Puede observarse que no aparece un valor binario en la columna de la geometría, sino que aparece el valor de un multipolígono, que es una de las primitivas de dibujo que soporta Postgis. Además se puede ver cómo se pueden leer fácilmente las coordenadas de los puntos que forman el polígono que se utiliza para describir a la ciudad.

Este formato de mostrar las geometrías se conoce como WKT  y es un estándar utilizado por muchas aplicaciones GIS. Con una geometría en WKT se podría abrir una aplicación GIS y visualizar de forma gráfica las geometrías obtenidas con el Geomeondrian.

En cuando a visualización ya hemos visto como Geomeondrian es capaz de convertir atributos geométricos en representaciones textuales legibles por un ser humano. Pero la extensión que se hace sobre el MDX permite además poder utilizar operadores espaciales en la consulta de la base de datos.

La figura siguiente muestra un ejemplo de consulta para mostrar las ventas en las tiendas de estados unidos y su área de venta en kilómetros cuadrados. Para ello se utilizan algunas funciones espaciales sobre un campo de tipo geometría.

La función ST_textunderscore Area devuelve el área de un polígono. La función ST_textunderscore Transform hace una conversión de sistema de coordenadas desde el sistema de coordenadas origen a un sistema de coordenadas destino que es el que necesita la función del área para operar. La función ST_textunderscore UnionAgg realiza la unión de varios campos de tipo geometría en una única geometría.

Conclusiones

El Geomondrian es un servidor SOLAP que tiene un proceso de  instalación muy sencillo (simplemente hay que desplegar un war en el Tomcat).

La aplicación extiende el MDX añadiendo soporte de operaciones espaciales que abren la posibilidad de mostrar cualquier dato que se pueda calcular mediante primitivas geométricas. Los límites son los que la semántica del GeoMDX ofrezca.

En la implementación actual se echa de menos un visor espacial para poder visualizar en un mapa las geometrías mostradas en WKT. Aunque esta tarea es aparentemente sencilla y se podría añadir con poco esfuerzo utilizando aplicaciones web de GIS existentes que ya son capaces de representar geometrías en WKT.

Os presento un artículo de investigación que algunos miembros de GIS&Chips publicamos el mes pasado en la revista “Environmental modelling & software”. Trata de la creación de un API libre y abierto para el cálculo de las métricas del paisaje.

La utilidad de las métricas del paisaje en distintas aplicaciones ha sido generalmente aceptada y ha provocado que existan muchos paquetes de software diseñados para proporcionar cálculos y análisis de los patrones estructurales del paisaje. Es posible obtener más información aqui

Tras examinar detenidamente las herramientas más utilizadas (Fragstats, V-Late, PA4, etc), se ha podido extraer una serie de puntos fuertes y débiles con la finalidad de crear una lista de las características deseables en este tipo de software. Tras dicho análisis se consideró necesario el diseño de un API sin limitaciones en los datos de entrada, capaz de calcular a partir de datos vectoriales o raster, etc. Este API debería facilitar no sólo la construcción de aplicaciones propias, sino que también debería permitir añadir nuevas métricas y la investigación de nuevos paradigmas relacionados con las métricas del paisaje. Con estas premisas se ha comenzado a desarrollar una propuesta, basada en estándares abiertos y software libre, que se ha denominado landmetrics-DIY (“Do It Yourself”). Podéis encontrar una versión alfa junto con un sencillo interfaz de usuario haciendo checkout en el siguiente repositorio subversión:

http://www.gisandchips.org/svn/landmetrics_diy

Tendréis que volver a añadir las referencias a NTS y os daréis cuenta de que el desarrollo está aún muy verde. Era solamente una primera propuesta, pero creemos que las directrices definidas en el artículo deberían marcar los criterios a seguir en el desarrollo de proyectos geoespaciales científicos de este tipo.

Por el momento, este API puede calcular unas 40 métricas del paisaje a partir de ficheros vectoriales Shapefile de ESRI, pero estamos trabajando para completar su contenido, siguiendo las líneas que se explica en el artículo.

Aquellos que estéis interesados en saber más, podéis ver el poster que hemos añadido al principio del post o podréis encontrar el artículo en ScienceDirect

Hasta pronto.

He probado diferentes métodos (sin transformación, transformación EPSG) y el resultado es muy dispar. Sólo un de los métodos funciona bien, aquel cuya transformación utiliza la rejilla en formato NTv2 del IGN, que fue creado precisamente para la transformación de datos desde ED50.

ogr2ogr -s_srs "EPSG:23030" -t_srs epsg_destino destino origen

Este comando tampoco solucionaría los problemas para ED50, por lo que habría que retocar el comando de esta forma para que utilice la rejilla del IGN:

ogr2ogr -s_srs '+init=epsg:23030 +nadgrids=./peninsula.gsb +wktext' -t_srs epsg_destino origen destino

Coincidiréis conmigo que este comando con todas sus opciones no es muy intuitivo para recordar. Por esta razón he decidido utilizar un interfaz gráfico para el terminal (sí, no habéis oído mal) que permite lanzar los típicos formularios gráficos para rellenar los datos. En definitiva, no se hace nada más que encapsular dicho comando con cuadros de diálogo que se lanzan al X-windows. Hay algunas utilidades para ello: Xdialog (x-windows), dialog (para el terminal), kdialog (para KDE). Yo voy a utilizar uno llamado “Zenity”, que está en todas las distribuciones linux.
Finalmente he creado un script BASH que empotra las peticiones (selección del EPSG, de ficheros, alertas, comprobaciones, etc.) en diálogos hechos con Zenity.
Este es el codigo:

#!/bin/sh
zenity --info --text="OGR2OGR Dialog \n
Autor: José Manuel Mira \n
GISandCHIPS 2011 \n
Aplicación basada en OGR2OGR para convertir \n
shapefile proyectado en UTM con datum ED50 en huso 30 (EPSG:23030) a  \n
Geodésica mundial con datum WGS84 (epsg:4326) o \n
UTM ETRS89 (epsg:25830) utilizando GDAL. \n
NOTA: Esta utilidad necesita la rejilla del IGN \n
descargable en: \n

http://www.gisandchips.org/demos/j3m/ogr/peninsula.gsb"

# Comprobar que está instalado OGR2OGR
if which ogr2ogr > /dev/null; then
echo "GDAL está instalado"
else
echo "No está instalado GDAL" | zenity --error --text="No está instalado GDAL. Salir \n ¿Ubuntu? sudo apt-get install gdal-bin"
exit
fi
#Comprobamos que tienes la rejilla instalada en el directorio del script
if [ -f ./peninsula.gsb ]
then
  echo "OK, tienes la rejilla" | zenity --info --text="OK, tienes la rejilla"
else
  echo "Descargando rejilla ..."
  wget http://www.gisandchips.org/demos/j3m/ogr/peninsula.gsb 2>=1 | sed -u 's/.*\ \([0-9]\+%\)\ \+\([0-9.]\+\ [KMB\/s]\+\)$/\1\n# Downloading \2/' | zenity --progress --title="Descargando rejilla IGN ..."
fi

origen=$(zenity --file-selection --title="Selecciona un shapefile en proyección ED50")
case $? in
  0)
    echo "\"$origen\" selected.";;
  1)
    exit;;
   -1)
    exit;;
esac

# obtener extensión del fichero
ext=${origen#*.}
echo $ext
if [ $ext != 'shp' ]
then
   echo "No es un shapefile. Salir" | zenity --error --text="NO es un shapefile. Salir"
   exit
fi

destino=$(zenity --entry --text "Nombre del Shapefile de salida (sin .shp)?" --entry-text "destino"); echo $destino
if [ $? = 1 ]
then
 exit
else
 echo $?
fi

epsg=$(zenity  --list  --text "Selecciona la proyección de salida" --radiolist  --column "Sel" --column "EPSG" --column "Nombre" TRUE 4326 "Geodésica WGS84" FALSE 25830 "UTM ETS89 HUSO 30"  FALSE 25831 "UTM ETS89 HUSO 31"  FALSE 25829 "UTM ETS89 HUSO 29"); 

echo "ogr2ogr -s_srs '+init=epsg:23030 +nadgrids=./peninsula.gsb +wktext' -t_srs EPSG:"${epsg} ${destino}-${epsg}.shp $origen
ogr2ogr -s_srs '+init=epsg:23030 +nadgrids=./peninsula.gsb +wktext' -t_srs EPSG:${epsg} ${destino}-${epsg}.shp $origen
zenity --info --text "Shapefile ${destino}-${epsg}.shp creado satisfactoriamente"

Guardamos el script como conversor.sh, le damos permisos de ejecución y lo lanzamos con

sh ./conversor.sh

Estos son algunos de los diálogos que aparecen:




Sí, ya se lo que estáis pensando (tanto rollo para una sola línea de comando). Es cierto, pero resulta útil y sencillo para aquellos que no se quieren complicar la vida. Además, con ligeras modificaciones lo puedo hacer extensible a las islas, otros husos, escribir un log con comandos para reutilizarlos en el futuro en modo BATCH, etc)
He realizado algunas pruebas para comprobar la validez de los datos reproyectados utilizando coordenadas de vértices geodésicos y visualizadas sobre un fondo con el WMS del PNOA. El resultado salta a la vista: “lo ha cuadrado”

Sí deseas probar el script aquí tienes unas shapes de prueba (vértices geodésicos en 23030) de la Comunidad Valenciana)

Más adelante, cuando tenga tiempo, colgaré una librería que he desarrollado en mi tesis y explicaré en detalle el código.

Hasta pronto.

Espero que os guste. Cuando tenga un rato subo más!!

El problema de las extensiones de SQLite

SpatiaLite es fundamentalmente una librería (libspatialite) que actúa como extensión de SQLite: explota el índice R-Tree de SQLite y añade el tipo Geometry y una serie de funciones que a su vez invocan a las de GEOS y PROJ. El resultado, al igual que el tándem PostgreSQL / PostGIS, es una implementación de la especificación “Simple Features for SQL” del OGC pero sin necesidad de una infraestructura cliente-servidor. Este tipo de librerías, que contienen funciones que extienden las propias de SQLite, se cargan dinámicamente, es decir, es el desarrollador el encargado de enlazarlas en tiempo de ejecución en el ámbito de una conexión de base de datos. El problema es que este mecanismo de extensión constituye a su vez un agujero de seguridad, un “coladero” de software malicioso. Así que, por defecto, las distros de Linux proporcionan el binario de SQLite compilado de forma que no se autoriza la carga dinámica de extensiones. Supongo que al final, quienes lo necesitan, optan por compilar su propio binario de SQLite o utilizar el que proporciona SpatiaLite, donde la extensión está estáticamente enlazada. Los que trabajáis con Mono tenéis ahora otra opción más.

Org.Gisandchips.Sqlite

Podéis descargar el código fuente de este proveedor ADO.NET para SQLite desde el repositorio SVN de Gis&Chips:

svn co http://www.gisandchips.org/svn/sqliteprovider

Para vuestra comodidad, os proporciono el fichero de proyecto MonoDevelop. Yo compilo sin problemas en una máquina Ubuntu 11.04, con Mono 2.10, MonoDevelop 2.6 y bajo el framework .NET 4.

Si echáis un vistazo al código comprobaréis que se trata del mismo proveedor Mono.Data.Sqlite (obtenido de la master branch del repositorio público de Mono en github) con los siguientes añadidos:

• Wrappers para las funciones enable_load_extension() and load_extension() del API de sqlite3.
• Método público SqliteConnection.LoadExtension() que permite cargar dinámicamente extensiones de sqlite3.

El siguiente código de ejemplo es suficientemente explicativo:

using (var cn = new SqliteConnection ("Data Source=path/file.sqlite")) {
  cn.Open ();
  cn.LoadExtension ("libspatialite.so.2");
  //your sql commands follow
}

Respecto a la licencia, me remito a la de Mono (MIT), que es la más permisiva y entiendo que compatible con las instrucciones del autor original de System.Data.SQLite, Robert Simpson, quien liberó este proveedor como de dominio público. Que lo disfrutéis.

En estos eventos es habitual el intercambio de tarjetas de presentación y como no teníamos una para darnos a conocer decidimos hacer una prueba en LaTeX, que no salió del todo mal. En este post compartimos esta primera versión para que la podáis personalizar o simplemente ver como se ha hecho.

Cara tarjeta G&C, generado con el fichero "G&C_businessCard_front.tex".

Características a considerar en el diseño

Hace mucho tiempo di un par de cursos para el diseño de marca. Ilustrator, Freehand y Photoshop para el diseño de tarjetas, trípticos, carteles, carátulas de cds, etc. Hay que reconocer que todo esto es una profesión en sí y hay muchos detalles por conocer. En tal caso, lo que hacemos nosotros aquí es un “apaño”, por lo que nadie espere grandes cosas

Las características más a tener en cuenta en el diseño son el tamaño de la tarjeta, el logo, el fondo y la gama de colores. Otra cuestión a tener en cuenta es la problemática vector-raster.

El tamaño de la tarjeta se puede determinar por una serie de estándares, supongo que según el tamaño de los billetes de cada país (la cartera donde guardamos las tarjetas también depende de esto). Nuestra tarjeta (90x45mm) no se ajusta a ninguno de estos estándares aunque la podríamos adaptar muy fácilmente. En http://en.wikipedia.org/wiki/Business_card encontramos los tamaños más utilizados.

El logo es el  que figura en el blog, diseñado por Miguel. No obstante, eliminamos el texto pues lo incluiremos entre los otros elementos de la tarjeta.

El fondo de la cara frontal es una nube de tags basada en la de esta página, después en Gimp le hemos añadido un degradado lineal de negro a blanco. El fondo del reverso es un mapa de http://www.openstreetmap.org/ que indica la ubicación de la sede. Para completar la ubicación se podría añadir la dirección en texto.

El nombre de la asociación, el lema y los datos de contacto (web y Email) aparecen en la cara frontal. El contraste blanco/negro es uno de los más recurridos y elegantes aunque esto dependerá del gusto de cada uno, además son los colores que aparecen en el blog.

Finalmente, la problemática raster-vector hace que no todos los programas sean igual de prácticos al combinar imágenes con texto en un marco tan pequeño. El texto podría verse mal al rasterizarse en Gimp o las imágenes verse mal por no estar bien escaladas al tamaño del papel. Para evitar problemas es recomendable preparar todos los componentes del diseño en un marco de 90x45mm. Además, conviene que el texto sea siempre vectorial.

Reverso de la tarjeta G&C, generado con el fichero "G&C_businessCard_back.tex".

¿Porqué LaTeX?

Seguramente, podríamos haber creado la tarjeta solamente con Gimp, Impress, Inkscape o cualquier otro software comercial… incluso en un editor de textos cualquiera. O con una máquina de escribir… No obstante, para crear la tarjeta hemos utilizado Gimp y LaTeX.

Básicamente, prefiero esta opción porqué es FOSS y estándar. No necesitamos conocer demasiados botones y LaTeX lleva una larga trayectoria de más de 20 años que hace pensar que la plantilla que creemos seguirá siendo funcional en el futuro. Además, ya utilizo LaTeX para hacer pósters, artículos, presentaciones, materiales docentes, el currículum y mi tesis. De este modo, la pregunta sería ¿porqué aprender otros programas?

Desgraciadamente, no puedo aquí hacer un tutorial de LaTeX para geógrafos pero, si necesitáis escribir un documento estructurado de grandes dimensiones (tesis) o muchos de otro tipo (informes) os animo a que le echéis un vistazo a este proyecto tan interesante pues a medio plazo os ahorrará muchos dolores de cabeza.

Como está hecha

Las plantillas que he encontrado en la web (buscando “Business card; LaTeX”) son más bien sosas con un estilo próximo al de la máquina de escribir y además son bastante complejas. Para evitar estos problemas hemos utilizado la clase “Beamer” de LaTeX, que permite entre otras cosas añadir una imagen de fondo y formatear el espacio de la tarjeta apropiadamente.

Primero, hemos utilizado Gimp para crear los fondos de la tarjeta como dos ficheros PNG. A continuación, LaTeX se ha utilizado para la composición, añadir textos y dibujos vectoriales (líneas). En el reverso de la tarjeta hemos añadido un mapa de http://www.openstreetmap.org/ para aportar una localización más detallada, aunque también se podría añadir la dirección postal.

La plantilla se compone de tres ficheros TeX, uno para cada cara de la tarjeta y, tras compilar las dos caras, el tercero para crear una composición de 12 tarjetas a dos caras.

Esta plantilla es altamente mejorable, es una versión 0.1 para salir del paso, pero os puede dar ideas. La idea de aprovechar OSM para la localización (idea de José) es interesante y evita tener que pensar en hacer un croquis como he visto en algunas tarjetas. Además, creo que los detalles de la tarjeta (la nube de tags como fondo, el mapa…) además de quedar bastante bien aportan información extra sobre nosotros si se saben leer.

Código en latex

Finalmente, os paso el código de los tres ficheros de LaTeX necesarios para crear esta tarjeta. Si leéis el código veréis donde podéis introducir vuestras propias imágenes de fondo o el logo que queráis. Espero que os sea de ayuda y os genere otras ideas.

G&C_businessCard_front.tex

% GIS&Chips nice business card
% By Benito M. Zaragozí
% Version 0.1 released 02/11/2011
% Further releases in: www.gisandchips.org
% This template is composed by three .tex files (businessCard.tex, businessCard_back.tex and businessCardx10.tex) for preparing a basic layout.

\documentclass{beamer}
\usepackage[utf8x]{inputenc}
\usepackage[spanish]{babel}
\usepackage{hyperref}
\usepackage[absolute,showboxes,verbose,overlay]{textpos}
\usepackage{geometry}
\geometry{paperwidth=90mm, paperheight=45mm, layoutwidth=90mm, layoutheight=45mm, left=0mm, top=0mm, right=0mm, bottom=0mm}

\usetheme{default}

\setbeamertemplate{navigation symbols}{}
\setbeamertemplate{background canvas}{\includegraphics[height=\paperheight]{background.png}}

\begin{document}

\frame{
\begin{columns}
	\begin{column}{35mm}
		\includegraphics[width=44mm]{gischips_logo.png}
	\end{column}

	\begin{column}{40mm}
		\color{white}{Asociación\\ \vspace{1mm} \huge{GIS \& Chips}}
		\normalsize
		\vspace{1mm}
		\line(1,0){80}
		\vspace{-1mm}
		\color{white}{Geografía útil para llevar}
	\end{column}
\end{columns}

\vspace{-1mm}
\normalsize
\color{black}{\hspace{5mm}\line(1,0){220}\\ \hspace{2mm} \url{www.gisandchips.org} \hspace{3mm} \url{info@gisandchips.org}}
}
\end{document}

G&C_businessCard_back.tex

% GIS&Chips nice business card
% By Benito M. Zaragozí
% Version 0.1 released 02/11/2011
% Further releases in: www.gisandchips.org
% This template is composed by three .tex files (businessCard.tex, businessCard_back.tex and businessCardx10.tex) for preparing a basic layout.

\documentclass{beamer}
\usepackage[utf8x]{inputenc}
\usepackage[spanish]{babel}
\usepackage{geometry}
\geometry{margin=2mm, paperwidth=90mm, paperheight=45mm}

\usetheme{default}
\setbeamertemplate{navigation symbols}{}
\setbeamertemplate{background canvas}{\hspace{2mm} \includegraphics[height=41.5mm]{mapbackground.png}}

\begin{document}
	\frame{
		\vspace{-1.5mm}
		\colorbox{white}{GIS\&Chips en la Universidad de Alicante}\\
		\vspace{31mm}
		\begin{flushright}
		\colorbox{white}{\tiny {\url{www.openstreetmap.org}}}
		\end{flushright}
	}
\end{document}

G&C_businessCard_x12.tex

% GIS&Chips nice business card
% By Benito M. Zaragozí
% Version 0.1 released 02/11/2011
% Further releases in: www.gisandchips.org
% This template is composed by three .tex files (businessCard.tex, businessCard_back.tex and businessCardx10.tex) for preparing a basic layout.

\documentclass[10pt,a4paper]{minimal}
\usepackage{graphicx}
\usepackage[margin=1cm]{geometry}

\begin{document}
\thispagestyle{empty}
\noindent

\begin{center}
\hspace{1mm}\includegraphics[scale=1]{businessCard.pdf} \includegraphics[scale=1]{businessCard.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard.pdf} \includegraphics[scale=1]{businessCard.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard.pdf} \includegraphics[scale=1]{businessCard.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard.pdf} \includegraphics[scale=1]{businessCard.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard.pdf} \includegraphics[scale=1]{businessCard.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard.pdf} \includegraphics[scale=1]{businessCard.pdf} \vspace{1mm}

\pagebreak
%%%%%%Back with map
\hspace{1mm}\includegraphics[scale=1]{businessCard_back.pdf} \includegraphics[scale=1]{businessCard_back.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard_back.pdf} \includegraphics[scale=1]{businessCard_back.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard_back.pdf} \includegraphics[scale=1]{businessCard_back.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard_back.pdf} \includegraphics[scale=1]{businessCard_back.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard_back.pdf} \includegraphics[scale=1]{businessCard_back.pdf} \vspace{1mm}

\hspace{1mm}\includegraphics[scale=1]{businessCard_back.pdf} \includegraphics[scale=1]{businessCard_back.pdf} \vspace{1mm}

\end{center}
\pagebreak
\end{document}

Test de rendimiento SIMD con Mono.Simd.Vector8s

Las técnicas de computación paralela encuentran en los sistemas de información geográfica una de sus aplicaciones más elocuentes dada la naturaleza masiva de los datos geográficos, en particular si hablamos de estructuras de datos raster. No en vano, allá por el año 1997, Richard Healey, junto con otros compañeros del Dpto. de Geografía de la Universidad de Edimburgo, editó “Parallel Processing Algorithms For GIS”. El profesor Healey es un geógrafo que ya a principios de los 90 planteaba a sus alumnos problemas de geoprocesamiento paralelo que éstos habían de resolver con ADA, el único lenguaje accesible que por aquel entonces prefiguraba el paradigma OOP y permitía definir tareas concurrentes. Desde entonces se han producido muchos avances en paralelización de procesos, tanto a nivel de software como de hardware. Sin embargo, existe un tipo de paralelización que no parece haber alcanzado demasiada notoriedad en el campo de los GIS: se trata del conjunto de instrucciones SIMD (Single Instruction, Multiple Data) que permite realizar operaciones aritméticas entre pares de vectores mediante una sola instrucción. No se trata por tanto de concurrencia de procesos, sino de paralelismo a nivel de datos. Lo más curioso es que, desde 1996, los fabricantes de microprocesadores han ido incorporando extensiones que soportan en mayor o menor medida este juego de instrucciones.

Cualquier PC o videoconsola actual permite ejecutar instrucciones SIMD. Evidentemente los desarrolladores de juegos hace tiempo que explotan este tipo de paralelización. Sin embargo, encontrar referencias recientes a SIMD en el ámbito de la Geomática resulta más complicado: aquí tenéis una interesante ponencia de la conferencia de usuarios de GRASS del 2004. Seguramente habrá más trabajos publicados, así que hacednos saber lo que encontréis.

Yo sinceramente no sabía nada de esto hasta que leí este post de Miguel de Icaza, anunciando la publicación de Mono.Simd. Este API es un paso en la estrategia de Mono por entrar de lleno en la industria del “gaming”. ¿Y para las desarrolladores de GIS? En mi opinión significa un abanico de posibilidades a la hora de paralelizar geoprocesos que queda pendiente de explorar.

En este artículo sólo quiero dejar constancia de un test para un caso extremo de optimización con SIMD. En la práctica carece de utilidad, mi única intención es medir diferencias de tiempo en un escenario absolutamente favorable a la paralelización mediante vectores. Dado un grid de elevaciones, donde cada valor de altitud puede representarse como un entero de 16 bits, se trata simplemente de doblar el valor de cada pixel para “exagerar” el modelo del terreno. Estas son las especificaciones de la máquina en que he realizado el test:

• Intel Q8300 @ 2.50GHz
• 3.9 GiB RAM
• Ubuntu 11.04 i686
• Mono 2.10.5

El test se ha realizado sobre 4 ficheros ARC/INFO ASCII GRID de distintos tamaños. Para cada fichero se han tomado 5 muestras del tiempo transcurrido en el cálculo aritmético sin aceleración de hardware, pixel a pixel, tomando finalmente la mediana. La misma metodología de muestreo se ha empleado para obtener los tiempos con aceleración de hardware, es decir, empleando Mono.Simd para multiplicar los pixels de 8 en 8 mediante una única instrucción. Para ello es necesario agrupar previamente los valores de altitud y almacenar cada vector en una lista de tipo Mono.Simd.Vector8s. Estos son los resultados:

Como veis, en todos los casos, con SIMD el proceso de cálculo se completa en torno a 3,5 veces más rápido. Esto es una buena noticia, pero no significa que este tipo de paralelismo sea la panacea. En mi opinión la conclusión es que, si con un escenario favorable se puede reducir el tiempo de proceso en más de un 70%, merece la pena explorar otros escenarios con verdadera utilidad práctica.

Si alguien se pregunta por la incidencia que puede tener esta mejora en el rendimiento si contabilizamos el proceso de lectura y carga en memoria de los datos, me atrevería a decir que escasa. En cualquier caso, los procesos de lectura y manejo de grandes volúmenes de información raster constituyen un problema distinto y existen diversas técnicas para su optimización.

Si os pica la curiosidad, podeis averiguar el conjunto de instrucciones SIMD que soporta vuestro procesador ejecutando el código de ejemplo que aparece en la documentación de la clase Mono.Simd.SimdRuntime (es posible que tengáis que hacer algunas modificaciones). Para aquellos que quieran reproducir un test de rendimiento con sus propios ficheros ASCII GRID, aquí os dejo el código C# que he usado (por cierto, ni es óptimo ni puedo garantizar que esté libre de bugs):

using System;
using System.IO;
using Mono.Simd;
using System.Collections.Generic;
using System.Globalization;
using System.Diagnostics;
using System.Linq;

namespace SIMDTest
{
 class MainClass
 {
  public static void Main (string[] args)
  {
   bool useSIMD = true;
   string path = string.Empty;
   CultureInfo formatProvider = null;
   try {
    if (args.Length == 0) {
     Console.WriteLine ("Usage: SIMDTest [ASC GRID file name] [specific culture name]");
     return;
    }
    if (args.Length > 0)
     path = args[0];
    if (args.Length > 1)
     formatProvider = CultureInfo.GetCultureInfo (args[1]);
    Console.WriteLine ("Use hardware acceleration? [y / n]? : ");
    ConsoleKeyInfo key = Console.ReadKey ();
    while (key.Key != ConsoleKey.Y && key.Key != ConsoleKey.N) {
     Console.WriteLine ("Press 'y' to use hardware acceleration, 'n' otherwise: ");
     key = Console.ReadKey ();
    }
    Console.WriteLine ();
    Console.WriteLine ("Processing, please wait ...");
    useSIMD = (key.Key == ConsoleKey.Y);
    ElevationModel dem = new ElevationModel(path, formatProvider);
    Stopwatch stopwatch = new Stopwatch ();
    if (useSIMD) {
     Vector8s scale = new Vector8s (2);
     foreach (IList vPage in dem.EnumerateVectors ()) {
      stopwatch.Start ();
      for (int i = 0; i < vPage.Count; i ++)
       vPage[i] = vPage[i] * scale;
      stopwatch.Stop ();
     }
    } else {
     foreach (IList page in dem.EnumerateData ()) {
      stopwatch.Start ();
      for (int i = 0; i < page.Count; i ++)
       page[i] = (short) (page[i] * 2);
      stopwatch.Stop ();
     }
    }
    Console.WriteLine ("Time elapsed: {0}", stopwatch.Elapsed);
   } catch (Exception ex) {
    Console.WriteLine (ex.Message);
   }
  }
 }

 class ElevationModel
 {
  string _path;
  IFormatProvider _formatProvider = CultureInfo.CurrentCulture;
  short _noData = -9999;
  int BOD = 0;

  public ElevationModel (string path, IFormatProvider formatProvider)
  {
   if (string.IsNullOrEmpty (path))
    throw new ArgumentNullException ();
   else if (!File.Exists (path))
    throw new FileNotFoundException ();
   else {
    _path = path;
    if (formatProvider != null)
     _formatProvider = formatProvider;
    this.SetMetadata ();
   }
  }

  public ElevationModel (string path) : this(path, null)
  {
  }

  void SetMetadata ()
  {
   using (StreamReader reader = new StreamReader (_path)) {
    //skip GRID header
    for (int i = 0; i < 5; i++)
     reader.ReadLine ();
    _noData = short.Parse (reader.ReadLine ().Split (' ')[1], _formatProvider);
    reader.BaseStream.Position = 0;
    reader.DiscardBufferedData ();
    char[] buffer = new char[512];
    reader.Read(buffer, 0, buffer.Length);
    BOD = this.GetBeginOfData (buffer);
   }
  }

  int GetBeginOfData (char[] source)
  {
   int i = 0;
   bool isNewHeaderLine = false;
   bool isNewLine = true;
   foreach (char c in source) {
    if (char.IsLetter (c) && isNewLine) {
     isNewHeaderLine = true;
     isNewLine = false;
    } else if ((char.IsDigit (c) || c == '-') && isNewLine && !isNewHeaderLine) {
     break;
    } else if (char.IsControl (c)) {
     isNewLine = true;
     isNewHeaderLine = false;
    }
    i++;
   }
   return i;
  }

  IList<short> GetPage (char[] source, int sourceIndex, int sourceLength, IList pendingDigits)
  {
   IList<short> page = new List<short> ();
   for (int i = sourceIndex; i < sourceLength; i++) {
    if (char.IsDigit (source[i]) || source[i] == '-')
     pendingDigits.Add (source[i]);
    else if (pendingDigits.Count > 0 && (char.IsWhiteSpace (source[i]) || char.IsControl (source[i]))) {
     short n = _noData;
     string s = new string(pendingDigits.ToArray ());
     if (!short.TryParse(s, NumberStyles.Integer, _formatProvider, out n))
      Console.WriteLine ("Error parsing {0}", s);
     page.Add (n);
     pendingDigits.Clear ();
    }
   }
   return page;
  }

  public IEnumerable<IList<short>> EnumerateData ()
  {
   using (StreamReader reader = new StreamReader (_path)) {
    char[] buffer = new char[25000 * 1024];
    int beginOfData = 0;
    IList pendingDigits = new List ();
    int charCount = 0;
    while ((charCount = reader.Read (buffer, 0, buffer.Length)) > 0) {
     if (beginOfData == 0) {
      beginOfData = this.BOD;
      yield return this.GetPage (buffer, beginOfData, charCount, pendingDigits);
     } else
      yield return this.GetPage (buffer, 0, charCount, pendingDigits);
    }
   }
  }

  IList<Vector8s> GetVectorPage(IList page, IList pendingValues)
  {
   List vPage = new List ();
   foreach (short n in page) {
    pendingValues.Add (n);
    if (pendingValues.Count == 8 ) {
     vPage.Add (new Vector8s (pendingValues[0], pendingValues[1],
      pendingValues[2], pendingValues[3], pendingValues[4],
      pendingValues[5], pendingValues[6], pendingValues[7]));
     pendingValues.Clear ();
    }
   }
   return vPage;
  }

  public IEnumerable<IList<Vector8s>> EnumerateVectors ()
  {
   int lastPageSize = -1;
   IList pendingValues = new List ();
   foreach (IList page in this.EnumerateData ()) {
    IList vPage = this.GetVectorPage (page, pendingValues);
    // complete last vector in last page
    if (lastPageSize != -1 && vPage.Count < lastPageSize) {
     if (pendingValues.Count > 0 && pendingValues.Count < 8 ) {
      for (int i = pendingValues.Count; i < 8; i++)
	pendingValues.Add(_noData);
      vPage.Add (new Vector8s (pendingValues[0], pendingValues[1],
       pendingValues[2], pendingValues[3], pendingValues[4],
       pendingValues[5], pendingValues[6], pendingValues[7]));
     }
    }
    lastPageSize = vPage.Count;
    yield return vPage;
   }
  }
 }
}

En algunas ocasiones la solución a estos problemas sería “trocear” el raster y procesarlo por partes. Esto lo podríamos hacer con varios software y en todos ellos sería interesante poder automatizar la tarea al máximo.

1. Modelo digital raster

En este post propongo realizar una prueba con RSAGA, que es un modulo de R que permite acceder a las funciones disponibles en la consola de SAGA GIS. Para este post he trabajado en Windows, con R 2.10.1 y SAGA 2.0.4, aunque supongo que no habrá problemas con usar otras versiones más recientes. Además, deberéis instalar en R el paquete RSAGA.

2. Shapefile de polígonos con que trocearemos el raster

Para esta práctica podéis descargar dos ficheros (Imagenes 1 y 2), uno raster y otro vectorial (polígonos). A efectos didácticos se utilizan ficheros pequeños, pero esto resultaría más útil en caso de necesitar fragmentar más el raster. Por ejemplo, en artículos anteriores he presentado código para trabajar imágenes a nivel de parcelas agrícolas. En aquellos artículos trabajábamos las imágenes de una en una y ya estaban recortadas, pero las parcelas podrían ser miles y no queremos hacer eso a mano.

Conociendo SAGA GIS

Lo principal, como en tantos otros casos, sería conocer como se trocea un fichero raster con SAGA GIS. Para ello abrimos la GUI de SAGA y vamos a la pestaña “modules” que tiene forma de TOC, generalmente se encuentra a la izquierda de la ventana principal. En esta pestaña podemos encontrar las distintas librerías que agrupan módulos que se relacionan por algún motivo. Explorar estas librerías y practicar con ellas es el mejor modo de estar seguro de lo que se quiere hacer.

3. TOC de SAGA GIS donde podemos explorar las librerías y los módulos disponibles.

Finalmente, he decidido que lo mejor es empezar por separar el shapefile original por polígonos (shapes-tools/separate shapes) y después realizar un recorte por cada uno de los nuevos shapefiles (shapes-grid/clip grid with polygon). Practicar con el interfaz gráfico es un bueno modo de estar seguros de los parámetros que nos pide cada módulo.

Antes de empezar a usar RSAGA he tenido que averiguar los nombres exactos de los módulos y los parámetros que aceptan. Esto lo he hecho con los métodos rsaga.get.modules() y rsaga.get.usage(). Aunque generalmente se puede averiguar también viendo los nombres en el SAGA GUI o mirando en la carpeta donde se contienen las *.dll.

Trabajando con RSAGA

## Abrimos la consola de R y cargamos la librería
library(RSAGA)

## Seleccionamos el fichero raster que queremos trocear y una capa vectorial que queramos usar como límites. También especificamos el directorio donde van los outputs. Para trabajar con Windows recomiendo rutas sin espacios.
raster<- file.choose()
poligonos<- file.choose()
directorio <-choose.dir()

## Consultamos hasta encontrar la herramienta que nos separa un shapefile en varios, obteniendo un shapefile por cada polígono, o lo que quisiéramos. Por ejemplo:
## rsaga.get.modules("shapes_grid")
## rsaga.get.usage("shapes_tools", 7)
## Se ejecuta el método con los parámetros necesarios.
rsaga.geoprocessor(lib="shapes_tools", module=7, param=list(SHAPES=poligonos, PATH=directorio, NAMING=0, FIELD=6))

Una vez ejecutado este proceso ya disponemos de un shapefile por cada polígono, esto significa que aunque el shapefile original tuviera miles de registros solamente cargaremos en memoria uno a la vez. Esto ya supone un ahorro importante de recursos.

# Listamos los nuevos shapefiles
shapefiles <- list.files(directorio, full.names=T, pattern="\\.shp")
## Por último utilizamos RSAGA para realizar tantos recortes del raster como polígonos habíamos extraído
for(i in 1:length(shapefiles)){
select <- shapefiles[i]</code>
rsaga.geoprocessor(lib="shapes_grid", module=7, param=list(OUTPUT=select, INPUT=raster, POLYGONS=select))
}#fin bucle

Conclusiones

Tras ver este ejemplo nos queda más claro el uso de RSAGA y se visualiza bien como se podrían automatizar gran número de tareas habituales en los GIS. Intentar realizar estas tareas herramientas básicas y generales podría ser algo trabajoso.

En nuestro fichero de prácticas había 18 polígonos por lo que nuestro resultado es una carpeta donde encontramos guardados 18 shapefiles y 18 ficheros raster. Los recortes quedan como se puede ver en la siguiente imagen:

4. Uno de los 18 raster resultantes de los recortes y al fondo el shapefile original

En caso de ser necesario, se podría realizar más tareas sobre todos estos ficheros, con este software o con otros programas. Incluso se plantea la posibilidad de usar multithreading en R para ejecutar tareas con RSAGA. Proximamente haremos alguna demostración de esto.

Antes de acabar, solo quiero recordar quando trabajamos con RSAGA en Windows es posible que nos afecten algunas particularidades sobre el modo en que se forman las rutas (espacios, barras de directorios…) y los mensajes de error que aparecen no son de gran ayuda.

——————————————————————

Si queréis contactar podéis enviarme un email (asunto: gisandchips):

Benito M. Zaragozí

benito.zaragozi@ua.es