基於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 voidfileRead(const char* 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】