Los sistemas de referencia de coordenadas (CRS) o sistemas de referencia espacial (SRS) son complicados, al menos para mí que no soy geógrafo y tengo un conocimiento superficial de GIS: existen más de 13.000 sistemas diferentes y siempre me cuesta saber en qué sistema están las ubicaciones de un conjunto de datos georreferenciados. Siempre he tenido la intuición de que es algo complejísimo, sin saber exactamente a qué nivel. Solo cuando leí este maravilloso post de 2011 sobre el tema (que por supuesto no hace falta leer para conseguir cambiar de CRS un conjunto de datos), empecé a tomar consciencia de la magnitud de la complejidad: el post explica la incapacidad de las distintas herramientas de transformación para conseguir precisión en las conversiones, dando resultados que pueden diferir en torno a 20km. Apasionante lectura también para personas interesadas en saber curiosidades como de dónde vienen los códigos EPSG, hoy uno de los estándares para sistemas de referencia de coordenadas.
Pero voy a tema. A mí que trabajo casi siempre desarrollando para web me interesa en general que los datos geográficos usen latitud y longitud. Es el sistema de coordenadas que se suele usar en web y en aplicaciones GPS: Google Maps, Bing Maps, OpenStreetMap usan latitud y longitud. Si los datos que recibo usan un sistema cartesiano, tengo que transformarlos.
Detectar el sistema de referencia de coordenadas de un conjunto de datos
Lo primero que hay que hacer es saber en qué sistema de referencia están. Esto puede no ser fácil. Para ello hay que mirar los meta-datos del conjunto de datos y no todos los formatos guardan información sobre el CRS. Cuando no sé en qué sistema están los datos, uso QGIS: abro el archivo con los datos y luego lo exporto a geojson. Este formato guarda meta-datos sobre el CRS. Por ejemplo:
{
"type": "FeatureCollection",
"name": "projects",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3857" } },
"features": []
}
Otra forma de detectar el CRS de un conjunto de datos es usar el programa de línea de comandos projinfo, que forma parte del paquete de software PROJ. Sin embargo projinfo únicamente detecta el CRS de archivos .prj (información geoespacial de los archivos shape) o .geotiff, hasta donde he podido leer en la documentación. Así que es QGIS el que me suele ayudar en esta tarea.
Cambiar el sistema de referencia de coordenadas de un conjunto de datos
Hay multitud de maneras de cambiar el CRS de un conjunto de datos. Yo busco herramientas que permitan automatizar la tarea dentro de los procesos de limpieza de datos que hago. Para este tipo de procesos, los programitas que integran el paquete de software PROJ me funcionan de maravilla.
PROJ integra librerías para varios lenguajes de programación. Yo he usado dos: cs2cs en la línea de comandos, y pyproj en scripts de Python.
Transformar EPSG:3857 (Web Mercator) a EPSG:4326 (WGS 84) con cs2cs
EPSG:4386 es el sistema que usan la mayoría de cartografías web como OpenStreetMap o Google Maps. Para transformar datos en otro CRS, por ejemplo EPSG:3857, a EPSG:4326 con cs2cs en la línea de comandos:
echo "20037508.34 20037508.34 100" | cs2cs -f "%.6f" EPSG:3857 EPSG:4326
El modificador -f permite especificar el formato de salida tras la conversión. Si no se usa -f, cs2cs devuelve por omisión grados, minutos y segundos. Para obtener grados decimales, que es lo que yo necesito para web, hay que añadir -f "%.6f" donde 6 es el número de cifras decimales.
cs2cs devuelve tres coordenadas: x y z. En proyecciones planas la tercera es siempre 0, por supuesto:
85.051129 180.000000 0.000000
cs2cs puede transformar un conjunto de coordenadas si se le suministra un archivo con pares por líneas. Por ejemplo, el siguiente listado en EPSG:3857 alojado en el archivo coords.list:
echo coords.list 231362.441622 5082306.625994 -97756.310090 5108535.370440 -674547.189064 4487249.462584 -220554.193010 5360852.515295 -36260.920048 5255629.363755 245097.478426 5071621.614563 -455214.864321 4964279.036953 -53978.135458 4628117.767552
Se puede transformar masivamente de la siguiente manera:
cs2cs -f "%.6f" EPSG:3857 EPSG:4326 coords.list 41.472631 2.078364 0.000000 41.648932 -0.878160 0.000000 37.344603 -6.059560 0.000000 43.320300 -1.981272 0.000000 42.628723 -0.325737 0.000000 41.400673 2.201748 0.000000 40.673346 -4.089265 0.000000 38.343871 -0.484894 0.000000
Transformar EPSG:3857 (Web Mercator) a EPSG:4326 (WGS 84) con pyproj en Python
Para transformaciones de CRS dentro de un proceso de limpieza o estructuración de datos más amplio, quizás tenga más sentido usar una librería como pyproj en Python. La transformación en Python se podría hacer de la siguiente manera, suponiendo que el conjunto de datos esté en formato CSV:
import pandas as pd
from pyproj import Transformer
# define la conversión de EPSG:3857 a EPSG:4326
transformer = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)
# variables para el dataset de origen y el transformado
input_csv = "input.csv"
output_csv = "output.csv"
# abre el archivo de origen
df = pd.read_csv(input_csv)
# variables para las columnas con la información geográfica
x_col = "longitude"
y_col = "latitude"
# hace la conversión
df[[x_col, y_col]] = df.apply(lambda row: transformer.transform(row[x_col], row[y_col]), axis=1, result_type="expand")
# guarda los datos transformados
df.to_csv(output_csv, index=False)