Seguramente muchos de vosotros habéis estado en la tesitura de tener que proyectar vuestras capas de datos gis que usan la vetusta y a extinguir proyección UTM con datum ED50. Lo de que dicha proyección “no esta de moda” no es algo causal, puesto que hay un Real Decreto, el 1071, que define como sistema geodésico oficial de referencia para España aquel con el datum ETRS89 y el elipsoide SGR80. Bien es cierto que el decreto permite la ambivalencia de ambos sistemas (ETRS89 y ED50) hasta el año 2015; fecha a partir del cual diremos adiós definitivamente a ED50 (Elipsoide Hayford).
Por otra parte, con el auge de los wikimaps es bastante habitual tener que realizar reproyecciones de las capas gis con ED50 “de toda la vida” al sistema geodésico mundial WGS84, que es lo que admiten casi todas las geowikis (OpenStreetMap, ikiMaps, geoNames, etc.).

Pero, ¿que es lo que ocurre cuando utilizamos de forma genérica las utilidades de reproyección de los programas de GIS?. Casi todos, generan datos erróneos, que de ninguna manera se corresponden con la realidad. He encontrado translaciones en la X de más de 100 metros, y unos 50 en la Y respecto a la coordenada donde debería estar.

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.

Para este artículo he creído conveniente proponer un método de reprojección que no haga uso de ningún programa GIS, y que podamos lanzarlo cómodamente desde la consola (eso de “cómodamente” a alguno le sonará a coña). Lo que quería era simplemente un programa que le indicase mi dato en UTM ED50 y me lo proyecte a lo que quiera. Finalmente he decidido crear un híbrido que esté a medio camino entre el “terminal” Linux y un interfaz gráfico (en este caso “GTK”).
Como no se trata de reinventar la rueda, he utilizado una herramienta con la que me siento muy cómodo y me ha sacado más de una vez de un apuro. Se trata de “ogr2ogr”, un comando que forma parte de la suite de GDAL. Este comando está diseñado para cambiar de formato nuestros datos, pero también para reproyectar. Si sintaxis está llena de modificadores que lo hacen a veces un poco críptico.
Entrando en harina, sí deseamos reproyectar un fichero en ED50 a otra proyección la sintaxis es:

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)