[gdal-dev] Coordinate transformation failing

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

[gdal-dev] Coordinate transformation failing

Stefan Moebius-2
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


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


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

smime.p7s (3K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate transformation failing

Even Rouault-3
Stefan,

Le ven. 22 mars 2019 à 10:44, Stefan Moebius <[hidden email]> a écrit :
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]],

--> OK, the TOWGS84 is wrong here. The rotation (4th,5th,6th) and scale factor difference (7th) terms of the TOWGS84 are obviously nonsense.
They should be respectively in arc-seconds and in PPM.
From what you mention, I believe this is due to
which is likely correct by itself.
The issue here is that as you mentionned the file was produced by GDAL 1.11.1 which incorrectly stored the values
of the rotation and sale factor difference as arc-seconds and PPM whereas Erdas expected radians and unity-based scale factor difference.
Thus when reading this with GDAL >=2.2, the file is wrongly interpretated.
You should override the TOWGS84 values of this WKT with the corrected ones.

Even


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Coordinate transformation failing

Stefan Moebius-2

Thanks, Even.

Is there a way to automatically detect and solve this when reading a file? Detection is obviously possible by checking for the failing transformation, but how about the solution? It seems that the towgs84 parameters aren’t strictly necessary – is simply removing them a reasonable (and generic) workaround?

 

--

Stefan

 

From: Even Rouault <[hidden email]>
Sent: Friday, March 22, 2019 12:05
To: Stefan Moebius <[hidden email]>
Cc: [hidden email]
Subject: Re: [gdal-dev] Coordinate transformation failing

 

Stefan,

 

Le ven. 22 mars 2019 à 10:44, Stefan Moebius <[hidden email]> a écrit :

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]],

 

--> OK, the TOWGS84 is wrong here. The rotation (4th,5th,6th) and scale factor difference (7th) terms of the TOWGS84 are obviously nonsense.

They should be respectively in arc-seconds and in PPM.

From what you mention, I believe this is due to

which is likely correct by itself.

The issue here is that as you mentionned the file was produced by GDAL 1.11.1 which incorrectly stored the values

of the rotation and sale factor difference as arc-seconds and PPM whereas Erdas expected radians and unity-based scale factor difference.

Thus when reading this with GDAL >=2.2, the file is wrongly interpretated.

You should override the TOWGS84 values of this WKT with the corrected ones.

 

Even

 


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


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev

smime.p7s (3K) Download Attachment