数据读写

本章我们就开始按顺序详细介绍GDAL读写数据的过程。最后将提供完整的读写流程代码。完整代码部分中有整个创建、两影像相加、写入的流程,如果已经大致了解GDAL的读写流程,可以直接参照。

GdalDriver

首先,GDAL对每种格式提供了一个驱动 GdalDriverGdalDriver将对对应格式的数据进行管理,例如读取、创建、删除、重命名、复制、从已有数据创建新数据集等。所以,我们所有程序开头,都将添加GDALAllRegister() 函数,注册所有GDAL支持的数据驱动。

再次强调,所有的程序都要首先调用GDALAllRegister() 函数,否则将无法打开任何数据。

数据读取

GDAL中数据读取的步骤如下:

  • 打开数据集

  • 打开数据集下所需的波段

  • 读取数据

下面分步讲解如何操作:

打开数据集

打开 multibands 步骤如下:使用 GDALOpen 或者 GDALOpenShared函数,传入文件名,打开数据集。 GDALOpenGDALOpenShared参数一样,区别在于,在同一线程,如果是相同的文件,多个 GDALOpenShared打开的其实是同一个 GDALDataset的引用(如果在不同线程使用,为了保证线程安全,返回的将是不同的对象)。

      #include "gdal_priv.h"
      #include "cpl_conv.h" 

      //...中间的程序

      //所有程序前先加上
      GDALAllRegister();
      //pszFilename代表文件名,GA_ReadOnly表示以只读方式打开
      //也可以使用GA_Update
      //GDALOpenShared和GDALOpen可以互换
      GDALDataset  *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
      if( poDataset == NULL )
      {
            ...;//打开失败处理
      }

打开 multidatasets方式与上面稍有区别,因为有多个子数据集,所以需要获取真实数据的话,实际上是获取子数据集中的波段,可以使用gdalinfo 获取如下的 subdatasetsdomain :

  Subdatasets:
        SUBDATASET_1_NAME=NETCDF:"example.nc":auditTrail
        SUBDATASET_1_DESC=[2x80] auditTrail (8-bit character)
        SUBDATASET_2_NAME=NETCDF:"example.nc":data
        SUBDATASET_2_DESC=[1x61x172] data (32-bit floating-point)
        SUBDATASET_3_NAME=NETCDF:"example.nc":lat
        SUBDATASET_3_DESC=[61x172] lat (32-bit floating-point)
        SUBDATASET_4_NAME=NETCDF:"example.nc":lon
        SUBDATASET_4_DESC=[61x172] lon (32-bit floating-point)

打开子数据集也使用 GDALOpen 或者 GDALOpenShared函数,但是传入的不是整个数据集的名称,而是需要打开的SUBDATASET_N_NAME= (_N_中N为第几个子数据集,即1 2 3 4...)后面的字符串,即子数据集名称。例如需要打开上例中data子数据集,使用如下代码:

      #include "gdal_priv.h"
      #include "cpl_conv.h" 

      //...中间的程序

      //所有程序前先加上
      GDALAllRegister();
      //也可以使用GA_Update
      //GDALOpenShared和GDALOpen可以互换
      //pszFilename代表子数据集的NAME,注意双引号需要转义,GA_ReadOnly表示以只读方式打开
      char * pszFilename = "NETCDF:\"example.nc\":data";
      GDALDataset  *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
      if( poDataset == NULL )
      {
            ...;//打开失败处理
      }

打开了数据集后,两种数据接下来的处理方式都是一样的。

打开波段

GDAL中, 波段起始索引为1 ,所以循环时需要特别注意。使用GDALDataset->GetRasterBand(int BandIndex) 函数,可以打开波段。

      #include "gdal_priv.h"
      #include "cpl_conv.h" 
      //...中间的程序

      //所有程序前先加上
      GDALAllRegister();
      GDALDataset  *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );//打开文件
      if( poDataset == NULL )
      {
            //...//打开失败处理
      }
      //打开数据集同上
      //获取影像相关信息
      int i,j,width,height,bandNum;
      width = inDS->GetRasterXSize();//影像宽度
      height = inDS->GetRasterYSize();//影像高度
      bandNum = inDS->GetRasterCount();//影像波段数量

      //将第一波段读入
      double *data = new double[width*height];
      GDALRasterBand *band = ds1->GetRasterBand(1);//获取第一波段,波段从1开始
      //如果是在循环内,需要注意,波段从1开始;
      //for(i=0;i < bandNum;i++){
      // GDALRasterBand *band = ds1->GetRasterBand(i+1);//这里注意下
      // ......对波段处理
      //}

读取内容

使用 RasterIO函数获取波段中的具体内容,该函数的功能强大,参数较多,我们先用代码演示,然后用图来表示。

函数说明:

      ///CPLErr GDALRasterBand::RasterIO (   
      ///@param eRWFlag,       //读取或者写入,GF_Read或GF_Write
      ///@param nXOff,         //起始点x坐标
      ///@param nYOff,         //起始点y坐标
      ///@param nXSize,        //所需读取(写入)块宽度
      ///@param nYSize,        //所读(写)块高度
      ///@param * pData,       //所读(写)数据,指针,
      ///@param nBufXSize,   //一般跟nXSize一致,用于缩放图像,
      ///                                图像将按nBufXSize/nXsize在x尺度缩放(会自动重采样)
      ///                               ,一般不需要调整
      ///@param nBufYSize,   //一般跟nYSize一致
      ///@param eBufType,      //与pData的实际类型一致,GDT_Float64
      ///                                代表double,其他的可以跳到定义查看  
      ///@param nPixelSpace,  
      ///                                设置为0为自动判断,一般设为0
      ///                                表示的是当前像素值和下一个像素值之间的间隔,单位是字节
      ///                                例:byte类型,就是1,double类型,就是8
      ///@param nLineSpace      
      ///                                设置为0为自动判断,一般设为0
      ///                                表示的是当前行和下一行的间隔,,单位是字节
      ///                                例如,一行300像素,类型为int,此时就是300*4 = 1200
      ///@return                   //是否成功,成功返回CE_None ,失败返回 CE_Failure 
      CPLErr GDALRasterBand::RasterIO (   
            GDALRWFlag eRWFlag, 
            int       nXOff, 
            int       nYOff, 
            int       nXSize,              
            int       nYSize,   
            void * pData,      
            int       nBufXSize,               
            int       nBufYSize,   
            GDALDataType      eBufType,
            int       nPixelSpace,        
            int       nLineSpace      
      ) 

示意图:

RasterIO.png

读取单幅影像第一波段原始数据代码示例:

      #include "gdal_priv.h"
      #include "cpl_conv.h" 
      //...中间的程序

      //所有程序前先加上
      GDALAllRegister();
      GDALDataset  *poDataset= (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );//打开文件
      if( poDataset == NULL )
      {
            //...打开失败处理
      }
      //打开数据集同上
      //获取影像相关信息
      int i,j,width,height,bandNum;
      width = inDS->GetRasterXSize();//影像宽度
      height = inDS->GetRasterYSize();//影像高度
      bandNum = inDS->GetRasterCount();//影像波段数量

      //将第一波段读入
      double *data = new double[width*height];
      GDALRasterBand *band = ds1->GetRasterBand(1);//获取第一波段,波段从1开始
      band->RasterIO(GF_Read,0,0,width,height,data,width,height,GDT_Float64,0,0);

warning

使用RasterIO函数时需要注意以下几点:

  • 波段索引从一开始计算

  • pData数组要有足够的大小

  • pData数组数据类型要和eBufType对应

  • 读取完成后必须要释放数组和关闭数据集,波段不需要关闭

关闭数据集

程序结束前,必须要关闭数据集,使用 GDALClose() 函数。

      GDALClose(ds2);

数据写入

数据写入分为几个步骤

  • 获取驱动

  • 创建数据集

  • 写入数据

  • 关闭数据集

获取驱动

gdaldriver 中已经说明,不论读取还是写入, GDALDriver都是很重要的存在。如果需要输出特定格式,我们必须首先获取该格式的驱动,才能创建该格式的数据集。获取方式如下:

      GDALDriver *poDriver;

      //一般使用tif作为输出,如果有特殊需求,请参阅gdal文档,一般修改这里改为其他驱动名称即可,例如PNG
      const char *pszFormat ="GTiff";
      poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);//获取特殊的驱动。
      if(poDriver == NULL) {
            return;
      }

创建数据集

根据不同输出需要,创建数据集有两种方式:

  • 输出的影像与输入的影像大小、投影都相同,只有数据不同,就可以从已经读入的数据集中创建

  • 投影或大小不同,需要创建新的数据集

下面分别说明两种方式如何创建数据集,下文都是以输出GeoTiff作为示例,不同的GdalDriver 中,可以设置的参数不同,具体请参考GDAL支持格式中对每种格式的说明。

从已读入的数据集中创建

从已经读入的数据集中创建数据集,可以使用 GDALDriverCreateCopy函数。

CreateCopy 函数原型如下:

      //从已有的dataset中创建
      //@param pszFilename, 输出文件名
      //@param poSrcDS, 已有的dataset
      //@param bStrict, TRUE表示严格等价,一般设置为FALSE,表示拷贝副本用作编辑
      //@param papszOptions, 参数设置,具体可以参考每个driver的文档
      //@param  pfnProgress 回调函数指针
      //@param *pProgressData 传入回调函数的数据
      GDALDataset * GDALDriver::CreateCopy       (
            const char *       pszFilename, 
            GDALDataset *       poSrcDS, 
            int       bStrict, 
            char **       papszOptions, 
            GDALProgressFunc       pfnProgress, 
            void *       pProgressData      
      )

使用 CreateCopy 函数创建数据集代码如下:

      GDALDataset *OutDs;
      char **papszOptions = NULL;//设置压缩、存储方式等,每种格式不同,可以直接设置为NULL
                                             //也可以根据需要设置,例如:
      //papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
      //papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );

      //创建与ds1相同的坐标信息、相同大小的文件
      OutDs = poDriver->CreateCopy(OutPut1,ds1,FALSE,papszOptions,NULL,NULL);

直接创建

直接创建数据集比从已有数据集中创建会多出创建投影的步骤:

创建数据集

直接集中创建数据集,可以使用 GDALDriverCreate 函数。

Create 函数原型如下:

      //Create文件
      //@param pszFilename, 输出文件名
      //@param nXSize, 文件宽度
      //@param nYSize, 文件高度
      //@param nBands, 波段数
      //@param eType, 数据类型
      //@param papszOptions       参数设置,具体可以设置的属性参考每个driver的文档
      //@return dataset
      GDALDataset * GDALDriver::Create ( 
            const char *       pszFilename, 
            int       nXSize, 
            int       nYSize, 
            int       nBands, 
            GDALDataType       eType, 
            char **       papszOptions 
      )

使用 Create 函数创建数据集代码如下:

      GDALDataset *OutDs;
      ///设置:
      char **papszOptions = NULL;
      papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
      papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
      ///...
      ///创建512*512,单波段的,byte类型的数据
      OutDs = poDriver->Create( pszDstFilename, 512, 512, 1, GDT_Byte, 
                                                papszOptions );
创建投影

在 coordinatesystem 和 affinegeotransform中已经介绍过GDAL中数据如何投影,下面我们将用代码示例:

      ///设置仿射变换参数
      double adfGeoTransform[6] = { 444720, 30, 0, 3751320, 0, -30 };
      OutDs->SetGeoTransform( adfGeoTransform );

      //设置投影
      OGRSpatialReference oSRS;
      char *pszSRS_WKT = NULL;
      oSRS.SetUTM( 11, TRUE );
      oSRS.SetWellKnownGeogCS( "NAD27" );
      oSRS.exportToWkt( &pszSRS_WKT );
      poDstDS->SetProjection( pszSRS_WKT );
      CPLFree( pszSRS_WKT );//使用完后释放

写入数据

写入数据与读取数据基本相同,都是使用 GDALRasterBandRasterIO函数,只是第一个参数变为 GF_Write ,下面就直接给出代码,不再另作说明。

      float * outData = new float [512*512];

      //对outData进行计算、处理
      /************************************************************************/
      /*********************************处理************************************/
      /************************************************************************/

      GDALRasterBand *outBand = OutDs->GetRasterBand(1);
      outBand->RasterIO(GF_Write,0,0,512,512,outData,512,512,GDT_Float32,0,0);

关闭数据集

数据集必须关闭后,数据才会写入到输出文件中,否则数据将在缓存中,使用GDALClose() 函数关闭数据集。

      GDALClose(OutDs);

完整代码

      #include "gdal_priv.h"
      int main(int argc, char* argv[])
      {
            GDALAllRegister();//注册类型,读取写入任何类型影像必须加入此句
            //以只读方式打开文件
            GDALDataset *ds1 = (GDALDataset *) GDALOpen("Input1.tif", GA_ReadOnly);
            if(ds1 == NULL) {
                  printf("不能打开第一个文件,请检查文件是否存在!");
                  return -1;
            }
            int WidthAll = ds1->GetRasterXSize();
            int HightAll = ds1->GetRasterYSize();
            int BandNum = ds1->GetRasterCount();
            float** InBuf1 = new float*[BandNum];
            int i,j,k;

            //分波段读取
            for(i=0; i<BandNum; i++) {
                  InBuf1[i] = new float[WidthAll*HightAll];
                  GDALRasterBand *band = ds1->GetRasterBand(i+1);//获取波段,波段从1开始
                  band->RasterIO(GF_Read,
                                       0,
                                       0,
                                       WidthAll,
                                       HightAll,
                                       InBuf1[i], 
                                       WidthAll,
                                       HightAll,
                                       GDT_Float32,
                                       0,
                                       0);
            }
            //打开第二个文件,假设第一个文件和第二个文件大小、波段数相同
            //以只读方式打开文件
            GDALDataset *ds2 = (GDALDataset *) GDALOpen("Input2.tif", GA_ReadOnly);

            if(ds2 == NULL) {
                  printf("不能打开第二个文件,请检查文件是否存在!");
                  return ;
            }

            //第二个文件中的数据
            float** InBuf2 = new float*[BandNum];

            for(i=0; i<BandNum; i++) {
                  InBuf2[i] = new float[WidthAll*HightAll];
                  GDALRasterBand *band = ds2->GetRasterBand(i+1);
                  band->RasterIO(GF_Read,
                                       0,
                                       0,
                                       WidthAll,
                                       HightAll,
                                       InBuf2[i], 
                                       WidthAll,
                                       HightAll,
                                       GDT_Float32,
                                       0,
                                       0);
            }
            //写入的数据(可以最后创建写入文件,需要先创建写入数据)
            float** outBuf = new float*[BandNum];

            for(i=0; i<BandNum; i++) {
                  outBuf[i] = new float[WidthAll*HightAll];
            }

            ////////////////////////////////////////////////////
            /////////////////////处理///////////////////////////
            ////////////////////////////////////////////////////
            //简单相加示例,实际上,这部分最好写入函数,或者类中单独封装起来,方便使用
            for(j=0; j<BandNum; j++) {
                  for(k=0; k<WidthAll*HightAll; k++) {
                        outBuf[j][k] = InBuf1[j][k]+InBuf2[j][k];
                  }
            }
            ////////////////////////////////////////////////////
            /////////////////////处理///////////////////////////
            ////////////////////////////////////////////////////
            //创建写入图像:
            GDALDriver *poDriver;
            const char *pszFormat ="GTiff";
            poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
            if(poDriver == NULL)
                return;
            GDALDataset *OutDs;
            char **papszOptions = NULL;
            char **papszMetadata;
            //这里的参数全是指tif格式的参数,如果是其他格式,请把这里所有注释掉,或者参照文档,自行设定
            
            //设置压缩类型,envi只认得packbits压缩.
            //papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
            
            //设置压缩比,可以不用,有些时候压缩反而更大,因为是无损的,除非图片很大
            //papszOptions = CSLSetNameValue( papszOptions, "ZLEVEL", "9" );
            
            //设置bsq或者BIP存储 bsq:BAND,bip:PIXEL
            papszOptions = CSLSetNameValue( papszOptions, "INTERLEAVE", "BAND" );

            OutDs = poDriver->Create( "output.tif", 
                                                  WidthAll,
                                                  HightAll, 
                                                  BandNum, 
                                                  GDT_Float32, 
                                                  papszOptions );
            double geos[6];
            
            //获取第一幅图的地理变化信息,这里如果与原图不同,请自己计算
            ds1->GetGeoTransform(geos);
            char *pszProjection = NULL;
            char *pszPrettyWkt = NULL;
            //设置为第一幅图的投影,如果与原图不同,请跳过这个if部分,到下个部分直接设置投影
            //设置投影,如果需要特殊投影,请找到wkt串,自行建立投影后传入
            OGRSpatialReferenceH  hSRS;
            if( GDALGetProjectionRef( ds1 ) != NULL ){
                pszProjection = (char *) GDALGetProjectionRef( ds1 );

                hSRS = OSRNewSpatialReference(NULL);
                if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )
                {
                    OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );
                    OutDs->SetProjection(pszPrettyWkt);
                }
            }
            //pszPrettyWkt实际上是ogc wkt串或者是proj4
            //如果指定投影的,可以在http://spatialreference.org/ 网站中搜索
            //找到所需要的投影后,例如西安80  3度带,中央经线117E,就是EPSG:2384,点开后
            //点击ogc wkt或者proj4,抄下来就行,proj4比较短,而且没有双引号不需要转义,比较简单
            //pszPrettyWkt = "+proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000  
            //                        +y_0=0 +a=6378140 +b=6356755.288157528 +units=m +no_defs";
            //OutDs->SetProjection(pszPrettyWkt);
            //设置坐标
            OutDs->SetGeoTransform(geos);
            CPLFree( pszPrettyWkt );
            //写入数据到outds中
            for(i=0; i<BandNum; i++) {
                  GDALRasterBand *outBand = OutDs->GetRasterBand(i+1);
                  outBand->RasterIO(GF_Write,
                                            0,
                                            0,
                                            WidthAll,
                                            HightAll,
                                            outBuf[i],
                                            WidthAll,
                                            HightAll,
                                            GDT_Float32,
                                            0,
                                            0);
            }

            //关闭dataset之后,数据才会真正写入到文件中,否则都是在缓存中!!!
            GDALClose(OutDs);
            GDALClose(ds1);
            GDALClose(ds2);
            
            //清理环境,删除new的数组、变量等
            for(i=0; i<BandNum; i++) {
                  delete[] InBuf1[i];
                  delete[] InBuf2[i];
                  delete[] outBuf[i];
            }

            delete[] InBuf1;
            delete[] InBuf2;
            delete[] outBuf;

            getchar();
            return 0;
      }

其他处理

可以使用 gdaldriver 对已有影像数据进行删除、复制、重命名等处理,可以参见GdalDriver类中的具体函数说明,下面使用一些例子来说明:

      GDALAllRegister();//注册类型,读取写入任何类型影像必须加入此句
      GDALDriver *poDriver;
      const char *pszFormat ="GTiff";
      poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
      poDriver->Delete("TempFile.tif");//删除
      poDirver->Rename("newName.tif","oldName.tif");//重命名
      poDirver->CopyFiles("Replica.tif","original.tif");//复制