1

How can I make sure that the coordinate transformations that I am applying using GDAL/OSR match transformations in ESRI/ArcGIS?

Specifically 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:

Original coordinate: -74.442799450000, 41.406811727778
ESRI point after transform: -74.4427951819612, 41.4068026527349
OGR point after transform: -74.44281012167299, 41.40671904426617

enter image description here

Code used to generate the ODR 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)

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")

1 Answer 1

1

After messaging the gdal-dev mailing list, I received this answer which resolved my issue:

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.