WGS84 vs ETRS89 vs ED50 vs Madrid 1870

En esta entrada voy a comparar los sistemas de coordenadas WSG84, ETRS89, ED50 y el vetustísimo Madrid 1870. Además, lo voy a hacer mal y luego voy a explicar no solo por qué sino por qué no es culpa mía.

Primero, las coordenadas de Sol (el Kilómetro 0, para ser más precisos) en WGS84 (EPSG:4326):

library(sf)
options(digits = 10)
sol_wsg84 <- st_sfc(st_point(
    c(40.416634493768065, -3.703811417868093)),
    crs = 4326)
st_coordinates(sol_wsg84)
#             X            Y
# 1 40.41663449 -3.703811418

Ahora, en ETRS89 (EPSG:4258):

sol_etrs <- st_transform(sol_wsg84, crs = 4258)
st_coordinates(sol_etrs)
#             X            Y
# 1 40.41663449 -3.703811418

Ahora en el caduco ED50 (EPSG:4230):

sol_ed50 <- st_transform(sol_wsg84, crs = 4230)
st_coordinates(sol_ed50)
#             X            Y
# 1 40.41663449 -3.703811418

Finalmente, en el vetusto Madrid 1870 (EPSG:4903):

sol_mad1870 <- st_transform(sol_wsg84, crs = 4903)
st_coordinates(sol_mad1870)
#             X            Y
# 1 44.10457338 -3.703811418

Como puede observarse, en los tres primeros sistemas de referencia, las coordenadas del punto son iguales (¡pero no, no es cierto; en realidad son distintos!) y es sorprendente lo muy al norte que se estuvo considerando Madrid durante muchos años (si es que sf_transform es de fiar).

¿Por qué coinciden las coordenadas? sf delega los cálculos en la librería PROJ (libproj en Linux) y…

$ projinfo -o PROJ -s EPSG:4326 -t EPSG:4258
Candidate operations found: 1
-------------------------------------
Operation No. 1:

INVERSE(EPSG):1149, Inverse of ETRS89 to WGS 84 (1), 1.0 m, Europe - ETRS89

PROJ string:
+proj=noop

En resumen, que PROJ decide que el cambio de sistema de referencia entre WSG84 y ETRS89 es la identidad. Lo mismo sucede con ED50. Lo cual no es ni cierto ni satisfactorio pero es la mentira piadosa en la que R + sf + PROJ nos hacen vivir.