#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)
Change History (17)
by , 8 years ago
| Attachment: | GDAL_Bug.ipynb added |
|---|
comment:1 by , 8 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 , 8 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 , 8 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 , 8 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 , 8 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 , 8 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 they are diverging by 30+ orders of magnitude.
comment:7 by , 8 years ago
| Priority: | normal → high |
|---|
comment:8 by , 8 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 , 8 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 , 8 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 , 8 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 , 8 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 , 8 years ago
| Summary: | Bilinear, Cubic, Cubic Spline, and Lanczos do not work on ReprojectImage() python bindings → OpenCL 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 , 8 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)
iPythonNB