Ho cercato di approfondire un pochino l'argomento quando mi è riuscito di mettere le mani sul cellulare di mia moglie (prima o poi mi stronca
), purtroppo devo dire che è difficilissimo realizzare una conversione "da se", troppi i dati occorrenti, complessi i metodi da applicarsi, tieni presente che nell'ambito dello stesso sistema bisogna eseguire delle roto-traslazioni per compensare le deformazioni del geoide di riferimento, Ti invito a leggere, giusto per informazione generale,
questo pdf, approfondendo ulteriormente si scopre che è un vero e proprio macello.
Allo stato suggerirei di tenerti le coordinate wgs84 che hai o di utilizzare il Vortex dello I.G.M. se Ti è proprio necessario, tieni presente che vi sarebbero comunque lievi errori di approssimazione rispetto alla cartografia.
Comunque, ho provato a fare una sbrigativa traduzione del codice javascript contenuto nello stralcio di Computer programming indicato in precedente post, non Ti serve perché funziona solo per conversione wgs84 => ed50/utm, il test indicato lo supera ma dubito molto che possa, comunque, essere utile a livello generalizzato per la cartografia europea (basata sul sistema ED50) ma solo per il caso specifico posto (per dette roto-traslazioni), pur se praticamente inutile, posto la traduzione
# -*- coding: utf-8 -*-
import math
class Orientation:
def __init__(self, x0, y0, z0, K, Ex, Ey, Ez):
self.z0 = z0
self.x0 = x0
self.y0 = y0
self.K = K
self.Ex = Ex
self.Ey = Ey
self.Ez = Ez
class RotationElipsoid:
def __init__(self, a, alpha, orientation):
self.orientation = orientation
self.a = a
self.alpha = alpha
semiminor = self.a - self.a * self.alpha
self.e2 = (self.a**2 - semiminor**2) / self.a**2
# RHO - Raggio di curvatura del meridiano
def RHO (self, LAMBDA):
return (self.a * (1 - self.e2) /
math.sqrt((1 - self.e2 * math.sin(LAMBDA))**2)**3)
# N_G - Gran Normale
def N_G(self, LAMBDA):
return self.a / math.sqrt(1 - self.e2 * math.sin(LAMBDA)**2);
def molodenskij(fromElipsoid, destElipsoid, PHI, LAMBDA, H):
alpha = fromElipsoid.alpha
a = fromElipsoid .a
Da = fromElipsoid.a - destElipsoid.a
Dalpha = fromElipsoid.alpha - destElipsoid.alpha
DEx = destElipsoid.orientation.Ex - fromElipsoid.orientation.Ex
DEy = destElipsoid.orientation.Ey - fromElipsoid.orientation.Ey
DEz = destElipsoid.orientation.Ez - fromElipsoid.orientation.Ez
Dx0 = destElipsoid.orientation.x0 - fromElipsoid.orientation.x0
Dy0 = destElipsoid.orientation.y0 - fromElipsoid.orientation.y0
Dz0 = destElipsoid.orientation.z0 - fromElipsoid.orientation.z0
DK = destElipsoid.orientation.K - destElipsoid.orientation.K
RHO = fromElipsoid.RHO(LAMBDA)
N_G = fromElipsoid.N_G(LAMBDA)
# scorciatoie
sinPHI = math.sin(PHI)
cosPHI = math.cos(PHI)
sinLAMBDA = math.sin(LAMBDA)
cosLAMBDA = math.cos(LAMBDA)
tanLAMBDA = math.tan(LAMBDA)
tempVar = RHO + H
tempVar2 = a**2 / N_G + H
tempVar6 = (1 - alpha)**2
tempVar5 = (1 - tempVar6) * N_G
tempVar3 = (N_G + H) * cosLAMBDA
tempVar4 = (tempVar6 * N_G + H) / (N_G + H)
o1 = (sinLAMBDA * cosPHI / tempVar) * Dx0
o2 = (sinLAMBDA * sinPHI / tempVar) * Dy0
o3 = (cosLAMBDA / tempVar) * Dz0
o4 = (tempVar2 / tempVar) * sinPHI * DEx
o5 = (tempVar2 / tempVar) * cosPHI * DEy
o6 = (tempVar5 / tempVar) * cosLAMBDA * sinLAMBDA * (DK + Da / a)
o7 = ((( RHO + tempVar6 * N_G) / tempVar) * cosLAMBDA * sinLAMBDA *
( Dalpha / (1 - alpha )))
dLAMBDA = o1 + o2 - o3 + o4 - o5 - o6 - o7
dPHI = (sinPHI / tempVar3 * Dx0 -cosPHI / tempVar3 * Dy0 - tempVar4 *
tanLAMBDA * cosPHI * DEx - tempVar4 * tanLAMBDA * sinPHI *
DEy + DEz)
dH = (-cosLAMBDA * cosPHI * Dx0 - cosLAMBDA * sinPHI * Dy0 - sinLAMBDA *
Dz0 + tempVar5 * sinLAMBDA * cosLAMBDA * sinPHI * DEx - tempVar5 *
sinLAMBDA * cosLAMBDA * cosPHI * DEy - tempVar2 * DK + a / N_G *
Da - tempVar6 * N_G * sinLAMBDA**2 *
(Dalpha / (1 - alpha)))
return {'LAMBDA' : LAMBDA + dLAMBDA, 'PHI': PHI + dPHI , 'H': H + dH}
# esegue l’integrale definito di di una funzione con
# il metodo dei trapezi in N parti
def Integral(myFunction, begin, end, parts):
step = (end - begin) / parts
temp = myFunction(begin) / 2
begin += step
for c in range(parts - 1):
temp += myFunction(begin)
begin += step
temp += myFunction(begin) / 2
return temp * step
def gauss(surface, m):
precision = 10000
PHI = m['PHI']
LAMBDA = m['LAMBDA']
H = m['H']
a = surface.a
e2 = surface.e2
# funzione da integrare
def Integral_Function(val):
tempVal = (a * (1 - e2))
return tempVal / math.sqrt((1 - e2 * math.sin(val)**2)**3)
# calcolo il fuso e sistemo la longitudine
fuse_length = 6 * math.pi / 180
if (PHI > math.pi):
fuse = math.floor((PHI - math.pi) / fuse_length)
else:
fuse = math.floor((PHI + math.pi) / fuse_length)
deltalength = (fuse - 30) * fuse_length + fuse_length / 2
dPHI = PHI - deltalength
# scorciatoie
sinLAMBDA = math.sin(LAMBDA)
cosLAMBDA = math.cos(LAMBDA)
tanLAMBDA = math.tan(LAMBDA)
tempVar1 = math.sqrt (1 - e2 * sinLAMBDA**2)
tempVar2 = e2 / math.sqrt(1 - e2)**2
valoreIntegrale = Integral(Integral_Function, 0, LAMBDA, precision)
x = (valoreIntegrale + (1/2 * a * sinLAMBDA * cosLAMBDA / tempVar1 ) *
dPHI**2 + (1/24 * a * sinLAMBDA * cosLAMBDA**3) /
tempVar1 * (5 - tanLAMBDA**2 + 9 * tempVar2 *
cosLAMBDA**2 + 4 * tempVar2**2 *
cosLAMBDA**4) * dPHI**4 + 1/720 * a *
sinLAMBDA * cosLAMBDA**5 / tempVar1 *
(61 - 58 * tanLAMBDA**2 + tanLAMBDA**4 + 270 * tempVar2 *
cosLAMBDA**2 - 330 * tempVar2 * cosLAMBDA**2 *
tanLAMBDA**2) * dPHI**6)
y = ((a * cosLAMBDA / tempVar1) * dPHI + (1/6) * a *
cosLAMBDA**3 / tempVar1 * (1 - tanLAMBDA**2 +
tempVar2 * cosLAMBDA**2) * dPHI**3 +
(1/120) * a * cosLAMBDA**5 / tempVar1 *
(5 - 18 * tanLAMBDA**2 + tanLAMBDA**4 + 14 * tempVar2 *
cosLAMBDA**2 - 58 * tempVar2 * cosLAMBDA**2 *
tanLAMBDA**2) * dPHI**5)
return [x,y,H];
# Ellissoidi di rotazione di riferimento
WGS84 = RotationElipsoid(6378137.0,
1 / 298.257222101,
Orientation(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
ED50 = RotationElipsoid(6378388,
1 / 297,
Orientation(-87.0, -98.0, -121.0, 0.0, 0.0, 0.0, 0.0))
def convert(longitudine, latitudine, quota):
''' Esegue la conversione completa da WGS84 a ED50 '''
m = molodenskij(WGS84, ED50, longitudine, latitudine, quota)
g = gauss(ED50, m)
# correzione coordinate per UTM
g[0] = g[0] * 0.9996
g[1] = g[1] * 0.9996 + 500000
return {'x': g[0], 'y': g[1], 'h': g[2]}
def to_radians(d, m, s):
''' converte gradi, primi, secondi in radianti '''
return math.pi / 180 * (d + m / 60 + s / 3600)
if __name__ == '__main__':
sel = 's'
while sel == 's':
ok = True
tmp = input('Inserire longitudine : ')
if not tmp:
ok = False
else:
values = tmp.split()
if len(values) == 1:
LAMBDA = to_radians(float(values[0]), 0.0, 0.0)
elif len(values) == 3:
LAMBDA = to_radians(*[float(x) for x in values])
else:
ok = False
tmp = input('Inserire latitudine : ')
if not tmp:
ok = False
else:
values = tmp.split()
if len(values) == 1:
PHI = to_radians(float(values[0]), 0.0, 0.0)
elif len(values) == 3:
PHI = to_radians(*[float(x) for x in values])
else:
ok = False
tmp = input('Inserire altitudine : ')
if not tmp:
ok = False
else:
H = float(tmp)
if ok:
converted = convert(PHI, LAMBDA, H)
print('Valori ED50-UTM :')
print(' - x :', round(converted['x'], 0))
print(' - y :', round(converted['y'], 0))
print(' - h :', round(converted['h'], 0))
else:
print('Dati non validi')
sel = input('\nPer altra conversione premere "s"')
sel = sel.lower()
Mi rincresce di non saper fare di più
Ciao
[Edit] Mi è venuto in mente adesso : dato che sei un user linux, potrebbe considerarsi l'utilizzo di "
cs2cs", sempre che Tu abbia le idee chiare in merito ai sistemi di riferimento da traslare (che poi sono le incognite più complesse)
