基于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(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】