Iker Perez de Albeniz

En el post Creando aplicación de Bizkaibus para Android (Parte 3) comentaba que había conseguido convertir coordenadas en formato UTM a Geográficas (latitud, longitud), pero que seguía teniendo el problema de la conversión de Batum, es decir, los datos geográficos los tenia en el sistema ED50 y tenia que pasarlos a WGS84 que es el sistema usado por Google Earth y actualmente, el estándar de facto. Tras muchos buscar e indagar, vi que la mejor opción era usar un sistema de “Rejilla” es decir, hacer la transformación usando una tabla que en función de la posición me da un coeficiente de compensación. Es decir, una tabla donde busco mi posición y me dice cuanto tengo que restar o sumar a mi posición para que sea el equivalente en otro sistema.

En mi aso necesitaba una rejilla de transformación de ED50 a WGS84.. y encontré la rejilla publicada por el Centro Nacional de Información Geográfica que permitía la transformación de ED50 a ETRS89 que es el actual sistema oficial en España y prácticamente igual que WGS84.

Ya solo necesitaba saber como utilizar dicha rejilla (la cual tiene un formato denominado NTv2). Primero analicé la rejilla y pude intuir su estructura interna ya que con cualquier editor hexadecimal (y si has hecho un poco de ingeniería inversa alguna vez) se podía ver que los datos estaban bastante acotados en 8 y 4 bytes. Así que lo primero que hice es una librería en c# que convertía ese archivo rejilla en una clase Rejilla manejable por código. Mas tarde encontré información sobre dicho formato, por lo que pude confirmar que mi aproximación era correcta.

Una vez convertida la rejilla a una clase solo necesitaba saber como hacer los cálculos a partir de la información de la rejilla para convertir las posiciones geográficas de ED50 a WGS84.

Googleando acabe en el blog referente en temas de GIS en castellano que es el blog de Gabriel Ortiz y en su foro pregunte por información acerca de el algoritmo de conversión. El propio Gabriel me recomendó documentación y ejemplos de código en VisualBasic y Matlab hechos por Javier González Matesanz perteneciente al Instituto Geográfico. A partir de dichos ejemplos pude implementar un método de conversión para mi librería .net.

Aun la librería para manejar archivo rejilla (NTv2) no esta completa, pero de cara a corregir los datos geográficos generados por la API de OpenBizkaibus hice una implementación en python de dicho algoritmo. Además de cara a no tener que convertir el archivo rejilla en una estructura de clases, directamente convertir la rejilla a dicha estructura y la guarde como un archivo sped2et.py como una variable sped2et_file importable, ahorrando un montón de tiempo de ejecución.

El archivo sped2et.py se puede encontrar aquí.

Y el ejemplo del algoritmo de conversión y como usar sped2et.py se muestra a continuación.

[code lang="python"]
from sped2et import sped2et_file
import math

def convertWithGrid(longitude, latitude):

numpages = len(sped2et_file['PAGES'])
M = []
n = []
i = []
j = []
p = []
sublist = []

msel = 1000000
longitude = -longitude * 3600
latitude = latitude * 3600
pagec = 0
selpage = 0
for page in sped2et_file['PAGES']:
M.append(1 + round((sped2et_file['PAGES'][pagec]['N_LAT'] - sped2et_file['PAGES'][pagec]['S_LAT'])/sped2et_file['PAGES'][pagec]['LAT_INC']))
n.append(1 + round((sped2et_file['PAGES'][pagec]['W_LONG'] - sped2et_file['PAGES'][pagec]['E_LONG']) / sped2et_file['PAGES'][pagec]['LONG_INC']))

aux1 = 1 +((latitude - sped2et_file['PAGES'][pagec]['S_LAT']) / sped2et_file['PAGES'][pagec]['LAT_INC'])
if aux1 >= 0:
i.append(math.floor(aux1))
else:
i.append(math.ceil(aux1))

aux2 = 1 +((longitude - sped2et_file['PAGES'][pagec]['E_LONG']) / sped2et_file['PAGES'][pagec]['LONG_INC'])
if aux2 >= 0:
j.append(math.floor(aux2))
else:
j.append(math.ceil(aux2))

if i[pagec] > 0 and j[pagec] > 0 and i[pagec] < M[pagec] and j[pagec] < n[pagec]:
if sped2et_file['PAGES'][pagec]['LAT_INC'] < msel:
msel = sped2et_file['PAGES'][pagec]['LAT_INC']
selpage = pagec
pagec = pagec + 1

p.append(n[selpage] * (i[selpage] - 1) + j[selpage])
p.append(p[0] + 1)
p.append(p[0] + n[selpage])
p.append(p[2] + 1)

for point in p:
sublist.append(sped2et_file['PAGES'][selpage]['NODES'][int(point)])

lataux = sped2et_file['PAGES'][selpage]['S_LAT'] + (i[selpage] - 1) * sped2et_file['PAGES'][selpage]['LAT_INC']
lonaux = sped2et_file['PAGES'][selpage]['E_LONG'] + (j[selpage] - 1) * sped2et_file['PAGES'][selpage]['LONG_INC']
y = (latitude - lataux) / sped2et_file['PAGES'][selpage]['LAT_INC']
x = (longitude - lonaux) / sped2et_file['PAGES'][selpage]['LONG_INC']

#COEFICIENTES
a0 = sublist[0]['iLat'];
a1 = sublist[1]['iLat']- sublist[0]['iLat']
a2 = sublist[2]['iLat'] - sublist[0]['iLat']
a3 = sublist[0]['iLat'] + sublist[3]['iLat'] - sublist[1]['iLat'] - sublist[2]['iLat']
ip = a0 + a1 * x + a2 * y + a3 * x * y;
olatitude = latitude + ip;

b0 = sublist[0]['iLong'];
b1 = sublist[1]['iLong'] - sublist[0]['iLong'];
b2 = sublist[2]['iLong'] - sublist[0]['iLong'];
b3 = sublist[0]['iLong'] + sublist[3]['iLong'] - sublist[1]['iLong'] - sublist[2]['iLong']
ib = b0 + b1 * x + b2 * y + b3 * x * y;
olongitude = longitude + ib;

olongitude = -olongitude / 3600;
olatitude = olatitude / 3600;
return (olongitude, olatitude)

longitude = -2.51087447913458
latitude = 42.8368613892524
(olongitude, olatitude) = convertWithGrid(longitude, latitude)
print "%s,%s --> %s,%s" % (latitude,longitude,olatitude,olongitude)
[/code]

En cuanto este lista también haré pública la librería .Net.