Archivo de la etiqueta: WGS 84

¿Latitud/Longitud o Longitud/Latitud? El caso de los Shapefiles

La verdad es que no me imaginaba la complejidad que tiene todo lo que tenga que ver con sistemas de coordenadas de referencia: a pesar de existir estándares, casi nadie los respeta.

Un caso muy curioso es el de los Shapefiles. Los archivos Shapefile llevan asociado un archivo .prj que define la proyección cartográfica en la que se almacenan las coordenadas de las entidades en el archivo.

El contenido de este archivo es una cadena Well Known Text, que forma parte del estándar Simple Feature Access del Open Geospatial Consortium.

Este estándar define que el orden en el que se almacenan las coordenadas lo define el sistema de coordenadas de referencia, luego si el sistema de coordenadas de referencia es un sistema de coordenadas geográfico con los ejes (Longitud, Latitud), las coordenadas a almacenar deben ir en el orden (Longitud, Latitud). Si por el contrario el sistema de coordenadas es también geográfico pero con los ejes (Latitud, Longitud), las coordenadas se deben almacenar en el orden (Latitud, Longitud).

Como el archivo .prj asociado a los Shapefiles es una cadena Well Kown Text todo hacía pensar que los programas que cargan archivos Shapefile y tienen asociado un sistema de coordenadas de referencia, deberían ser lo suficientemente inteligentes como para saber el orden de los ejes.

Digi3D 2011 cumple al 100% con el estándar. Además obtiene los parámetros de los sistemas de coordenadas de referencia de la base de datos EPSG Geodetic Parameter Dataset, y en esta base de datos todos los sistemas de coordenadas geográficos tienen sus ejes ordenados de la forma (Latitud, Longitud).

Por lo tanto, como Digi3D 2011 es estricto cumpliendo el estándar, si seleccionamos como sistema de coordenadas de salida uno geográfico (basado en EPSG), el archivo generado tendrá las coordenadas ordenadas (Latitud, Longitud).

En teoría, y confiando en que todos los programas cumplieran con el estándar tal y como lo hace Digi3D 2011, decidí exportar el famoso modelo de Bronchales a Shapefile para cargarlo con ArcGIS Explorer Desktop y cual es mi sorpresa que el modelo aparece con los ejes cambiados y además en Kenia.

Captura de pantalla de ArcGIS Explorer mostrando el modelo de bronchales con los ejes cambiados y en Kenia

El contenido del archivo .prj del Shapefile cargado con ArcGIS Explorer en la captura de pantalla anterior es el siguiente:

GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree (supplier to define representation)",0.01745329251994328,AUTHORITY["EPSG","9122"]],AXIS["Lat", North],AXIS["Long", East],AUTHORITY["EPSG","4326"]]

que es una cadena Well Known Text en la que se comprueba que los parámetros se han obtenido de la base de datos EPSG (el sistema de coordenadas geográfico WGS 84 tiene asociado el número 4326).

Si te fijas en la cadena Well Kown Text se están definiendo los ejes del sistema geográfico, se indica que el primero es Latitud y que el segundo es Longitud.

De todos modos el estándar define que si una cadena Well Kown Text incluye una autoridad (EPSG en este caso), y el sofware que carga el archivo en cuestión implementa esa autoridad, debe hacer caso omiso a la información que aparece en la georeferenciación y crear el sistema de coordenadas basándose en el código (4326 en este caso). De forma que pensé: Quizás ArcGIS explorer trate de forma especial el sistema de coordenadas 4326, de modo que quité todas las referencias a la autoridad EPSG del archivo .prj creando lo siguiente:

GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree (supplier to define representation)",0.01745329251994328],AXIS["Lat",North],AXIS["Long",East]]

Pero desafortunadamente el resultado fue el mismo: modelo representado con los ejes cambiados y ubicado en Kenia.

Parece que que ArcGIS Explorer espera siempre (Longitud/Latitud) independientemente de lo que se indique en el sistema de coordenadas.

Por último, y ya para confirmar mis sospechas, cambié el sentido de los ejes para ver si afectaba en algo a la representación/ubicación con la siguiente cadena:

GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree (supplier to define representation)",0.01745329251994328],AXIS["Long",East], AXIS["Lat",North]]

pero tampoco.

El siguiente conjunto de pruebas lo hice con Global Mapper con idénticos resultados: seleccionando como sistema de coordenadas de visualización Geographic (Latitude/Longitude) y con los tres archivos .prj asociados el resultado es el mismo.

GlobalMapperMostrandoBronchalesMal

Es posible que Digi3D 2011 no estuviera generando bien el archivo .prj, así que decidí estropearlo, y al intentar cargarlo en Global Mapper este mostró un cuadro de diálogo indicando que el Shapefile no tiene información de proyección (por lo tanto si que se estaba generando bien, obviamente).

Si se acepta el cuadro de diálogo que indica que el Shapefile no tiene asociado un sistema de coordenadas, se muestra otro que nos permite seleccionar el sistema de coordenadas en el que están las coordenadas del archivo Shapefile, y ¡sorpresa! existe un botón titulado: Init From EPSG…. Genial, tan solo tenemos que pulsarlo y seleccionar el sistema de coordenadas 4326, que según EPSG, tiene los ejes ordenados de la forma (Latitud, Longitud) por lo tanto debe mostrar correctamente el modelo de Bronchales

Pues no. El resultado es el mismo, lo que me hizo pensar que nadie sigue el estándar (al menos en el caso de Shapefiles).

Investigando un poco, llegué a una discusión muy interesante entre usuarios de Microsoft SQL Server Spatial y los ingenieros que han desarrollado ese producto (dentro de muy poco hablaré mucho sobre SQL Server espacial y Digi3D). Ahí se dice que SQL Server Espacial cumple al 100% con el estándar (igual que Digi3D) y que por lo tanto las coordenadas para el caso del sistema de coordenadas 4326 se almacenan como Latitud/Longitud. Un determinado usuario sin embargo me dió la pista sobre cómo se comportan el resto de programas:

Coordinates in SHP (ESRI “shapefiles”), MIF (MapInfo) and other files produced by current GIS packages may be in many different coordinate systems. For lat/lon coordinate systems the order of coordinates is lon/lat in those systems.

y un poco más adelante dice lo siguiente (lo que te hace pensar que realmente sabe del tema)

Whew! … that’s what I mean by “universal” usage, as the above formats comprise, what? 99.99%? of the data out there. If anyone can point to even a single repository of “Latitude / Longitude” data that uses (lat, lon) ordering instead of (lon, lat) ordering, I’d be very grateful to see the URL to that repository.

Lo que significa que el resto de programas hacen caso omiso del estándar. En realidad me imagino que será porque existen archivos Shapefile en coordenadas geográficas desde mucho antes que este estándar. Piensa que el estándar es de 2006 (pero aún no he encontrado nadie que lo respete al 100% en lo referente a sistemas de coordenadas de referencia, quitando Digi3D 2011) y sin embargo el formato Shapefile data de 1998.

Aquí llegamos a un compromiso: Digi3D 2011 cumple al 100% el estándar y creo personalmente que si existe un estándar es para cumplirlo, pero no puedes dar la espalda a lo que en teoría es el estándar de facto: toda esa cartografía que ya existe, que está almacenada en coordenadas geográficas y generada antes de que existiera el estándar, por programas que posiblemente no tengan ni idea de que existe un estándar, y que símplemente asocian un archivo .prj a cada Shapefile generado copiando una plantilla con el mismo nombre que el archivo Shapefile y con extensión .prj.

Así que la única solución que se me ha ocurrido ha sido permitir configurar en Herramientas/Configuración/Importador Exportador de Shapefiles una opción para indicar el orden de las coordenadas en caso de sistemas de coordenadas geográficos. Por defecto (y para seguir el estándar de facto, que no es el que se debería seguir) la opción es Longitud/Latitud. De esta manera podemos saltarnos el estándar en el caso de Shapefiles. También podrás encontrar estas opciones para DWG y DGN.

El resto de exportadores importantes .bin, .bind, .kml, … no permiten esta configuración pues los binarios de Digi3D (tanto .bin como .bind) cumplen al 100% con el estándar y los .kml tienen asociado de forma implícita el sistema de coordenadas geográfico WGS 84 y además lo hace bien: primero se almacena la Latitud y luego la Longitud.

Actualizando la base de datos Epsg v7.9 para añadir transformaciones de cambio de datum ED50 -> ETRS89 y ED50 -> WGS84

Este post es avanzado, y supone que acabas de descargar una versión nueva de la base de datos EPSG y la has exportado al formato Microsoft SQL Server Compact tal y como se explica en el post Cómo se ha creado la base de datos geográfica EPSG en formato SQL Server Compact que se instala con Digi3D 2011.

La última versión de los archivos de rejilla de mínima curvatura NTv2 suministrados por el Instituto Geográfico Nacional data del 2009, sin embargo no existe ninguna transformación en la base de datos EPSG v7.9 (la mas moderna a día de hoy 17/04/2012) que utilice este archivo de mínima curvatura.

En este post te explico cómo modificar la base de datos para añadir una transformación que utilice el archivo de mínima curvatura actualizado.

La base de datos EPSG define 12 transformaciones para transformar entre los sistemas de coordenadas geográficos ED50 y ETRS89 y 40 transformaciones para transformar entre ED50 y WGS 84.

Puedes comprobarlo mediante las siguientes consultas SQL:

SELECT COUNT(*)
FROM Coordinate_Operation
WHERE SOURCE_CRS_CODE=4230 and TARGET_CRS_CODE=4258 and COORD_OP_TYPE='transformation';

SELECT COUNT(*)
FROM [Coordinate_Operation]
WHERE SOURCE_CRS_CODE=4230 and TARGET_CRS_CODE=4326 AND COORD_OP_TYPE='transformation';

Muchas de ellas consisten simplemente en desplazamientos geográficos (es decir, transformar de coordenadas geográficas a geocéntricas y realizar un desplazamiento en geocéntricas para luego volver a transformar a geográficas), por lo que tienen poca precisión.

De entre todos los algoritmos para transformar entre sistemas geográficos, el más preciso es el método interpolación mediante una rejilla de mínima curvatura. Esta rejilla le indica al programa: “Para este par de coordenadas latitud, longitud en el sistema de coordenadas de referencia ED50, le corresponde este otro par de coordenadas latitud, longitud en el sistema de coordenadas de referencia ETRS89”. Luego el programa tan solo tiene que interpolar mediante una interpolación bilineal. El resultado es muy preciso (20 cm para el caso de España).

Es responsabilidad de cada país/organismo crear su archivo de rejilla y de publicarlo. Cualquiera podría inventarse su propio formato de rejilla (lo que obligaría a que los desarrolladores de programas de cartografía y sistemas de información geográfica tuvieran que implementar un importador específico para cada formato de rejilla) pero existe uno estandarizado denominado NTv2.

En el caso de España, el organismo encargado de crear y actualizar ese archivo de rejilla es el Instituto Geográfico Nacional que publica este archivo a través de su página en internet: www.ign.es.

Descarga los archivos de rejilla entrando en www.ign.es: en el menú Herramientas de la página web, puedes seleccionar la opción: Rejilla para cambio de Datum entre ED50 y ETRS89 (en formato NTV2). Verás que puedes descargar dos rejillas: una para Península y otra para Baleares. Descarga el archivo Península, comprobarás que el nombre del archivo descargado es PENR2009.gsb.

Como la creación y mantenimiento de los archivos NTv2 depende del centro geográfico nacional de cada país, la base de datos EPSG no distribuye estos datos, sin embargo, Digi3D 2011 incorpora una copia de estos archivos en la carpeta %ProgramData%Digi3D 2011OpenGis, de modo que en ese directorio encontrarás el archivo PENR2009.gsb.

El único inconveniente es que la base de datos original EPSG (la que acabas de descargar de http://www.epsg.org/Geodetic.html no dispone de ninguna transformación que utilice este archivo, únicamente tiene transformaciones que utilizan versiones anteriores de ese archivo.

Vamos a comprobarlo:

La base de datos EPSG dispone de una tabla denominada ‘Coordinate_Operation Method’ que registra los códigos, nombres y descripciones de las distintas transformaciones.

Consultemos el código asociado al algoritmo NTv2 mediante la siguiente consulta SQL:

SELECT * FROM [Coordinate_Operation Method]
WHERE COORD_OP_METHOD_NAME='NTv2';

Puedes comprobar que el algoritmo NTv2 tiene asociado el código 9615.

Podemos consultar también los parámetros que recibe un determinado algoritmo. Comprobemos los parámetros que recibe el algoritmo 9615:

La base de datos EPSG dispone de una tabla ‘Coordinate_Operation Parameter Usage’ que nos indica los parámetros de cada algoritmo. Estos parámetros (como todo en la base de datos EPSG) tienen asignado un código. Disponemos de otra tabla ‘Coordinate_Operation Parameter’ que nos muestra una descripción de cada tipo de parámetro.

Para averiguar los parámetros que recibe la transformación con código 9615 podemos ejecutar la siguiente consulta SQL:

SELECT *
FROM [Coordinate_Operation Parameter Usage]
WHERE COORD_OP_METHOD_CODE=9615;

Eso nos devolverá un único registro indicando que el único parámetro que recibe el algoritmo 9615 es el parámetro cuyo código es el 8656.

Podemos consultar el significado de ese parámetro con la siguiente consulta SQL:

SELECT *
FROM [Coordinate_Operation Parameter]
WHERE PARAMETER_CODE=8656;

Por último, si quieres hacerlo bien, puedes hacerlo con una consulta JOIN de la siguiente manera:

SELECT *
FROM [Coordinate_Operation Parameter Usage] AS A
JOIN [Coordinate_Operation Parameter] AS B
ON A.PARAMETER_CODE = B.PARAMETER_CODE
WHERE A.COORD_OP_METHOD_CODE=9615;

De entre todas las posibles transformaciones (te recuerdo que había 12) entre ED50 y ETRS89 que define la base de datos EPSG, únicamente hay dos que utilicen el algoritmo 9615: La transformación 15895: ED50 to ETRS89 (11) y la operación 15932: ED50 to ETRS89 (12).

Puedes comprobarlo mediante la siguiente consulta SQL:

SELECT *
FROM Coordinate_Operation
WHERE SOURCE_CRS_CODE=4230 and TARGET_CRS_CODE=4258 and COORD_OP_TYPE='transformation' AND COORD_OP_METHOD_CODE=9615;

Si nos leemos los comentarios de estas dos transformaciones, rápidamente averiguamos que la transformación 15895 ha sido sustituida por la 15932, de modo que la que nos interesará utilizar es la 15932, pero tal y como vamos a comprobar muy pronto, esta transformación no nos sirve, pues el archivo de rejilla a que utiliza esta transformación no es “PENR2009.gsb”.

Los parámetros de las transformaciones se almacenan en la tabla ‘Coordinate_Operation Parameter Value’.

Esta tabla almacena los parámetros de todos los algoritmos utilizados en todas las posibles transformaciones especificadas en la base de datos EPSG. Actualmente tiene 15743 registros.

Cada registro tiene varios campos:

COORD_OP_CODE COORD_OP_METHOD_CODE PARAMETER_CODE PARAMETER_VALUE PARAM_VALUE_FILE_REF UOM_CODE

Si comprobamos los registros de esta tabla para el código de operación 15932, comprobaremos que el nombre del archivo de rejilla asignado para esa operación es “SPED2ETV2.gsb”.

SELECT * FROM [Coordinate_Operation Parameter Value] WHERE COORD_OP_CODE=15932;

Por lo tanto esta operación no nos sirve. El archivo “SPED2ETV2.gsb” es el nombre que tenía la versión antigua de los archivos de rejilla, de hecho, la operación 15932 se añadió a la base de datos EPSG con fecha 27-03-2007, y el nombre de archivo “PENR2009” nos hace pensar que este archivo se creó en el año 2009.

Podrías pensar que una solución sería modificar este registro de la base de datos y cambiar la palabra “SPED2ETV2.gsb” por “PENR2009.gsb”, pero esto no se puede hacer porque la base de datos EPSG es inmutable. No se permiten cambios. Podemos añadir información, transformaciones, sistemas de coordenadas, pero nunca cambiar nada.

Así que no nos queda más remedio que crear una transformación nueva.

En la documentación de EPSG te indican que si quieres añadir información, utilices siempre códigos con valores superiores a 32767, por lo tanto vamos a crear una transformación con código 32768 que nos servirá para transformar entre los sistemas de coordenadas geográficos ED50 (4230) y ETRS89 (4326), para la zona de España con el algoritmo 9615 y utilizando como archivo de rejilla NTv2 el archivo “PENR2009.gsb”.

Para ello ejecutaremos la consulta:

INSERT INTO Coordinate_Operation (COORD_OP_CODE, COORD_OP_NAME, COORD_OP_TYPE, SOURCE_CRS_CODE, TARGET_CRS_CODE, COORD_TFM_VERSION, COORD_OP_VARIANT, AREA_OF_USE_CODE, COORD_OP_SCOPE, COORD_OP_ACCURACY, COORD_OP_METHOD_CODE, REMARKS, INFORMATION_SOURCE, DATA_SOURCE, REVISION_DATE, SHOW_OPERATION, DEPRECATED)
VALUES (32768, 'ED50 a ETRS89 (14)', 'transformation', 4230, 4258, 'IGN-Esp v3', 14, 3429, 'For applications to an accuracy of 10-15cm (95% confidence) for Spain mainland and about 4cm (95%) for Balearic Islands.',0.2,9615,'Reemplaza a ED50 to ETRS89(12)','Instituto Geográfico Nacional, www.ign.es', 'OGP', 2009-01-01, 1, 0);

Y ya puestos, y considerando que el sistema de coordenadas geográfico ETRS89 es idéntico al WGS84, vamos a crear una transformación para transformar entre los sistemas de coordenadas geográficos ED50 (4230) y WGS84 (4326), esta vez con código 32769 con la siguiente consulta:

INSERT INTO Coordinate_Operation (COORD_OP_CODE, COORD_OP_NAME, COORD_OP_TYPE, SOURCE_CRS_CODE, TARGET_CRS_CODE, COORD_TFM_VERSION, COORD_OP_VARIANT, AREA_OF_USE_CODE, COORD_OP_SCOPE, COORD_OP_ACCURACY, COORD_OP_METHOD_CODE, REMARKS, INFORMATION_SOURCE, DATA_SOURCE, REVISION_DATE, SHOW_OPERATION, DEPRECATED)
VALUES (32769, 'ED50 a WGS84 (43)', 'transformation', 4230, 4326, 'IGN-Esp v3', 43, 3429, 'For applications to an accuracy of 10-15cm (95% confidence) for Spain mainland and about 4cm (95%) for Balearic Islands.',0.2,9615,'Valor de los parámetros de la transformación ED50 a ETRS89 (14). Asume que ETRS89 y WGS 84 se consideran iguales dentro de la precisión de la transformación. Reemplaza ED50 to WGS 84 (41)','Instituto Geográfico Nacional, www.ign.es', 'OGP', 2009-01-01, 1, 0);

Tan solo nos queda añadir en la tabla ‘Coordinate_Operation Parameter Value’ los parámetros para estas dos transformaciones.

Para ello tan solo tenemos que ejecutar las consultas:

INSERT INTO [Coordinate_Operation Parameter Value]
VALUES (32768, 9615, 8656, NULL, 'PENR2009.gsb', NULL);

INSERT INTO [Coordinate_Operation Parameter Value]
VALUES (32769, 9615, 8656, NULL, 'PENR2009.gsb', NULL);

y ya lo tenemos. Ahora cuando Digi3D detecte que tiene que transformar entre los sistemas de coordenadas de referencia ED50 y ETRS89 o ED50 y WGS 84 nos mostrará esta transformacion en el cuadro de diálogo de posibles transformaciones.