[gdal-dev] Coordinate transformation failing

[Stefan Moebius]

Hi again,

We recently encountered an issue with failing coordinate transformations.
The scenario is that a current (Java-based) software using GDAL 2.4.0/PROJ
5.2.0 reads HFA files generated a long time ago by another software using
GDAL 1.11.1/PROJ 4.8.0. The HFA files happen to come in a DHDN projection
system. The current software then transforms the files into a UTM projection
system and fails.

Investigating this further, we've come up with the following minimal test
showing the issue:

The current version of gdalinfo returns this coordinate system:
PROJCS["DHDN_VfD2_Germany_zone_3",
    GEOGCS["GCS_DHDN",
        DATUM["Deutsche_Hauptdreiecksnetz",
            SPHEROID["Bessel_1841",6377397.16,299.15281],
            TOWGS84[583,68,399.5,-0,-0,578614.3160276701,11300000]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",3500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
The old version of gdalinfo returns the same WKT string, but without the
TOWGS84 attribute.

The following code shows the issue when using GDAL 2.4.0/PROJ 5.2.0, where
the first test produces correct results, and the second one returns
infinity:

public class GdalTest
{
    @Test
    public void withoutTOWGS84()
    {
        final SpatialReference utm32 = new SpatialReference();
        utm32.ImportFromProj4("+proj=utm +zone=32 +datum=WGS84 +units=m
+no_defs");

        final SpatialReference dhdn3 = new SpatialReference();
        dhdn3.ImportFromProj4("+proj=tmerc +lat_0=0 +lon_0=9 +k=1
+x_0=3500000 +y_0=0 +ellps=bessel +units=m +no_defs");

        final CoordinateTransformation trafo = new
CoordinateTransformation(utm32, dhdn3);
        System.out.println(Arrays.toString(trafo.TransformPoint(473600.0,
5552000.0)));  // [3473592.6555052917, 5553652.669486293, 0.0]
    }

    @Test
    public void withTOWGS84()
    {
        final SpatialReference utm32 = new SpatialReference();
        utm32.ImportFromProj4("+proj=utm +zone=32 +datum=WGS84 +units=m
+no_defs");

        final SpatialReference dhdn3 = new SpatialReference();
        dhdn3.ImportFromProj4("+proj=tmerc +lat_0=0 +lon_0=9 +k=1
+x_0=3500000 +y_0=0 +ellps=bessel
+towgs84=583,68,399.5,-0,-0,578614.3160276701,11300000 +units=m +no_defs");

        final CoordinateTransformation trafo = new
CoordinateTransformation(utm32, dhdn3);
        System.out.println(Arrays.toString(trafo.TransformPoint(473600.0,
5552000.0)));  // [Infinity, Infinity, Infinity]
    }
}

Attempting the same things using gdaltransform also shows unexpected
behavior: although the results are not infinite, they are still wrong when
using the towgs84 parameters.

gdaltransform -s_srs "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
-t_srs "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel
+units=m +no_defs"
473600 5552000
3473592.65550529 5553652.66948488 0

gdaltransform -s_srs "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
-t_srs "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +ellps=bessel
+towgs84=583,68,399.5,-0,-0,578614.3160276701,11300000 +units=m +no_defs"
473600 5552000
-5235124.00823006 5717391.73687511 -5306247.65833537

So in summary, we have files produced by an old version of GDAL which we can
no longer use with the current version of GDAL. Are we doing something
wrong? Is this a bug? Is there a work-around?

Regards,
Stefan

Stefan Möbius
Development Manager
Amdocs Open Network

+49 351 652 47617 (internal 2044 47617)

amdocs open network – service acceleration

-------------- next part --------------
A non-text attachment was scrubbed...
Name: smime.p7s
Type: application/pkcs7-signature
Size: 2526 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20190322/f7d64df4/attachment-0001.bin>
-------------- next part --------------
This email and the information contained herein is proprietary and confidential and subject to the Amdocs Email Terms of Service, which you may review at https://www.amdocs.com/about/email-terms-of-service <https://www.amdocs.com/about/email-terms-of-service>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20190322/f7d64df4/attachment-0001.html>