[gdal-dev] Matching ESRI projection/transformations using GDAL/OSR

[adamgutonski]

Thank you so much! This appears to have resolved the issue.

Sent with [ProtonMail](https://protonmail.com/) Secure Email.

------- Original Message -------
On Thursday, February 17th, 2022 at 10:36 AM, Even Rouault <even.rouault at spatialys.com> wrote:

> Hi,
>
>> I am attempting to project and transform geometry using OSR. The projection/transformation is from ITRF2014 via Proj4 string to NAD83(CSRS)v7 (EPSG:8255). My problem is that I need the transformation output to match the output I receive when transforming the geometry in ESRI using the geometry.projectAs function. It seems that the OSR transformation creates a point that is always shifted from the ESRI geometry. Any advice on how to match projections and transformations between open source and ESRI would be extremely helpful.
>>
>> Original coordinate: -74.442799450000, 41.406811727778
>>
>> ESRI point after transform: -74.4427951819612, 41.4068026527349
>>
>> OGR point after transform: -74.44281012167299, 41.40671904426617
>>
>> Code used to generate the OGR point:
>
> You've likely used the coordinate order wrong. GDAL 3 honours by default authority (EPSG) axis order, so for EPSG this means lat/lon for geographic CRS.
>
> So either switch your coordinate, or call itrf.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) and to_nad.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
>
> With that, you'll get the same results as Esri software
>
>>>> print(outpoint)
> POINT (-74.442810121673 41.4067190442662 0)
>
> Even
>
>> itrf = osr.SpatialReference()
>> itrf.ImportFromProj4(
>>
>> 'GEOGCS["ITRF2014",DATUM["International_Terrestrial_Reference_Frame_2014",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433],AUTHORITY["EPSG",9000]]'
>>
>> )
>>
>> to_nad = osr.SpatialReference()
>> to_nad.ImportFromEPSG(
>>
>> 8255
>>
>> )
>>
>> opts = osr.CoordinateTransformationOptions()
>> opts.SetOperation(
>>
>> "ITRF2014_To_NAD_1983_CSRS_v7_7par"
>>
>> )
>> transformer = osr.CreateCoordinateTransformation(itrf, to_nad, opts)
>>
>> outpoint = ogr.Geometry(ogr.wkbPoint)
>> outpoint.AddPoint(*test_coords)
>> outpoint.AssignSpatialReference(itrf)
>> outpoint.Transform(transformer)
>>
>> With ESRI I am using the ArcPy geometry function, projectAs when taking in the name of a transformation. When ArcGIS reprojects from ITRF2014 to NAD83(CSRS)v7 (EPSG:8255) with the transformation "ITRF2014_To_NAD_1983_CSRS_v7_7par", a point is generated ~30 ft away from the point generated by GDAL/OSR
>>
>> Code used to generate the Esri point:
>>
>> ITRF_SPATIAL_REF = arcpy.SpatialReference(
>>
>> 'ITRF2014.prj'
>>
>> )
>> NAD_SPATIAL_REF= arcpy.SpatialReference(
>>
>> 8255
>>
>> )
>>
>> longitude = test_coords[
>>
>> 0
>>
>> ]
>> latitude = test_coords[
>>
>> 1
>>
>> ]
>> point = arcpy.Point(
>>     longitude,
>>     latitude
>> )
>> return_point = arcpy.PointGeometry(point, ITRF_SPATIAL_REF).projectAs(NAD_SPATIAL_REF,
>>
>> "ITRF2014_To_NAD_1983_CSRS_v7_7par"
>>
>> )
>>
>> Thank you for your time,
>> adamg
>>
>> Sent with [ProtonMail](https://protonmail.com/) Secure Email.
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>>
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>
> --
> http://www.spatialys.com
> My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20220217/ff673d75/attachment.html>