matrici di trasformazione
Locked
 
|
ciao,
se si volesse indagare (nel senso di confrontare) sulle matrici di trasformazione esiste qualche possibilità in qgis, mediante pyqgis, di visualizzarle? eventualmente da linea di comando con proj? grazie, giuliano PS: scrivo qui sotto, per i più curiosi, il motivo della domanda; Salvo C. alias elyparker aveva sollevato tempo fa un problema di mancato allineamento fra mappe di origine (e SR) diversa; ora siamo riusciti, tramite un plugin sperimentale di cui parlo in altro post, a costruire una matrice di trasformazione che sembra dare un risultato migliore: il confronto fra le due matrici poteva dare indicazioni utili sulle cause del mancato allineamento e sulle eventuali contromisure da adottare; _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
Re: matrici di trasformazione
|
Ciao Giuliano,
ti chiedo scusa, probabilmente non ho capito bene la tua domanda, ma intendi i parametri di rototraslazione tra i diversi SR? Quelli li trovi direttamente all'interno delle definizioni delle stringhe proj, in corrispondenza del parametro +towgs84. Ad esempio, per quanto riguarda il sistema Gauss-Boaga fuso Ovest (EPSG 3003), hai i seguenti valori di rototraslazione: "+towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68" Ovvero: Traslazione X (Dx_BF): -104.1 m Traslazione Y (Dy_BF): -49.1 m Traslazione Z (Dz_BF): -9.9 m Rotazione X (Rx_BF): 0.971 secondi Rotazione Y (Ry_BF): -2.917 secondi Rotazione Z (Rz_BF): 0.714 secondi Fattore di scala (M_BF): -11.68 ppm Le coordinate nel sistema WGS84, partendo dal sistema Gauss-Boaga, sono le seguenti (tratto da: https://trac.osgeo.org/proj/wiki/GenParms): x_out = M_BF*( x_in - Rz_BF*y_in + Ry_BF*z_in) + Dx_BF y_out = M_BF*( Rz_BF*x_in + y_in - Rx_BF*z_in) + Dy_BF z_out = M_BF*(-Ry_BF*x_in + Rx_BF*y_in + z_in) + Dz_BF Nel caso volessi combinare due sistemi di riferimento che non siano WGS84, secondo la logica di proj dovresti comunque passare da quel sistema come "base di riferimento". Quindi, volendo passare dal sistema Gauss-Boaga al sistema ED50, occorrerà combinare due trasformazioni di Helmert: 1) da Gauss-Boaga a WGS84 2) da WGS84 a ED50 (invertendo i parametri contenuti nella stringa proj). QGIS effettua al volo questa sequenza di trasformazioni. Se ho capito male la domanda...beh, almeno ho messo in "bella" queste informazioni! ;-) M. |
|
|
Il giorno Thu, 13 Feb 2014 00:08:24 -0800 (PST)
Mattia De Agostino <[hidden email]> ha scritto: > Ciao Giuliano, ciao Mattia, > ti chiedo scusa, probabilmente non ho capito bene la tua domanda, ma > intendi i parametri di rototraslazione tra i diversi SR? > Quelli li trovi direttamente all'interno delle definizioni delle > stringhe proj, in corrispondenza del parametro +towgs84. Ad esempio, > ...... hai capito benissimo e ti ringrazio della chiara esposizione; li avevo visti ma non avendo un background da geografo non avevo ben chiaro il quadro; ora è certamente meglio :-) > Le coordinate nel sistema WGS84, partendo dal sistema Gauss-Boaga, > sono le seguenti (tratto da: > https://trac.osgeo.org/proj/wiki/GenParms): x_out = M_BF*( x_in > - Rz_BF*y_in + Ry_BF*z_in) + Dx_BF y_out = M_BF*( Rz_BF*x_in + > y_in - Rx_BF*z_in) + Dy_BF z_out = M_BF*(-Ry_BF*x_in + Rx_BF*y_in > + z_in) + Dz_BF a parte la considerazione che la scala è quindi isotropa, non mi tornano le rotazioni: Rx_BF, Ry_BF e Rz_BF sono angoli, nella matrice ci dovrebbero andare i coseni (nella diagonale principale) ed i seni; i coseni, trattandosi di angoli molto piccoli, potrebbero essere approssimati ad 1.00, ma seni dell'ordine di 0.917 non sono riferiti ad angoli piccoli (e -2.917 addirittura non esiste); ovviamente possono essere prodotti di più fattori, devo approfondire :-( grazie comunque della spiegazioen molto utile (e scusa la mia pedanteria :-( > Nel caso volessi combinare due sistemi di riferimento che non siano > WGS84, secondo la logica di proj dovresti comunque passare da quel > sistema come "base di riferimento". ..... io devo lavorare su un Cassini/Soldner, però la struttura è la stessa, no? semmai ti rompo ancora; > .... > Se ho capito male la domanda...beh, almeno ho messo in "bella" queste > informazioni! ;-) beh, credo di avere acquisito meriti imperituri per averti obbligato a farlo :-) :-) > M. grazie, ciao, giuliano PS: grazie anche per l'altra splendida spiegazione sul geoide; magari più avanti ti farò una domanda su gps, coordinate geocentriche, topocentriche, ecc. ma lo farò restando nel thread che hai aperto :-) _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
Re: matrici di trasformazione
|
|
La tipologia di trasformazione che viene applicata da PROJ - e quindi da QGIS - è una trasformazione di Helmert a 3 o a 7 parametri (a seconda del numero di valori del parametro +towgs84). Puoi trovare l'indicazione di come funziona la trasformazione di Helmert qui (da pag. 19 in poi): http://labtopo.ing.unipg.it/files_sito/compiti/georef_d.pdf Come vedi dal documento allegato, i termini Rx_BF, Ry_BF e Rz_BF sono gli elementi della matrice di rotazione linearizzata. Per piccoli angoli (e qui parliamo di secondi sessagesimali!) l'approssimazione è accettabile. Ciao Mattia |
|
|
Il giorno Fri, 14 Feb 2014 01:41:52 -0800 (PST)
Mattia De Agostino <[hidden email]> ha scritto: ciao Mattia, > giulianc51 wrote > > a parte la considerazione che la scala è quindi isotropa, non mi > > tornano le rotazioni: Rx_BF, Ry_BF e Rz_BF sono angoli, ... > > La tipologia di trasformazione che viene applicata da ... > http://labtopo.ing.unipg.it/files_sito/compiti/georef_d.pdf > > Come vedi dal documento allegato, i termini Rx_BF, Ry_BF e Rz_BF sono > gli elementi della matrice di rotazione linearizzata. Per piccoli > angoli (e qui parliamo di secondi sessagesimali!) l'approssimazione è > accettabile. ti ringrazio; ovviamente non volevo spaccare il capello, era per rendermi conto dell'origine dei parametri; > Ciao > Mattia ciao, giulian _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
|
|
Il giorno Fri, 14 Feb 2014 12:43:00 +0100
giulianc51 <[hidden email]> ha scritto: > .... in merito al problema delle proiezione su cui mi sto applicando (si fa per dire :-) trovo questo problema su cui chiedo luce a qualche (paziente e gentile) esperto; sto provando la trasformazione di un punto in coordinate proiettate GaussBoaga fuso Est (epsg 3003) in coordinate geografiche Wgs84 (epsg 4326); il punto (a caso) è: E = 1525480.000 N = 5023760.000 H = 0.000 (messe nel file proj_dati-1.txt); eseguendo il comando: ~$ invproj +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000 +y_0=0 +ellps=intl +towgs=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68,+units=m +no_defs proj_dati-1.txt ottengo le coordinate: 9d19'31.263"E 45d21'57.756"N 0.000; eseguo la stessa operazione sul PCN fra Roma 1940 zona 1 e ETRS89-ETRF89 ed ottengo: 9.32502553;45.36671369 equivalenti a: 9d19'30.092” 45d22'00.169”; _ipotizzando io_ che "Roma 1940 zona 1" ~ "GaussBoaga fuso Est" e "ETRS89-ETRF89" ~ "WGS84" (*) trovo i dati discordanti: è dovuto all'uso dei grigliati? è sbagliata l'ipotesi? sto sbagliando qualcos'altro? grazie infinite, ciao, giuliano (*) se non ho frainteso l'ho visto nel sito PCN; _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
|
|
Giuliano,
gauss boaga zona 1 è ovest e non est. :) E Quindi ... Ciao Roberto Inviato da iPhone Il giorno 21/feb/2014, alle ore 12:34, giulianc51 <[hidden email]> ha scritto: > Il giorno Fri, 14 Feb 2014 12:43:00 +0100 > giulianc51 <[hidden email]> ha scritto: > >> .... > > in merito al problema delle proiezione su cui mi sto applicando (si fa > per dire :-) trovo questo problema su cui chiedo luce a qualche > (paziente e gentile) esperto; > > sto provando la trasformazione di un punto in coordinate proiettate > GaussBoaga fuso Est (epsg 3003) in coordinate geografiche Wgs84 (epsg > 4326); > > il punto (a caso) è: E = 1525480.000 N = 5023760.000 H = 0.000 (messe > nel file proj_dati-1.txt); > > eseguendo il comando: > ~$ invproj +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000 > +y_0=0 +ellps=intl > +towgs=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68,+units=m > +no_defs proj_dati-1.txt > ottengo le coordinate: > 9d19'31.263"E 45d21'57.756"N 0.000; > > eseguo la stessa operazione sul PCN fra Roma 1940 zona 1 e > ETRS89-ETRF89 ed ottengo: > 9.32502553;45.36671369 > equivalenti a: > 9d19'30.092” 45d22'00.169”; > > _ipotizzando io_ che "Roma 1940 zona 1" ~ "GaussBoaga fuso Est" e > "ETRS89-ETRF89" ~ "WGS84" (*) trovo i dati discordanti: > > è dovuto all'uso dei grigliati? > > è sbagliata l'ipotesi? > > sto sbagliando qualcos'altro? > > > grazie infinite, ciao, > giuliano > > (*) se non ho frainteso l'ho visto nel sito PCN; > > _______________________________________________ > [hidden email] > http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss > Questa e' una lista di discussione pubblica aperta a tutti. > I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. > 666 iscritti al 22.7.2013 [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
|
|
Il giorno Fri, 21 Feb 2014 17:31:49 +0100
Geodrinx <[hidden email]> ha scritto: > Giuliano, > > gauss boaga zona 1 è ovest e non est. sì, certo, chiedo scusa dell'errore; grazie Roberto; > :) > > E Quindi ... > > Ciao > Roberto ciao, giuliano _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
Re: matrici di trasformazione
|
|
In reply to this post by Giuliano Curti
Ciao Giuliano,
nessun errore (non ho rifatto in realtà la procedura con PROJ, ma sono sicuro che tu l'abbia applicata correttamente!), le discordanze che trovi sono assolutamente ragionevoli...purtroppo! La differenza è legata proprio alla combinazione di tutti i problemi che hai già individuato tu: 1) il PCN utilizza i grigliati IGM, ovvero una maglia di punti abbastanza fitta di cui sono note le coordinate in entrambi i sistemi di riferimento, mentre PROJ utilizza una trasformazione di Helmert a 7 parametri tale da ridurre gli scarti su tutto il fuso Ovest; 2) "ETRS89-ETRF89" ~ "WGS84". Si tratta di un'approssimazione, visto che WGS84 è un datum globale, valido quindi per tutto il globo terrestre, mentre ETRS89 è un datum locale, europeo, frutto dell'eliminazione del moto medio della placca euroasiatica (per ridurre gli spostamenti legati alla deriva dei continenti, che sono di entità non trascurabile purtroppo). Queste due approssimazioni, combinate insieme, generano differenze di quegli ordini di grandezza purtroppo. Mattia |
|
|
Il giorno Sun, 23 Feb 2014 23:58:39 -0800 (PST)
Mattia De Agostino <[hidden email]> ha scritto: > Ciao Giuliano, ciao Mattia, > nessun errore (non ho rifatto in realtà la procedura con PROJ, ma sono > sicuro che tu l'abbia applicata correttamente!), ti ringrazio della fiducia, ma credo di averci messo del "mio" nel casino, infatti ho lasciato il parametro +ellps=intl, quindi la mia chiamata (ma chiedo conferma) non prevedeva trasformazione (modifica di datum) ma solo conversione dal sistema piano al sistema geografico sottostante (e dalle mie prove che rinvio ad una prossima mail) dovrebbero tornare; quindi ricapitolando, chiamata: $ echo '1525480.00 5023760.00 0.00' | invproj +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000 +y_0=0 +ellps=intl risultato: 9d19'31.263"E 45d21'57.756"N 0.00 (penso corretto); mi mettono invece in difficoltà le chiamate: ~$ echo '1525480.00 5023760.00 0.00' | invproj +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000 +y_0=0 +ellps=WGS84 e: $ echo '1525480.00 5023760.00 0.00' | invproj +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000 +y_0=0 +ellps=WGS84 +towgs=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68,+units=m +no_defs che danno entrambe: 9d19'31.335"E 45d22'0.803"N 0.00 da cui sembra che: a) il parametro +towgs=.... non influisca b) il risultato è comunque diverso da quello del PCN (9d19'30.092” 45d22'00.169”); > Mattia in una prossima mail cerco di dare un quadro più completo all'argomento; grazie infinite, ciao, giuliano _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
|
|
In reply to this post by Mattia De Agostino
Il giorno Sun, 23 Feb 2014 23:58:39 -0800 (PST)
Mattia De Agostino <[hidden email]> ha scritto: > Ciao Giuliano, ciao Mattia, eccomi con una quadro più ampio, spero, della situazione; ho provato a farmi qualche nozione più precisa delle "operazioni", come le chiama Epsg, sulle coordinate; le posto nella speraza che qualche esperto abbia voglia di darmi indicazioni per proseguire su una fruttuosa strada di apprendimento dei sistemi di riferimento geografici (pro domo mea :-); spero possano essere utili anche ad altri per risparmiare loro qualche fatica in questa materia tutt'altro che semplice anche se, come diceva Menecmo al grande Alessandro: non vi sono scorciatoie "regie"....... provo a contestualizzare e completare una risposta datami in precedenza da Mattia: > Quelli li trovi direttamente all'interno delle definizioni delle > stringhe proj, in corrispondenza del parametro +towgs84. Ad esempio, > ........ > "+towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68" > Ovvero: > Traslazione X (Dx_BF): -104.1 m > ...... > Le coordinate nel sistema WGS84, partendo dal sistema Gauss-Boaga, > sono le seguenti (tratto da: > https://trac.osgeo.org/proj/wiki/GenParms): x_out = M_BF*( x_in > - Rz_BF*y_in + Ry_BF*z_in) + Dx_BF y_out = M_BF*( Rz_BF*x_in + > y_in - Rx_BF*z_in) + Dy_BF z_out = M_BF*(-Ry_BF*x_in + Rx_BF*y_in > + z_in) + Dz_BF che rischia di essere parziale e fuorviante; NON intendo assolutamente criticare Mattia, ANZI lo ringrazio caldamente perchè lo sbattimento di testa provocato da quella risposta mi ha stimolato a procedere; bando alle ciance adesso ed entriamo in argomento; 1) l'applicazione diretta, da me tentata in prima istanza, delle formule indicate a coordinate epsg 3003 non ha senso; non c'è modo di ricavare le coordinate geografiche cercate (wgs84 epsg 4326); 2) probabilmente molti altri, ma certamente l'Handbook di Gloeckler per la DoD americana è molto chiaro: la trasformazione si applica alle coordinate geocentriche, non a quelle piane; 3) i prametri di rotazione (Rx_BF, Ry_BF e Rz_BF) non si applicano come indicati, ma vanno convertiti in radianti; per cui la catena, per come l'ho finora capita, dovrebbe essere questa (almeno per la trasformazione GB->WBS64, per le altre situazioni occorre fare le aggiunte o sottrazioni del caso): a) coordinate proiettate -> coordinate geografiche (relative all'elissoide di riferimento, ad es. GB/1 -> INTL24); b) coordinate geografiche -> coordinate geocentriche (sempre riferite allo stesso ellissoide); c) coordinate geocentriche (INTL24) -> coordinate geocentriche (WGS84): è solo a questo punto che interviene la trasformazione con i parametri +towgs di cui parla Mattia nella sua risposta; d) coordinate geocentriche -> coordinate geografiche (sull'elissoide selezionato; nell'esempio WGS84); ammesso che questo sia corretto, mi chiedo: 1) se non esista un metodo per passare da coordinate proiettate a coordinate geocentriche, fondendo i passi a + b, (proverò a vedere se è riconducibile a prodotto di matrici); 2) il perchè di una strana assenza: le coordinate topocentriche che non fanno capolino nelle formule finora incontrate :-) ho riportato tutte le operazioni indicate sopra in fogli di calcolo che sono a disposizione di chiunque volesse esaminarli, correggerli, ampliarli, ecc.; per quanto mi riguarda continuerò la ricerca postando gli eventuali risultati, nella speranza di non rompere troppo; grazie, ciao, giuliano _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
Re: matrici di trasformazione
|
|
Ciao Giuliano,
nessun problema! Anzi, mi fa molto piacere parlare di questi argomenti! Confermo che la prima operazione è quella corretta: Se ho capito bene la logica di PROJ (non si finisce mai di imparare), anche il comportamento delle ultime due stringhe è corretto. Nelle ultime due operazioni, stai chiedendo a PROJ di trasformare in coordinate sessagesimali delle coordinate cartografiche che hanno: - una proiezione cartografica trasversa di mercatore (OK) - un'origine per le latitudini di 0° (OK) - un'origine per le longitudini di 9° (OK) - coefficiente di contrazione del cilindro secante di 0.9996 (OK) - falsa origine su Est di 1500000 (OK) - ellissoide di riferimento WGS84 (mentre Roma40-MonteMario utilizza l'ellissoide Internazionale "itnl"). I parametri "towgs" in questo caso non vengono applicati proprio perché sei già nel sistema WGS84, quindi non devi fare alcun cambio di datum. I risultati ti vengono coincidenti con quelli del PCN, se effettui la trasformazione da UTM32N a geografiche per il sistema ETRF89, a patto che tu sottragga 1000000 alla tua Est (perché il sistema UTM ha falsa origine 500000). Di fatto, con quelle operazioni hai fatto creare a PROJ un sistema di riferimento nuovo, del tutto simile al UTM 32N a meno di una falsa origine di 1500000. Direi che la catena di operazioni è assolutamente corretta (o meglio, anche io avrei fatto così!). ;-) In effetti le trasformazioni di Helmert vanno applicate alle coordinate geocentriche, sono stato troppo speditivo l'altra volta (pardon!). Pronto ad essere smentito da altri utenti, ma temo che la risposta sia negativa... il fatto è che le coordinate proiettate dipendono dalla longitudine. Quindi temo che la successione dei passaggi a+b sia obbligata. Riuscire a passare attraverso un prodotto di matrici secondo me è possibile, ma solo come adattamento locale in una zona molto ristretta (qualche km?)... se vuoi una procedura "standard", valida indipendentemente dalla zona, il passaggio tra i diversi sistemi di coordinate è obbligatorio. Ci mancherebbe! Siamo qui per supportarci a vicenda, il dibattito aiuta tutti (spero!). Anzi, spero che la mia risposta possa esserti stata d'aiuto! Ciao! |
|
|
Il giorno Mon, 24 Feb 2014 02:58:20 -0800 (PST)
Mattia De Agostino <[hidden email]> ha scritto: > Ciao Giuliano, > nessun problema! Anzi, mi fa molto piacere parlare di questi > argomenti! ciao Mattia, grazie :-) > ....... > Anzi, spero che la mia risposta possa esserti stata d'aiuto! mi sento come un pitone che ha appena ingerito un facocero..... mi ci vorrà qualche tempo per riacquisire un pò di lucidità :-) a più tardi per una qualsiasi risposta sensata (si fa per dire :-) :-) > Ciao! grazie, ciao, giuliano _______________________________________________ [hidden email] http://lists.gfoss.it/cgi-bin/mailman/listinfo/gfoss Questa e' una lista di discussione pubblica aperta a tutti. I messaggi di questa lista non hanno relazione diretta con le posizioni dell'Associazione GFOSS.it. 666 iscritti al 22.7.2013 |
«
Return to Gfoss -- Geographic Free and Open Source Software - Italian mailing list
|
1 view|%1 views
| Free forum by Nabble | Edit this page |