<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Steve,</p>
<p>you may try PR <a class="moz-txt-link-freetext" href="https://github.com/OSGeo/PROJ/pull/3126">https://github.com/OSGeo/PROJ/pull/3126</a> or
<a class="moz-txt-link-freetext" href="https://github.com/OSGeo/PROJ/pull/3127">https://github.com/OSGeo/PROJ/pull/3127</a>, 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. <br>
</p>
<p>Even<br>
</p>
<div class="moz-cite-prefix">Le 23/03/2022 à 23:42, Steve Riddell a
écrit :<br>
</div>
<blockquote type="cite"
cite="mid:CY4PR01MB2406A2CE96E342DD7E5D795FA6189@CY4PR01MB2406.prod.exchangelabs.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta name="Generator" content="Microsoft Word 15 (filtered
medium)">
<style>@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}@font-face
{font-family:Consolas;
panose-1:2 11 6 9 2 2 4 3 2 4;}@font-face
{font-family:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}@font-face
{font-family:"Cascadia Mono";
panose-1:2 11 6 9 2 0 0 2 0 4;}p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}pre
{mso-style-priority:99;
mso-style-link:"HTML Preformatted Char";
margin:0in;
font-size:10.0pt;
font-family:"Courier New";}span.HTMLPreformattedChar
{mso-style-name:"HTML Preformatted Char";
mso-style-priority:99;
mso-style-link:"HTML Preformatted";
font-family:Consolas;}span.EmailStyle23
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:windowtext;}.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}div.WordSection1
{page:WordSection1;}</style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
<div class="WordSection1">
<p class="MsoNormal">Hi Even,<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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.
<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">COMPD_CS["TempTM + CGVD28 height - HT2_0",<o:p></o:p></p>
<p class="MsoNormal"> PROJCS["Custom",<o:p></o:p></p>
<p class="MsoNormal"> GEOGCS["NAD83(CSRS)",<o:p></o:p></p>
<p class="MsoNormal">
DATUM["NAD83_Canadian_Spatial_Reference_System",<o:p></o:p></p>
<p class="MsoNormal"> SPHEROID["GRS
1980",6378137,298.257222101,<o:p></o:p></p>
<p class="MsoNormal">
AUTHORITY["EPSG","7019"]],<o:p></o:p></p>
<p class="MsoNormal"> TOWGS84[0,0,0,0,0,0,0],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","6140"]],<o:p></o:p></p>
<p class="MsoNormal"> PRIMEM["Greenwich",0,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","8901"]],<o:p></o:p></p>
<p class="MsoNormal">
UNIT["degree",0.0174532925199433,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","9122"]],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","4617"]],<o:p></o:p></p>
<p class="MsoNormal"> PROJECTION["Transverse_Mercator"],<o:p></o:p></p>
<p class="MsoNormal">
PARAMETER["latitude_of_origin",49.351346659616],<o:p></o:p></p>
<p class="MsoNormal">
PARAMETER["central_meridian",-123.20266499149],<o:p></o:p></p>
<p class="MsoNormal"> PARAMETER["scale_factor",1],<o:p></o:p></p>
<p class="MsoNormal">
PARAMETER["false_easting",15307.188],<o:p></o:p></p>
<p class="MsoNormal">
PARAMETER["false_northing",6540.975],<o:p></o:p></p>
<p class="MsoNormal"> UNIT["Meters",1],<o:p></o:p></p>
<p class="MsoNormal"> AXIS["Easting",EAST],<o:p></o:p></p>
<p class="MsoNormal"> AXIS["Northing",NORTH]],<o:p></o:p></p>
<p class="MsoNormal"> VERT_CS["CGVD28 height - HT2_0",<o:p></o:p></p>
<p class="MsoNormal"> VERT_DATUM["Canadian Geodetic
Vertical Datum of 1928",2005,<o:p></o:p></p>
<p class="MsoNormal">
EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","5114"]],<o:p></o:p></p>
<p class="MsoNormal"> UNIT["metre",1,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","9001"]],<o:p></o:p></p>
<p class="MsoNormal"> AXIS["Gravity-related height",UP],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","5713"]]]<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Experimenting, I removed the “TOWGS84”
directive from BOTH the source and target systems, and I do
get an answer instead of an error:<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">For test point: (<span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black">49.351346658, -123.202664992,
289.809)<o:p></o:p></span></p>
<p class="MsoNormal"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black">I get:
(15307.1880,6540.9748,307.744)</span><o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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:<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal" style="text-autospace:none"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black">strPipe
</span><span style="font-size:9.5pt;font-family:"Cascadia
Mono";color:teal">=</span><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black">
</span><span style="font-size:9.5pt;font-family:"Cascadia
Mono";color:#A31515">"+proj=pipeline\</span><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black"><o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:#A31515"> +step +proj=axisswap
+order=2, 1\</span><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black"><o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:#A31515"> +step +proj=unitconvert
+xy_in=deg +xy_out=rad\</span><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black"><o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:#A31515"> +step +inv
+proj=vgridshift +grids=HT2_0.gtx +multiplier=1\</span><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black"><o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:#A31515"> +step +proj=tmerc
+lat_0=49.351346659616 +lon_0=-123.20266499149 +k=1\</span><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black"><o:p></o:p></span></p>
<p class="MsoNormal"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:#A31515"> +x_0=15307.188
+y_0=6540.975 +ellps=GRS80"</span><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black">;<o:p></o:p></span></p>
<p class="MsoNormal"><span
style="font-size:9.5pt;font-family:"Cascadia
Mono";color:black"><o:p> </o:p></span></p>
<p class="MsoNormal">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.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Any arrows left in your quiver?<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Thanks,<o:p></o:p></p>
<p class="MsoNormal">Steve<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div style="border:none;border-top:solid #E1E1E1
1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal"><b>From:</b> Steve Riddell <br>
<b>Sent:</b> Friday, March 18, 2022 6:17 PM<br>
<b>To:</b> Even Rouault
<a class="moz-txt-link-rfc2396E" href="mailto:even.rouault@spatialys.com"><even.rouault@spatialys.com></a>;
<a class="moz-txt-link-abbreviated" href="mailto:gdal-dev@lists.osgeo.org">gdal-dev@lists.osgeo.org</a><br>
<b>Subject:</b> RE: OGRCreateCoordinateTransformation()<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">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.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div style="border:none;border-top:solid #E1E1E1
1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal"><b>From:</b> Even Rouault <<a
href="mailto:even.rouault@spatialys.com"
moz-do-not-send="true" class="moz-txt-link-freetext">even.rouault@spatialys.com</a>>
<br>
<b>Sent:</b> Friday, March 18, 2022 5:05 PM<br>
<b>To:</b> Steve Riddell <<a
href="mailto:sriddell@geocue.com" moz-do-not-send="true"
class="moz-txt-link-freetext">sriddell@geocue.com</a>>;
<a href="mailto:gdal-dev@lists.osgeo.org"
moz-do-not-send="true" class="moz-txt-link-freetext">gdal-dev@lists.osgeo.org</a><br>
<b>Subject:</b> Re: OGRCreateCoordinateTransformation()<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<p>Steve,<o:p></o:p></p>
<p>Fix in PROJ queued in <a
href="https://github.com/OSGeo/PROJ/pull/3123"
moz-do-not-send="true" class="moz-txt-link-freetext">https://github.com/OSGeo/PROJ/pull/3123</a><o:p></o:p></p>
<p>The complexity of dealing with legacy features TOWGS84[] and
PROJ4_GRIDS is highly stressing PROJ pipeline computation
engine.<o:p></o:p></p>
<p>Using the plain EPSG codes EPSG:4955 -> EPSG:4617+5713
would be much preferable here<o:p></o:p></p>
<p>Even<o:p></o:p></p>
<div>
<p class="MsoNormal">Le 18/03/2022 à 21:14, Steve Riddell a
écrit :<o:p></o:p></p>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<p class="MsoNormal">Hi,<o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal">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.<o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal">COMPD_CS["NAD83(CSRS) + Ellipsoid
(Meters)",<o:p></o:p></p>
<p class="MsoNormal"> GEOGCS["NAD83(CSRS)",<o:p></o:p></p>
<p class="MsoNormal">
DATUM["NAD83_Canadian_Spatial_Reference_System",<o:p></o:p></p>
<p class="MsoNormal"> SPHEROID["GRS
1980",6378137,298.257222101,<o:p></o:p></p>
<p class="MsoNormal">
AUTHORITY["EPSG","7019"]],<o:p></o:p></p>
<p class="MsoNormal"> TOWGS84[0,0,0,0,0,0,0],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","6140"]],<o:p></o:p></p>
<p class="MsoNormal"> PRIMEM["Greenwich",0,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","8901"]],<o:p></o:p></p>
<p class="MsoNormal"> UNIT["degree",0.0174532925199433,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","9122"]],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","4617"]],<o:p></o:p></p>
<p class="MsoNormal"> VERT_CS["Ellipsoid (Meters)",<o:p></o:p></p>
<p class="MsoNormal"> VERT_DATUM["Ellipsoid",2002],<o:p></o:p></p>
<p class="MsoNormal"> UNIT["metre",1,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","9001"]],<o:p></o:p></p>
<p class="MsoNormal"> AXIS["ellipsoidal height",UP]]]<o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal">COMPD_CS["NAD83(CSRS) + CGVD28 height -
HT2_0",<o:p></o:p></p>
<p class="MsoNormal"> GEOGCS["NAD83(CSRS)",<o:p></o:p></p>
<p class="MsoNormal">
DATUM["NAD83_Canadian_Spatial_Reference_System",<o:p></o:p></p>
<p class="MsoNormal"> SPHEROID["GRS
1980",6378137,298.257222101,<o:p></o:p></p>
<p class="MsoNormal">
AUTHORITY["EPSG","7019"]],<o:p></o:p></p>
<p class="MsoNormal"> TOWGS84[0,0,0,0,0,0,0],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","6140"]],<o:p></o:p></p>
<p class="MsoNormal"> PRIMEM["Greenwich",0,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","8901"]],<o:p></o:p></p>
<p class="MsoNormal"> UNIT["degree",0.0174532925199433,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","9122"]],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","4617"]],<o:p></o:p></p>
<p class="MsoNormal"> VERT_CS["CGVD28 height - HT2_0",<o:p></o:p></p>
<p class="MsoNormal"> VERT_DATUM["Canadian Geodetic
Vertical Datum of 1928",2005,<o:p></o:p></p>
<p class="MsoNormal">
EXTENSION["PROJ4_GRIDS","HT2_0.gtx"],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","5114"]],<o:p></o:p></p>
<p class="MsoNormal"> UNIT["metre",1,<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","9001"]],<o:p></o:p></p>
<p class="MsoNormal"> AXIS["Gravity-related
height",UP],<o:p></o:p></p>
<p class="MsoNormal"> AUTHORITY["EPSG","5713"]]]<o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal">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]]]]]'<o:p></o:p></p>
<pre> <o:p></o:p></pre>
<pre><span style="font-family:"Calibri",sans-serif">Thanks in advance for any assistance!<o:p></o:p></span></pre>
<pre><span style="font-family:"Calibri",sans-serif">Best regards,</span><o:p></o:p></pre>
<pre><span style="font-family:"Calibri",sans-serif">Steve</span><o:p></o:p></pre>
<pre><span style="font-family:"Calibri",sans-serif"> </span><o:p></o:p></pre>
<p class="MsoNormal"> <o:p></o:p></p>
<p class="MsoNormal"><span
style="font-size:13.5pt;font-family:"Tahoma",sans-serif">Steve
Riddell</span><o:p></o:p></p>
<p class="MsoNormal"><span
style="font-family:"Tahoma",sans-serif">Software
Engineering</span><o:p></o:p></p>
<p class="MsoNormal"><b><span style="color:#525252">GeoCue
Group, Inc.</span></b><o:p></o:p></p>
<p class="MsoNormal"><span style="font-size:10.0pt">520 6th
Street | Madison, AL 35756 USA</span><o:p></o:p></p>
<p class="MsoNormal"><b><span style="font-size:10.0pt">Phone:</span></b><span
style="font-size:10.0pt"> 256.461.8289 |
<b>Fax:</b> 256.461.8249 </span><o:p></o:p></p>
<p class="MsoNormal"><i><span style="color:#0070C0">LIDAR/Drone
Mapping Software & Services
</span></i><i><span style="font-size:9.0pt">– </span><span
style="color:#538135">True View® 3D Imaging Sensors
</span></i><i><span style="font-size:9.0pt">–</span><span
style="color:#C45911"> Image/LIDAR Data Management
Solutions</span></i><o:p></o:p></p>
<p class="MsoNormal"><span style="font-size:10.0pt"> </span><a
href="http://www.geocue.com/" moz-do-not-send="true"><span
style="font-size:10.0pt;color:#0070C0">www.geocue.com</span></a><span
style="font-size:10.0pt;color:#0070C0">
<a href="https://support.geocue.com/"
moz-do-not-send="true">support.geocue.com</a> </span><o:p></o:p></p>
<p class="MsoNormal"><span style="font-size:10.0pt"> </span><o:p></o:p></p>
<p class="MsoNormal"> <o:p></o:p></p>
<pre><span style="font-size:11.0pt;font-family:"Calibri",sans-serif"> </span><o:p></o:p></pre>
</blockquote>
<pre>-- <o:p></o:p></pre>
<pre><a href="http://www.spatialys.com" moz-do-not-send="true" class="moz-txt-link-freetext">http://www.spatialys.com</a><o:p></o:p></pre>
<pre>My software is free, but my time generally not.<o:p></o:p></pre>
</div>
</blockquote>
<pre class="moz-signature" cols="72">--
<a class="moz-txt-link-freetext" href="http://www.spatialys.com">http://www.spatialys.com</a>
My software is free, but my time generally not.</pre>
</body>
</html>