Opened 7 years ago

Closed 7 years ago

Last modified 7 years ago

#6669 closed defect (fixed)

OpenCL Warping: Bilinear, Cubic, Cubic Spline, and Lanczos do not work with Float32 datatype + nodata

Reported by: espg Owned by: warmerdam
Priority: high Milestone:
Component: default Version: 1.11.5
Severity: major Keywords:
Cc:

Description

Here is a short script for showing the issue in an iPython notebook:

%pylab inline

# Fixing figure output for all figures
rcParams['axes.labelsize'] = 12
rcParams['xtick.labelsize'] = 12
rcParams['ytick.labelsize'] = 12
rcParams['legend.fontsize'] = 12
rcParams['figure.dpi'] = 160
rcParams['savefig.dpi'] = 200

from osgeo import gdal
import numpy as np
import osr

def resample_tif(tif_file, pixel_spacing=0.5):
    source = gdal.Open(tif_file)
    source.GetRasterBand(1).SetNoDataValue(-32767)
    wgs84 = osr.SpatialReference() # slopppy 
    wgs84.ImportFromEPSG(3413)     # only for polar stereographic!
    # Get the Geotransform vector
    geo_t = source.GetGeoTransform()
    x_size = source.RasterXSize # Raster xsize
    y_size = source.RasterYSize # Raster ysize
    # Work out the boundaries of the new dataset in the target projection
    ulx, uly = geo_t[0], geo_t[3]
    lrx = geo_t[0] + geo_t[1]*x_size
    lry = geo_t[3] + geo_t[5]*y_size
    # Now, we create an in-memory raster
    mem_drv = gdal.GetDriverByName('MEM')
    # The size of the raster is given the new projection and pixel spacing
    # Using the values we calculated above. Also, setting it to store one band
    # and to use Float32 data type.
    dest = mem_drv.Create('', int((lrx - ulx)/pixel_spacing), \
            int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)

    # Calculate the new geotransform
    new_geo = ( ulx, pixel_spacing, geo_t[2], \
                uly, geo_t[4], -pixel_spacing )
    # Set the geotransform
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(source.GetProjection())
    # Perform the projection/resampling 

    return dest, source

# Change this for your system
input_file = "/Users/grigsbye/Desktop/April_21st_DEMs/IODMS3_20120321_15202432_09868_DEM.tif"

testD, testS = resample_tif(input_file)

gdal.ReprojectImage(testS, testD, None, None, gdal.GRA_Bilinear)

# What it looks like
imshow(testD.GetRasterBand(1).ReadAsArray(),vmax = 20, vmin=0)
colorbar()

# What the source image looks like
imshow(testS.GetRasterBand(1).ReadAsArray(),vmax = 20, vmin=0)
colorbar()

# Things working with Nearest Neighbor
gdal.ReprojectImage(testS, testD, None, None, gdal.GRA_NearestNeighbour)

imshow(testD.GetRasterBand(1).ReadAsArray(),vmax = 20, vmin=0)
colorbar()

Attachments (1)

GDAL_Bug.ipynb (200.2 KB ) - added by espg 7 years ago.
iPythonNB

Download all attachments as: .zip

Change History (17)

by espg, 7 years ago

Attachment: GDAL_Bug.ipynb added

iPythonNB

comment:1 by espg, 7 years ago

Image is too big to attach (6 mb), you can download it here:

http://cires1.colorado.edu/shane/IODMS3_20120321_15202432_09868_DEM.tif

comment:2 by espg, 7 years ago

Note that this also doesn't work in GDAL 2.1, and that it did previously work on GDAL 1.11.3

comment:3 by Even Rouault, 7 years ago

Could you be more specific about you think that does not "work" ? I do get different results in the nearest and bilinear cases.

comment:4 by espg, 7 years ago

on the source dataset, I get the following:

In [4]: testD, testS = resample_tif(input_file)

In [5]: gdal.ReprojectImage(testS, testD, None, None, gdal.GRA_Bilinear)
Out[5]: 0

In [6]: np.max(testS.GetRasterBand(1).ReadAsArray())
Out[6]: 15.672032

In [7]: np.min(testS.GetRasterBand(1).ReadAsArray())
Out[7]: -32767.0

Looking at the output from bilinear interpolation, I get the following:

In [8]: np.max(testD.GetRasterBand(1).ReadAsArray())
Out[8]: 3.4028235e+38

In [9]: np.min(testD.GetRasterBand(1).ReadAsArray())
Out[9]: 0.0

Bilinear interpolation is interpolating to 3.4028235e+38, which is a pretty big problem.

The non-masked data mean is about a half meter:

In [11]: np.mean(testS.GetRasterBand(1).ReadAsArray() > -32767.0)
Out[11]: 0.58823436976198462

Taking from the center of the bilinearly interpolated data scene, I get:

In [12]: testD.GetRasterBand(1).ReadAsArray()[600:605,600:605]
Out[12]: 
array([[  3.40282347e+38,   3.40282347e+38,   3.40282347e+38,
          3.40282347e+38,   3.40282347e+38],
       [  3.40282347e+38,   3.40282347e+38,   3.40282347e+38,
          3.40282347e+38,   3.40282347e+38],
       [  3.40282347e+38,   3.40282347e+38,   3.40282347e+38,
          3.40282347e+38,   3.40282347e+38],
       [  3.40282347e+38,   3.40282347e+38,   3.40282347e+38,
          3.40282347e+38,   3.40282347e+38],
       [  3.40282347e+38,   3.40282347e+38,   3.40282347e+38,
          3.40282347e+38,   3.40282347e+38]], dtype=float32)


…which is clearly garbage.

comment:5 by espg, 7 years ago

By the way, here's the target output when using nearest neighbor instead:

In [13]: gdal.ReprojectImage(testS, testD, None, None)  # Default is nearest
Out[13]: 0

In [14]: testD.GetRasterBand(1).ReadAsArray()[600:605,600:605]
Out[14]: 
array([[ 12.03216171,  12.04551601,  11.97155762,  11.88038635,
         11.90100574],
       [ 12.02534962,  12.05619907,  12.03511333,  11.96802902,
         11.93831539],
       [ 12.12434006,  12.13599014,  12.1388731 ,  12.02012444,
         11.90720558],
       [ 12.06089115,  12.08328438,  12.03147507,  11.96336269,
         11.88910961],
       [ 12.05824184,  12.08551979,  12.01619148,  11.9562254 ,
         11.87326145]], dtype=float32)

comment:6 by espg, 7 years ago

Anymore information needed for this bug? Behavior is correct when calling nearest neighbors or calling 'average' for the kernel window. Bilinear and other methods are expected to be different, but still within an order of magnitude... here the are diverging by 30+ orders of magnitude.

Version 0, edited 7 years ago by espg (next)

comment:7 by espg, 7 years ago

Priority: normalhigh

comment:8 by espg, 7 years ago

Went back and checked pervious versions of GDAL... this is broken on 1.11.3, but works on GDAL 1.9.2.

comment:9 by Even Rouault, 7 years ago

Even with your latest precisions, I cannot reproduce. I've testing the following script (resample_tif() being exactly what you provided):

from osgeo import gdal
import numpy as np
import osr

def resample_tif(tif_file, pixel_spacing=0.5):
    source = gdal.Open(tif_file)
    source.GetRasterBand(1).SetNoDataValue(-32767)
    wgs84 = osr.SpatialReference() # slopppy 
    wgs84.ImportFromEPSG(3413)     # only for polar stereographic!
    # Get the Geotransform vector
    geo_t = source.GetGeoTransform()
    x_size = source.RasterXSize # Raster xsize
    y_size = source.RasterYSize # Raster ysize
    # Work out the boundaries of the new dataset in the target projection
    ulx, uly = geo_t[0], geo_t[3]
    lrx = geo_t[0] + geo_t[1]*x_size
    lry = geo_t[3] + geo_t[5]*y_size
    # Now, we create an in-memory raster
    mem_drv = gdal.GetDriverByName('MEM')
    # The size of the raster is given the new projection and pixel spacing
    # Using the values we calculated above. Also, setting it to store one band
    # and to use Float32 data type.
    dest = mem_drv.Create('', int((lrx - ulx)/pixel_spacing), \
            int((uly - lry)/pixel_spacing), 1, gdal.GDT_Float32)

    # Calculate the new geotransform
    new_geo = ( ulx, pixel_spacing, geo_t[2], \
                uly, geo_t[4], -pixel_spacing )
    # Set the geotransform
    dest.SetGeoTransform(new_geo)
    dest.SetProjection(source.GetProjection())
    # Perform the projection/resampling 

    return dest, source

# Change this for your system
input_file = "IODMS3_20120321_15202432_09868_DEM.tif"

testD, testS = resample_tif(input_file)

gdal.ReprojectImage(testS, testD, None, None, gdal.GRA_Bilinear)

print np.min(testS.GetRasterBand(1).ReadAsArray())

print np.max(testS.GetRasterBand(1).ReadAsArray())

print np.min(testD.GetRasterBand(1).ReadAsArray())

print np.max(testD.GetRasterBand(1).ReadAsArray())

print testD.GetRasterBand(1).ReadAsArray()[600:605,600:605]

With latest trunk and 2.1 branch it returns:

-32767.0
15.672
0.0
15.6537
[[ 12.0340271   12.04712486  12.0134573   11.9291172   11.90681553]
 [ 12.06609249  12.08703804  12.08900166  12.02989483  11.9312067 ]
 [ 12.11827087  12.13494396  12.1448555   12.06794453  11.92060471]
 [ 12.08259392  12.10066986  12.07983112  12.0018177   11.89653778]
 [ 12.06880093  12.08003616  12.03974342  11.97593212  11.88527966]]

And with latest state of 1.11:

-32767.0
15.672
0.0
15.6621
[[ 12.03403282  12.045331    12.00694466  11.91388798  11.90785027]
 [ 12.06357765  12.08479118  12.08582973  12.02124023  11.93400478]
 [ 12.12120438  12.13504219  12.14491177  12.05905056  11.91540623]
 [ 12.07710361  12.09484959  12.0701704   11.98960972  11.89205551]
 [ 12.06372452  12.07659531  12.0347271   11.9680872   11.87340927]]

The differences between 1.11 and later are expected since in 2.0 fixes have been made in various resampling algorithms, but in all cases, the results don't show the huge values you get.

(Those tests have been made on Ubuntu 14.04 64 bit)

Do you get the same issue if you do "gdalwarp IODMS3_20120321_15202432_09868_DEM.tif out.tif -s_srs EPSG:3413 -t_srs EPSG:3413 -srcnodata -32767 -overwrite" ?

Is there anything particular in your build? Hum, do you have opencl enabled ?

comment:10 by espg, 7 years ago

running "gdalwarp IODMS3_20120321_15202432_09868_DEM.tif out.tif -s_srs EPSG:3413 -t_srs EPSG:3413 -srcnodata -32767 -overwrite" gives an image of the same size/dimensions, with valid data.

Thanks for having a look. I do have OpenCL and armadillo enabled (as well as libxml), but I'll recompile without OpenCL and see if that changes things...

comment:11 by Even Rouault, 7 years ago

I forgot the "-r bilinear" in the above gdalwarp command line.

I'm pretty sure this must be due to OpenCL then (OpenCL is used for all resampling methods, except nearest neighbour). I guess it is somehow experimental, especially for non byte data types. Or perhaps this is due to the nodata value. Or perhaps an issue with your particular OpenCL implementation (could you run with CPL_DEBUG=ON defined as environment variable ? This should tell more about it)

You should normally be able to disable OpenCL at runtime by passing options = [ 'USE_OPENCL=NO' ] to gdal.ReprojectImage (untested though)

comment:12 by espg, 7 years ago

Well, recompiling without OpenCL fixes things. Thanks again for the help tracking down the problem!

I had tried testing without nodata values on the other version of GDAL compiled with OpenCL, and the same problem occurred, so the issue in OpenCL is likely with the 32 Float data type (or elsewhere). Is it worth renaming this bug to an OpenCL resampling bug and leaving it an open issue in GDAL, or should I close this bug (i.e., assume it's a problem in the upstream OpenCL)? My guess is that if OpenCL is installed and can't handle 32 floats types, that the expected behavior would be to return the operation to the CPU on non-OpenCL code...

comment:13 by Even Rouault, 7 years ago

Summary: Bilinear, Cubic, Cubic Spline, and Lanczos do not work on ReprojectImage() python bindingsOpenCL Warping: Bilinear, Cubic, Cubic Spline, and Lanczos do not work with Float32 datatype + nodata

I've renamed the ticket. Would you mind trying with opencl again with CPL_DEBUG=ON enabled, so to get more info on your OpenCL runtime ? The issue might be generic to any implementation, but in case this would be specific to your implementation. That said, nobody is currently maintaining the opencl backend so I can't promise this will be fixed any time soon.

comment:14 by espg, 7 years ago

Here's the result:

~/Desktop ➤ gdalwarp IODMS3_20120321_15202432_09868_DEM.tif out.tif -r bilinear -s_srs EPSG:3413 -t_srs EPSG:3413 -srcnodata -32767 -overwrite --config CPL_DEBUG ON

GDAL: GDALOpen(IODMS3_20120321_15202432_09868_DEM.tif, this=0x7fca63427410) succeeds as GTiff.
Creating output file that is 1413P x 1216L.
GDAL: GDALDriver::Create(GTiff,out.tif,1413,1216,1,Float32,0x0)
Processing input file IODMS3_20120321_15202432_09868_DEM.tif.
WARP: Copying metadata from first source to destination dataset
Copying nodata values from source IODMS3_20120321_15202432_09868_DEM.tif to destination out.tif.
WARP: srcNoData=-32767.000000 dstNoData=-32767.000000
WARP: calling GDALSetRasterNoDataValue() for band#0
OpenCL: Connected to NVIDIA GeForce GTX 780M.
GDAL: GDALWarpKernel()::GWKOpenCLCase()
Src=0,0,1413x1216 Dst=0,0,1413x1216
0...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: GDALClose(IODMS3_20120321_15202432_09868_DEM.tif, this=0x7fca63427410)
GDAL: GDALClose(out.tif, this=0x7fca63428660)

I had tried on a subset that didn't have any nodata values, and the same problem occurred there (bad output)... similar debug trace:

~/Desktop ➤ gdalwarp test_in.tif out2.tif -r bilinear -s_srs EPSG:3413 -t_srs EPSG:3413 -overwrite --config CPL_DEBUG ON

GDAL: GDALOpen(test_in.tif, this=0x7fc73bd168a0) succeeds as GTiff.
Creating output file that is 100P x 100L.
GDAL: GDALDriver::Create(GTiff,out2.tif,100,100,1,Float32,0x0)
Processing input file test_in.tif.
WARP: Copying metadata from first source to destination dataset
Using internal nodata values (e.g. 0) for image test_in.tif.
Copying nodata values from source test_in.tif to destination out2.tif.
WARP: srcNoData=0.000000 dstNoData=0.000000
WARP: calling GDALSetRasterNoDataValue() for band#0
OpenCL: Connected to NVIDIA GeForce GTX 780M.
GDAL: GDALWarpKernel()::GWKOpenCLCase()
Src=0,0,100x100 Dst=0,0,100x100
0...10...20...30...40...50...60...70...80...90...100 - done.
GDAL: GDALClose(test_in.tif, this=0x7fc73bd168a0)
GDAL: GDALClose(out2.tif, this=0x7fc73bc102a0)

comment:15 by Even Rouault, 7 years ago

Resolution: fixed
Status: newclosed

In 36104:

OpenCL wrapper: fix wrapping of Float32 data type with NVIDIA OpenCL (fixes #6669)

comment:16 by Even Rouault, 7 years ago

In 36105:

OpenCL wrapper: fix wrapping of Float32 data type with NVIDIA OpenCL (fixes #6669)

Note: See TracTickets for help on using tickets.