[gdal-dev] problem with interpolation on GDAL RasterIO ?

[Andrea Battisti]

Hi all,


I think I found a glitch in how bilinear interpolation is applied when 
reading a sub block of a dataset.

I have put together few lines of code to help showing the issue (see 
below) in case any one can help me on this. (tested on GDAL 2.4.2).

The region being read has values in ~ -2800..-2300 range, and there is a 
no_data block with value -340282265508890445205022487695511781376.

Now if I do not apply any interpolation, I get the expected range of 
values back (not masking nodata here for simplicity), eg:

computed min, max: -340282265508890445205022487695511781376.000000, 
-2360.786377

While if I enable _bilinear_ interpolation, I get:
computed min, max: -340282265508890445205022487695511781376.000000, 0.000000

A whole block within the nodata area is set to 0 while should just stay 
as nodata (eg -340282265508890445205022487695511781376) and the maximum 
should stay around -2360.


To add some fun, I tried to extract the same block using 
'gdal_translate' utility, both with and without interpolation, and in 
that case I get consistent results:

gdal_translate -srcwin 0 9880 2525 4442  -outsize 704 400 
NAC_DTM_APOLLO17_3.TIF -r near nsub_n.tif
gdalinfo -mm nsub_n.tif
--> Computed Min/Max=-2867.847,-2360.786

gdal_translate -srcwin 0 9880 2525 4442  -outsize 704 400 
NAC_DTM_APOLLO17_3.TIF -r bilinear nsub_bl.tif
gdalinfo -mm nsub_bl.tif
--> Computed Min/Max=-2867.411,-2360.956

I am not familiar enough with the GDAL internals to dig more but I guess 
gdal_translate might add more care for edge cases.


The dataset that is referenced is an elevation model and can be 
downloaded from:
http://lroc.sese.asu.edu/data/LRO-L-LROC-5-RDR-V1.0/LROLRC_2001/DATA/SDP/NAC_DTM/APOLLO17_3/NAC_DTM_APOLLO17_3.TIF


Thanks in advance!
Andrea

==================================

#define SNC 2525
#define SNR 4442
#define XOFF 0
#define YOFF 9880

int main() {

         GDALRasterIOExtraArg psExtraArg;
         GDALDatasetH hDataset;
         GDALDriverH hDriver;
         CPLErr err;
         double *buf, bmin, bmax;
         void *bp;
         unsigned int i;

         printf( "GDAL rasterio test\n"  );

         GDALAllRegister();
         hDataset = GDALOpenShared( "NAC_DTM_APOLLO17_3.TIF", GA_ReadOnly );
         hDriver = GDALGetDatasetDriver( hDataset );

         printf( "Driver: %s - %s\n", GDALGetDriverShortName( hDriver ), 
GDALGetDriverLongName( hDriver )  );
         printf( "b,r,c: %d, %d, %d\n", GDALGetRasterCount( hDataset ), 
GDALGetRasterYSize( hDataset ), GDALGetRasterXSize( hDataset ) );

         // read sub region
         buf = malloc( BNC * BNR * sizeof(double) );
         bp = (void *)buf;


         // use interpolation while reading
         INIT_RASTERIO_EXTRA_ARG(psExtraArg);
         psExtraArg.eResampleAlg = GRIORA_NearestNeighbour; // 
GRIORA_NearestNeighbour | GRIORA_Bilinear

         err = GDALDatasetRasterIOEx(
                 hDataset,
                 GF_Read, // GDALRWFlag eRWFlag,
                 XOFF, // int nDSXOff,
                 YOFF, // int nDSYOff,
                 SNC, // int nDSXSize,
                 SNR, // int nDSYSize,
                 bp, // void *pBuffer,
                 BNC, // int nBXSize,
                 BNR, // int nBYSize,
                 GDT_Float64, // GDALDataType eBDataType,
                 1, // int nBandCount,
                 NULL, // int *panBandCount,
                 0, // GSpacing nPixelSpace,
                 0, // GSpacing nLineSpace,
                 0, // GSpacing nBandSpace,
                 &psExtraArg // GDALRasterIOExtraArg *psExtraArg
         );

         bmin = 1e50;
         bmax = -bmin;
         for (i=0; i < (BNC*BNR); ++i) {
                 if (buf[i] > bmax) bmax = buf[i];
                 if (buf[i] < bmin) bmin = buf[i];
         }
         printf( "computed min, max: %f, %f\n", bmin, bmax );

         GDALClose(hDataset);

         return 0;
}


==================================



-- 
Andrea Battisti
Computer Science Engineer