short integer overflow with r.in.gdal

[dylan]

Noticed that importing Arc ASCII grid files with r.in.gdal can lead to data corruption when the input file contains integer values larger than what can be stored in a short integer (>32767). See the test case and attached file below.

test case:

# setup in spearfish60 location
g.region res=1000 -pa

# make a map with values larger than what a short integer can hold
r.mapcalc "large_int_map = rand(10000, 78000)"

# check raster info:
r.info large_int_map -r
min=10000
max=77594

# export to Arc ASCII format
r.out.arc in=large_int_map out=big_int_arc_file.asc



# import with GDAL:
r.in.gdal -o in=big_int_arc_file.asc out=big_int_arc_file

# we are stuck with short integer values and data corruption:
r.info big_int_arc_file -r
min=-31293
max=32420

It looks like the culprit may be associated with the following lines from the r.out.gdal source (main.c:677 | revision 31323):

/* -------------------------------------------------------------------- */
/*      Select a cell type for the new cell.                            */
/* -------------------------------------------------------------------- */
    eRawGDT = GDALGetRasterDataType( hBand );
    
    switch(eRawGDT) {
      case GDT_Float32:
      case GDT_Float64:
        data_type = FCELL_TYPE;
        eGDT = GDT_Float32;
        complex = FALSE;
        break;

      case GDT_Byte:
        data_type = CELL_TYPE;
        eGDT = GDT_Int32;
        complex = FALSE;
        G_set_cell_format(0);
        /* raster_open_new_func = G_open_raster_new_uncompressed;*/ /* ?? */
        break;

      case GDT_Int16:
      case GDT_UInt16:
        data_type = CELL_TYPE;
        eGDT = GDT_Int32;
        complex = FALSE;
        G_set_cell_format(1);
        /* raster_open_new_func = G_open_raster_new_uncompressed;*/ /* ?? */
        break;

      default:
        data_type = CELL_TYPE;
        eGDT = GDT_Int32;
        complex = FALSE;
	G_set_cell_format(3);
        break;
    }

[dylan]

example arc info ascii grid

[dylan]

forgot to mention: one possible work around is to only import arc info ascii grid with r.in.arc -- which does not have this problem.

[dylan]

This appears to be a problem with GDAL not recognizing the ASCII grid data type, as GDAL-specific operations such as gdal_translate produce erroneous results. GDAL is interpreting the attached input file as Int16.

Using the attached file:

gdal_translate  big_int_arc_file.asc test.tiff

the resulting file is produced:

gdalinfo -stats test.tiff 
Driver: GTiff/GeoTIFF
Files: test.tiff
Size is 7, 7
Coordinate System is `'
Origin = (592000.000000000000000,4926000.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  592000.000, 4926000.000) 
Lower Left  (  592000.000, 4919000.000) 
Upper Right (  599000.000, 4926000.000) 
Lower Right (  599000.000, 4919000.000) 
Center      (  595500.000, 4922500.000) 
Band 1 Block=7x7 Type=Int16, ColorInterp=Gray
  Minimum=-31293.000, Maximum=32420.000, Mean=216.837, StdDev=18087.589
  NoData Value=-9999
  Metadata:
    STATISTICS_MINIMUM=-31293
    STATISTICS_MAXIMUM=32420
    STATISTICS_MEAN=216.83673469388
    STATISTICS_STDDEV=18087.589140406

[dylan]

Since this is a GDAL bug, I will change the type to "task", and make a new ticket on the GDAL tracker. The fix for now is to use r.in.arc for any AA Grid files which contain large (32bit) integer data.

[dylan]

Replying to dylan:

Since this is a GDAL bug, I will change the type to "task", and make a new ticket on the GDAL tracker. The fix for now is to use r.in.arc for any AA Grid files which contain large (32bit) integer data.

​GDAL ticket created.

[neteler]

The related GDAL ticket was closed. Can we close here, too?

Markus

[dylan]

Right. Closed as the source bug was identified and fixed in gdal_trunk.