Floating Point Exception with geodetic to mercator reprojection

[jasonbeverage]

What is the bug?

My application uses the GDAL python bindings to generate map tiles. I've noticed that occasionally it was throwing a "Floating Point Exception" and crashing when reprojecting Geodetic (epsg:4326) to Mercator (epsg:3857). I'm using an older version of GDAL in my application (3.4.1) but this problem still happens in GDAL 3.8.5 with proj 9.4.0 as I've shown in the Steps to Reproduce section.

I don't think this is anything to do with Python as it seems the crash is deep in proj_trans.

I've ran the application through gdb and the stack trace that I get when it crashes is:

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

Steps to reproduce the issue

warp_issue.zip contains all the data needed to reproduce the issue. The vrt_test.py application opens the given file and loads band 1 into memory using band.ReadAsArray() over and over and again, opening and closing the dataset each time.

The mt_rainier_90m_merc.vrt file was created with the command:

gdalwarp -t_srs epsg:3857 -r bilinear -of VRT mt_rainier_90m.tif mt_rainier_90m_merc.vrt

docker run -it -v .:/code ghcr.io/osgeo/gdal:ubuntu-full-3.8.5 bash
cd /code
python3 vrt_test.py mt_rainier_90m_merc.vrt

Sometimes it will crash on the first run, sometimes it will take quite a few runs before it crashes, but it's always a Floating Point Exception in proj_trans.

Read data from mt_rainier_90m_merc.vrt length=1399
finished 182
Starting attempt 183
Read data from mt_rainier_90m_merc.vrt length=1399
finished 183
Starting attempt 184
Read data from mt_rainier_90m_merc.vrt length=1399
finished 184
Starting attempt 185
Read data from mt_rainier_90m_merc.vrt length=1399
finished 185
Starting attempt 186
Read data from mt_rainier_90m_merc.vrt length=1399
finished 186
Starting attempt 187
Read data from mt_rainier_90m_merc.vrt length=1399
finished 187
Starting attempt 188
Read data from mt_rainier_90m_merc.vrt length=1399
finished 188
Starting attempt 189
Read data from mt_rainier_90m_merc.vrt length=1399
finished 189
Starting attempt 190
Floating point exception

Versions and provenance

The latest version I tried was gdal 3.8.5 with proj 9.4.0 via the ghcr.io/osgeo/gdal:ubuntu-full-3.8.5 docker image.

Additional context

No response