My application uses GDAL 3.4.1 and proj 8.2.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()
I've been able to reproduce this with multiple files, but the one I've been testing the most has a gdalinfo like this:
Driver: GTiff/GeoTIFF
Files: mt_rainier_90m.tif
Size is 1201, 1201
Coordinate System is:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-122.000416666666666,47.000416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-122.0004167, 47.0004167) (122d 0' 1.50"W, 47d 0' 1.50"N)
Lower Left (-122.0004167, 45.9995833) (122d 0' 1.50"W, 45d59'58.50"N)
Upper Right (-120.9995833, 47.0004167) (120d59'58.50"W, 47d 0' 1.50"N)
Lower Right (-120.9995833, 45.9995833) (120d59'58.50"W, 45d59'58.50"N)
Center (-121.5000000, 46.5000000) (121d30' 0.00"W, 46d30' 0.00"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Gray
NoData Value=-32768
I've ran it through GDB and I get this stack trace where the app crashes
rogram received signal SIGFPE, Arithmetic exception.
__atan_fma (x=1.0486944322247436) at ../sysdeps/ieee754/dbl-64/s_atan.c:116
116 ../sysdeps/ieee754/dbl-64/s_atan.c: No such file or directory.
(gdb) backtrace
#0 __atan_fma (x=1.0486944322247436) at ../sysdeps/ieee754/dbl-64/s_atan.c:116
#1 0x00007ffff42b912d in ?? () from /lib/x86_64-linux-gnu/libproj.so.22
#2 0x00007ffff4301a8f in ?? () from /lib/x86_64-linux-gnu/libproj.so.22
#3 0x00007ffff430db08 in proj_trans () from /lib/x86_64-linux-gnu/libproj.so.22
#4 0x00007ffff4311311 in ?? () from /lib/x86_64-linux-gnu/libproj.so.22
#5 0x00007ffff4301918 in ?? () from /lib/x86_64-linux-gnu/libproj.so.22
#6 0x00007ffff430db08 in proj_trans () from /lib/x86_64-linux-gnu/libproj.so.22
#7 0x00007ffff6338c72 in OGRProjCT::TransformWithErrorCodes (this=this@entry=0x55555612d8a0, nCount=nCount@entry=84, x=<optimized out>, y=<optimized out>, z=<optimized out>, t=0x0, panErrorCodes=0x55555637c230)
at ogrct.cpp:2465
#8 0x00007ffff633b788 in OGRProjCT::TransformWithErrorCodes (panErrorCodes=<optimized out>, t=<optimized out>, z=<optimized out>, y=<optimized out>, x=<optimized out>, nCount=<optimized out>,
this=0x55555612d8a0) at /usr/include/c++/11/bits/stl_vector.h:363
#9 OGRProjCT::Transform (this=0x55555612d8a0, nCount=84, x=<optimized out>, y=<optimized out>, z=<optimized out>, t=<optimized out>, pabSuccess=0x55555645b990) at ogrct.cpp:2038
#10 0x00007ffff63eb473 in GDALReprojectionTransform (pTransformArg=0x55555611a0a0, bDstToSrc=1, nPointCount=84, padfX=0x55555663b830, padfY=0x55555663bad0, padfZ=0x55555663bd70, panSuccess=0x55555645b990)
at gdaltransformer.cpp:3180
#11 0x00007ffff63e8a7f in GDALGenImgProjTransform (pTransformArgIn=0x555555d9fb50, bDstToSrc=1, nPointCount=84, padfX=0x55555663b830, padfY=0x55555663bad0, padfZ=<optimized out>, panSuccess=0x55555645b990)
at gdaltransformer.cpp:2529
#12 0x00007ffff63e99f9 in GDALApproxTransform (pCBData=0x55555611d800, bDstToSrc=<optimized out>, nPoints=84, x=0x55555663b830, y=0x55555663bad0, z=0x55555663bd70, panSuccess=0x55555645b990)
--Type <RET> for more, q to quit, c to continue without paging--
p:3883
#13 0x00007ffff640db95 in GDALWarpOperation::ComputeSourceWindow (this=this@entry=0x55555618e6b0, nDstXOff=nDstXOff@entry=0, nDstYOff=nDstYOff@entry=896, nDstXSize=nDstXSize@entry=512,
nDstYSize=nDstYSize@entry=128, pnSrcXOff=pnSrcXOff@entry=0x7fffffffcd78, pnSrcYOff=0x7fffffffcd80, pnSrcXSize=0x7fffffffcd88, pnSrcYSize=0x7fffffffcd90, pdfSrcXExtraSize=0x7fffffffcbe8,
pdfSrcYExtraSize=0x7fffffffcbe0, pdfSrcFillRatio=0x0) at gdalwarpoperation.cpp:2754
#14 0x00007ffff640f644 in GDALWarpOperation::WarpRegionToBuffer (this=0x55555618e6b0, nDstXOff=0, nDstYOff=896, nDstXSize=512, nDstYSize=128, pDataBuf=0x555556a24fb0, eBufDataType=GDT_Int16,
nSrcXOff=<optimized out>, nSrcYOff=<optimized out>, nSrcXSize=<optimized out>, nSrcYSize=<optimized out>, dfSrcXExtraSize=<optimized out>, dfSrcYExtraSize=<optimized out>, dfProgressBase=0
I've tried this using GDAL 3.8.4 and Proj 9.3.1 and it happens much less often but it still fails with a Floating Point Exception sometimes.
Does anyone see something I'm obviously doing wrong and have some tips?
