[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