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....@baesystems.com
10920 Technology Pl, San Diego, CA
From: Even Rouault <even.roua...@spatialys.com>
Sent: Monday, September 9, 2024 6:13 PM
To: Fox, Shawn D (US) <shawn....@baesystems.us>; gdal-dev@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.
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
