3

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.

1 Answer 1

0

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]);

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.