GDAL Coordinate Transformation – Using NTv2 and CS2CS Successfully in Java Console

gdaljavaproj

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 Infinitys 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.

SpatialReference dhdnGk3 = new SpatialReference();
dhdnGk3.ImportFromProj4(
    "+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");

SpatialReference etrs89Utm32 = new SpatialReference();
etrs89Utm32.ImportFromProj4(
    "+proj=utm +ellps=GRS80 +zone=32 +nadgrids=@null");

CoordinateTransformation transformation = new CoordinateTransformation(dhdnGk3, etrs89Utm32);

double[] res = transformation.TransformPoint(3399371.190396, 5930724.531323);

System.out.printf("East, North: %s, %s%n", res[0], res[1]);
Related Question