[gdal-dev] Upgrading to GDAL 3.4.2 from 1.10.0 shows different outputs from gdaltransform and corresponding C++ API
Fox, Shawn D (US)
shawn.fox at baesystems.us
Mon Sep 9 17:55:32 PDT 2024
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 end ---
I had to search the web for some of the strings such as " NAD27 to WGS 84 (41)" to understand what that meant. Evidently there are EPSG codes for the various sets of steps for transformations which I did not know.
https://epsg.io/8594 (NAD27 to WGS 84 (41))
More importantly, I still don't understand why the behavior of GDAL would have changed for the exact same inputs, regardless of the advertised precision. I suppose that perhaps the transformation chosen by GDAL for that operation changed at some point. I have read in the tutorial at https://gdal.org/en/latest/tutorials/osr_api_tut.html that GDAL picks the most appropriate steps based on the coordinates, but that one can specify a proj pipeline to make it pick a specific set of steps. Do you know how to print the proj steps that GDAL will use when transforming? I did not see an obvious way to do that within the class API declarations of OGRCoordinateTransformationH.
Given the limits of this type of transformation, perhaps the differences I am seeing are a non-issue but I still need to be able to explain this to others why are regression tests fail. At a minimum, I am trying to understand why the results are different. For now I will try to use the proj pipelies to test the C++ Apis with various inputs to see if I can reproduce the results observed with previous versions of GDAL. I've also tested this on multiple platforms, windows 10 and Red Hat 8 so I don't believe that this is due to floating point limitations of CPUs because I observe the exact same results with GDAL 3.1.0 on windows 10 as I observe with GDAL 3.4.2 on Red Hat 8. Therefore, I don't believe I can simply blame calculations with floating point for the differences.
Thanks,
Shawn Fox
-------------------------------------------------------------------
Date: Mon, 26 Aug 2024 18:26:47 +0200
From: Even Rouault <even.rouault at spatialys.com>
To: "Fox, Shawn D (US)" <shawn.fox at baesystems.us>,
"'gdal-dev at lists.osgeo.org'" <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
Message-ID: <e2bbdd85-c464-4fb6-bfb8-19bcc3676b51 at spatialys.com>
Content-Type: text/plain; charset="utf-8"; Format="flowed"
Shawn,
The results depends on which PROJ grids are available on your system (or if you allowed network access).
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
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 NAD27 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
3) You would get the """"most precise result"""" (with lots of quotes, because this is highly subjective, and assumes that NAD83(HARN) is equivalent to WGS 84, according to the comment in the EPSG database) with the "NAD27 to WGS 84 (41)" EPSG transformation, with an advertized accuracy of 1.5m, which combines the CONUS grid and us_noaa_cshpgn:
$ echo 32.804191 -117.258911 | PROJ_NETWORK=ON cs2cs -d 8 NAD27 WGS84 32.80424065??? -117.25978105 0.00000000
or explicitly with the pipeline show by "projinfo NAD27 WGS84 --bbox
-118,32,-117,33 --spatial-test intersects --single-line":
$ 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=hgridshift +grids=us_noaa_cshpgn.tif +step +proj=unitconvert
+xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1
32.80424065?? -117.25978105??? 0.00000000?????????? inf
==> All those results are valid. Some are more precise than others, under certain assumptions... But basically speaking about WGS 84 introduces an intrinsic inaccuracy of about 2 meters (cf https://www.youtube.com/watch?v=M2ck3cAGvhg). I don't think there are any published time-dependent transformations for NAD27 to WGS84 or any of its realization.
Even
Le 24/08/2024 ? 03:14, Fox, Shawn D (US) via gdal-dev a ?crit?:
>
> Hello GDAL users.
>
> I have noticed that from GDAL 1.10 ? GDAL 3.1.0 that gdaltransform is
> producing different results for NAD27 GCS to WGS84 GCS
> transformations.? The differences in longitude in the output is more
> significant.? Having tested the differences in feet between the points
> output from each version I can see a difference of ~5 feet primarily
> due to the change in the longitude value.
>
> Would anyone look at the following results and share your thoughts
> about why the result of this transformation is different in a newer
> version?? Could this be related to the concept of WGS 84 realization
> being different in newer versions of GDAL?? I am also linking GDAL
> 3.4.2 with PROJ 8.1.1 at this time.? These results were produced with
> GDAL 3.4.2, but the same results occur with GDAL 3.1.0 as well.? GDAL
> 3.4.2 is what I am currently trying to upgrade to.? The results from
> the C++ programming API are similar in difference.? I just used
> gdaltransform for the purposes of this message because it was a more
> convenient way of producing something that others might be able to
> reproduce.
>
> *GDAL 1.10.0 results*
>
> ./gdaltransform -s_srs NAD27 -t_srs WGS84
>
> -117.258911 32.804191
>
> -117.25978144845 32.8042393819884 1.32247805595398e-07
>
> *GDAL 3.4.2 results*
>
> ./gdaltransform -s_srs NAD27 -t_srs WGS84
>
> -117.258911 32.804191
>
> -117.259764476223 32.8042388409449 0
>
> I have read that the support for specifying epochs related to
> realizations is not available until GDAL 3.8.0.
>
> Here are some references I?ve read so far.
>
> Transforming between WGS84 Realizations | Journal of Surveying
> Engineering | Vol 148, No 2 (ascelibrary.org)
> <https://ascelibrary.org/doi/10.1061/%28ASCE%29SU.1943-5428.0000389>
>
> RFC 81: Support for coordinate epochs in geospatial formats ? GDAL
> documentation
> <https://gdal.org/development/rfc/rfc81_coordinate_epoch.html>
>
> Thanks,
>
> Shawn Fox
>
> San Diego, CA
>
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/gdal-dev
--
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/20240826/6677cc17/attachment.htm>
------------------------------
Subject: Digest Footer
_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/gdal-dev
------------------------------
End of gdal-dev Digest, Vol 243, Issue 20
*****************************************
More information about the gdal-dev
mailing list