search

OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.8k stars 2.51k forks source link

Different results on different computers in GDAL 3.7.0 #8080

Closed andreas-wagner-ptv closed 1 year ago

andreas-wagner-ptv commented 1 year ago

Expected behavior and actual behavior.

GDAL 3.6.4 produced the same results on all our test systems.

With the update to GDAL 3.7.0 on some systems the results were differing. On these systems with differing results, the results were differing in the same way each run. The deviation was in the 4. -5. significant digit. The deviations were reproducible. Proj was used in Version 9.2.0 both with the test with GDAL 3.6.4 and GDAL 3.7.0.

CPUs on systems with result differing in GDAL 3.7.0 compared to GDAL 3.6.4:

CPUs on systems with same results in GDAL 3.7.0 and GDAL 3.6.4:

Operating systems on systems result differing in GDAL 3.7.0 compared to GDAL 3.6.4:

Operating systems on systems with same results in GDAL 3.7.0 and GDAL 3.6.4:

I expect the result to be the same regardless of the system it is run on.

Steps to reproduce the problem.

Unfortunately, I can't provide an isolated test for this problem and I can't share our code that uses GDAL.

Operating system

Windows

GDAL version and provenance

GDAL 3.7.0 and Proj 9.2.0 obtained via vcpkg

rouault commented 1 year ago

Which part of GDAL do you use exactly ? Without a reproducer or very detailed indications, there is no way this ticket is actionable. Or you'll need to run a "git bisect" session to figure out which commit causes the different in results. I'd note that for some parts of GDAL, small differences in numeric output might be expected.

andreas-wagner-ptv commented 1 year ago

In the test with the deviations we use 'OGRCoordinateTransformation' to do three consecutive coordinate transformations and calculate the direct distance between two points.

To calculate the distances we transform to WGS-84 and use GeographicLib::Geodesic::WGS84.Inverse.

from_point: longitude: 8.349849881622 lattitude: 49.0053751692096 to_point: longitude: 8.34508155270866 lattitude: 49.0015850416239

As distance we calculated 0.547157161km and that stayed the same after each transformation with GDAL 1.6.4. With GDAL 1.7.0 we also get this up to the coordinate system "ETRS_1989_LCC" but after the last transformation from "ETRS_1989_LCC" to "DHDN_3_Degree_Gauss_Zone_3_TOWGS" we get 0.547092630km on some systems instead.

We use the following projection definitions: 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]'

'PROJCS["World_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",0],UNIT["Meter",1]]'

'PROJCS["ETRS_1989_LCC",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",4000000.0],PARAMETER["False_Northing",2800000.0],PARAMETER["Central_Meridian",10.0],PARAMETER["Standard_Parallel_1",35.0],PARAMETER["Standard_Parallel_2",65.0],PARAMETER["Latitude_Of_Origin",52.0],UNIT["Meter",1.0]]'

'PROJCS["DHDN_3_Degree_Gauss_Zone_3_TOWGS",GEOGCS["GCS_Deutsche_Hauptdreiecksnetz",DATUM["D_Deutsche_Hauptdreiecksnetz",SPHEROID["Bessel_1841",6377397.155,299.1528128],TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,0.0000067]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",3500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]'

rouault commented 1 year ago

Are you using the same binary builds of GDAL and PROJ on all those platforms ? Can you for example try to run the following on the different platforms (on Windows, to avoid quoting issues, you may need to move the WKT definitions in text file and reference those text files as the value of -s_srs and -t_srs)

echo 3883369.85816388 2479594.16053989 | gdaltransform  -s_srs 'PROJCS["ETRS_1989_LCC",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",4000000.0],PARAMETER["False_Northing",2800000.0],PARAMETER["Central_Meridian",10.0],PARAMETER["Standard_Parallel_1",35.0],PARAMETER["Standard_Parallel_2",65.0],PARAMETER["Latitude_Of_Origin",52.0],UNIT["Meter",1.0]]' -t_srs 'PROJCS["DHDN_3_Degree_Gauss_Zone_3_TOWGS",GEOGCS["GCS_Deutsche_Hauptdreiecksnetz",DATUM["D_Deutsche_Hauptdreiecksnetz",SPHEROID["Bessel_1841",6377397.155,299.1528128],TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,0.0000067]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",3500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]'

I get "3452508.00005929 5429988.50422708" on my system

andreas-wagner-ptv commented 1 year ago

The binaries of GDAL, PROJ and our test executable are the same for all the platforms. They are build on a separate build system. We link dynamically against GDAL and PROJ and there is no GDAL and PROJ installed globally.

On the systems with no deviations, I get the same result as on your system ("3452508.00005929 5429988.50422708 0"). On systems with a deviation, I get ("3452438.40012114 5429874.11158339 0").

rouault commented 1 year ago

@andreas-wagner-ptv

Can you paste the output of

set PROJ_DEBUG=3  # or whatever way of setting environment variables on your system
echo 3883369.85816388 2479594.16053989 | gdaltransform  -s_srs 'PROJCS["ETRS_1989_LCC",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",4000000.0],PARAMETER["False_Northing",2800000.0],PARAMETER["Central_Meridian",10.0],PARAMETER["Standard_Parallel_1",35.0],PARAMETER["Standard_Parallel_2",65.0],PARAMETER["Latitude_Of_Origin",52.0],UNIT["Meter",1.0]]' -t_srs 'PROJCS["DHDN_3_Degree_Gauss_Zone_3_TOWGS",GEOGCS["GCS_Deutsche_Hauptdreiecksnetz",DATUM["D_Deutsche_Hauptdreiecksnetz",SPHEROID["Bessel_1841",6377397.155,299.1528128],TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,0.0000067]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False_Easting",3500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]' --debug on

and

echo  3883369.85816388 2479594.16053989 | gdaltransform -ct +proj=pipeline +step +inv +proj=lcc +lat_0=52 +lon_0=10 +lat_1=35 +lat_2=65 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +step +proj=push +v_3 +step +proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=598.1 +y=73.7 +z=418.2 +rx=0.202 +ry=0.045 +rz=-2.455 +s=6.7e-06 +convention=position_vector +step +inv +proj=cart +ellps=bessel +step +proj=pop +v_3 +step +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel
andreas-wagner-ptv commented 1 year ago

Thank you for your quick response.

The coordinates deviate if proj.db cannot be found. Setting the environment variable PROJ_DATA gives the correct result. In our code, we do not want to rely on this environment variable to be set, but want to set this folder programmatically. Currently we call the the Proj-API function proj_context_set_search_paths for the default context at program startup. My guess is that gdal 3.7.0 might create additional proj-contexts that then don't know where the proj.db is located. What is the correct way to set the proj dir for use in gdal programatically? I had a look at CPLSetConfigOption() and https://gdal.org/user/configoptions.html but did not find anything. In the gdal-source code, I found a change between 3.6.4 and 3.7.0 in ogr/ogr_proj_p.cpp that refers to the config option "PROJ_DATA". Do I have to set PROJ_DATA via CPLSetConfigOption()?

Edit: I guess, I can also use OSRSetPROJSearchPaths() instead that is also available in gdal 3.6.4

rouault commented 1 year ago

I can also use OSRSetPROJSearchPaths() instead that is also available in gdal 3.6.4

yes this is the recommended way. Starting with GDAL 3.7, CPLSetConfigOption("PROJ_DATA", ...) should also do a similar job since there's a specific behaviour for catching the setting of this config option and configure PROJ appropriately.

Currently we call the the Proj-API function proj_context_set_search_paths for the default context at program startup.

That's supposed to work as well.

andreas-wagner-ptv commented 1 year ago

I just verified that calling OSRSetPROJSearchPaths(...) solved our problem. Just calling proj_context_set_search_paths(NULL, ...) was not enough on some sytems. Thank you for your help.

andreas-wagner-ptv commented 1 year ago

I just found out that on the test systems with deviations, the environment variable PROJ_LIB was set.