[Topografia] Realizzare un programma di passaggio da coordinate geografiche WGS84 a coordinate cartografiche UTM

[Jack871]

Ciao.

Come da titolo devo scrivere il sorgente di un programma che converta delle coordinate geografiche riferite al sistema WGS84 in coordinate cartesiane nella rappresentazione UTM-WGS84.

Una volta completato, dovrò testare il programma convertendo una serie di punti in coordinate geografiche e confrontare i risultati ottenuti con quelli restituiti dal software CartLab.


Per prima cosa devo avere ben chiaro il procedimento con cui avviene tale trasposizione.


Le formule per la conversione sono queste (correggetemi se sbaglio):

${ (E = mu*x(varphi, lambda-lambda_0) + E_0),(N = mu*y(varphi, lambda-lambda_0)) :}$

dove:

[*:1h9fyo6j]$(varphi, lambda)$ sono le coordinate geografiche del punto (note)[/*:m:1h9fyo6j]
[*:1h9fyo6j]$(E, N)$ sono le coordinate cartografiche del punto[/*:m:1h9fyo6j]
[*:1h9fyo6j]$mu = 0,9996$ è il fattore di contrazione[/*:m:1h9fyo6j]
[*:1h9fyo6j]$x(varphi, lambda), y(varphi, lambda)$ sono le formule della proiezione di Gauss[/*:m:1h9fyo6j]
[*:1h9fyo6j]$lambda_0$ è la longitudine del meridiano di riferimento del fuso[/*:m:1h9fyo6j]
[*:1h9fyo6j]$E_0 = 500$ è la falsa origine (espressa in km)[/*:m:1h9fyo6j][/list:u:1h9fyo6j]

Le formule di Gauss sono queste:




Prima di imbarcarmi nella stesura del programma vero e proprio, vorrei convertire "a mano" le coordinate di almeno un punto.


Considero allora il primo punto dato in coordinate geografiche (WGS84):


[size=150]45°48'49.58500"N 13°25'54.27733"E[/size]


Dalla longitudine deduco che il punto si trova nel fuso 33 della rappresentazione UTM, il cui meridiano centrale ha longitudine $lambda_0 = 15°$.


Inizio convertendo gli angoli dati da sessagesimali a sessadecimali e da sessadecimali a radianti:

$45°+(48/60)°+(49.58500/3600)° = 45°.81377$

$13°+(25/60)°+(54.27733/3600)° = 13°.43174$

$varphi = {45.81377 * pi}/180 = 0.79960$

$lambda = {13.43174 * pi}/180 = 0.46886$


Riporto i valori dell'ellissoide WGS84 che mi occorrono:

$a = 6378137$ semisasse maggiore (espresso in metri)
$e^2 = 0.006694379990$ eccentricità prima al quadrato


Calcolo il raggio di curvatura della sezione normale "gran-normale":

$N = a / sqrt(1-e^2*sin^2 varphi) = 6378141.15764$


Calcolo gli altri valori ricorrenti nelle formule di Gauss:

$eta^2 = e^2/{1-e^2}*cos^2 varphi = 0.0067395$

$t = tan varphi = 0.013957$

$lambda - lambda_0 = 0.26180$


La formula della coordinata y nella proiezione di Gauss richiede il calcolo di $l_varphi = int_0^varphi rho dvarphi$ (lunghezza dell'arco di meridiano dall'equatore alla latitudine $varphi$), ossia il calcolo del seguente integrale:

$int {a(1-e^2)}/(1-e^2 sin^2 varphi)^{3/2} dvarphi$

che mi mette in difficoltà... Per il momento allora mi limito a calcolare la coordinata Est.


Ricavo il valore della coordinata x dalle formule di Gauss (ho troncato lo sviluppo della serie al termine $lambda^3$):

$x = (lambda-lambda_0) N cos varphi + 1/6 (lambda-lambda_0)^3 N cos^3 varphi (1-t^2+eta^2) = 1669634.75264$


Infine ricavo la coordinata Est con le formule di trasformazione:

$E = mu*x(varphi, lambda-lambda_0) + E_0 = 2168966,89874$


Adesso è arrivato il momento di confrontare il risultato trovato con quello prodotto da CartLab:



Ahimè siamo ben oltre i 5 metri di margine concessi al programma che devo realizzare relativamente ai risultati restituiti da CartLab.


Molto probabilmente avrò commesso degli errori nei calcoli, spero almeno che il procedimento sia a grosse linee corretto.


Sarò eternamente grato a chiunque vorrà aiutarmi. Grazie!

[Jack871]

Ciao.

Vi ragguaglio sui progressi fatti. Purtroppo il programma non raggiunge ancora il suo scopo, ossia ci sono ancora delle discrepanze importanti tra i valori calcolati da me e quelli restituiti da CartLab.

Questo è il codice in JavaScript della funzione di trasformazione delle coordinate:

//semiasse maggiore ed eccentricita' al quadrato dell'ellissoide WGS84
a = 6378137;
e2 = 0.00669437999;

//fattore di contrazione e false origini
var k = 0.9996;
var est0 = 500000;
var nrd0 = 10000000;


function WGS84toUTM(p)
{
	//da angoli sessagesimali ad angoli sessadecimali
	lat = Number(p[0][0]) + Number(p[0][1])/60 + Number(p[0][2])/3600;
	lon = Number(p[1][0]) + Number(p[1][1])/60 + Number(p[1][2])/3600;

	//calcolo del fuso e della longitudine del rispettivo meridiano centrale
	fus = 1 + Math.floor((lon+180)/6);
	mer = 6*(fus-1) - 177;

	//calcolo di phi e lambda
	//phi = latitudine in radiandi
	//lmb = longitudine meno longitudine meridiano centrale (in radianti)
	phi = lat*Math.PI/180;
	lmb = (lon-mer)*Math.PI/180;

	//calcolo del raggio di curvatura della sezione normale "gran-normale"
	N = a/Math.sqrt(1-e2*Math.pow(Math.sin(phi),2));

	//calcolo di valori ricorrenti negli sviluppi
	//eta al quadrato e tangente di phi al quadrato
	eta2 = e2/(1-e2)*Math.pow(Math.cos(phi),2);
	t2 = Math.pow(Math.tan(phi),2);

	//calcolo della lunghezza dell'arco di meridiano
	b0 = phi*(1 - 1/4*e2 - 3/64*Math.pow(e2,2) - 5/256*Math.pow(e2,3));
	b0 = b0 - Math.sin(2*phi)*(3/8*e2 - 3/32*Math.pow(e2,2) - 45/1024*Math.pow(e2,3));
	b0 = b0 + Math.sin(4*phi)*(15/256*Math.pow(e2,2) + 45/1024*Math.pow(e2,3));
	b0 = b0 - Math.sin(6*phi)*(35/3072*Math.pow(e2,3));
	b0 = a*b0;

	//calcolo degli altri coefficienti degli sviluppi
	b1 = N*Math.cos(phi);
	b2 = 1/2*N*Math.sin(phi)*Math.cos(phi);
	b3 = 1/6*N*Math.pow(Math.cos(phi),3)*(1-t2+eta2);
	b4 = 1/24*N*Math.sin(phi)*Math.pow(Math.cos(phi),3)*(5-t2+9*eta2+4*Math.pow(eta2,2));
	b5 = 1/120*N*Math.pow(Math.cos(phi),5)*(5-18*t2+Math.pow(t2,2)+14*eta2-58*t2*eta2);
	b6 = 1/720*N*Math.sin(phi)*Math.pow(Math.cos(phi),5)*(61-58*t2+Math.pow(t2,2)+270*eta2-330*t2*eta2);

	//calcolo delle coordinate cartografiche
	est = k*(b1*lmb+b3*Math.pow(lmb,3)+b5*Math.pow(lmb,5))+est0;
	nrd = k*(b0+b2*Math.pow(lmb,2)+b4*Math.pow(lmb,4)+b6*Math.pow(lmb,6));
	if(nrd<0) nrd = nrd+nrd0;

	return {est:est, nrd:nrd};
}

Qui ci sono i risultati prodotti dal codice rapportati con quelli generati da CartLab.



Come potete osservare per la coordinata Nord lo scarto si attesta sui 5.8 km (ancora troppo elevato), mentre per la coordinata Est lo scarto è di 673 km.


Per il calcolo dell'arco di meridiano ho usato il seguente sviluppo:




Grazie.