Geographic vs geocentric coordinates-2015-02-03-ita

Confronto fra Coordinate Geogra�che e

Geocentriche. Un nuovo algoritmo per il calcolo

inverso in forma chiusa basato su alcune

proprietà dell'ellisse.

Giuliano Curti

03/02/2015

Indice

1 Premessa 3

2 La geometria dell'ellisse 5

2.1 Una imprimitura . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1.1 Primo approccio (algoritmo del giardiniere) . . . . . . . . 52.1.2 Secondo approccio . . . . . . . . . . . . . . . . . . . . . . 72.1.3 Terzo approccio . . . . . . . . . . . . . . . . . . . . . . . . 72.1.4 Quarto approccio . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Dati canonici dell'ellisse . . . . . . . . . . . . . . . . . . . . . . . 82.3 Tangente e Normale . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3.1 Notazione algebrica . . . . . . . . . . . . . . . . . . . . . 112.3.2 Notazione vettoriale . . . . . . . . . . . . . . . . . . . . . 13

2.4 Interessanti proprietà geometriche . . . . . . . . . . . . . . . . . 142.4.1 Relazioni fra gli angoli . . . . . . . . . . . . . . . . . . . . 152.4.2 Relazione fra le ordinate . . . . . . . . . . . . . . . . . . . 162.4.3 La normale principale . . . . . . . . . . . . . . . . . . . . 162.4.4 La normale bisettrice dei raggi vettori in Q dai fuochi . . 16

3 Trasformazione di Coordinate 20

3.1 Da Coordinate Geogra�che a Coordinate Geocentriche . . . . . . 213.2 Da Coordinate Geocentriche a Coordinate Geogra�che . . . . . . 22

3.2.1 Soluzione approssimata proposta . . . . . . . . . . . . . . 223.2.1.1 Test numerico della funzione interpolante . . . . 25

3.2.2 Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Referimenti bibliogra�ci 30

1

Aggiornamenti

� 2015-02-03 estensione della mia congettura (�3.2): DA COMPLETARE

� 2014-12-13 aggiunto capitolo sulla normale bisettrice dei raggi vettori inQ dai fuochi

� 2014-12-10 impianto iniziale

2

Capitolo 1

Premessa

In un personale percorso di avvicinamento e comprensione della manipolazionedelle coordinate fra le diverse notazioni, mi sono imbattuto nel problema del pas-saggio fra coordinate geogra�che (geodetiche, polari ellissoidiche) e geocentriche(cartesiane rettangolari).

Come noto in letteratura il passaggio diretto dalle prime alle seconde èsemplice e immediato, basato su semplici funzioni trigonometriche come saràripreso in seguito; il passaggio inverso, da coordinate rettangolari a coordinategeogra�che, risulta considerevolmente più complesso, basato sostanzialmente sudue metodi1:

� metodo iterativo

� metodo diretto.

Un esempio di quest'ultimo, forse il più noto, è quello elaborato da K. M.Borkowsky nei due testi del 1987 e 1989 di cui riporto sotto i riferimenti; l'el-eganza sintetica del procedimento è parzialmente o�uscata dalla necessità dirisolvere un'equazione di quarto grado, pertanto è elevato l'interesse per unasoluzione, anche solo approssimata, che abbini alla comodità della forma chiusala semplicità di calcolo.

Nel mio studio dell'ellisse credo di essermi accorto di proprietà non riscon-trate in letteratura che volevo segnalare al lettore; alcune di queste propri-età (�2.1.3, �2.1.4, �2.4) poi mi hanno suggerito una congettura (�3.2) su unasoluzione approssimata in forma chiusa che presento appena sotto, dopo averdescritto le proprietà dell'ellisse di cui dicevo e richiamato alcune altre proprietànote, in modo da dare al lettore un quadro completo della situazione.

Il mio comprensibile entusiamo per la �scoperta� è molto scemato quandoho rivisto la soluzione proposta da B. Bowring: nonostante Gloeckler, op. cit.pag. 31, la annoveri fra le soluzioni iterative, già al primo passo questa soluzione

1Si rinvia alla copiosa documentazione contenuta nei riferimenti bibliogra�ci citati perun'ampia rassegna sull'argomento.

3

CAPITOLO 1. PREMESSA 4

o�re un'ottima approssimazione; pazienza, mi sono comunque divertito un saccosentendomi, per qualche giorno, nell'M32.

Tornando con i piedi per terra, ricordo che, al contrario, non sono un matem-atico professionista, pertanto con�do nella bontà del lettore che spero accoglieràle mie elucubrazioni come un gioco dell'intelletto e mi segnali con pazienzafraterna eventuali errori e sviste.

Con molto azzardo a�anco una versione in lingua inglese che ovviamenterichiede al lettore una pazienza ancora maggiore.

Tutto il contenuto del presente documento è rilasciato sotto licenza CC-BY

.

2MMM: Mondo delle Menti Matematiche, sinonimo per me del monte Olimpo, sacro aiGreci

Capitolo 2

La geometria dell'ellisse

2.1 Una imprimitura

L'ellisse può essere vista come un cerchio scalato in modo anisotropico, cioècon due fattori di scala Sx e Sy diversi nelle due direzioni ortogonali; nel casoparticolare che studieremo avremo Sy < Sx, situazione che darà luogo alla tipicaforma con l'asse a, lungo la X, maggiore dell'asse b, lungo la Y , senza per questoperdere generalità in quanto i fattori di scala possono essere semplicementeinvertiti o, ancor più semplicemente, l'ellisse può essere ruotata di 90°.

La completa uniformità del genitore (cerchio) con la sua totale simmetriaradiale non lascia trasparire le molteplici sfaccettature di cui dispone l'ellisse:è proprio su questo che voglio so�ermarmi in questa prima parte dove cercheròdi evidenziare i suoi molteplici aspetti passando da quelli più consueti ai due�nali (�2.1.3, �2.1.4) che mi sembrano sfuggiti, come dicevo in premessa, allaletteratura sull'argomento.

2.1.1 Primo approccio (algoritmo del giardiniere)

5

CAPITOLO 2. LA GEOMETRIA DELL'ELLISSE 6

Immaginiamo di avere i due punti F1 = (−c, 0) e F2 = (c, 0)1 dispostisull'asse X e simmetrici rispetto l'origine e di disegnare i punti Q trattenuti dauna cordicella di lunghezza �ssa ancorata ai punti F1 e F2, cioè F1Q + F2Q =d1 + d2 = d.

Quando Q è sull'asse orizzontale abbiamo:

� Q = (a, 0)

� F1Q = d1 = a+ c

� F2Q = d2 = a− c

� F1Q+ F2Q = d1 + d2 = a+ c+ a− c = 2a = d.

Quando Q è sull'asse verticale abbiamo:

� Q = (0, b)

� F1Q = F2Q = d1 = d2 =√c2 + b2

� F1Q+ F2Q = 2F2Q = d1 + d2 = 2√c2 + b2 = 2a, da cui

� c2 + b2 = a2 e c = ±√a2 − b2.

Nel caso generale abbiamo:F1Q+ F2Q = d1 + d2 =

=

√(x+ c)

2+ y2 +

√(x− c)2 + y2 =

= 2a;con alcuni arrangiamenti otteniamo√(x+ c)

2+ y2 = 2a−

√(x− c)2 + y2

ed elevando al quadrato una prima volta otteniamo

(x+ c)2+ y2 =

(2a−

√(x− c)2 + y2

)2

=

= 4a2 − 4a

√(x− c)2 + y2 + (x− c)2 + y2

da cui

x2 + 2cx+ c2 + y2 = 4a2 − 4a

√(x− c)2 + y2 + x2 − 2cx+ c2 + y2

e

cx = a2 − a√(x− c)2 + y2

risistemato in

a

√(x− c)2 + y2 = a2 − cx;

elevando al quadrato una seconda volta otteniamo

a2x2 − 2a2cx+ a2c2 + a2y2 =(a2 − cx

)2= a4 − 2a2cx+ c2x2

da cui(a2 − c2

)x2 + a2y2 =

(a2 − c2

)a2;

sostituendo b2 = a2 − c2, abbiamob2x2 + a2y2 = b2a2

e �nalmente la forma canonicax2

a2 + y2

b2 = 1.

1I due punti F1e F2 sono detti fuochi.


2.1.2 Secondo approccio

Dalla forma canonica (cfr. Grewal et al, op. cit., Appendix C.3.5)x2

a2 + y2

b2 = 1 = cos2 t+ sin2 t = a2 cos2 ta2 + b2 sin2 t

b2 = (a cos t)2

a2 + (b sin t)2

b2

otteniamox = a cos t, y = b sin t.

2.1.3 Terzo approccio

Immaginiamo una interpolazione trigonometrica del punto Q fra A = (a, 0)a B = (0, b) data da

Q = (x, y) = (a cos t, b sin t)da cui(xa ,

yb

)= (cos t, sin t);

per una nota identità trigonometrica abbiamo

cos2 t+ sin2 t =(xa

)2+(yb

)2= x2

a2 + y2

b2 = 1,la forma canonica implicita dell'equazione dell'ellisse (essendo la forma es-

plicita y2 = b2(1− x2

a2

)= b2

a2

(a2 − x2

)→ y = ±r

√a2 − x2).

Possiamo porrey = b sin t = a

ab sin t = ar sin t


eQ = (x, y) = a (cos t, r sin t)dove r = b

a .

2.1.4 Quarto approccio

Dato il triangolo A1A2B, per note proprietà dei triangoli simili abbiamoQ1H

P1H= b

a = r.

Vogliamo applicare una logica analoga per cui dato un cerchio di raggio a,cerchiamo la curva i cui punti Q, identi�cati dalla verticale per i punti P del

cerchio, abbiano ordinata QH nel rapporto QH

PH= r rispetto all'ordinata PH di

P ; detti t l'angolo ̂A1OP e u l'angolo ̂A1OQ, abbiamo

� xQ = xP = a cos t

� yQ = ryP = ra sin t = baa sin t = b sin t

che rappresentano le coordinate dei punti dell'ellisse di semiassi a e b.

2.2 Dati canonici dell'ellisse

Nel nostro studio ci riferiamo alla �gura sottostante


in cui abbiamo

� a semiasse maggiore

� b semiasse minore

� r = ba rapporto dei semiassi

� c = ae ascissa del fuoco, da cui

� e = ca eccentricità;

� H la proiezione di Q sull'asse X

� P l'intersezione della verticale da Q con il cerchio esterno (Re = A)

� angolo parametrico t = P̂OH

� angolo u = Q̂OH;

considerando l'uguaglianza c2 + b2 = a2 abbiamo

� c2 = a2e2 = a2 − b2

→ e2 = a2−b2a2


� f = a−ba schiacciamento, da cui

� f = a−ba = 1− b

a = 1− r ⇒ r = ba = 1− f rapporto dei semiassi

� b2 = a2 − a2e2 =(1− e2

)a2 → b2

a2 = r2 = 1− e2 e �nalmente

� r = ba = 1− f =

√1− e2 da cui

� (1− f)2 = 1− 2f + f2 = 1− e2 → 2f − f2 = (2− f) f = e2

� d = OQ =√x2 + y2 =

=√x2 + b2

a2 (a2 − x2) =

=√

a2−b2a2 x2 + b2 =

=√e2x2 + b2

o d =√

a2−b2a2 a2 cos2 t+ b2 =

=√(a2 − b2) cos2 t+ b2 =

=√a2 cos2 t+ b2 (1− cos2 t) =

=√a2 cos2 t+ b2 (1− cos2 t) =

=√a2 cos2 t+ b2 sin2 t risultato cui potevamo arrivare direttamente da

d = OQ =√x2 + y2 e x = a cos t, y = b sin t;

dopo di che possiamo riscrivere

� OH = x = d cosu =

=√a2 cos2 t+ b2 sin2 t cosu =

= a cos t

da cui cosu = a cos t√a2 cos2 t+b2 sin2 t

� QH = y = d sinu =

=√a2 cos2 t+ b2 sin2 t sinu =

= b sin t

da cui sinu = b sin t√a2 cos2 t+b2 sin2 t

.

�yx = d sinu

d cosu =

= tanu =

= b sin ta cos t =

= ba tan t

= r tan t =

= (1− f) tan t


� tan t = tanur =

= 11−f tanu.

2.3 Tangente e Normale

Calcoliamo la tangente e la normale in Q = (x0, y0).

2.3.1 Notazione algebrica

� equazione dell'ellisse

x2

a2 + y2

b2 = 1

→ y2

b2 = 1− x2

a2 = a2−x2

a2

→ y =√

b2

a2 (a2 − x2) =

= ba

√a2 − x2 =

= r√a2 − x2 =

= r(a2 − x2

) 12


� derivata prima

y′

x =[r(a2 − x2

) 12

]′x

= 12r(a2 − x2

)− 12 (−2x) =

= − 2rx2√a2−x2

=

= − r2xy =

= − r2

tanu =

= − rtan t

� coe�ciente angolare della normale N

mN = tan v =

= − 1mT

=

= a2yb2x =

= tanur2 =

= tan tr

da cui

� tanu = r2 tan v2 and tan t = r tan v3

� sin t = tan t√1+tan2 t

= r tan v√1+r2 tan2 v

=ba

sin vcos v√

1+ b2

a2sin2 vcos2 v

= b sin v√a2 cos2 v+b2 sin2 v

4

� cos t = 1√1+tan2 t

= 1√1+r2 tan2 v

= 1√1+ b2

a2sin2 vcos2 v

= a cos v√a2 cos2 v+b2 sin2 v

per i punti sul contorno dell'allisse:

� x = a cos t = a2 cos v√a2 cos2 v+b2 sin2 v

� z = b sin t = b2 sin v√a2 cos2 v+b2 sin2 v

e nel caso generale:

� x = a cos t+ h cos v =

(a2√

a2 cos2 v+b2 sin2 v+ h

)cos v 5

2Risultato coincidente la formula fornita da wikipedia ψ (ϕ) = arctan(1− e2

)tanϕ con le

sostituzioni ψ = u, ϕ = v, 1− e2 = r2 [http://en.wikipedia.org/wiki/Latitude].3Risultato coincidente la formula fornita da wikipediaβ(φ) = arctan

[√1− e2 tanφ

]=

arctan [r tanφ] con le sostituzioni β = t, ϕ = v,√1− e2 = r

[http://en.wikipedia.org/wiki/Latitude].4Questo risultato ed il seguente possono essere derivati dalle relazioni tanu = r tan t e

tan t = r tan v trovate in questo paragrafo e dai risultati ottenuti nella sezione 2.2.5Borkowsky, Accurate Algorithms to Transform Geocentric to Geodetic Coordinates, 1989,

fornisce dr = a cos t+ h cosϕ.


� z = b sin t+ h sin v =

(b2√

a2 cos2 v+b2 sin2 v+ h

)sin v6.

Ponendo N2 = a2 cos2 v + b2 sin2 v == a2

(1− sin2 v

)+ b2 sin2 v =

= a2 −(a2 − b2

)sin2 v =

= a2(1− a2−b2

a2 sin2 v)

= a2(1− e2 sin2 v

)otteniamo le espressioni:

� x =(a2

N + h)cos v

� z =(b2

N + h)sin v.

2.3.2 Notazione vettoriale

� tangente T = (x′t, y′t) = (−a sin tdt, b cos tdt)

� normale N = (b cos tdt, a sin tdt)

� coe�cient angolare di N : mn = a sin tdtb cos tdt =

tan tr =

tanur

r = tanur2 .

Quindi abbiamo:

� equazione della normale y − y0 = mN (x− x0) = tan v (x− x0) → y =tan v (x− x0) + y0

� punto S = (xS , yS):

� xS = 0

� yS = tan v (xS − x0) + y0 = y0 − x0 tan v

� SQ2= (x0 − xS)2 + (y0 − yS)2 =

= x20 + [y0 − y0 + x0 tan v]2=

= x20(1 + tan2 v

)=

= x20

(1 + sin2 v

cos2 v

)=

= x20cos2 v+sin2 v

cos2 v =

=x20

cos2 v =

= a2 cos2 tcos2 v =

= a2

cos2 v 1cos2 t

=

6Borkowsky, Accurate Algorithms to transform Geodetic coordinates, 1989, fornisce z =b sin t+ h sinϕ.


= a2

cos2 v(1+tan2 t) =

= a2

cos2 v(1+r2 tan2 v) =

= a2

cos2 v+r2 cos2 v tan2 v =

= a2

cos2 v+r2 sin2 v=

= a2

1−sin2 v+(1−e2) sin2 v=

= a2

1−sin2 v+sin2 v−e2 sin2 v=

= a2

1−e2 sin2 v=

= a4

N27

� SQ = a√1−e2 sin2 v

=

= a2

N8

� point U = (xU , yU ):

� yU = 0

� xU = x0 +yu−y0tan v = x0 − y0

tan v

� UQ2= (x0 − xU )2+(y0 − yU )2 =

(x0 − x0 + y0

tan v

)2+y20 = y20

(1 + 1

tan2 v

)=

y20

(1 + cos2 v

sin2 v

)= y20

(sin2 v+cos2 v

sin2 v

)=

y20sin2 v

= = b2 sin2 tsin2 v

= b2

sin2 v 1sin2 t

=

b2

sin2 v(1+ 1tan2 t

)= b2

sin2 v(1+ 1r2 tan2 v

)= b2

sin2 v+ sin2 vr2 tan2 v

= b2

sin2 v+ cos2 v1−e2

= b2

sin2 v(1−e2)+cos2 v

1−e2

=

=(1−e2)(ar)2

sin2 v−e2 sin2 v+1−sin2 v=

(1−e2)r2a2

1−e2 sin2 v=

(1−e2)2a2

1−e2 sin2 v=(1− e2

)2SQ

2

� UQ =(1− e2

)SQ =

(1−e2)a√1−e2 sin2 v

.

2.4 Interessanti proprietà geometriche

Quì presento, ai paragra� 1, 2 e 3, le osservazioni, forse originali di cui dicevoall'inizio; al paragrafo 4 presento invece una proprietà nota e che vale la penaricordare.

7Borkowsky, Accurate Algorithms to Transform Geocentric to Geodetic Coordinates, 1989,

calcola SQ2= a2

cos2 v+r2 sin2 v= a4

a2 cos2 v+a2r2 sin2 v= a4

(a cos v)2+(b sin v)2.

8Borkowsky, Accurate Algorithms to Transform Geocentric to Geodetic Coordinates, 1989,

calcola SQ =√

a4

(a cos v)2+(b sin v)2= a2√

(a cos v)2+(b sin v)2.


2.4.1 Relazioni fra gli angoli

Da risultati precedenti sappiamo che

1. tanu = r tan t

2. tan t = r tan v;

poichè gli angoli u e t, uniti dalla prima relazione, rappresentano gli angoli alcentro O individuati rispettivamente dalle linee OQ e OP , dove P è l'intercettadella verticale da Q al cerchio esterno, possiamo interpretare gli angoli t e v, chegodono della seconda analoga relazione, come gli angoli al centro O individuatidalle linee OP e OP1 dove P1 è l'intersezione della verticale da Q1 al cerchioesterno e Q1 l'intersezione di OP con l'ellisse.

Pertanto risulta facile trovare per via geometrica la normale (latitudine v)in un punto Q dell'ellisse:

� si trova il punto P innalzando la verticale da Q �no al cerchio di raggio a


� il raggio vettore OP individua la latitudine parametrica t del punto Q

� si trova il punto Q1 intersezione di OP con l'ellisse

� si trova il punto P1 innnalzando la verticale da Q1 �no al cerchio

� il raggio vettore OP1 individua l'angolo v che è latitudine parametrica diQ1 e latitudine geodetica di Q.

2.4.2 Relazione fra le ordinate

Dal precedente paragrafo sappiamo che QH

OH= tanu = r tan t = r PH

OHda cui

QH

PH= r9.

2.4.3 La normale principale

Precedentemente abbiamo usato il parametroN =√a2 cos2 v + b2 sin2 v; poichè

(a cos v, b sin v) sono le coordinate del punto Q1 abbiamo che N è la lunghezzadel vettore OQ1.

Pertanto la normale principale SQ = a2

N trovata al paragrafo �2.3.2 puòessere vista come rapporto fra il quadrato del raggio del cerchio circoscritto e lalunghezza del vettore OQ1.

2.4.4 La normale bisettrice dei raggi vettori in Q dai fuochi

Avevo dimenticato un'altra proprietà, questa volta molto nota e anch'essa degnadi essere ricordata, dell'ellisse; come si vede dalla �gura

9Questa proprietà è già stata usata nel quarto approccio (�2.1.4).


quando il punto Q1 è sulla verticale, si ha completa simmetria: F1O = F2Oper de�nizione, Q1O comune, quindi F1Q1 = F2Q1; la normale coincide conl'asse verticale e, essendo gli angoli ̂F1Q1O = ̂F2Q1O per note proprietà deitriangoli, la normale è bisettrice dell'angolo ̂F1Q1F2 formato dai raggi vettoriF1Q1 e F2Q1: vogliamo analizzare se questa proprietà vale anche per un puntoQ qualsiasi.

Mi sembra ragionevole impiegare le proprietà di spazio euclideo del piano,sfruttando in particolare:

� il prodotto interno: (u,v) =∑i uivi

10

� la norma euclidea: |u| =√(u,u) =

√∑i uiui =

√∑i u

2i

da cui deriva la nota formula per il calcolo dell'angolo compreso fra due vettori

cosα = (u,v)|u||v| .

Data l'ellisse in forma vettoriale abbiamo11

� punto generico Q = (x, y) = (a cosα, b sinα) da cui

� cosα = xa , sinα = y

b

� vettore tangente in Q: t = (−a sinα, b cosα)

� vettore normale in Q: n = (b cosα, a sinα) =

=(bax,

ab y)

10I simboli in grassetto rappresentano i vettori, quelli normali rappresentano i valori scalari.11In questo paragrafo usiamo il simbolo α per indicare la latitudine parametrica.


� vettore F1Q = (x+ c, y)

� modulo |F1Q| =√(x+ c)

2+ y2 =

=√x2 + 2cx+ c2 + y2 =

=√x2 + 2cx+ c2 + b2

a2 (a2 − x2) =

=√x2 + 2cx+ a2 − b2 + b2 − b2

a2x2 =

=√

a2−b2a2 x2 + 2cx+ a2 =

=√

c2

a2x2 + 2cx+ a2 =

=

√(a+ c

ax)2

=

= a+ cax

� vettore F2Q = (x− c, y)

� modulo |F2Q| =√(x− c)2 + y2 =

=√x2 − 2cx+ c2 + y2 =

=√x2 − 2cx+ c2 + b2

a2 (a2 − x2) =

=√x2 − 2cx+ a2 − b2 + b2 − b2

a2x2 =

=√

a2−b2a2 x2 − 2cx+ a2 =

=√

c2

a2x2 − 2cx+ a2 =

=

√(a− c

ax)2

=

= a− cax

� cos γ =(F1Q,n)|F1Q||n| =

=((x+c,y),( bax,

ab y))

( cax+a)|n|=

=(x+c) bax+

ab y

2

( cax+a)|n|=

=bax

2+ bca x+

aba2−x2

a2b2

( cax+a)|n|=

=bax

2+ bca x+ab−

bax

2

( cax+a)|n|=

=(a+ c

ax)b(a+ c

ax)|n|=

= b|n|


� cos δ =(F2Q,n)|F2Q||n| =

=((x−c,y),( bax,

ab y))

( cax−a)|n|=

=(x−c) bax+

ab y

2

( cax−a)|n|=

=bax

2− bca x+aba2−x2

a2b2

( cax−a)|n|=

=bax

2− bca x+ab−bax

2

( cax−a)|n|=

=(a− cax)b(a− cax)|n|

=

= b|n|

quindi cos γ = cos δ → γ = δ, cioè gli angoli formati dai raggi vettori con la nor-male sono uguali o, detto altrimenti, la normale è bisettrice dell'angolo formatodai raggi vettori in Q, quindi la proprietà che avevamo visto nel caso particolarein cui il punto Q era situato sull'asse verticale rappresenta una proprietà validaper ogni punto dell'ellisse. Questa proprietà è nota da tempo e sfruttata daiprogettisti a �ni acustici.

Capitolo 3

Trasformazione di Coordinate

In ambito geogra�co il sistema di riferimento è forzatamente tridimensionale; idue sistemi che stiamo considerando sono:

� sistema geodetico (polare) basato sui vettori mutuamente indipendenti

� longitudine λ

� latitudine ϕ

� altitudine (elissoidica) h

� sistema cartesiano rettangolare basato sui vettori mutuamente indipen-denti

� asse X

� asse Y

� asse Z.

I due sistemi condividono l'origine O e l'asse Z; gli assi X e Y sono disposti nelpiano equatoriale e l'asse X coincide con la traccia del piano verticale da cui sicalcolano le longitudini.

Pertanto avremo tipicamente la relazione f : PG (λ, ϕ, h) → PC (x, y, z);tuttavia la simmetria radiale dell'ellissoide permette una comoda sempli�cazioneche adotteremo senza indugio: adotteremo come spazio il piano contenente l'asseverticale e passante per il punto P dove avremo la relazione sempli�cata g :PG (ϕ, h)→ PC (r, z) dopo la sostituzione r =

√x2 + y2.

Nel percorso inverso avremo g−1 : PC (r, z) → PG (ϕ, h), l'ulteriore oper-azione λ = arctan y

x consente di completare la relazione inversa f−1 : PC (x, y, z)→PG (λ, ϕ, h).

20

CAPITOLO 3. TRASFORMAZIONE DI COORDINATE 21

3.1 Da Coordinate Geogra�che a Coordinate Geo-

centriche

Con riferimento alla �gura, dato il punto P nelle coordinate geogra�che(λ, ϕ, h) abbiamo:

� v = SQ = a√1−e2 sin2 ϕ

� dr = (v + h) cosϕ

� x = dr cosλ

� y = dr sinλ

� q = UQ =(1− e2

)SQ =

(1− e2

)v = r2v

� z = (q + h) sinϕ =[(1− e2

)v + h

]sinϕ =

[r2v + h

]sinϕ.


3.2 Da Coordinate Geocentriche a Coordinate Ge-

ogra�che

Come anticipato la relazione inversa risulta notevolmente più complessa e siriduce in sostanza alla scrittura di un'equazione di quarto grado1; per quantola soluzione di questa equazione, e quindi della trasformazione inversa, sia allaportata di qualsiasi moderno strumento di calcolo, mi sembra di aver osservatoalcune caratteristiche che consentono di congetturare una soluzione approssima-ta in forma chiusa; nel seguito descrivo queste osservazioni e poi la soluzioneproposta.

3.2.1 Soluzione approssimata proposta

Con riferimento alla �gura abbiamo gli elementi:

1. A = (xA, yA) punto in esame

2. tanw = yAxA

coe�ciente angolare della retta OA

3. tan t1 = 1r tanw tangente della latitudine parametrica 1

4. Q1 = (x1, y1) = (a cos t1, b sin t1) = a (cos t1, r sin t1) punto sulla super�ciedell'ellisse di latitudine parametrica t1

5. d = Q1A distanza Q1A

6. tan t2 = tanw tangente della latitudine parametrica 2

1Mi sono intestardito più volte, nonostante le mie insu�cienti conoscenze matematiche,a cercare qualche espressione più semplice, ma tutte le strade �nora tentate (punto Q comepiede della perpendicolare da P all'ellisse, punto Q come punto più vicino a P , ecc.) si sonorivelate illusorie: lo sbocco è sempre e solo un'equazione di quarto grado.


7. Q2 = (x2, y2) = (a cos t2, b sin t2) = a (cos t2, r sin t2) punto sulla super�ciedell'ellisse di latitudine parametrica t2

8. tan v2 = 1r tan t2 = 1

r tanw = tan t1 tangente della normale in Q2

9. A0 = (x0, y0) intersezione della retta OA con la verticale per Q2

w = (cosw, sinw)

v2 = (cos v2, sin v2)

A0 = Q1 + αw = Q2 + βv2 da cui[x0y0

]=

[x1y1

]+ α

[coswsinw

]=

[x2y2

]+ β

[cos v2sin v2

]α

[coswsinw

]−β[

cos v2sin v2

]=

[cosw − cos v2sinw − sin v2

] [αβ

]=

[x2y2

]−[x1y1

]=[

x2 − x1y2 − y1

]e quindi[αβ

]=


]−1 [x2 − x1y2 − y1

];

sostituendo a = cosw, b = cos v2, c = sinw, d = sin v2 abbiamo[a −b 1 0c −d 0 1

]∼[

1 − ba

1a 0

1 −dc 0 1c

]∼[

1 − ba

1a 0

0 ba −

dc − 1

a1c

]∼[

1 − ba

1a 0

0 bc−adac − 1

a1c

]∼[

1 − ba

1a 0

0 1 − cbc−ad

abc−ad

]∼[

1 − ba

1a 0

0 ba − bc

a(bc−ad)b

bc−ad

]∼

[1 0 1

a −bc

a(bc−ad)b

bc−ad0 b

a − bca(bc−ad)

bbc−ad

]

∼

[1 0 bc−ad−bc

a(bc−ad)b

bc−ad0 1 −c

bc−ada

bc−ad

]∼[

1 0 −dbc−ad

bbc−ad

0 1 −cbc−ad

abc−ad

]∼[a −bc −d

]−1=[ −d

bc−adb

bc−ad−c

bc−ada

bc−ad

]= 1

bc−ad

[−d b−c a

]= 1

D

[− sin v2 cos v2− sinw cosw

]dove ab-

biamo posto D = bc− ad = sinw cos v2 − cosw sin v2; check[cosw − cos v2sinw − sin v2

] [cosw − cos v2sinw − sin v2

]−1=


]1D


]=

= 1D

[− cosw sin v2 + sinw cos v2 cosw cos v2 − cosw cos v2− sinw sin v2 + sinw sin v2 sinw cos v2 − cosw sin v2

]= 1

D

[D 00 D

]=[

1 00 1

]c.v.d.;

pertanto[αβ

]=


]−1 [x2 − x1y2 − y1

]= 1

D


] [x2 − x1y2 − y1

],

cioè{α = −(x2−x1) sin v2+(y2−y1) cos v2

D

β = −(x2−x1) sinw+(y2−y1) coswD


e quindi

A0 =

[x0y0

]=

[x1y1

]+ α

[coswsinw

]=

=

[x1 +

−(x2−x1) sin v2+(y2−y1) cos v2D cosw

y1 +−(x2−x1) sin v2+(y2−y1) cos v2

D sinw

]=

=

[x1(sinw cos v2−cosw sin v2)−(x2−x1) cosw sin v2+(y2−y1) cosw cos v2

Dy1(sinw cos v2−cosw sin v2)−(x2−x1) sinw sin v2+(y2−y1) sinw cos v2

D

]=

=

[x1 sinw cos v2−x1 cosw sin v2−x2 cosw sin v2+x1 cosw sin v2+y2 cosw cos v2−y1 cosw cos v2

Dy1 sinw cos v2−y1 cosw sin v2−x2 sinw sin v2+x1 sinw sin v2+y2 sinw cos v2−y1 sinw cos v2

D

]=

=

[x1 sinw cos v2−x2 cosw sin v2+y2 cosw cos v2−y1 cosw cos v2

D−y1 cosw sin v2−x2 sinw sin v2+x1 sinw sin v2+y2 sinw cos v2

D

]=

=

[x2 sinw cos v2−x2 sinw cos v2+x1 sinw cos v2−x2 cosw sin v2+y2 cosw cos v2−y1 cosw cos v2

Dy2 cosw sin v2−y2 cosw sin v2−y1 cosw sin v2−x2 sinw sin v2+x1 sinw sin v2+y2 sinw cos v2

D

]=

=

[x2(sinw cos v2−cosw sin v2)+x1 sinw cos v2−x2 sinw cos v2+y2 cosw cos v2−y1 cosw cos v2

Dy2(sinw cos v2−cosw sin v2)+y2 cosw sin v2−y1 cosw sin v2−x2 sinw sin v2+x1 sinw sin v2

D

]=

=

[x2(sinw cos v2−cosw sin v2)

D + −(x2−x1) sinw+(y2−y1) coswD cos v2

y2(sinw cos v2−cosw sin v2)D + −(x2−x1) sinw+(y2−y1) cosw

D sin v2

]=

=

[x2 +

−(x2−x1) sinw+(y2−y1) coswD cos v2

y2 +−(x2−x1) sinw+(y2−y1) cosw

D sin v2

]=

=

[x2y2

]+ −(x2−x1) sinw+(y2−y1) cosw

D

[cos v2sin v2

]=

=

[x2y2

]+ β

[cos v2sin v2

]c.v.d.

10. d0 = Q1A0 distanza Q1A0

Spostando il punto A lungo la retta OQ1 abbiamo

1. A ≡ Q1 → dd0

= 0 → u = w → tan t = 1r tanu = 1

r tanw → tan v =1r tan t =

1r2 tanw

2. A ≡ A0 → dd0

= 1 → t = w → tan t = tanw → tan v = 1r tan t =

1r tanw

3. A→∞→ dd0→∞→ v = w → tan v = tanw → tan t = r tan v = r tanw

e possiamo esprimere le due disequazioni

1. tan t = st tanw con r ≤ st ≤ 1r e

2. tan v = sv tanw con 1 ≤ sv ≤ 1r2 .


Usando la prima cerchiamo la funzione di interpolazione st = f(dd0

)tale da

fornire

1. f (0) = 1r

2. f (1) = 1

3. f (∞) = r;

proviamo la funzione s = α arctanβq + γ, dopo aver sostituito dd0

= q:

1. st=0 = α arctan 0 + γ = γ = 1r

2. st=1 = α arctanβ + γ = 1 da cui β = tan 1−γα = tan

1− 1r

α = tan r−1αr

3. st→∞ = α arctan∞+γ = απ2 +γ = r da cui α = 2π (r − γ) = 2

π

(r − 1

r

)=

2πr2−1r

4. e quindi β = tan r−1αr = tan r−1

2πr2−1r r

= tan r−1r2−1

π2 = tan 1

1+rπ2

pertanto adottiamo la funzione interpolante s = 2πr2−1r arctan

(tan 1

1+rπ2

)q+ 1

r .

Nota: il rapporto q = dd0

comporta il calcolo degli incrementi dx e dy,l'elevamento al quadrato e l'estrazione della radice della somma per il calcolodelle due distanze d e d0; è più e�ciente applicare il teorema di Talete ai triangolirettangoli simili di ipotenusa Q1A e Q1A0 con cui q = d

d0= dx

dx0dove abbiamo

posto dx = xA − x1 e dx0 = x0 − x12.

3.2.1.1 Test numerico della funzione interpolante

1. r = 12 = 0.5

(a) α = 2πr2−1r = 2

π

14−1

12

= − 3π

(b) β = tan 11+r

π2 = tan 1

32

π2 = tan π

3 = 1.732

(c) γ = 1r = 2

(d) s = α arctanβt+ γ = − 3π arctan 1.732t+ 2

(e) st=0 = − 3π arctan 0 + 2 = 2

(f) st=1 = − 3π arctan 1.732 + 2 = − 3

ππ3 + 2 = 1

(g) st→∞ = − 3ππ2 + 2 = 1

2

2. r = 0.8

(a) α = 2πr2−1r = 2

π0.64−1

0.8 = − 910π

2Questo metodo ha anche il vantaggio di individuare con un valore negativo il caso di puntidi altitudine negativa, al di sotto della super�cie terrestre, che qui non sono trattati.


(b) β = tan 11+r

π2 = tan π

3.6 = 1.19175

(c) γ = 1r = 1.25

(d) s = α arctanβt+ γ = − 910π arctan 1.19175t+ 1.25

(e) st=0 = − 910π arctan 0 + 1.25 = 1.25

(f) st=1 = − 910π arctan 1.19175+1.25 = − 9

10π10π36 +1.25 = − 1

4+1.25 = 1

(g) st→∞ = − 910π

π2 + 1.25 = − 9

20 + 1.25 = 0.8.

Usando la seconda e la funzione di interpolazione sv = f (t) = α arctanβq + γavremmo

1. f (0) = 1r2

2. f (1) = 1r

3. f (∞) = 1;

e quindi

1. st=0 = α arctan 0 + γ = γ = 1r2

2. st=1 = α arctanβ + γ = 1r da cui β = tan

1r−γα = tan

1r−

1r2

α = tan r−1αr2 =

− tan 1−rαr2

3. st→∞ = απ2 + γ = 1 da cui α = 2π (1− γ) = 2

π

(1− 1

r2

)= 2

πr2−1r2 =

− 2π

1−r2r2

4. e quindi β = tan r−1αr2 = tan r−1

2 r2−1

r2r2π = tan 1

1+rπ2

pertanto adottiamo la funzione interpolante s = 2πr2−1r2 arctan

(tan 1

1+rπ2

)q +

1r2 .

Test numerici della funzione interpolante:

1. r = 12 = 0.5

(a) α = 2πr2−1r2 = 2

π

14−1

14

= − 2π

3441 = − 6

π

(b) β = tan 11+r

π2 = tan 1

1+ 12

π2 = tan 2

3π2 = tan π

3 = 1.732

(c) γ = 1r2 = 1

14

= 4

(d) s = 2πr2−1r2 arctan

(tan 1

1+rπ2

)t+ 1

r2 = − 6π arctan 1.732t+ 4

(e) st=0 = − 6π arctan 0 + 4 = 4

(f) st=1 = − 6π arctan 1.732 + 4 = − 6

ππ3 + 4 = −2 + 4 = 2

(g) st→∞ = − 6ππ2 + 4 = −3 + 4 = 1

2. r = 0.8


(a) α = 2πr2−1r2 = 2

π0.64−10.64 = − 9

8π

(b) β = tan 11+r

π2 = tan 1

1+0.8π2 = tan π

3.6 = 1.19175

(c) γ = 1r2 = 1

0.64 = 1.5625

(d) s = 2πr2−1r2 arctan

(tan 1

1+rπ2

)t+ 1

r2 = − 98π arctan 1.19175t+ 1.5625

(e) st=0 = − 98π arctan 0 + 1.5625 = 1.5625

(f) st=1 = − 98π arctan 1.19175 + 1.5625 = − 9

8π10π36 + 1.5625 = − 5

16 +1.5625 = 1.25

(g) st→∞ = − 98π

π2 + 1.5625 = − 9

16 + 1.5625 = 1.

3.2.2 Test

Ho condotto alcuni test mediante un piccolo script python3 confrontando i risul-tati con esempi trovati in letteratura:

1. Esempio da Università di Pisa, Geodesia, pag. 39-40)

(a) dati ellissoide

Ellissoide INTL24 (Hayford)semiasse maggiore 6378388.0semiasse minore 6356911.94613

rapporto 0.996632996633schiacciamento 298.257

eccentricità al quadrato 0.00672267002233

(b) calcolo diretto

coordinate geogra�che coordinate cartesianelon 1d 16m 51.012s X 4616659.048lat 43d 42m 45.418s Y 103221.763h 0.000000 Z 4385144.185

(c) calcolo inverso

coordinate cartesiane Bowring Gloeckler mio metodoX 4616659.048 lon 1d 16' 51.012" 1d 16' 51.012" 1d 16' 51.012"Y 103221.763 lat 43d 42' 45.418" 43d 42' 45.418" 43d 42' 45.418"Z 4385144.185 h 0.000000 0.000000 0.000000

2. Esempio da Università di Pisa, Geodesia, pag. 41

(a) dati ellissoide

3Lo script è disponibile per chiunque voglia prenderne visione.


Ellissoide INTL24 (Hayford)semiasse maggiore 6378388.0semiasse minore 6356911.94613



(b) calcolo diretto

coordinate geogra�che coordinate cartesianelon 30d 0m 0.0s X 3912744.523lat 45d 0m 0.00s Y 2259024.104h h=350.000 Z 4487676.524

(c) calcolo inverso



(a) dati ellissoide

Ellissoide WGS84semiasse maggiore 6378137.0semiasse minore 6356752.29822



(b) calcolo diretto


(c) calcolo inverso



(a) dati ellissoide


Ellissoide WGS84semiasse maggiore 6378137.0semiasse minore 6356752.29822



(b) calcolo diretto


(c) calcolo inverso


Capitolo 4

Referimenti bibliogra�ci

Per interessanti annotazioni e alcuni metodi costruttivi dell'ellisse si vedano:

� Grewal et al, Global Positioning Systems, Inertial Navigation and Inte-gration, Appendix C

� Borkowsky, Transformation of Geocentric to Geodetic Coordinates with-out Approximations, 1987

� Borkowsky, Accurate Algorithms to Transform Geocentric to GeodeticCoordinates, 1989

� Gloeckler et al, Handbook for Transformation of Datums, Projections,Grids and common Coordinate Systems (TEC-SR-7), 1986

� Richard H. Rapp, Geometric Geodesy, 2 voll, 1991

� la pagina web http://it.wikipedia.org/wiki/Ellisse.

30

Geographic vs geocentric coordinates-2015-02-03-ita

Science

Transcript of Geographic vs geocentric coordinates-2015-02-03-ita