Script en Python para convertir Lat Long a ITM (Irish Transverse Mercator)

Question

Actualmente estoy trabajando en un script de python que convierte Lat Long (WGS84) a ITM. Las fórmulas utilizadas se basan en un documento publicado por OSGB (Apéndice B):

http://www.leica-geosystems.co.uk/downloads123/gb/gps/gps1200/other/OS%20Transformations%20and%20OSGM02%20user%20guide_en.pdf .

Obtengo buenos resultados en los Este pero no en los Norte (desde unos pocos mm hasta aproximadamente 1,5m) y no puedo entender por qué. ¿Alguien sabe las razones? Gracias por vuestra ayuda.

Andrea

# this script convert Latitude and Longitude to ITM

# import modules
import math

# GRS80 Ellipsoid constants
a = 6378137.0
b = 6356752.3141
e2 = (a**2-b**2)/a**2

# ITM projection constants
F = 0.99982
Lat0 = 0.933751149817
Long0 = -0.13962634016
E0 = 600000.0
N0 = 750000.0

# input Lat degrees, minutes and seconds for latitude
dLat = float(raw_input('Enter latitude degrees: '))
mLat = float(raw_input('Enter latitude minutes: '))
sLat = float(raw_input('Enter latitude seconds: '))

# convert Lat to decimal degrees
ddLat = dLat + (mLat/60) + (sLat/3600)
print 'Latitude dd:', ddLat

# convert Lat to radians
rLat = (ddLat*math.pi)/180.0
Lat = rLat
print 'Latitude r:', Lat
print

# input Long degrees, minutes and seconds for longitude
dLong = float(raw_input('Enter longitude degrees: '))  #The degrees have to be entered with a posite sign
mLong = float(raw_input('Enter longitude minutes: '))
sLong = float(raw_input('Enter longitude seconds: '))

# convert Long to decimal degrees
ddLong = dLong + (mLong/60) + (sLong/3600)
print 'Longitude dd:', ddLong

# convert Long to radians
rLong = -(ddLong*math.pi)/180.0 # the negative sign account for west direction
Long = rLong
print 'Longitude r:', Long
print

# calculate constants for transformation
n = (a - b)/(a + b)
v = a*F*(1-e2*(math.sin(Lat)**2))**-0.5
p = a*F*(1-e2)*(1-e2*(math.sin(Lat)**2))**-1.5
n2 = v/p-1
print 'n:', n
print 'v:', v
print 'p:', p
print 'n2:', n2
print

M1 = (1+n+5.0/4.0*n**2+5.0/4.0*n**3)*(Lat-Lat0)
M2 = (3*n+3*n**2+21.0/8.0*n**3)*math.sin(Lat-Lat0)*math.cos(Lat+Lat0)
M3 = (15.0/8.0*n**2+15.0/8.0*n**3)*math.sin(2*(Lat-Lat0))*math.cos(2*(Lat+Lat0))
M4 = 35.0/24.0*n**3*math.sin(3*(Lat-Lat0))*math.cos(3*(Lat+Lat0))
M = b*F*(M1-M2+M3-M4)
print 'M1:', M1
print 'M2:', M2
print 'M3:', M3
print 'M4:', M4
print 'M:', M
print

I = M+N0
print 'I:', I
II = v/2*math.sin(Lat)*math.cos(Lat)
print 'II:', II
III = v/24*math.sin(Lat)*(math.cos(Lat))**3*(5-(math.tan(Lat)**2)+9*n2)
print 'III:', III
IIIA = v/720*math.sin(Lat)*(math.cos(Lat)**5)*(61-58*(math.tan(Lat)**2)+(math.tan(Lat)**4))
print 'IIIA:', IIIA
IV = v*math.cos(Lat)
print 'IV:', IV
V = v/6*(math.cos(Lat)**3)*(v/p-(math.tan(Lat)**2))
print 'V:', V
VI = v/120*(math.cos(Lat)**5)*(5-18*(math.tan(Lat)**2)+(math.tan(Lat)**4)+14*n2-58*(math.tan(Lat)**2)*n2)
print 'VI:', VI
print

# calculate Eastings and Northings
N = I+II*(Long-Long0)**2+III*(Long-Long0)**4+IIIA*(Long-Long0)**6
E = E0+IV*(Long-Long0)+V*(Long-Long0)**3+VI*(Long-Long0)**5
print 'Easting: ', E
print 'Northing:', N

Answer 1

Como om_henners aconsejó que es mejor utilizar la biblioteca disponible para este propósito, ya que está implementado y probado por muchas personas ...

Por lo tanto, eche un vistazo a pyproj Librería Python.

A continuación se muestra un código de ejemplo para reproyectar WGS-84 long/lat a ITM ( EPSG:2157 ) x,y:

from pyproj import Proj, transform

def reproject_wgs_to_itm(longitude, latitude):
    prj_wgs = Proj(init='epsg:4326')
    prj_itm = Proj(init='epsg:2157')
    x, y = transform(prj_wgs, prj_itm, longitude, latitude)
    return x, y

print reproject_wgs_to_itm(-7.748108, 53.431628)

Answer 2

Sólo soy un novato así que no puedo comentar nada pero espero que mi respuesta le sirva a alguien. Necesitaba convertir WGS84 Latitud, Longitud a ITM para un datalogger inteligente que estaba construyendo así que convertí su código de python a C++. Usé Arduino IDE para compilar y probar el código que contenía 8 ejemplos alrededor y dentro de Irlanda. Funcionó bien; la orientación estaba dentro de 0,06 m y el norte estaba dentro de 0,41 m de los valores ITM obtenidos por GridInQuest. Aquí está el código C++ con un ejemplo para ahorrar espacio.

#include <math.h>
#define pi 3.14159265358979323846

// GRS80 Ellipsoid constants;
double a = 6378137.0;
double b = 6356752.3141;
double e2 = (pow(a, 2) - pow(b, 2)) / pow(a, 2);
// ITM projection constants;
double F = 0.99982;
double Lat0 = 0.933751149817;
double Long0 = -0.13962634016;
double E0 = 600000.0;
double N0 = 750000.0;
double dLat, mLat, dLong, mLong = 0;
double N, E = 0;

void setup() {
  // put your setup code here, to run once:;
  Serial.begin(57600);
}

void loop() {
  // Example with degrees, minutes and decimal minutes for latitude and longitude;
  dLat = 54; mLat = 37;  //latitude degrees, minutes
  dLong = 6; mLong = 11; //longitude degrees, minutes (use positive for degrees west of Greenwich)
  latLong2ITM(dLat, mLat, dLong, mLong); // Convert to ITM

  while (-1); //Wait here forever
}

// ==============================================================================
// FUNCTIONS
// This function converts Latitude and Longitude to ITM;
void latLong2ITM (double dLat, double mLat, double dLong, double mLong) {

  // convert Lat to decimal degrees;
  double ddLat = dLat + (mLat / 60);
  // convert Lat to radians;
  double rLat = (ddLat * pi) / 180.0;
  double Lat = rLat;
  // convert Long to decimal degrees;
  double ddLong = dLong + (mLong / 60);
  // convert Long to radians;
  double rLong = -(ddLong * pi) / 180.0; // the negative sign account for west direction;
  double Long = rLong;

  // calculate constants for transformation;
  double n = (a - b) / (a + b);
  double v = a * F * pow((1 - e2 * (pow(sin(Lat), 2))), -0.5);
  double p = a * F * (1 - e2) * pow((1 - e2 * pow(sin(Lat), 2)), -1.5);
  double n2 = v / p - 1;

  double M1 = (1 + n + 5.0 / 4.0 * pow(n, 2) + 5.0 / 4.0 * pow(n, 3)) * (Lat - Lat0);
  double M2 = (3 * n + 3 * pow(n, 2) + 21.0 / 8.0 * pow(n, 3)) * sin(Lat - Lat0) * cos(Lat + Lat0);
  double M3 = (15.0 / 8.0 * pow(n, 2) + 15.0 / 8.0 * pow(n, 3)) * sin(2 * (Lat - Lat0)) * cos(2 * (Lat + Lat0));
  double M4 = 35.0 / 24.0 * pow(n, 3) * sin(3 * (Lat - Lat0)) * cos(3 * (Lat + Lat0));
  double M = b * F * (M1 - M2 + M3 - M4);

  double I = M + N0;
  double II = v / 2 * sin(Lat) * cos(Lat);
  double III = v / 24 * sin(Lat) * pow(cos(Lat), 3) * (5 - pow(tan(Lat), 2) + 9 * n2);
  double IIIA = v / 720 * sin(Lat) * pow(cos(Lat), 5) * (61 - 58 * pow(tan(Lat), 2) + pow(tan(Lat), 4));
  double IV = v * cos(Lat);
  double V = v / 6 * pow(cos(Lat), 3) * (v / p - pow(tan(Lat), 2));
  double VI = v / 120 * pow(cos(Lat), 5) * (5 - 18 * pow(tan(Lat), 2) + pow(tan(Lat), 4) + 14 * n2 - 58 * pow(tan(Lat), 2) * n2);

  // calculate Eastings and Northings;
  N = I + II * pow((Long - Long0), 2) + III * pow((Long - Long0), 4) + IIIA * pow((Long - Long0), 6);
  E = E0 + IV * (Long - Long0) + V * pow((Long - Long0), 3) + VI * pow((Long - Long0), 5);
  Serial.println();
  Serial.print("Easting: "); Serial.println(E);
  Serial.print("Northing:"); Serial.println(N);
}

Answer 3

En caso de que alguien todavía esté interesado en una conversión pura de Python Irish Transverse Mercator hacia adelante,

from math import *

# Meridian Arc    
def arcmer(a,equad,lat1,lat2):
    b=a*sqrt(1-equad)
    n=(a-b)/(a+b)
    a0=1.+((n**2)/4.)+((n**4)/64.)
    a2=(3./2.)*(n-((n**3)/8.))
    a4=(15./16.)*((n**2)-((n**4)/4.))
    a6=(35./48.)*(n**3)

    s1=a/(1+n)*(a0*lat1-a2*sin(2.*lat1)+a4*sin(4.*lat1)-a6*sin(6.*lat1))
    s2=a/(1+n)*(a0*lat2-a2*sin(2.*lat2)+a4*sin(4.*lat2)-a6*sin(6.*lat2))
    return s2-s1

# Gauss-Kruger Projection
def geogauss(lat,lon,a,equad,lat0,lon0):

    lat0=radians(lat0)
    lon0=radians(lon0)

    lat=radians(lat)
    lon=radians(lon)

    lon=lon-lon0

    N=a/sqrt(1-equad*(sin(lat))**2)
    RO=a*(1-equad)/((1-equad*(sin(lat)**2))**(3./2.))

    k1=(N/RO)+(4.*(N**2)/(RO**2))-((tan(lat))**2)

    k2=(N/RO)-((tan(lat))**2)

    k3=N/RO*(14.-58.*((tan(lat)))**2)+40.*((tan(lat))**2)+((tan(lat))**4)-9.

    x=lon*N*cos(lat)+(lon**3)/6.*N*((cos(lat))**3)*k2+(lon**5)/120.*N*((cos(lat))**5)*k3

    y=arcmer(a,equad,lat0,lat)+(lon**2)/2.*N*sin(lat)*cos(lat)+((lon**4)/24.)*N*sin(lat)*((cos(lat))**3)*k1

    return x,y

# Irish Transverse Mercator
def itm(lat,lon):
    # GRS-80
    a = 6378137.
    equad =0.00669437999

    # Natural Origin 
    lat0=53.5
    lon0=-8.

    coords = geogauss(lat,lon,a,equad,lat0,lon0)

    k0=0.999820
    x = (k0*coords[0])+600000.
    y = (k0*coords[1])+750000.
    return x,y

# Test values    
lat=53.5
lon=-8.

xy = itm(lat,lon)

print "x= %.16f" %xy[0]
print "y= %.16f" %xy[1]