Discussione:
[Gfoss] Concetti di base, conversione 3003 -> 4326
unknown
17 years ago
Permalink
Salve a tutti,
sto lavorando con dei dati geografici provenienti da alcune schede
ministeriali, e ho un po' di confusione in testa, spero che qualcuno mi
possa aiutare :).

I dati geografici delle schede dovrebbero essere, in realtà ne sono
quasi certo, in formato Gauss-Boaga fuso Ovest, EPSG:3003; il mio
problema è che convertendole in WGS84 EPSG:4326 ed utilizzando le
coordinate con Google Earth ottengo uno scostamento dal punto reale.

Ad esempio, la chiesa di S. Pietro a Tula in provincia di SS:

Coordinate: 1499340 4515144
Conversione in WGS84*: 8.992178 40.785764

*Per la conversione ho utilizzato PostGIS:
SELECT AsText(Transform(GeomFromText('POINT(1499340 4515144)', 3003),
4326));

Se utilizzate Google Earth e zoomate a sufficienza vedete come le
coordinate distino notevolmente dalla chiesa di S. Pietro che è
evidenziata con un punto rosso.

Non fidandomi di Google Earth ho inserito le coordinate in una tabella
Postgres (+PostGIS) e utilizzando QGIS le ho sovrapposte alle ortofoto
del WMS della regione Sardegna; anche in questo caso ho riscontrato lo
stesso scostamento rilevato con Google Earth (nell'ordine dei 100m).

WMS:
http://webgis.regione.sardegna.it/wmsconnector/com.esri.wms.Esrimap/ras_wms?request=getcapabilities&service=WMS&version=1.1.1
Layer per le Ortofoto: tutti gli Ortofoto IT2006

La cosa che mi lascia molto perplesso, e che mi fa sospettare di non
aver capito niente :), è che se utilizzo CartLab e converto le
coordinate Gauss-Boaga da piane a geografiche rimanendo nel sistema di
riferimento ROMA40 e poi converto queste nel sistema WGS84 ottengo delle
coordinate che se viste in Google Earth o sovrapposte al WMS della
Sardegna combaciano (a parte un errore che per me potrebbe essere
trascurabile) con il punto cercato.

Esempio:

Coordinate: 1499340 4515144
Coordinate geografiche in ROMA40: 8.992178 40.78657
Coordinate geografiche in WGS84: 8.99181 40.78720

Ammetto di avere un po' di confusione in testa, se l'esperimento con
CartLab fosse fallito avrei pensato ad un errore di rilevamento (anche
se un po' grossolano), ma in questo caso non so bene cosa pensare;
credevo che con i dati georiferiti accompagnati sempre dall'indicazione
del sistema di riferimento, QGIS (o un altro qualsiasi software GIS) mi
avrebbe sovrapposto correttamente i punti convertendo le coordinate da
un sistema all'altro, cosa mi sfugge?

Grazie a tutti per le delucidazioni
david
unknown
17 years ago
Permalink
Ciao David,

L'SRID 3003 non sa niente su come trasformasi in WGS. Ne ci sono towgs
(3 o 7 parametri) ne riferimento a una griglia.

Se ho capito bene, ci sono alcuni parametri in giro (o per la Sardegna
forse si trovano oppure nella data base di EPGS?) ma al meno sul
"mainland" italiano non sono precisi. (Ma penso che isole potrebbero
essere messi meglio ;-)

Per una conversione precisa, occorre conoscere le deformazioni locali
della rete geodetica del IGM95 come incorporato nel programma traspunto
(prodotto dal MATTM e dovrebbe essere disponibili su richiesta per enti
pubblici--o illegalmente dal web (cerca con google)) oppure piu'
ufficiale con le griglie necessarie dell'IGMI usato da software come
Verto (del IGMI) o CartaLab (e forse altri).

Per poter usare le griglie IGMI da proj.4 e di consequenza tutti tipo
di sw open source (incl PostGIS), ho scritto un primo prototipo di un
convertitore del formato del IGMI usato per le loro griglie al format
"standard de fatto" NTv2 che puo essere letto da proj.4.

Il software dovrebbe funzionare, il formato e' corretto (verificato
con una applicazione Australiana), ma non ho avuto tempo di fare prove
con proj.4 e comparare al output del Verto. Per questo ancora non ho
pubblicato il codice (che dovrebbe essere GPL).

Se hai le griglie IGMI del tuo territorio e vuoi provarlo, te lo mando
molto volentiermente.

saluti
-b



On Thu, 13 Mar 2008 15:48:40 +0100
...
--
Bud P. Bruegger, Ph.D. +39-0564-488577 (voice), -21139 (fax)
European Chair, Global Collaboration Forum on eID
Chair, Porvoo Subgroup on collab. govs/operating systems
Leader of the Permanent eID Status Observatory (PESO) project
Servizio Elaborazione Dati e-mail: bud at comune.grosseto.it
Comune di Grosseto jabber: bud at jabber.no
Via Ginori, 43 http://www.comune.grosseto.it/
58100 Grosseto (Tuscany, Italy)
http://www.comune.grosseto.it/interopEID/
unknown
17 years ago
Permalink
Ciao,
prova a creare un nuovo SRID in PostGIS usando questi parametri:

+proj=tmerc +lat_0=0 +lon_0=9 +k=0.999600 +x_0=1500000 +y_0=0

+ellps=intl +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68


Che corrispondono alla proiezione 3003 con i fattori di correzione

Noi abbiamo avuto tempo fa lo stesso problema ed abbiamo risolto con
questa procedura.

ciao a tutti
Leonardo
...
unknown
17 years ago
Permalink
Post by unknown
Ciao,
+proj=tmerc +lat_0=0 +lon_0=9 +k=0.999600 +x_0=1500000 +y_0=0
+ellps=intl +units=m +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68
Che corrispondono alla proiezione 3003 con i fattori di correzione
Leo, questi però sono i parametri di correzione validi per l'Italia
peninsulare...
Io userei invece i parametri specifici per la Sardegna (EPSG code: 1662):
+towgs84=-168.6,-34,38.6,-0.374,-0.679,-1.379,-9.48

(http://www.epsg-registry.org/ --> Retrieve by code --> 1662)

ciao
Antonio
unknown
17 years ago
Permalink
[snip]
Post by unknown
Leo, questi però sono i parametri di correzione validi per l'Italia
peninsulare...
Questi parametri di correzione, come sono ottenuti?
I parametri della trasformazione indicati da Leo sono reperibili da
http://www.epsg-registry.org/ --> Retrieve by code --> 1660
Post by unknown
+towgs84=-168.6,-34,38.6,-0.374,-0.679,-1.379,-9.48
(http://www.epsg-registry.org/ --> Retrieve by code --> 1662)
Qua l'accuratezza dichiarata è sui 4m.
Interessante, io utilizzavo http://spatialreference.org/ e il 1662 non è
presente.
Semplicemente perchè in Spatial Reference sono catalogati solo i sistemi
di coordinate e non le trasformazioni, a differenza del database EPSG
originale che identifica in maniera univoca tutti gli elementi
caratterizzanti tali sistemi (ellissoidi, area geografica di interesse,
trasformazioni, ecc.).
Domanda su PostGIS, mi sono creato l'SRS 93003, in pratica un 3003 con
il campo proj4text modificato secondo la mail di Leonardo.
Ho appena finito l'inserimento dei dati (GeomFromText('XXXX', 93003)) in
una colonna che ha come SRSID di default il 4326, se adesso modifico
l'SRS 93003 con le tue correzioni, PostGIS mi mostra i cambiamenti al
volo?
Dovrebbe... non ho mai sperimentato le trasformazioni in PostGIS,
poichè generalmente preferisco controllarle a livello client e non di
DBMS, ma dato che i parametri sono quelli la trasformazione andrà bene
nei limiti della sua accuratezza.
Facci sapere se ottieni dei risultati accettabili ai tuoi fini.

Ciao
Antonio
unknown
17 years ago
Permalink
On Fri, 14 Mar 2008 09:57:36 +0100
...
Se fai una view con ST_Transform, e' buono di definire anche un indice:

<quote src=
"http://postgis.refractions.net/pipermail/postgis-users/2007-September/017197.html"

PostGIS is pretty
good in doing transforms on the fly, so using a view is a good
approach. Keep in mind that you may have to create a functional index
on your multipoint geometry column if you want to use an index in
Albers. CREATE INDEX new_pt_tbl_idx ON new_pt_tbl USING GIST ON
(ST_Transform(my_multi_pt, <Albers SRID>));
</quote>

o questo:

<quote
src="http://www.postgresonline.com/journal/index.php?/archives/9-How-to-create-an-index-based-on-a-function.html">
Here is another example taken from PostGIS land. Often times you
provide your data in various transformation, but for space savings and
row seek reasons, you want to only transform your data to the less used
projections as needed. One way to do this is to create functional
indexes on the commonly used transformations and create views or just
write raw SQL that uses these alternative transformations.


CREATE INDEX parcels_idx_geom_nadm
ON parcels USING gist(ST_Transform(the_geom, 26986))

So now when I do a select like this that lists all buildings within 100
meters of my NAD 83 MA Meter State Plane point of interest:


SELECT bldg_name, ST_Transform(the_geom,26986) As newgeom
FROM buildings
WHERE ST_DWithin(ST_Transform(the_geom, 26986) ,ST_GeomFromText
('POINT(235675.754215375 894022.495855985)', 26986), 100)

it will use indexes
</quote>
--
Bud P. Bruegger, Ph.D. +39-0564-488577 (voice), -21139 (fax)
European Chair, Global Collaboration Forum on eID
Chair, Porvoo Subgroup on collab. govs/operating systems
Leader of the Permanent eID Status Observatory (PESO) project
Servizio Elaborazione Dati e-mail: bud at comune.grosseto.it
Comune di Grosseto jabber: bud at jabber.no
Via Ginori, 43 http://www.comune.grosseto.it/
58100 Grosseto (Tuscany, Italy)
http://www.comune.grosseto.it/interopEID/
unknown
17 years ago
Permalink
Sto cercando di sintetizzare in un micro-howto questo thread:
http://jgrasstechtips.blogspot.com/2008/03/how-to-find-and-use-bursa-wolf.html
Commenti e correzioni ben accette.

Mi manca un link: come siete arrivati a trovare 1660, 1662 a partire dal 3003?
C'e' una procedura?

Andrea
...
Continua a leggere su narkive:
Loading...