基於GDAL庫，讀取.grd文件（以海洋地形數據為例）C++版
技術背景
海洋地形數據主要是通過美國全球地形起伏數據（GMT）獲得，數據格式為grd（GSBG）二進制數據，打開軟件通過是Surfer軟件，surfer軟件可進行數據的編輯處理，以及進一步的可視化表達等功能操作；由於Surfer軟件不支持二次開發，沒有提供相應的SDK供開發者進行使用，所以這一切只能通過相應類似的技術進行實現，首先，數據的讀取，如何通過編程實現數據的讀取操作呢？這裏就要說一下GIS軟件所使用的一個開源庫-GDAL，GDAL庫的具體解釋資料，請查閱官方網站【https://www.gdal.org/index.html】，由於後期要進行數據入庫的步驟，所以本文提供的是一種采用C++語言（Qt平臺）進行讀取的方法，前面的GDAL庫的編譯方法請參考博客4【gdal庫編譯並適配至vs2017】，下面具體講一下環境配置。
環境配置有兩種方法，這個看個人習慣，主要是看自己使用時是否是配置環境變量，編譯好的庫文件主要有一下文件夾
其中bin文件夾主要是編譯後的運行文件要使用的，也就是Qt編譯出來的debug文件或者是release文件會在運行的時候調用這個文件夾裏面的程序，而lib文件夾是你編程時候要使用的庫文件，編譯是否正確跟這個庫文件有關。方法一：配置環境變量法，將bin文件夾添加到系統環境變量path內，然後在新建的項目文件添加庫文件（外部庫），添加完成後可見*.pro文件裏面多了
這樣幾個引用語句，這就表明已經成功導入項目了，直接開始後面的編輯工作了；方法二：直接導入法，直接將bin文件夾裏面的*.dll文件導入到項目的debug或release文件裏面，然後在新建的項目文件添加庫文件（外部庫），添加完成後*.pro文件和上文提到的一樣，多了那三行引用語句。集體來說，我推薦第一種方法，因為一勞永逸，不用每個項目都要去導一下。
數據讀取
下面進行數據的讀取，數據的讀取步驟和正常GDAL庫讀取差不多，沒什麽太大的區別，具體理論步驟請大家參考博客3【gdal讀寫圖像分塊處理（精華版）】，這裏面唯一要註意的就是數據的緩沖區大小，因為grd為二進制數據，數據的值是否正確和你緩沖區大小設置有關系，這裏我是用的是float，也就是代碼中的GDT_Float32，不說了，直接看代碼吧。
讀取-頭文件
1 #ifndef GRDFILEREAD_H 2 #define GRDFILEREAD_H 3 4 class grdFileRead 5 { 6 public: 7 void fileRead(constchar* pszFile); 8 }; 9 10 #endif // GRDFILEREAD_H
讀取-源文件
1 #include "grdfileread.h" 2 3 #include <QtDebug> 4 #include <gdal_priv.h> 5 6 void grdFileRead::fileRead(const char* pszFile) 7 { 8 GDALAllRegister(); 9 GDALDataset *poDataset; 10 //使用只讀方式打開圖像 11 poDataset = (GDALDataset*) GDALOpen( pszFile,GA_ReadOnly ); 12 if( poDataset == NULL ){ 13 printf( "File: %s不能打開！\n",pszFile); 14 return; 15 } 16 17 printf( "Driver:%s/%s\n", 18 poDataset->GetDriver()->GetDescription(), 19 poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME) ); 20 21 //輸出圖像的大小和波段個數 22 printf( "Size is%dx%dx%d\n", 23 poDataset->GetRasterXSize(),poDataset->GetRasterYSize(), 24 poDataset->GetRasterCount()); 25 26 27 //輸出圖像的投影信息 28 if( poDataset->GetProjectionRef() != NULL ) 29 printf( "Projectionis `%s‘\n", poDataset->GetProjectionRef() ); 30 31 //輸出圖像的坐標和分辨率信息 32 double adfGeoTransform[6]; 33 if( poDataset->GetGeoTransform( adfGeoTransform) == CE_None ){ 34 printf( "Origin =(%.6f,%.6f)\n",adfGeoTransform[0], adfGeoTransform[3]); 35 printf( "PixelSize = (%.6f,%.6f)\n",adfGeoTransform[1], adfGeoTransform[5]); 36 } 37 38 //讀取第一個波段 39 GDALRasterBand *poBand = poDataset->GetRasterBand( 1 ); 40 41 //獲取該波段的最大值最小值，如果獲取失敗，則進行統計 42 int bGotMin, bGotMax; 43 double adfMinMax[2]; 44 adfMinMax[0] = poBand->GetMinimum( &bGotMin); 45 adfMinMax[1] = poBand->GetMaximum( &bGotMax); 46 47 if( ! (bGotMin&& bGotMax) ) 48 GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax); 49 50 printf( "Min=%.3fd,Max=%.3f\n", adfMinMax[0], adfMinMax[1] ); 51 52 int nXSize = poBand->GetXSize(); 53 int nYSize = poBand->GetYSize(); 54 float *pafScanline = new float[nXSize]; 55 56 //讀取圖像數據 57 for(int i = 0; i< 10/*nYSize*/;i++){ 58 poBand->RasterIO(GF_Read, 0, i, nXSize,1, pafScanline, nXSize,1, GDT_Float32, 0, 0 ); 59 60 QString LineDataInfo = ""; 61 for(int j = 0; j< 10/*nXSize*/;j++){ 62 if(j == 0){ 63 LineDataInfo = QString("%1").arg(pafScanline[j], 0, ‘f‘, 0); 64 }else{ 65 LineDataInfo = LineDataInfo + ", " + QString("%1").arg(pafScanline[j], 0, ‘f‘, 1); 66 } 67 68 } 69 qDebug() << LineDataInfo << endl; 70 } 71 72 delete []pafScanline; 73 74 //關閉文件 75 GDALClose((GDALDatasetH)poDataset); 76 }
主函數調用
1 #include <QCoreApplication> 2 3 #include "grdfileread.h" 4 #include <QWidget> 5 6 int main(int argc, char *argv[]) 7 { 8 QCoreApplication a(argc, argv); 9 10 grdFileRead gfr; 11 gfr.fileRead("F:/Data File/test/E135N30_sf.grd"); 12 13 return a.exec(); 14 }
運行結果
至此，文件讀取完成。
致謝
感謝李民錄老師的指導，以及相關技術博主的技術分享，謝謝你們！
參考博客
1、GDAL遙感影像讀取與顯示【https://blog.csdn.net/zhouxuguang236/article/details/8070090】
2、Qt配置GDAL（Qt 5.6.1+MSVC 2013+64 bit）【https://blog.csdn.net/u010670734/article/details/53106786?locationNum=13&fps=1】
3、gdal讀寫圖像分塊處理（精華版）【https://blog.csdn.net/elfxwt_study/article/details/9071261】
4、gdal庫編譯並適配至vs2017【https://blog.csdn.net/Dragonzxc/article/details/80356883】
基於GDAL庫，讀取.grd文件（以海洋地形數據為例）C++版