C++ GDALRasterBand::GetRasterDataType方法代码示例

int GDAL_EDBFile::ReadBlock( int channel,
                             int block_index, void *buffer,
                             int win_xoff, int win_yoff,
                             int win_xsize, int win_ysize )

{
    GDALRasterBand *poBand = poDS->GetRasterBand(channel);
    int nBlockXSize, nBlockYSize;
    int nBlockX, nBlockY;
    int nWidthInBlocks;
    int nPixelOffset;
    int nLineOffset;

    if( GetType(channel) == CHN_UNKNOWN )
    {
        ThrowPCIDSKException("%s channel type not supported for PCIDSK access.",
                             GDALGetDataTypeName(poBand->GetRasterDataType()) );
    }

    poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );

    nWidthInBlocks = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
    
    nBlockX = block_index % nWidthInBlocks;
    nBlockY = block_index / nWidthInBlocks;

    nPixelOffset = GDALGetDataTypeSize(poBand->GetRasterDataType()) / 8;
    nLineOffset = win_xsize * nPixelOffset;

/* -------------------------------------------------------------------- */
/*      Are we reading a partial block at the edge of the database?     */
/*      If so, ensure we don't read off the database.                   */
/* -------------------------------------------------------------------- */
    if( nBlockX * nBlockXSize + win_xoff + win_xsize > poBand->GetXSize() )
        win_xsize = poBand->GetXSize() - nBlockX * nBlockXSize - win_xoff;

    if( nBlockY * nBlockYSize + win_yoff + win_ysize > poBand->GetYSize() )
        win_ysize = poBand->GetYSize() - nBlockY * nBlockYSize - win_yoff;

    CPLErr eErr = poBand->RasterIO( GF_Read, 
                                    nBlockX * nBlockXSize + win_xoff, 
                                    nBlockY * nBlockYSize + win_yoff,
                                    win_xsize, win_ysize, 
                                    buffer, win_xsize, win_ysize, 
                                    poBand->GetRasterDataType(),
                                    nPixelOffset, nLineOffset, NULL );

    if( eErr != CE_None )
    {
        ThrowPCIDSKException( "%s", CPLGetLastErrorMsg() );
    }

    return 1;
}

int GDAL_EDBFile::WriteBlock( int channel, int block_index, void *buffer)

{
    GDALRasterBand *poBand = poDS->GetRasterBand(channel);
    int nBlockXSize, nBlockYSize;
    int nBlockX, nBlockY;
    int nWinXSize, nWinYSize;
    int nWidthInBlocks;

    if( GetType(channel) == CHN_UNKNOWN )
    {
        ThrowPCIDSKException("%s channel type not supported for PCIDSK access.",
                             GDALGetDataTypeName(poBand->GetRasterDataType()) );
    }

    poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );

    nWidthInBlocks = (poBand->GetXSize() + nBlockXSize - 1) / nBlockXSize;
    
    nBlockX = block_index % nWidthInBlocks;
    nBlockY = block_index / nWidthInBlocks;

/* -------------------------------------------------------------------- */
/*      Are we reading a partial block at the edge of the database?     */
/*      If so, ensure we don't read off the database.                   */
/* -------------------------------------------------------------------- */
    if( nBlockX * nBlockXSize + nBlockXSize > poBand->GetXSize() )
        nWinXSize = poBand->GetXSize() - nBlockX * nBlockXSize;
    else
        nWinXSize = nBlockXSize;

    if( nBlockY * nBlockYSize + nBlockYSize > poBand->GetYSize() )
        nWinYSize = poBand->GetYSize() - nBlockY * nBlockYSize;
    else
        nWinYSize = nBlockYSize;

    CPLErr eErr = poBand->RasterIO( GF_Write, 
                                    nBlockX * nBlockXSize, 
                                    nBlockY * nBlockYSize,
                                    nWinXSize, nWinYSize,
                                    buffer, nWinXSize, nWinYSize,
                                    poBand->GetRasterDataType(), 0, 0, NULL );

    if( eErr != CE_None )
    {
        ThrowPCIDSKException( "%s", CPLGetLastErrorMsg() );
    }

    return 1;
}

Raster* import_raster(string raster_filename, int band_number) {
    GDALAllRegister();
    GDALDataset*  poDataset = (GDALDataset *) GDALOpen( raster_filename.c_str(), GA_ReadOnly );
    if( poDataset == NULL ) { 
        cerr << "Error: Could not open raster data file" << endl;
        exit(1); 
    }

    fprintf(stderr, "Driver: %s/%s\n", poDataset->GetDriver()->GetDescription(), poDataset->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) );
    fprintf(stderr, "Size is %dx%dx%d\n", poDataset->GetRasterXSize(), poDataset->GetRasterYSize(), poDataset->GetRasterCount() );

    if( poDataset->GetProjectionRef()  != NULL ) cerr << "Projection is `" << poDataset->GetProjectionRef() << "'" << endl;;

    
    GDALRasterBand* poBand = poDataset->GetRasterBand( band_number );
    int nBlockXSize, nBlockYSize;
    poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
    fprintf(stderr, "Block=%dx%d Type=%s, ColorInterp=%s\n",
            nBlockXSize, nBlockYSize,
            GDALGetDataTypeName(poBand->GetRasterDataType()),
            GDALGetColorInterpretationName( poBand->GetColorInterpretation()) );

    Raster* raster = extract_raster_attributes( poDataset );
    raster->band = poBand;
    return raster;
}

void DemReader::MetaDataInfo(GDALDataset* gdalDataset, MapControllerPtr mapController)
{
	FileMetaData* header = mapController->GetMapModel()->GetMetaData();

	header->driverDesc = gdalDataset->GetDriver()->GetDescription();
	header->driverMetaData = gdalDataset->GetDriver()->GetMetadataItem(GDAL_DMD_LONGNAME);
	Log::Inst().Write("Driver: " + header->driverDesc + " \\ " + header->driverMetaData);

	header->xSize = gdalDataset->GetRasterXSize();
	header->ySize = gdalDataset->GetRasterYSize();
	header->bands = gdalDataset->GetRasterCount();
	Log::Inst().Write("Size (x,y): " + 
		StringTools::ToString(header->xSize) + ", " +
		StringTools::ToString(header->ySize));
	Log::Inst().Write("Bands: " + StringTools::ToString(header->bands));

	header->projection = std::string(gdalDataset->GetProjectionRef());
	Log::Inst().Write("Projection: " + header->projection);

	double adfGeoTransform[6];
	gdalDataset->GetGeoTransform( adfGeoTransform );

	header->originX = adfGeoTransform[0];
	header->originY = adfGeoTransform[3];
	Log::Inst().Write("Origin: " +
		StringTools::ToString(adfGeoTransform[0]) + ", " +
		StringTools::ToString(adfGeoTransform[3]));

	header->pixelSizeX = adfGeoTransform[1];
	header->pixelSizeY = adfGeoTransform[5];
	Log::Inst().Write("Pixel Size: " +
		StringTools::ToString(adfGeoTransform[1]) + ", " +
		StringTools::ToString(adfGeoTransform[5]));

	GDALRasterBand* gdalBand = gdalDataset->GetRasterBand(1);
	int nBlockXSize, nBlockYSize;
	gdalBand->GetBlockSize(&nBlockXSize, &nBlockYSize);

	header->dataType = std::string(GDALGetDataTypeName(gdalBand->GetRasterDataType()));
	Log::Inst().Write("Block, Type: " + StringTools::ToString(nBlockXSize) + ", " +
		StringTools::ToString(nBlockYSize) + ", " +
		std::string(header->dataType));

	header->colourInterpretation = std::string(GDALGetColorInterpretationName(gdalBand->GetColorInterpretation()));	
	Log::Inst().Write("Color Interpretation: " + header->colourInterpretation);

	header->extents.x = adfGeoTransform[0];
	header->extents.y = adfGeoTransform[3] + gdalBand->GetYSize()*adfGeoTransform[5];
	Log::Inst().Write("Lower, Left (x,y): " +
		StringTools::ToString(header->extents.x) + ", " +
		StringTools::ToString(header->extents.y));

	header->extents.dx = adfGeoTransform[0]+gdalBand->GetXSize()*adfGeoTransform[1];
	header->extents.dy = adfGeoTransform[3];
	Log::Inst().Write("Upper, Right (x,y): " +
		StringTools::ToString(header->extents.dx) + ", " +
		StringTools::ToString(header->extents.dy));

	Log::Inst().Write("");
}

bool getRawValuesFromFile(string fname,vector<vector<float>>& vecs)
{
   
   //vector<float> temp = vector<float>()
   GDALDataset *poDataset;
   GDALAllRegister();
   poDataset= (GDALDataset*) GDALOpen(fname.c_str(),GA_ReadOnly);
   if(poDataset == NULL)
   {
      cout << "OUCH!" << endl;
      return false;
   }
   cout << "Data size: " << GDALGetRasterXSize(poDataset) << " " << GDALGetRasterYSize(poDataset) << endl;

   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);
   int width = poBand->GetXSize(); 
   int height = poBand->GetYSize();
   int bands = poDataset->GetRasterCount();
   float *pafScanline;
   std::cout << "Before allocation" << adfMinMax[0] << " " << adfMinMax[1] << endl;
   int dsize = 256;
   pafScanline = (float *) CPLMalloc(sizeof(float)*width*height);
   vector<vector<float>> out = vector<vector<float>>(height,vector<float> (width,0));
   poBand->RasterIO(GF_Read,0,0,width,height,pafScanline,width,height,GDT_Float32,0,0);
   cout << "After allocation" << endl;
   for(int i = 0; i < height; i++)
   {
    for(int j = 0; j < width; j++)
    {
      //cout << i << j << endl << pafS;
      out[i][j] = pafScanline[i*width+j];
    }
   }
   CPLFree(pafScanline);
   //for(auto i : out)
   //for(auto j : i)
   //		cout << j << endl;
  cout << "After allocation" << endl;
  vecs = out;
   return true;
}

int main()
{
    GDALDataset  *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset *) GDALOpen( "GE01.tif", GA_ReadOnly );
    printf("Working! \n");

    if( poDataset != NULL ){
    	//Get 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] );
		}

		//Fetch 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() );

		//Close Dataset
		GDALClose(poDataset);
		//Exit      
		return 0;
	}


}

void bin2ascii(char *binFile, char *asciiFile){
    GDALDataset *layer;
    GDALAllRegister();
    layer = (GDALDataset *)GDALOpen(binFile, GA_ReadOnly);
    
    int bandNumber = 1;
    GDALRasterBand *band = layer->GetRasterBand(bandNumber);
    GDALDataType type    = band->GetRasterDataType();

    double ranges[6];
    layer->GetGeoTransform(ranges);

    ofstream outFile;
    outFile.open(asciiFile);
     
    outFile<<"ncols "<<layer->GetRasterXSize()<<"\n";
    outFile<<"nrows "<<layer->GetRasterYSize()<<"\n";
    outFile<<"xllcorner "<<ranges[0]<<"\n";
    outFile<<"yllcorner "<<(ranges[3] + layer->GetRasterXSize() * ranges[4] + layer->GetRasterYSize() * ranges[5])<<"\n";
    outFile<<"cellsize " <<ranges[1]<<"\n";
    outFile<<"nodata_value "<<"-9999"<<"\n";
    //getchar(); getchar();

    int cols = layer->GetRasterXSize();
    int rows = layer->GetRasterYSize();
    double NODATA_VAL = band->GetNoDataValue();

    cout<<"NODATA_VALUE= "<<NODATA_VAL<<"\n";

    int size = GDALGetDataTypeSize(type) / 8;
    void *data = CPLMalloc(size);
    //CPLErr err = band->RasterIO(GF_Read, 0, 0, cols, rows, data, cols, rows, type, 0, 0);
	//getchar(); getchar();
    //for(int j=0; j<rows*cols; j++){
	//int col, row;
    for(int row=0; row<rows; row++){
	for(int col=0; col<cols; col++){	
	 
	CPLErr err = band->RasterIO(GF_Read, col, row, 1, 1, data, 1, 1, type, 0, 0);
	double tempVal = readValueB2A(data, type, 0);
	outFile<<( tempVal != NODATA_VAL ? tempVal : -9999)<<" "; 
	
    	    //if((j+1)%cols == 0){
	//	cout<<"\n";
		//getchar(); getchar();
	//	}
	}
	outFile<<"\n"; 
	//getchar();
    }
	cout<<"Bin2Ascii.. Done!\n";
	outFile.close();
	//getchar(); getchar();
}

GDALDataType peekGDALType(const std::string &filename) {
  GDALAllRegister();
  GDALDataset *fin = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
  assert(fin!=NULL);

  GDALRasterBand *band   = fin->GetRasterBand(1);
  GDALDataType data_type = band->GetRasterDataType();

  GDALClose(fin);

  return data_type;
}

int
NBHeightMapper::loadTiff(const std::string& file) {
#ifdef HAVE_GDAL
    GDALAllRegister();
    GDALDataset* poDataset = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
    if (poDataset == 0) {
        WRITE_ERROR("Cannot load GeoTIFF file.");
        return 0;
    }
    const int xSize = poDataset->GetRasterXSize();
    const int ySize = poDataset->GetRasterYSize();
    double adfGeoTransform[6];
    if (poDataset->GetGeoTransform(adfGeoTransform) == CE_None) {
        Position topLeft(adfGeoTransform[0], adfGeoTransform[3]);
        mySizeOfPixel.set(adfGeoTransform[1], adfGeoTransform[5]);
        const double horizontalSize = xSize * mySizeOfPixel.x();
        const double verticalSize = ySize * mySizeOfPixel.y();
        myBoundary.add(topLeft);
        myBoundary.add(topLeft.x() + horizontalSize, topLeft.y() + verticalSize);
    } else {
        WRITE_ERROR("Could not parse geo information from " + file + ".");
        return 0;
    }
    const int picSize = xSize * ySize;
    myRaster = (int16_t*)CPLMalloc(sizeof(int16_t) * picSize);
    for (int i = 1; i <= poDataset->GetRasterCount(); i++) {
        GDALRasterBand* poBand = poDataset->GetRasterBand(i);
        if (poBand->GetColorInterpretation() != GCI_GrayIndex) {
            WRITE_ERROR("Unknown color band in " + file + ".");
            clearData();
            break;
        }
        if (poBand->GetRasterDataType() != GDT_Int16) {
            WRITE_ERROR("Unknown data type in " + file + ".");
            clearData();
            break;
        }
        assert(xSize == poBand->GetXSize() && ySize == poBand->GetYSize());
        if (poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, myRaster, xSize, ySize, GDT_Int16, 0, 0) == CE_Failure) {
            WRITE_ERROR("Failure in reading " + file + ".");
            clearData();
            break;
        }
    }
    GDALClose(poDataset);
    return picSize;
#else
    UNUSED_PARAMETER(file);
    WRITE_ERROR("Cannot load GeoTIFF file since SUMO was compiled without GDAL support.");
    return 0;
#endif
}

/**
  @brief  Determine data type of a GDAL file's first layer
  @author Richard Barnes ([email protected])

  @param[in]  filename   Filename of file whose type should be determined
*/
GDALDataType peekGDALType(const std::string &filename) {
  GDALAllRegister();
  GDALDataset *fin = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
  if(fin==NULL)
    throw std::runtime_error("Unable to open file '"+filename+"'!");

  GDALRasterBand *band   = fin->GetRasterBand(1);
  GDALDataType data_type = band->GetRasterDataType();

  GDALClose(fin);

  return data_type;
}

void HazePerfection::invertImage(GDALDataset *pDataset, GDALDataset *dstDataset, float a)
{
	GDALRasterBand *band = pDataset->GetRasterBand(1);
	GDALDataType dataType = band->GetRasterDataType();
	float *pixelData = new float[nXSize*nYSize];
	band->RasterIO(GF_Read, 0, 0, nXSize, nYSize, pixelData, nXSize, nYSize, GDT_Float32, 0, 0);
	
	GDALRasterBand *dstBand = dstDataset->GetRasterBand(1);
	for (int j = 0; j < nYSize; j++)
	{
		for (int i = 0; i < nXSize; i++)
		{
			pixelData[j*nXSize + i] = a - pixelData[j*nXSize + i];
		}
	}
	dstBand->RasterIO(GF_Write, 0, 0, nXSize, nYSize, pixelData, nXSize, nYSize, GDT_Float32, 0, 0);
	delete[]pixelData;
}

wxGISRasterLayer::wxGISRasterLayer(const wxString &sName, wxGISDataset* pwxGISDataset) : wxGISLayer(sName, pwxGISDataset)
{
    wxGISRasterDataset* pwxGISRasterDataset = wxDynamicCast(pwxGISDataset, wxGISRasterDataset);
	if(pwxGISRasterDataset)
	{
        if(m_sName.IsEmpty())
            m_sName = pwxGISRasterDataset->GetName();

        m_SpatialReference = pwxGISRasterDataset->GetSpatialReference();
		m_FullEnvelope = pwxGISRasterDataset->GetEnvelope();

        //TODO: load or get all renderers and check if i render can draw this dataset. If yes - set it as current
		if(pwxGISRasterDataset->GetBandCount() >= 3)
        {
            m_pRenderer = new wxGISRasterRGBARenderer(pwxGISDataset);
        }
        else
        {
		    GDALDataset* poGDALDataset = pwxGISRasterDataset->GetMainRaster();
		    if(!poGDALDataset)
			    poGDALDataset = pwxGISRasterDataset->GetRaster();
		    if(!poGDALDataset)
			    return;

		    GDALRasterBand* pBand = poGDALDataset->GetRasterBand(1);
            GDALColorInterp eColorInterpretation = pBand->GetColorInterpretation();
            if( eColorInterpretation == GCI_PaletteIndex )
            {
                m_pRenderer = new wxGISRasterRasterColormapRenderer(pwxGISDataset);
            }
            else if(pBand->GetRasterDataType() == GDT_Int32)
            {
                m_pRenderer = new wxGISRasterPackedRGBARenderer(pwxGISDataset);
            }
            else// if( eColorInterpretation == GCI_GrayIndex )
            {
    		    //TODO: else RasterStretchColorRampRenderer
                m_pRenderer = new wxGISRasterGreyScaleRenderer(pwxGISDataset);
            }
        }
	}
}

  void loadGDAL(const std::string &filename, int xOffset=0, int yOffset=0, int part_width=0, int part_height=0){
    assert(empty());
    assert(xOffset>=0);
    assert(yOffset>=0);

    GDALDataset *fin = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
    assert(fin!=NULL);

    GDALRasterBand *band = fin->GetRasterBand(1);
    auto data_type       = band->GetRasterDataType();

    total_width  = band->GetXSize();
    total_height = band->GetYSize();
    no_data      = band->GetNoDataValue();

    if(xOffset+part_width>=total_width)
      part_width  = total_width-xOffset;
    if(yOffset+part_height>=total_height)
      part_height = total_height-yOffset;

    if(part_width==0)
      part_width = total_width;
    view_width = part_width;

    if(part_height==0)
      part_height = total_height;
    view_height = part_height;

    view_xoff = xOffset;
    view_yoff = yOffset;

    std::cerr<<"Allocating: "<<view_height<<" rows by "<<view_width<<" columns"<<std::endl;
    data = InternalArray(view_height, Row(view_width));
    for(int y=yOffset;y<yOffset+view_height;y++)
      band->RasterIO( GF_Read, xOffset, y, view_width, 1, data[y-yOffset].data(), view_width, 1, data_type, 0, 0 );

    GDALClose(fin);
  }

GDALDatasetH CPL_STDCALL
GDALCreateWarpedVRT( GDALDatasetH hSrcDS, 
                     int nPixels, int nLines, double *padfGeoTransform,
                     GDALWarpOptions *psOptions )
    
{
    VALIDATE_POINTER1( hSrcDS, "GDALCreateWarpedVRT", NULL );

/* -------------------------------------------------------------------- */
/*      Create the VRTDataset and populate it with bands.               */
/* -------------------------------------------------------------------- */
    VRTWarpedDataset *poDS = new VRTWarpedDataset( nPixels, nLines );
    int i;

    psOptions->hDstDS = (GDALDatasetH) poDS;

    poDS->SetGeoTransform( padfGeoTransform );

    for( i = 0; i < psOptions->nBandCount; i++ )
    {
        VRTWarpedRasterBand *poBand;
        GDALRasterBand *poSrcBand = (GDALRasterBand *) 
            GDALGetRasterBand( hSrcDS, i+1 );

        poDS->AddBand( poSrcBand->GetRasterDataType(), NULL );

        poBand = (VRTWarpedRasterBand *) poDS->GetRasterBand( i+1 );
        poBand->CopyCommonInfoFrom( poSrcBand );
    }

/* -------------------------------------------------------------------- */
/*      Initialize the warp on the VRTWarpedDataset.                    */
/* -------------------------------------------------------------------- */
    poDS->Initialize( psOptions );
    
    return (GDALDatasetH) poDS;
}

/**
  @brief  Retrieve height, width, data type, and geotransform from a GDAL file
  @author Richard Barnes ([email protected])

  @param[in]  filename   GDAL file to peek at
  @param[out] height     Height of the raster in cells
  @param[out] width      Width of the raster in cells
  @param[out] dtype      Data type of the file in question
  @param[out] geo_trans  Returns the SIX elements of the raster's geotransform
*/
void getGDALDimensions(
  const   std::string &filename,
  int32_t &height,
  int32_t &width,
  GDALDataType &dtype,
  double geotransform[6]
){
  GDALAllRegister();
  GDALDataset *fin = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
  if(fin==NULL)
    throw std::runtime_error("Could not open file '"+filename+"' to get dimensions.");

  GDALRasterBand *band = fin->GetRasterBand(1);
  
  dtype = band->GetRasterDataType();

  if(geotransform!=NULL && fin->GetGeoTransform(geotransform)!=CE_None)
    throw std::runtime_error("Error getting geotransform from '"+filename+"'.");

  height  = band->GetYSize();
  width   = band->GetXSize();

  GDALClose(fin);
}