<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<style><!--
/* Font Definitions */
@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:Tahoma;
panose-1:2 11 6 4 3 5 4 4 2 4;}
@font-face
{font-family:Consolas;
panose-1:2 11 6 9 2 2 4 3 2 4;}
@font-face
{font-family:"Lucida Sans Unicode";
panose-1:2 11 6 2 3 5 4 2 2 4;}
/* Style Definitions */
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-compose;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;
mso-ligatures:none;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
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]-->
</head>
<body lang="EN-US" link="blue" vlink="purple" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal">We’ve already set it back to the traditional axis order but that code snip wasn’t included. I had forgotten about that change we had to make to keep our existing code as is. I think the differences would be drastically different otherwise,
but we are seeing differences that are mostly in the 5-7 foot range. I have seen only one example of a difference greater than 7 feet. Based on the fact that the NAD27 to WGS84 transforms indicate a precision of only 1.5 meters even with the most precise
method, I could probably make a case that the differences I’m seeing are within that margin but I just want to find a way to prove what method GDAL is selecting under the covers. I need to be able to explain the differences after the upgrade of GDAL and PROJ.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">Going back to some of your earlier suggestions I did have some follow up questions. Here is a reference to some of your earlier comments with my questions initialed.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">1) The result you get with GDAL 1.10 can be reproduced with modern PROJ versions by allowing (and restricting) to using the CONUS / us_noaa_conus.tif grid which corresponds to the "NAD27 to WGS 84 (79)"
<o:p></o:p></p>
<p class="MsoNormal">EPSG transformation, with an advertized accuracy of 5.0 m:<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">$ echo 32.804191 -117.258911 0 | PROJ_NETWORK=ON cct -d 8 +proj=pipeline
<o:p></o:p></p>
<p class="MsoNormal">+step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg
<o:p></o:p></p>
<p class="MsoNormal">+xy_out=rad +step +proj=hgridshift +grids=us_noaa_conus.tif +step
<o:p></o:p></p>
<p class="MsoNormal">+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1<o:p></o:p></p>
<p class="MsoNormal">32.80423938?? -117.25978145??? 0.00000000?????????? inf<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">sfox – I am already obtaining that result with GDAL 1.10 and Proj 4.4.9 by specifying only the source and target CRS and the coordinates as shown in the sample code I showed. This is the result I’d like to get with GDAL 3.4.2, but I cannot
figure out how to get that.<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">sfox – the other two methods require this step, “+grids=us_noaa_conus.tif” and I try to run any of those pipelines through the cct program that is part of proj I get a file not found error so that is perplexing. I see the following output,
and I assume that setting PROJ_NETWORK=ON is supposed to allow the correct grids or files to be found. Why would a file not found or invalid error be presented?<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">PROJ_NETWORK=ON<o:p></o:p></p>
<p class="MsoNormal">echo $PROJ_NETWORK<o:p></o:p></p>
<p class="MsoNormal">ON<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">echo 32.804191 -117.258911 0 | cct -d 8 +proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=hgridshift +grids=us_noaa_conus.tif +step +proj=hgridshift +grids=us_noaa_cshpgn.tif +step
+proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">proj_create: Error 1029 (File not found or invalid): pipeline: Pipeline: Bad step definition: proj=hgridshift (File not found or invalid)<o:p></o:p></p>
<p class="MsoNormal">cct: Bad transformation arguments - (File not found or invalid)<o:p></o:p></p>
<p class="MsoNormal"> 'cct -h' for help<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">2) The result you get with GDAL 3.4.2 / PROJ 8.1.1 can be actually reproduced by using only the "NAD27 to WGS 84 (6)" EPSG Helmert transformation,?with an advertized accuracy of 7.0 m, which is the most precise one in the absence of any
grid:<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">$ echo 32.804191 -117.258911 | cs2cs -d 8 NAD27 WGS84 32.80423884??? -117.25976448 0.00000000<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">or explicitly with the pipeline show by <o:p></o:p></p>
<p class="MsoNormal">projinfo -s NAD27 -t WGS84 --bbox -118,32,-117,33 --spatial-test intersects --single-line<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">$ echo 32.804191 -117.258911 0 |? cct -d 8 +proj=pipeline +step
<o:p></o:p></p>
<p class="MsoNormal">+proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad
<o:p></o:p></p>
<p class="MsoNormal">+step +proj=push +v_3 +step +proj=cart +ellps=clrk66 +step +proj=helmert<o:p></o:p></p>
<p class="MsoNormal">+x=-8 +y=159 +z=175 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop<o:p></o:p></p>
<p class="MsoNormal">+v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap<o:p></o:p></p>
<p class="MsoNormal">+order=2,1<o:p></o:p></p>
<p class="MsoNormal">32.80423884?? -117.25976448??? 0.00000000?????????? inf<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<p class="MsoNormal">sfox - This is the result that I get with GDAL 3.4.2 by specifying only the source and target CRS and the coordinates. Am I missing something or is the default behavior in the newer version of GDAL 3.4.2 with Proj 8.1.1 giving me a result
that is based on a less precise method? What am I missing?<o:p></o:p></p>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:10.0pt;font-family:"Tahoma",sans-serif;color:black">Shawn Fox<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:10.0pt;font-family:"Tahoma",sans-serif;color:black">Sr Principal SW Engineer<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:10.0pt;font-family:"Tahoma",sans-serif;color:black">BAE Systems, Inc.<o:p></o:p></span></p>
<p class="MsoNormal" style="text-autospace:none"><span style="font-size:10.0pt;font-family:"Tahoma",sans-serif;color:black">Geospatial eXploitation Products<o:p></o:p></span></p>
<p class="MsoNormal"><b><span style="font-size:8.0pt">T: </span></b><span style="font-size:8.0pt">858-592-5310<b>
</b></span><span style="font-size:8.0pt">office phone number | <b>M: </b>858-337-2380 |
<b>E: </b>shawn.fox@baesystems.com <o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:8.0pt">10920 Technology Pl, San Diego, CA
</span><o:p></o:p></p>
</div>
<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 <even.rouault@spatialys.com> <br>
<b>Sent:</b> Monday, September 9, 2024 6:13 PM<br>
<b>To:</b> Fox, Shawn D (US) <shawn.fox@baesystems.us>; gdal-dev@lists.osgeo.org<br>
<b>Subject:</b> Re: [gdal-dev] Upgrading to GDAL 3.4.2 from 1.10.0 shows different outputs from gdaltransform and corresponding C++ API<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<table class="MsoNormalTable" border="1" cellspacing="0" cellpadding="0" width="97%" style="width:97.0%;margin-left:5.4pt;border-collapse:collapse;border:none">
<tbody>
<tr style="height:21.2pt">
<td width="97%" valign="top" style="width:97.0%;border:solid red 1.0pt;padding:0in 5.4pt 0in 5.4pt;height:21.2pt">
<p class="MsoNormal" align="center" style="mso-margin-top-alt:auto;margin-bottom:4.0pt;text-align:center;background:white">
<strong><u><span style="font-size:13.5pt;font-family:"Calibri",sans-serif;color:red">External Email Alert</span></u></strong><o:p></o:p></p>
</td>
</tr>
<tr style="height:21.2pt">
<td width="1440" valign="top" style="width:15.0in;border:solid red 1.0pt;border-top:none;padding:0in 5.4pt 0in 5.4pt;height:21.2pt">
<p class="MsoNormal" align="center" style="mso-margin-top-alt:3.0pt;margin-right:0in;margin-bottom:4.0pt;margin-left:0in;text-align:center;background:white">
<strong><span style="font-size:10.0pt;font-family:"Calibri",sans-serif;color:black">This email has been sent from an account outside of the BAE Systems network.</span></strong><o:p></o:p></p>
<p class="MsoNormal" align="center" style="mso-margin-top-alt:auto;margin-bottom:4.0pt;text-align:center">
<span style="font-size:7.5pt">Please treat the email with caution, especially if you are requested to click on a link, decrypt/open an attachment, or enable macros. For further information on how to spot phishing, access “Cybersecurity OneSpace Page” and report
phishing by clicking the button “Report Phishing” on the Outlook toolbar.</span><o:p></o:p></p>
</td>
</tr>
</tbody>
</table>
<p class="MsoNormal"><o:p> </o:p></p>
<p><o:p> </o:p></p>
<div>
<p class="MsoNormal">Le 10/09/2024 à 02:55, Fox, Shawn D (US) via gdal-dev a écrit :<o:p></o:p></p>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<pre>Evan,<o:p></o:p></pre>
<pre><o:p> </o:p></pre>
<pre>Thanks for putting some time into that response and showing some examples. Perhaps I should clarify exactly how our software is using GDAL. I thought that it was easier to show the problem using gdaltransform rather than source code, but that does lack some context. Here is crude, but simple example of how the C Apis are being used.<o:p></o:p></pre>
</blockquote>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<pre><o:p> </o:p></pre>
<pre>--- source code begin ---<o:p></o:p></pre>
<pre>OGRSpatialReferenceH sourceSRS = OSRNewSpatialReference(NULL);<o:p></o:p></pre>
<pre>OGRSpatialReferenceH targetSRS = OSRNewSpatialReference(NULL);<o:p></o:p></pre>
<pre><o:p> </o:p></pre>
<pre>errorSource = OSRSetWellKnownGeogCS( sourceSRS, "NAD27");<o:p></o:p></pre>
<pre>errorTarget = OSRSetWellKnownGeogCS( targetSRS, "WGS84");<o:p></o:p></pre>
<pre><o:p> </o:p></pre>
<pre>OGRCoordinateTransformationH ptrCT = OCTNewCoordinateTransformation( sourceSRS, targetSRS );<o:p></o:p></pre>
<pre>status = OCTTransform( ptrCT, 1, x_lon, y_lat, NULL );<o:p></o:p></pre>
</blockquote>
<p>Source code speaks for itself :-) This is a CRS axis order issue. Cf <a href="https://gdal.org/en/latest/tutorials/osr_api_tut.html#crs-and-axis-order">
https://gdal.org/en/latest/tutorials/osr_api_tut.html#crs-and-axis-order</a> for the whole discussion<o:p></o:p></p>
<p>With your current code and GDAL >= 3, the input axis order should be (lat, lon) since NAD27 in EPSG is declared as (lat, lon) ordered and GDAL 3 honours that by default<o:p></o:p></p>
<p>If you want to keep long, lat axis order, add OSRSetAxisMappingStrategy(errorSource, OAMS_TRADITIONAL_GIS_ORDER); and OSRSetAxisMappingStrategy(errorTarget, OAMS_TRADITIONAL_GIS_ORDER); before creating the coordinate transformation<o:p></o:p></p>
<p>Even<o:p></o:p></p>
<pre>-- <o:p></o:p></pre>
<pre><a href="http://www.spatialys.com">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>
</div>
</body>
</html>