使用GDAL进行栅格数据读取示例

来源:互联网 发布:如何辨别淘宝代购真假 编辑:程序博客网 时间:2024/04/30 11:43
 

当然,使用GDAL可以对其支持的数据进行操作,不局限在栅格数据。但目前主要解决的问题是栅格数据,更具体一点是针对ENVI的标准数据格式。

以下程序仅作为一个简单的应用示例,主要参考了http://www.gdal.org/gdal_tutorial.html中的内容。使用时需要用到一下头文件:
cpl_config.h     
cpl_conv.h       
cpl_error.h      
cpl_minixml.h    
cpl_port.h       
cpl_string.h     
cpl_vsi.h        
gdal.h           
gdal_frmts.h     
gdal_priv.h      
gdal_version.h   
另外需要连接gdal.lib(静态链接)或者gdal_i.lib(动态链接)。程序没有特别的功能,仅作示例而已。

gdal_sample.cpp
///////////////////////////////////
#include "./gdal/gdal_priv.h"

int main(int argc, char **argv)
{
if(argc != 2)
{
    printf( "Usage: gdal_sample <input_file>/n" );
    return 0;
}

char *pszFilename;
pszFilename = argv[1];

// Opening the File
      GDALDataset    *poDataset;
      GDALAllRegister();

      poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );
      if( poDataset != NULL )
      {
    // Getting Dataset Information
    double          adfGeoTransform[6];

    printf( "Driver: %s/%s/n",
      poDataset->GetDriver()->GetDescription(),
      poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );

    printf( "Size is %dx%dx%d/n",
      poDataset->GetRasterXSize(), poDataset->GetRasterYSize(),
      poDataset->GetRasterCount() );

    if( poDataset->GetProjectionRef()    != NULL )
     printf( "Projection is `%s'/n", poDataset->GetProjectionRef() );

    if( poDataset->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
     printf( "Origin = (%.6f,%.6f)/n",
       adfGeoTransform[0], adfGeoTransform[3] );

     printf( "Pixel Size = (%.6f,%.6f)/n",
       adfGeoTransform[1], adfGeoTransform[5] );
    }

    // Fetching a Raster Band
    GDALRasterBand    *poBand;
          int               nBlockXSize, nBlockYSize;
          int               bGotMin, bGotMax;
          double            adfMinMax[2];
        
          poBand = poDataset->GetRasterBand( 1 );
          poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
          printf( "Block=%dx%d Type=%s, ColorInterp=%s/n",
                  nBlockXSize, nBlockYSize,
                  GDALGetDataTypeName(poBand->GetRasterDataType()),
                  GDALGetColorInterpretationName(
                      poBand->GetColorInterpretation()) );

          adfMinMax[0] = poBand->GetMinimum( &bGotMin );
          adfMinMax[1] = poBand->GetMaximum( &bGotMax );
          if( ! (bGotMin && bGotMax) )
              GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);

          printf( "Min=%.3fd, Max=%.3f/n", adfMinMax[0], adfMinMax[1] );
        
          if( poBand->GetOverviewCount() > 0 )
              printf( "Band has %d overviews./n", poBand->GetOverviewCount() );

          if( poBand->GetColorTable() != NULL )
              printf( "Band has a color table with %d entries./n",
                       poBand->GetColorTable()->GetColorEntryCount() );

    // Reading Raster Data
    float *pafScanline;
          int     nXSize = poBand->GetXSize();

          pafScanline = (float *) CPLMalloc(sizeof(float)*nXSize);
          poBand->RasterIO( GF_Read, 0, 0, nXSize, 1,
                            pafScanline, nXSize, 1, GDT_Float32,
                            0, 0 );

    // Closing the Dataset
    GDALClose(poDataset);
      }

return 0;
}