<div dir="ltr"><div>Hi Even,</div><div><br></div><div>Thank you for your detailed reply.</div><div>So if I want to use 1169 or also 1660, that is the first candidate listed by</div><div><br></div><div>
<span style="font-family:monospace">projinfo -s EPSG:3003 -t EPSG:4326 --spatial-test intersects</span>

</div><div><br></div><div>I can transform the coordinates from 3003 to 4265, that is the geodetic SRS used by 3003</div><div><br></div><div><span style="font-family:monospace">sourceProj = osr.SpatialReference()<br>sourceProj.ImportFromEPSG(3003)<br><br>sourceGeo = osr.SpatialReference()<br>sourceGeo.ImportFromEPSG(4265)<br><br>proj2geo = osr.CoordinateTransformation(sourcePro,sourceGeo)</span></div><div><span style="font-family:monospace"><br></span></div><div><span style="font-family:monospace">point.Transform(proj2geo)<br></span></div><div><br></div><div>then I can use EPSG:1169. I tried the code snippet and it seems to be working properly.</div><div>Do you think this is a correct approach?</div><div><br></div><div>All the best</div><div><br></div><div>Stefano<br></div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno ven 5 mar 2021 alle ore 12:11 Even Rouault <<a href="mailto:even.rouault@spatialys.com">even.rouault@spatialys.com</a>> ha scritto:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
  
    
  
  <div>
    <p>Stefano,</p>
    <p><span style="font-family:monospace">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. <br>
      </span></p>
    <p>You could instead use</p>
    <p><span style="font-family:monospace"><span style="font-family:monospace">options.SetOperation(</span>"""+proj=pipeline<br>
          +step +inv +proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996
        +x_0=1500000 +y_0=0<br>
                +ellps=intl<br>
          +step +proj=push +v_3<br>
          +step +proj=cart +ellps=intl<br>
          +step +proj=helmert +x=-225 +y=-65 +z=9<br>
          +step +inv +proj=cart +ellps=WGS84<br>
          +step +proj=pop +v_3<br>
          +step +proj=unitconvert +xy_in=rad +xy_out=deg<br>
          +step +proj=axisswap +order=2,1""")<br>
      </span></p>
    <p><span style="font-family:monospace">as suggested by the output of
        projinfo -s EPSG:3003 -t EPSG:4326 --spatial-test intersects</span></p>
    <p><span style="font-family:monospace">Even<br>
      </span></p>
    <div>Le 05/03/2021 à 09:42, Stefano
      Iacovella a écrit :<br>
    </div>
    <blockquote type="cite">
      
      <div dir="ltr">
        <div>Hi all,</div>
        <div><br>
        </div>
        <div>I am trying to use GDAL, via Python bindings, to convert
          coordinates.</div>
        <div>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.</div>
        <div>I am using binaries from OSGEO4WIN, the release are 3.1.4
          for GDAL and 6.3.2 for PROJ</div>
        <div>This is a code snippet I am using to debug:</div>
        <div><br>
        </div>
        <div><span style="font-family:monospace">from osgeo import ogr<br>
            from osgeo import osr<br>
            <br>
            source = osr.SpatialReference()<br>
            source.ImportFromEPSG(3003)<br>
            <br>
            target = osr.SpatialReference()<br>
            target.ImportFromEPSG(4326)<br>
            <br>
            options = osr.CoordinateTransformationOptions()<br>
options.SetOperation("urn:ogc:def:coordinateOperation:EPSG::1169")<br>
            transform = osr.CoordinateTransformation(source, target,
            options)<br>
            <br>
            point = ogr.CreateGeometryFromWkt("POINT (1622987 4756952)")<br>
            point.Transform(transform)<br>
            <br>
            print point.ExportToWkt()</span></div>
        <div><br>
        </div>
        <div>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:</div>
        <div><br>
        </div>
        <div>
          <span style="font-family:monospace">transform =
            osr.CoordinateTransformation(source, target)</span><br>
        </div>
        <div><br>
        </div>
        <div>The point coordinates after transformation are ok but of
          course the conversion is performed using the Coordinate
          Operation automatically selected by the system.</div>
        <div>Am I doing something wrong?</div>
        <div><br>
        </div>
        <div>Thank you in advance for your help.</div>
        <div><br>
        </div>
        <div>All the best,</div>
        <div><br>
        </div>
        <div>Stefano<br>
        </div>
        <div><br>
        </div>
      </div>
      <br>
      <fieldset></fieldset>
      <pre>_______________________________________________
gdal-dev mailing list
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a>
<a href="https://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">https://lists.osgeo.org/mailman/listinfo/gdal-dev</a>
</pre>
    </blockquote>
    <pre cols="72">-- 
<a href="http://www.spatialys.com" target="_blank">http://www.spatialys.com</a></pre>
  </div>

</blockquote></div>