[gdal-dev] OGRCreateCoordinateTransformation()
Steve Riddell
sriddell at geocue.com
Wed Mar 23 15:42:05 PDT 2022
Hi Even,
I tried the fix with the CRS’s I sent, and it worked. But as I mentioned, my real target CRS is a custom projected system (below), and it still fails with what seems to be the same error as in my original email.
COMPD_CS["TempTM + CGVD28 height - HT2_0",
PROJCS["Custom",
GEOGCS["NAD83(CSRS)",
DATUM["NAD83_Canadian_Spatial_Reference_System",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6140"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4617"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",49.351346659616],
PARAMETER["central_meridian",-123.20266499149],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",15307.188],
PARAMETER["false_northing",6540.975],
UNIT["Meters",1],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]],
VERT_CS["CGVD28 height - HT2_0",
VERT_DATUM["Canadian Geodetic Vertical Datum of 1928",2005,
EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],
AUTHORITY["EPSG","5114"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Gravity-related height",UP],
AUTHORITY["EPSG","5713"]]]
Experimenting, I removed the “TOWGS84” directive from BOTH the source and target systems, and I do get an answer instead of an error:
For test point: (49.351346658, -123.202664992, 289.809)
I get: (15307.1880,6540.9748,307.744)
But since geoid separation at the test point is -18.238, Z should be 308.047. In fact, if I code up the following pipeline and initialize a transform object:
strPipe = "+proj=pipeline\
+step +proj=axisswap +order=2, 1\
+step +proj=unitconvert +xy_in=deg +xy_out=rad\
+step +inv +proj=vgridshift +grids=HT2_0.gtx +multiplier=1\
+step +proj=tmerc +lat_0=49.351346659616 +lon_0=-123.20266499149 +k=1\
+x_0=15307.188 +y_0=6540.975 +ellps=GRS80";
The call to Transform() gives: (15307.1880, 6540.9748, 308.042) -- XY matches above, but now I get the expected Z (-0.005). However, I’d like to avoid assembling pipelines for every conversion/transform.
Any arrows left in your quiver?
Thanks,
Steve
From: Steve Riddell
Sent: Friday, March 18, 2022 6:17 PM
To: Even Rouault <even.rouault at spatialys.com>; gdal-dev at lists.osgeo.org
Subject: RE: OGRCreateCoordinateTransformation()
Thanks very much, Even. For reporting the problem, I “packaged up” a case with CRS’s with EPGS codes. For the actual problem I’m working, the target CRS is a custom projection –no EPSG.
From: Even Rouault <even.rouault at spatialys.com<mailto:even.rouault at spatialys.com>>
Sent: Friday, March 18, 2022 5:05 PM
To: Steve Riddell <sriddell at geocue.com<mailto:sriddell at geocue.com>>; gdal-dev at lists.osgeo.org<mailto:gdal-dev at lists.osgeo.org>
Subject: Re: OGRCreateCoordinateTransformation()
Steve,
Fix in PROJ queued in https://github.com/OSGeo/PROJ/pull/3123
The complexity of dealing with legacy features TOWGS84[] and PROJ4_GRIDS is highly stressing PROJ pipeline computation engine.
Using the plain EPSG codes EPSG:4955 -> EPSG:4617+5713 would be much preferable here
Even
Le 18/03/2022 à 21:14, Steve Riddell a écrit :
Hi,
I’m getting an error (below) using GDAL trying to set up a transformation between NAD83(CSRS) geodetic, w/ ellipsoidal heights, and NAD83(CSRS) geodetic, CGVD28 heights. Similar transformations work fine on some datums, but fail on others. Hope someone can offer some insight. Note the if I don't use the PROJ4_GRIDS extension, I don't get the error, but I don't get the geoid adjustment either.
COMPD_CS["NAD83(CSRS) + Ellipsoid (Meters)",
GEOGCS["NAD83(CSRS)",
DATUM["NAD83_Canadian_Spatial_Reference_System",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6140"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4617"]],
VERT_CS["Ellipsoid (Meters)",
VERT_DATUM["Ellipsoid",2002],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["ellipsoidal height",UP]]]
COMPD_CS["NAD83(CSRS) + CGVD28 height - HT2_0",
GEOGCS["NAD83(CSRS)",
DATUM["NAD83_Canadian_Spatial_Reference_System",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6140"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4617"]],
VERT_CS["CGVD28 height - HT2_0",
VERT_DATUM["Canadian Geodetic Vertical Datum of 1928",2005,
EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],
AUTHORITY["EPSG","5114"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Gravity-related height",UP],
AUTHORITY["EPSG","5713"]]]
ERROR 6: Cannot find coordinate operations from `EPSG:4955' to `COMPOUNDCRS["NAD83(CSRS) + CGVD28 height - HT2_0",BOUNDCRS[SOURCECRS[GEOGCRS["NAD83(CSRS)",DATUM["NAD83 Canadian Spatial Reference System",ELLIPSOID["GRS 1980",6378137,298.257222101,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4617]]],TARGETCRS[GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]]],ABRIDGEDTRANSFORMATION["Transformation from NAD83(CSRS) to WGS84",METHOD["Position Vector transformation (geog2D domain)",ID["EPSG",9606]],PARAMETER["X-axis translation",0,ID["EPSG",8605]],PARAMETER["Y-axis translation",0,ID["EPSG",8606]],PARAMETER["Z-axis translation",0,ID["EPSG",8607]],PARAMETER["X-axis rotation",0,ID["EPSG",8608]],PARAMETER["Y-axis rotation",0,ID["EPSG",8609]],PARAMETER["Z-axis rotation",0,ID["EPSG",8610]],PARAMETER["Scale difference",1,ID["EPSG",8611]]]],BOUNDCRS[SOURCECRS[VERTCRS["CGVD28 height - HT2_0",VDATUM["Canadian Geodetic Vertical Datum of 1928"],CS[vertical,1],AXIS["gravity-related height",up,LENGTHUNIT["metre",1]],ID["EPSG",5713]]],TARGETCRS[GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,3],AXIS["latitude",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["longitude",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],AXIS["ellipsoidal height",up,ORDER[3],LENGTHUNIT["metre",1]],ID["EPSG",4979]]],ABRIDGEDTRANSFORMATION["CGVD28 height - HT2_0 to WGS84 ellipsoidal height",METHOD["GravityRelatedHeight to Geographic3D"],PARAMETERFILE["Geoid (height correction) model file","HT2_0.gtx",ID["EPSG",8666]]]]]'
Thanks in advance for any assistance!
Best regards,
Steve
Steve Riddell
Software Engineering
GeoCue Group, Inc.
520 6th Street | Madison, AL 35756 USA
Phone: 256.461.8289 | Fax: 256.461.8249
LIDAR/Drone Mapping Software & Services – True View® 3D Imaging Sensors – Image/LIDAR Data Management Solutions
www.geocue.com<http://www.geocue.com/> support.geocue.com<https://support.geocue.com/>
--
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/20220323/ead67550/attachment-0001.html>
More information about the gdal-dev
mailing list