Le 05/03/2021 à 12:54, Stefano Iacovella a écrit :
Hi Even,
Thank you for your detailed reply.
So if I want to use 1169 or also 1660, that is the first candidate
listed by
projinfo -s EPSG:3003 -t EPSG:4326 --spatial-test intersects
I can transform the coordinates from 3003 to 4265, that is the
geodetic SRS used by 3003
sourceProj = osr.SpatialReference()
sourceProj.ImportFromEPSG(3003)
sourceGeo = osr.SpatialReference()
sourceGeo.ImportFromEPSG(4265)
proj2geo = osr.CoordinateTransformation(sourcePro,sourceGeo)
point.Transform(proj2geo)
then I can use EPSG:1169. I tried the code snippet and it seems to be
working properly.
Do you think this is a correct approach?
Yes that's probably the best approach if you don't want to use PROJ
pipelines to do it in one step.
All the best
Stefano
Il giorno ven 5 mar 2021 alle ore 12:11 Even Rouault
<[email protected] <mailto:[email protected]>> ha
scritto:
Stefano,
urn:ogc:def:coordinateOperation:EPSG::1169 is a transformation
between EPSG:4265 (the geographic CRS of Monte Mario) and
EPSG:4326. You cannot use it directly when the source CRS is
EPSG:3003, which is a projected CRS based on EPSG:4265.
You could instead use
options.SetOperation("""+proj=pipeline
+step +inv +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=1500000
+y_0=0
+ellps=intl
+step +proj=push +v_3
+step +proj=cart +ellps=intl
+step +proj=helmert +x=-225 +y=-65 +z=9
+step +inv +proj=cart +ellps=WGS84
+step +proj=pop +v_3
+step +proj=unitconvert +xy_in=rad +xy_out=deg
+step +proj=axisswap +order=2,1""")
as suggested by the output of projinfo -s EPSG:3003 -t EPSG:4326
--spatial-test intersects
Even
Le 05/03/2021 à 09:42, Stefano Iacovella a écrit :
Hi all,
I am trying to use GDAL, via Python bindings, to convert coordinates.
I would like to select a specific transformation when converting
among CRS with different Datums. Reading the documentation my
understanding is that I need to create a
osr.CoordinateTransformation with options parameter valorized
properly.
I am using binaries from OSGEO4WIN, the release are 3.1.4 for
GDAL and 6.3.2 for PROJ
This is a code snippet I am using to debug:
from osgeo import ogr
from osgeo import osr
source = osr.SpatialReference()
source.ImportFromEPSG(3003)
target = osr.SpatialReference()
target.ImportFromEPSG(4326)
options = osr.CoordinateTransformationOptions()
options.SetOperation("urn:ogc:def:coordinateOperation:EPSG::1169")
transform = osr.CoordinateTransformation(source, target, options)
point = ogr.CreateGeometryFromWkt("POINT (1622987 4756952)")
point.Transform(transform)
print point.ExportToWkt()
I have no errors running it but the coordinates of the point,
after the transformation, are infinite. If I create the
CoordinateTransformation without options like this:
transform = osr.CoordinateTransformation(source, target)
The point coordinates after transformation are ok but of course
the conversion is performed using the Coordinate Operation
automatically selected by the system.
Am I doing something wrong?
Thank you in advance for your help.
All the best,
Stefano
_______________________________________________
gdal-dev mailing list
[email protected] <mailto:[email protected]>
https://lists.osgeo.org/mailman/listinfo/gdal-dev
<https://lists.osgeo.org/mailman/listinfo/gdal-dev>
--
http://www.spatialys.com <http://www.spatialys.com>
--
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
