[gdal-dev] OGRCreateCoordinateTransformation()
Even Rouault
even.rouault at spatialys.com
Wed Mar 23 16:06:21 PDT 2022
Steve,
you may try PR https://github.com/OSGeo/PROJ/pull/3126 or
https://github.com/OSGeo/PROJ/pull/3127, that might potentially help,
but they don't make consensus (see discussions of those PRs). Basically
it is a nightmare to deal with the old TOWGS84 / PROJ4_GRIDS approach,
that poorly specifies some elements of transformations (especially the
vertical part), and doesn't fit well in the ISO-19111 model used by
newer PROJ (we use and abuse the BoundCRS approach in complex ways). One
of the issue too is that we've realized late in the process that we
didn't emulate some aspects of the old transformation pipeline, and
there isn't general enthusiasm to change that now.
Even
Le 23/03/2022 à 23:42, Steve Riddell a écrit :
>
> 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>
> *Sent:* Friday, March 18, 2022 5:05 PM
> *To:* Steve Riddell <sriddell at geocue.com>; 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.
--
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/20220324/a47b70c8/attachment-0001.html>
More information about the gdal-dev
mailing list