GDAL對於raw資料支援的一個bug

下面的GDAL版本針對最新的1.10.0而言，對於GDAL1.9.2問題依舊。

最近在處理衛星的長條帶資料，映像的行數差不多200w行左右，格式的raw格式。可以使用VRT檔案來很方便的對raw資料進行管理。比如要提取raw資料中的某一部分等等。由於處理的資料是1536列的16bit資料，並對raw資料進行提取。比如指定從多少行開始，到多少行結束。關於VRT資料進行管理raw資料可以參考《使用GDAL開啟裸資料（RAW）》http://blog.csdn.net/liminlu0314/article/details/7753377和《使用GDAL下載並轉換SRTM的DEM資料（二）》http://blog.csdn.net/liminlu0314/article/details/8272394這兩篇文章。

本來功能已經完成，今天測試說起始位置超過70w行（最後定位到如果大於等於699051行就會出錯）提取出來的資料就不正確了，開始覺得可能是我的程式寫的有問題，因為我的程式映像的寬高位移量都是用int，最大值應該是2147483647，對於映像的寬高來說，這個數完全夠用，可是對於映像的位移量就不一定了。

映像的位移量就是指映像的第一個像元值先對於raw檔案的起始位置的位元組數，比如上面的問題，就是出在這個映像位移量，由於之前用的是int來進行儲存，對於最大值來說就是2147483647。

由於每行映像是1536列的16bit資料，所以每行就是3072位元組，699051行就是2147484672，這個值已經明顯大於int所能表示的最大值2147483647了，看來問題就是出在這裡了，按照這個思路將程式中的映像位移量全部用__int64進行表示，按理說不會出現問題，可是問題依舊。

初步查看了一下GDAL的原始碼，發現映像的位移量是用的一個名為vsi_l_offset的類型進行表示，查看定義會發現下面這樣的重定義：

/* ==================================================================== *//*      64bit stdio file access functions.  If we have a big size       *//*      defined, then provide protypes for the large file API,          *//*      otherwise redefine to use the regular api.                      *//* ==================================================================== */typedef GUIntBig vsi_l_offset;

以及

/* -------------------------------------------------------------------- *//*      64bit support                                                   *//* -------------------------------------------------------------------- */#if defined(WIN32) && defined(_MSC_VER)#define VSI_LARGE_API_SUPPORTEDtypedef __int64          GIntBig;typedef unsigned __int64 GUIntBig;#elif HAVE_LONG_LONGtypedef long long        GIntBig;typedef unsigned long long GUIntBig;#elsetypedef long             GIntBig;typedef unsigned long    GUIntBig;#endif

從上面的定義來看，vsi_l_offset實質上就是一個64為的無符號整型，按理說不會出錯啊，可是為什麼還是會出錯呢？繼續調試GDAL的代碼，在檔案f:\Work\3rdPart\gdal-1.10.0\frmts\vrt\vrtdataset.cpp的926行左右，會發現下面的代碼：

/* -------------------------------------------------------------------- *//*      Collect required information.                                   *//* -------------------------------------------------------------------- */        if( CSLFetchNameValue(papszOptions, "ImageOffset") != NULL )            nImageOffset = atoi(CSLFetchNameValue(papszOptions, "ImageOffset"));

最後發現是函數atoi出問題了，查看MSDN的說明，會發現atoi函數，對於超出int儲存範圍時會返回int儲存的最大值或者最小值。按照上面的位移量是2147484672，通過上面這句代碼轉換的結果就成2147483647，從而導致後續定位錯誤，導致提取的映像失敗。知道原因，我們就可以使用一個函數來替換atoi這個函數。下面是我寫的一個函數用來替換atoi來支援int64，具體如下：

// define atoi64 transform string to 64bit intGIntBig atoi64(const char * str);GUIntBig atoi64u(const char * str);// function implementations/* -------------------------------------------------------------------- *//*      string to 64bit int support(liml 2013-5-17)                     *//* -------------------------------------------------------------------- */#if defined(WIN32) && defined(_MSC_VER)GIntBig atoi64(const char * str) { return _atoi64(str); }GUIntBig atoi64u(const char * str) { return _strtoui64(str, NULL, 10); }#elif HAVE_LONG_LONGGIntBig atoi64(const char * str) { return strtoll(str, NULL, 10); }GUIntBig atoi64u(const char * str) { return strtoull(str, NULL, 10); }#elseGIntBig atoi64(const char * str) { return _atoi64(str); }GUIntBig atoi64u(const char * str) { return _strtoui64(str, NULL, 10); }#endif

上面的代碼在Windows平台下測試通過，Linux下沒有進行測試，正確性未知。有了上面的函數，就可以用函數atoi64來代替上面的atoi函數即可，然後重新編譯產生GDAL，替換原來的dll即可。

修改上面的代碼後，這個提取的問題解決，可是使用GDAL開啟之後關閉發現又出問題了，顯示的映像不正確。繼續調試代碼，還發現兩處代碼，分別是：檔案f:\Work\3rdPart\gdal-1.10.0\frmts\vrt\vrtrawrasterband.cpp的334行和403行左右：

    nImageOffset = atoi(CPLGetXMLValue( psTree, "ImageOffset", "0") );

/* -------------------------------------------------------------------- *//*      Set other layout information.                                   *//* -------------------------------------------------------------------- */    CPLCreateXMLElementAndValue(         psTree, "ImageOffset",         CPLSPrintf("%d",(int) poRawRaster->GetImgOffset()) );

對於上面第一處，將atoi函數替換為atoi64即可，對於第二處將%d替換為%llu並將後面的強制類型轉換(int)去掉即可。

通過上面的修改，對於raw格式的大映像位移量應該不會出現什麼別的問題了。一個無符號的int64最大值應該是18446744073709551615，換成位元組的話，差不多是2048的PB，這下應該足夠大了。應該不會有一個檔案大小能夠變態到這種地步。

PS：已經將這個bug提交到GDAL開發組，不出所料的話應該下個版本會修複，GDAL開發組修複問題還是非常及時的。