[gdal-dev] Upgrading to GDAL 3.4.2 from 1.10.0 shows different outputs from gdaltransform and corresponding C++ API

Tue Sep 10 17:41:31 PDT 2024

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.

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.

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)"
EPSG transformation, with an advertized accuracy of 5.0 m:

$ echo 32.804191 -117.258911 0 | PROJ_NETWORK=ON 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=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
32.80423938?? -117.25978145??? 0.00000000?????????? inf

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.

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?

PROJ_NETWORK=ON
echo  $PROJ_NETWORK
ON

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

proj_create: Error 1029 (File not found or invalid): pipeline: Pipeline: Bad step definition: proj=hgridshift (File not found or invalid)
cct: Bad transformation arguments - (File not found or invalid)
    'cct -h' for help

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:

$ echo 32.804191 -117.258911 | cs2cs -d 8 NAD27 WGS84 32.80423884??? -117.25976448 0.00000000

or explicitly with the pipeline show by
projinfo -s NAD27 -t WGS84 --bbox -118,32,-117,33 --spatial-test intersects --single-line

$ 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=push +v_3 +step +proj=cart +ellps=clrk66 +step +proj=helmert
+x=-8 +y=159 +z=175 +step +inv +proj=cart +ellps=WGS84 +step +proj=pop
+v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap
+order=2,1
32.80423884?? -117.25976448??? 0.00000000?????????? inf

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?

Shawn Fox
Sr Principal SW Engineer
BAE Systems, Inc.
Geospatial eXploitation Products
T: 858-592-5310 office phone number | M: 858-337-2380 | E: shawn.fox at baesystems.com
10920 Technology Pl, San Diego, CA

From: Even Rouault <even.rouault at spatialys.com>
Sent: Monday, September 9, 2024 6:13 PM
To: Fox, Shawn D (US) <shawn.fox at baesystems.us>; gdal-dev at lists.osgeo.org
Subject: Re: [gdal-dev] Upgrading to GDAL 3.4.2 from 1.10.0 shows different outputs from gdaltransform and corresponding C++ API

External Email Alert
This email has been sent from an account outside of the BAE Systems network.
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.

Le 10/09/2024 à 02:55, Fox, Shawn D (US) via gdal-dev a écrit :

Evan,

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.

--- source code begin ---

OGRSpatialReferenceH sourceSRS = OSRNewSpatialReference(NULL);

OGRSpatialReferenceH targetSRS = OSRNewSpatialReference(NULL);

errorSource = OSRSetWellKnownGeogCS( sourceSRS, "NAD27");

errorTarget = OSRSetWellKnownGeogCS( targetSRS, "WGS84");

OGRCoordinateTransformationH ptrCT = OCTNewCoordinateTransformation( sourceSRS, targetSRS );

status = OCTTransform( ptrCT, 1, x_lon, y_lat, NULL );

Source code speaks for itself :-) This is a CRS axis order issue. Cf https://gdal.org/en/latest/tutorials/osr_api_tut.html#crs-and-axis-order for the whole discussion

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

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

Even

--

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/20240911/672066cb/attachment.htm>