[Gfoss] Problema conversione coordinate

Salve,

ho un problema nella conversione di coordinate da cartografiche Gauss
Boaga a geografiche WGS84 (sto lavorando con la CTR della regione
Friuli Venezia Giulia). Volevo effettuare la conversione utilizzando
la libreria OSR sotto Python. Ecco il code fragment:

def tolatlong(points):
    gb = osgeo.osr.SpatialReference()
    gb.ImportFromEPSG(3004)
    wgs84 = osgeo.osr.SpatialReference()
    wgs84.ImportFromEPSG(4326)
    tr = osgeo.osr.CoordinateTransformation(gb, wgs84)
    for tp in points:
        p = tp[1]
        z = p['z']
        if z == 999999:
            z = 0
        else:
            z = z / 100.0
        (lon, lat, h) = tr.TransformPoint((p['x'] + 200000000) / 100.0,
                                          (p['y'] + 500000000) / 100.0,
                                          z)
        p['lat'] = lat
        p['lon'] = lon
        p['h'] = h

Ma il risultato è impreciso (dell'ordine dei 50 metri se poi vado a
vedere rispetto ad una carta correttamente georeferenziata). Guardando
il WKT per l'EPSG 3004 non mi è chiaro perché questa si riferisca ad
una proiezione di Mercatore. Volevo passare per una proiezione di UTM
(utilizzando la conversione indicata qua:
http://www.gpscomefare.com/cartografia/cartesucd/ctr.htm) ma poi mi
sono arenato nel spiegare a OSR il WKT per le coordinate cartografiche
UTM33 con datum Roma1940.

Probabilmente da buon neofita sto facendo qualche stupidata immane per
questo vi chiedo qualche suggerimento in merito sul come procedere.
Inoltre anche se c'è un modo diretto mi piacerebbe capire come
impostare il WKT per la proiezione di cui sopra (anche un RTFM mi va
bene).

Grazie!

--
Christian Pellegrin, see http://www.evolware.org/chri/
"Real Programmers don't play tennis, or any other sport which requires
you to change clothes. Mountain climbing is OK, and Real Programmers
wear their climbing boots to work in case a mountain should suddenly
spring up in the middle of the computer room."

On Sat, Nov 1, 2008 at 5:30 PM, chri <chripell@gmail.com> wrote:

Probabilmente da buon neofita sto facendo qualche stupidata immane per
questo vi chiedo qualche suggerimento in merito sul come procedere.
Inoltre anche se c'è un modo diretto mi piacerebbe capire come
impostare il WKT per la proiezione di cui sopra (anche un RTFM mi va
bene).

Ho provato a salvare lo shapefile con le coordinate Gauss Boaga e
convertirlo con:

ogr2ogr -f "ESRI Shapefile" -s_srs "roma40.srs" -t_srs "EPSG:4326"
wgs84_shp IP000GM.shp

dove roma40.srs è

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.999600 +x_0=2520000 +y_0=0
+ellps=intl +units=m
+towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

e tutto va alla grande. Non capisco perché il seguente frammento di codice:

import osgeo.osr
import osgeo.ogr
gb = osgeo.osr.SpatialReference()
gb.ImportFromEPSG(3004)
wgs84 = osgeo.osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
tr = osgeo.osr.CoordinateTransformation(gb, wgs84)
tr.TransformPoint(2365622, 5118141, 0)

non è equivalente. Qualche suggerimento?

--
Christian Pellegrin, see http://www.evolware.org/chri/
"Real Programmers don't play tennis, or any other sport which requires
you to change clothes. Mountain climbing is OK, and Real Programmers
wear their climbing boots to work in case a mountain should suddenly
spring up in the middle of the computer room."

On Sat, Nov 1, 2008 at 7:05 PM, chri <chripell@gmail.com> wrote:

gb.ImportFromEPSG(3004)

se sostituisco questo con:

gb.ImportFromProj4("+proj=tmerc +lat_0=0 +lon_0=15 +k=0.999600
+x_0=2520000 +y_0=0 +ellps=intl +units=m
+towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68")

va tutto a posto

--
Christian Pellegrin, see http://www.evolware.org/chri/
"Real Programmers don't play tennis, or any other sport which requires
you to change clothes. Mountain climbing is OK, and Real Programmers
wear their climbing boots to work in case a mountain should suddenly
spring up in the middle of the computer room."

chri ha scritto:

Ho provato a salvare lo shapefile con le coordinate Gauss Boaga e
convertirlo con:

ogr2ogr -f "ESRI Shapefile" -s_srs "roma40.srs" -t_srs "EPSG:4326"
wgs84_shp IP000GM.shp

dove roma40.srs è

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.999600 +x_0=2520000 +y_0=0
+ellps=intl +units=m
+towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

e tutto va alla grande. Non capisco perché il seguente frammento di codice:

import osgeo.osr
import osgeo.ogr
gb = osgeo.osr.SpatialReference()
gb.ImportFromEPSG(3004)
wgs84 = osgeo.osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
tr = osgeo.osr.CoordinateTransformation(gb, wgs84)
tr.TransformPoint(2365622, 5118141, 0)

non è equivalente. Qualche suggerimento?

Non è equivalente perchè EPSG:3004 è pari a

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.999600 +x_0=2520000 +y_0=0
+ellps=intl +units=m
+towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

meno +towgs84=-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68

cioè

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.999600 +x_0=2520000 +y_0=0
+ellps=intl +units=m

Ciao
Antonio

PS: Scusa la risposta tardiva.

On Mon, Nov 24, 2008 at 7:21 PM, Antonio Falciano <afalciano@yahoo.it> wrote:

Non è equivalente perchè EPSG:3004 è pari a

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.999600 +x_0=2520000 +y_0=0
+ellps=intl +units=m

Quindi quel EPSG si riferisce alla sola proiezione. Senza nessuna
informazione sul map datum (che sarà in un altro EPSG suppongo). Ho
capito bene?

PS: Scusa la risposta tardiva.

Grazie! Meglio tardi che mai!

--
Christian Pellegrin, see http://www.evolware.org/chri/
"Real Programmers don't play tennis, or any other sport which requires
you to change clothes. Mountain climbing is OK, and Real Programmers
wear their climbing boots to work in case a mountain should suddenly
spring up in the middle of the computer room."

christian pellegrin ha scritto:

On Mon, Nov 24, 2008 at 7:21 PM, Antonio Falciano wrote:

Non è equivalente perchè EPSG:3004 è pari a

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.999600 +x_0=2520000 +y_0=0
+ellps=intl +units=m

Quindi quel EPSG si riferisce alla sola proiezione. Senza nessuna
informazione sul map datum (che sarà in un altro EPSG suppongo). Ho
capito bene?

No, EPSG:3004 comprende al suo interno due informazioni: il *datum*
EPSG:4265 (+ellps=intl +pm=greenwich, ovvero l'ellissoide Internazionale
con primo meridiano passante per Greenwich) e la *proiezione* Gauss-Boaga,
data dall'insieme della proiezione Trasversa di Mercatore (+proj=tmerc),
coordinate del centro di emanazione (+lat_0=0 +lon_0=15), fattore di
contrazione (+k=0.999600) e false origini (+x_0=2520000 +y_0=0).

ciao
Antonio