AutoCreateWarpedVRT function creates incorrect VRT Dataset

[georges]

I noticed different behavior between gdal.AutoCreateWarpedVRT function from gdal_csharp.dll and gdalwarp command ran from command prompt. VRT dataset created with gdal.AutoCreateWarpedVRT function definitely outputs incorrect values, while those from gdalwarp are correct.

I am using Rhino5 application to generate a 3d mesh of the SRTM .tif topography model of Vesuvius. The .tif file is originally in WGS84 spatial reference system. What i do is create a new .tif file with projected reference system which uses Azimuthal Equidistant projection.

a) If I use the following gdalwarp command:

C:/gdalwin32-1.4.1/bin/gdalwarp.exe -s_srs EPSG:4326 -t_srs "+proj=aeqd +lat_0=40.81266 +lon_0=14.414252 +datum=WGS84 +units=m +no_defs" C:/vesuvius_wgs84.tif C:/vesuvius_aeqd.tif

And then open the new "vesuvius_aeqd.tif" file in Rhino5 application, everything is correct, and ​3d mesh representation of the .tif file looks like this.

b) However when I try to do the same thing, but with gdal C# libraries by calling the following functions through Rhino5 ironpython interpreter:

...
sourceSRS = osrc.SpatialReference("")
sourceSRS.ImportFromEPSG(4326)
num, datasetsourceSRS = inputCRS.ExportToWkt()

targetSRS = osrc.SpatialReference("")
targetSRS.ImportFromProj4("+proj=aeqd +lat_0=40.81266 +lon_0=14.414252 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
num, datasetOutputSRS = outputSRS.ExportToWkt()
datasetWarped = gdal.AutoCreateWarpedVRT(dataset, datasetsourceSRS, datasetOutputSRS, gdal.ResampleAlg.GRA_NearestNeighbour, 0)
...

​this is what I get.

I also tried to replace the ImportFromProj4 function with SetFromUserInput function, and I still get the upper incorrect result. Looks like for some reason the AutoCreateWarpedVRT function creates an incorrect VRT dataset. Just to mention that if I use EPSG codes, then the newly created VRT dataset is valid, and looks the same as the one generated with gdalwarp command. For example, if I use the EPSG:32633, instead of the upper one which uses Azimuthal Equidistant projection:

...
targetSRS.ImportFromEPSG(32633)
...

in that case, the AutoCreateWarpedVRT function creates a correct VRT dataset, the same one that gdalwarp would create with:

C:/gdalwin32-1.4.1/bin/gdalwarp.exe -s_srs EPSG:4326 -t_srs EPSG:32633 C:/vesuvius_wgs84.tif C:/vesuvius_32633.tif

So something is definitely wrong with the way AutoCreateWarpedVRT function uses Azimuthal Equidistant projection.

The .tif file with WGS84 SRS can be ​downloaded from here:

Just to mention that the ​40.81266, 14.414252 coordinates which appear in upper PROJ4 string, represent the center of the raster.

To call the gdalwarp from the command prompt, I am using ​gdal 1.4.1. As for the C# gdal libraries, I am using the ones from the release_1500_gdal_1.11.1_mapserver_6.4.1 as these were the last ones which supported Windows XP (I am currently using Windows XP SP3 and dotNET Framework 4.0).

I would welcome any kind of reply on this issue.

Thank you.

[Even Rouault]

The link you provided for vesuvius_wgs84.tif gives a raster whose origin is (138.662916666594271,35.407916666672250). So definitely not in Italy. There must be something wrong...

[georges]

Thank you for the reply. I uploaded a wrong .tif file. I corrected the upper link.

[Even Rouault]

Trying the equivalent with Python bindings (with both 1.11.4 and trunk) gives a checksum of 2450 with gdal.AutoCreateWarpedVRT(), that is the same as the one gotten with gdalwarp, so I'm not sure what might be wrong in your case.

from osgeo import gdal, osr
src_ds = gdal.Open('vesuvius2_wgs84.tif')
targetSRS = osr.SpatialReference()
targetSRS.ImportFromProj4("+proj=aeqd +lat_0=40.81266 +lon_0=14.414252 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
warped= gdal.AutoCreateWarpedVRT(src_ds, src_ds.GetProjectionRef(), targetSRS.ExportToWkt(), gdal.GRA_NearestNeighbour, 0)
print(warped.GetRasterBand(1).Checksum())

[georges]

Thank you for the quick reply once again Rouault. I can only guess that there is a slight difference between the gdal_charp.dll file and gdal python .py files. In the very beginning I tried calling functions from gdal python, but they did not work on ironpython 2.7 (which is what Rhino5 is using), so I turned up using c# ones (and successfully have been using them, up until this issue with Azimuthal Equidistant projection).

Can you check the upper .tif file by opening it in c# please? Thank you.

[Even Rouault]

Would you mind providing a standalone C# code that is the port of the above ? I'm not a C# dev.

[georges]

@rouault: It looks the issue was caused by my mistake in the code. Thank you for your support and time, and please accept my apology. Just this more, if I have any kind of issue in future: what does the ​Band.Checksum function do? Returns the number of elevation values in the raster? Wouldn't it be more convenient to actually sum up all those elevation values (in order to check if the two bands are equal), rather than just counting the number of elevation values?

[Even Rouault]

Checksum() adds values with some modulo logic. It isn't however a strong CRC or signature algorithm so you could potentially have bands with same Checksum() but different content. It is mainly used by the regression test suite to catch changes in results.

[georges]

Thank you for the help once again.