简单学习了一下GDAL中的GDALRasterBand中API的用法, 以前只用过GetRasterDateType(), RasterIO(), GetXSize(), GetYSize()这些简单的用法, 今天简单学习GDALRasterBand中更深入的用法, 记录如下:
1 获得波段统计数据
代码如下:
//影像统计数据-GetStatistics()
int bApproxOK = false;
int bForce=true;
double pdfMin, pdfMax, pdfMean, pdfStdDev;
pBand->GetStatistics(bApproxOK, bForce, &pdfMin, &pdfMax, &pdfMean, &pdfStdDev);
cout << "pdfMin: " << pdfMin <<" " << "pdfMax: " << pdfMax << " " << "pdfMean: " << pdfMean
<< " " << "pdfStdDev: " << pdfStdDev << " " << endl;
根据API文档说明, GetStatistics()函数原型如下:
CPLErr GDALRasterBand::GetStatistics (int bApproxOK,
int bForce,
double * pdfMin,
double * pdfMax,
double * pdfMean,
double * pdfStdDev
)
其中, 参数列表所代表的含义为:
- bApproxOK [in] bool型, 如果为TRUE则通过波段的overview(概览?这是啥东西?目前没看懂…)或者a subset of all tiles(所有瓦片的子集?这也不太懂…)得到统计数据. 反正我用的false.
- bForce [in] bool型, 如果为FALSE的话, 统计数据只在不重新扫描波段的情况下得到. 也就是说, FALSE效率高, 但是结果不能保证. 经过试验, 确实FALSE的结果不太稳定, 还是用TRUE为好.
- pdfMin [out] 得到波段的最小DN值.
- pdfMax [out] 得到波段的最大DN值
- pdfMean [out] 得到波段的DN值均值
- pdfStdDev [out] 得到波段的DN值标准差
经过实验, 可以得到输入波段的DN最大最小值, 均值和方差的数据, 还是很方便的.
2 获得波段的DN值的最大最小值
有时候, 不需要那么多统计数据, 只需要一个最大最小值就可以, 那么可以用到这三个方法: GetMaximum(), GetMinimum()以及GDALComputeRasterMinMax(). 代码如下:
//波段中最大值-GetMaximum()
double max = pBand->GetMaximum(&bGotMax);
cout << "波段最大值: " << max << endl;
//波段中的最小值-GetMinimum()
double min = pBand->GetMinimum(&bGotMin);
cout << "波段最小值: " << min << endl;
double maxmin[2];
if (!(bGotMax&&bGotMin))
{
GDALComputeRasterMinMax((GDALRasterBandH)pBand, true, maxmin);
}
cout << "再次计算的最小值为: " << maxmin[0] << endl;
cout << "再次计算的最大值为: " << maxmin[1] << endl;
其中, 他们的函数原型分别为:
double GDALRasterBand::GetMaximum(int *pbSuccess = NULL )
double GDALRasterBand::GetMinimum(int *pbSuccess = NULL )
void GDALComputeRasterMinMax(GDALRasterBandH hBand,
int bApproxOK,
double adfMinMax[2])
对于GetMaximum()和GetMinimum()的参数列表,pbSucess是bool型指针, 代表着得到最大最小值是否成功, 成功返回1, 失败返回0. 原文中描述如下:”对于不知道文件格式的波段, 最大(最小)值将会返回波段数据类型的最大(最小)值.”
也就是说,如果GDAL不能识别你输入影像的格式,这个GetMaximum()和GetMinimum()基本就废了…而我输入的影像是TIFF格式的,而且是BSQ排列的三个波段的影像,然而GDAL通过这两个函数给我返回的结果并不正确….所以没办法,我只好用了GDAL中面向C语言的一个方法,GDALComputeRasterStatistics().
GDALComputeRasterStatistics()参数列表说明如下:
- hBand [in] 输入波段(记得如果是GDALRasterBand类型的要强制转换到GDALRasterBandH类型).
- bApproxOK [in] 与上面GetStatistics()中的一样,反正我用的TRUE.
adfMinMax[2] [out] 输出得到的最大最小值,adfMinMax[1]为最小,adfMinMax[2]为最大.
用GDALComputeRasterMinMax()就得到了最大最小值,还是比较靠谱的.
3 获得波段的直方图
之前都不知道GDAL中还有这个功能,都是自己写统计直方图的方法,结果今天才发现有这么一个利器.
//波段直方图-GetHistogram()
cout << "输出直方图..." << endl;
const int DS = 1024;
GUIntBig iHistogram[DS];
pBand->GetHistogram(0-0.5, DS-0.5, DS, iHistogram, false, false, NULL, NULL);
//直方图输出到CSV中
FILE *pFile = fopen("BandHistogram.csv", "w");
for (int i = 0; i < DS; i++)
{
fprintf(pFile,"%d,%d\r",i,(int)iHistogram[i]);
}
fclose(pFile);
cout << "输出完毕!"<<endl;
其中,GetHistogram()函数原型如下:
CPLErr GDALRasterBand::GetHistogram(double dfMin,
double dfMax,
int nBuckets,
GUIntBig *panHistogram,
int bIncludeOutOfRange,
int bApproxOK,
GDALProgressFunc pfnProgress,
void *pProgressData)
参数列表解释:
- dfMin [in] 直方图x轴最小DN值,在我的代码中,我用的-0,5,这样可以去除边缘算不算这种蛋疼的问题.
- dfMax [in] 直方图x轴最大DN值,与上面类似.
- nBuckets [in] 分多少个块,如果要是一个灰度一级的话,nBuckets=dfMax-dfMin,以此类推.
- panHistogram [out] 得到直方图,注意:GUIntBig类型是GDAL中自己define的一个东西,其实就是usigned long long.
- bIncludeOutOfRange [in] bool型,如果是TRUE的话,低于直方图统计范围的值全部都算在panHistogram[0]中,高于直方图统计范围都算在panHistogram[nBuckets-1]中;如果为FALSE,则忽略超出直方图范围的值.
- bApproxOK [in] TRUE的话,返回的直方图也可以是一个近似的,不完全的.
- pfnProgress [in] 报告进度一个函数句柄
- pProgressData [in] 向pfnProgress传递的数据.
通过这个方法,我得到了波段的直方图,而且效率也还不错.
参考资料:GDAL的API文档 http://www.gdal.org/classGDALRasterBand.html#aaabf419931d0f505428f91cff54085cc