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

adamgutonski adamgutonski at protonmail.com
Thu Feb 17 06:54:59 PST 2022


Hello,

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:

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20220217/67062b98/attachment.html>


More information about the gdal-dev mailing list