下面的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_SUPPORTED typedef __int64 GIntBig; typedef unsigned __int64 GUIntBig; #elif HAVE_LONG_LONG typedef long long GIntBig; typedef unsigned long long GUIntBig; #else typedef 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 int GIntBig 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_LONG GIntBig atoi64(const char * str) { return strtoll(str, NULL, 10); } GUIntBig atoi64u(const char * str) { return strtoull(str, NULL, 10); } #else GIntBig 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开发组修复问题还是非常及时的。