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
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
http://www.spatialys.com
My software is free, but my time generally not.
_______________________________________________
gdal-dev mailing list
[email protected]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
