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

[Even Rouault]

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/6abcdde4/attachment-0001.html>