Our Java application uses GDAL (v3.2.0) and I am trying to convert some coordinates from EPSG:31467 (DHDN GK3) to EPSG:25832 (ETRS89 UTM32). For this I downloaded the grid shift file BETA2007.gsb and tested it using console version of PROJ (v.8.2.0) like that:
echo 3399371.190396 5930724.531323 |\
cs2cs -f "%.6f"\
+proj=tmerc +ellps=bessel +lat_0=0 +lon_0=9 +k_0=1 +x_0=3500000 +y_0=0 +units=m +nadgrids=BETA2007.gsb +to\
+proj=utm +ellps=GRS80 +zone=32 +nadgrids=@null
I get 399340.601862 5928794.177992 0.000000
as a result, which is correct. I have a CSV with test coordinates and transformation results for validation.
This works fine until I try to do the same using GDAL. I have the following Java code which always prints Infinity Infinity Infinity
:
SpatialReference dhdnGk3 = new SpatialReference();
dhdnGk3.SetProjCS("GK3 (DHDN)");
dhdnGk3.SetWellKnownGeogCS("EPSG:4314");
dhdnGk3.SetTM(0, 9, 1, 3500000, 0);
SpatialReference etrs89Utm32 = new SpatialReference();
etrs89Utm32.SetProjCS("UTM32N (ETRS89)");
etrs89Utm32.SetWellKnownGeogCS("EPSG:4258");
etrs89Utm32.SetUTM(32, 1);
CoordinateTransformationOptions options;
CoordinateTransformation transformation;
double[] res;
options = new CoordinateTransformationOptions();
options.SetOperation("" +
"+proj=tmerc +ellps=bessel +lat_0=0 +lon_0=9 +k_0=1 +x_0=3500000 +y_0=0 +units=m +nadgrids=./src/test/resources/BETA2007.gsb " +
"+to " +
"+proj=utm +ellps=GRS80 +zone=32 +nadgrids=@null"
);
transformation = new CoordinateTransformation(dhdnGk3, etrs89Utm32, options);
res = transformation.TransformPoint(3399371.190396, 5930724.531323, 0);
System.out.printf("X, Y, Z: %s, %s, %s%n", res[0], res[1], res[2]);
It starts working as soon as I drop the nadgrids
parameter. The result is however different from the reference value (which is to be expected since no grid shift file is used).
I tried transforming geographic coordinates instead of projected ones using the following PROJ String
echo 7.483333333333 53.500000000000 |\
cs2cs -f "%.12f"\
+proj=latlong +ellps=bessel +nadgrids=BETA2007.gsb +to\
+proj=latlong +ellps=GRS80 +nadgrids=@null
The story is the same: it works perfectly fine in the console but returns Infinity
s when I use it with GDAL.
Then I tried the whole thing with another grid shift file (CHENYX06_etrs.gsb) and EPSG:21781 (CH1903 LV03) as the source CRS but got the same result.
My idea was that in Java example it doesn't find the .gsb-file but this is not the case. I tried to give him a wrong path and he returned an error which implies that the file is found successfully.
Best Answer
After some time I found the answer. The problem was that I was defining SpacialReferences (strictly speaking, their geographic CS part) using EPSG codes but their definitions in the EPSG database doesn't contain any information about the grid shift files. I removed CoordinateTransformationOptions and rewrote the CRS definitons using the PROJ-Strings.