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
smime.p7s
Description: S/MIME cryptographic signature
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>
_______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org https://lists.osgeo.org/mailman/listinfo/gdal-dev
