My application uses GDAL 4.3.1 with the Python bindings and I'm using gdal.AutoCreateWarpedVRT to reproject a Geodetic (epsg:4326) to Mercator (epsg:3857) dataset.
Sometimes, but not always, I will get a "Floating point exception" in the and the application will crash in the warp kernel. I've turned on CPL_DEBUG=ON and I see an output like this
GDAL: GDALWarpKernel()::GWKRealCase() Src=0,545,641x113 Dst=0,640,512x128 GDAL: GDALWarpKernel()::GWKRealCase() Src=636,545,565x113 Dst=512,640,451x128 GDAL: GDALWarpKernel()::GWKRealCase() Src=0,655,641x113 Dst=0,768,512x128 GDAL: GDALWarpKernel()::GWKRealCase() Src=636,655,565x113 Dst=512,768,451x128 GDAL: GDALWarpKernel()::GWKRealCase() Src=0,765,641x113 Dst=0,896,512x128 Floating point exception
I've made a simple sample application that reproduces the issue:
from osgeo import gdal, osr
import argparse
def renderTile(filename):
# Open the dataset
src_ds = gdal.Open(filename)
# Warp the dataset to mercator using a WarpedVRT
dst_srs = osr.SpatialReference()
dst_srs.ImportFromEPSG(3857)
dst_wkt = dst_srs.ExportToWkt()
error_threshold = 0.125
resampling = gdal.GRA_Bilinear
warped_dataset = gdal.AutoCreateWarpedVRT(
src_ds,
None,
dst_wkt,
resampling,
error_threshold
)
# Read the first band from the warped dataset
band = warped_dataset.GetRasterBand(1)
data = band.ReadAsArray()
print("Read data from %s length=%s" % (filename, len(data)))
# Close the dataset
band = None
warped_dataset = None
src_ds = None
def main():
args = argparse.ArgumentParser()
args.add_argument("filename", help="filename to render")
args = args.parse_args()
filename = args.filename
for i in range(0, 100):
print("Starting attempt %s" % i)
renderTile(filename)
print("finished %s" % i)
print("done")
if __name__ == "__main__":
main()
This code appears to work correctly with the latest version of GDAL, so I may just have to upgrade, but does anyone see something I'm obviously doing wrong and have some tips?