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

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;
	}


}

SEXP
RGDAL_GetPaletteInterp(SEXP sxpRasterBand) {

  GDALRasterBand *pRasterBand = getGDALRasterPtr(sxpRasterBand);

  GDALPaletteInterp ePI = 
    pRasterBand->GetColorTable()->GetPaletteInterpretation();

  return(mkString_safe(GDALGetPaletteInterpretationName(ePI)));

}

GDALColorTable *GDALProxyPoolRasterBand::GetColorTable()
{
    GDALRasterBand* poUnderlyingRasterBand = RefUnderlyingRasterBand();
    if (poUnderlyingRasterBand == NULL)
        return NULL;

    if (poColorTable)
        delete poColorTable;
    poColorTable = NULL;

    GDALColorTable* poUnderlyingColorTable = poUnderlyingRasterBand->GetColorTable();
    if (poUnderlyingColorTable)
        poColorTable = poUnderlyingColorTable->Clone();

    UnrefUnderlyingRasterBand(poUnderlyingRasterBand);

    return poColorTable;
}

static GDALDataset *
BSBCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
               int bStrict, char ** papszOptions, 
               GDALProgressFunc pfnProgress, void * pProgressData )

{
    int  nBands = poSrcDS->GetRasterCount();
    int  nXSize = poSrcDS->GetRasterXSize();
    int  nYSize = poSrcDS->GetRasterYSize();

/* -------------------------------------------------------------------- */
/*      Some some rudimentary checks                                    */
/* -------------------------------------------------------------------- */
    if( nBands != 1 )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "BSB driver only supports one band images.\n" );

        return NULL;
    }

    if( poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte 
        && bStrict )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "BSB driver doesn't support data type %s. "
                  "Only eight bit bands supported.\n", 
                  GDALGetDataTypeName( 
                      poSrcDS->GetRasterBand(1)->GetRasterDataType()) );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Open the output file.                                           */
/* -------------------------------------------------------------------- */
    BSBInfo *psBSB;

    psBSB = BSBCreate( pszFilename, 0, 200, nXSize, nYSize );
    if( psBSB == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Prepare initial color table.colortable.                         */
/* -------------------------------------------------------------------- */
    GDALRasterBand	*poBand = poSrcDS->GetRasterBand(1);
    int			iColor;
    unsigned char       abyPCT[771];
    int                 nPCTSize;
    int                 anRemap[256];

    abyPCT[0] = 0;
    abyPCT[1] = 0;
    abyPCT[2] = 0;

    if( poBand->GetColorTable() == NULL )
    {
        /* map greyscale down to 63 grey levels. */
        for( iColor = 0; iColor < 256; iColor++ )
        {
            int nOutValue = (int) (iColor / 4.1) + 1;

            anRemap[iColor] = nOutValue;
            abyPCT[nOutValue*3 + 0] = (unsigned char) iColor;
            abyPCT[nOutValue*3 + 1] = (unsigned char) iColor;
            abyPCT[nOutValue*3 + 2] = (unsigned char) iColor;
        }
        nPCTSize = 64;
    }
    else
    {
        GDALColorTable	*poCT = poBand->GetColorTable();
        int nColorTableSize = poCT->GetColorEntryCount();
        if (nColorTableSize > 255)
            nColorTableSize = 255;

        for( iColor = 0; iColor < nColorTableSize; iColor++ )
        {
            GDALColorEntry	sEntry;

            poCT->GetColorEntryAsRGB( iColor, &sEntry );

            anRemap[iColor] = iColor + 1;
            abyPCT[(iColor+1)*3 + 0] = (unsigned char) sEntry.c1;
            abyPCT[(iColor+1)*3 + 1] = (unsigned char) sEntry.c2;
            abyPCT[(iColor+1)*3 + 2] = (unsigned char) sEntry.c3;
        }

        nPCTSize = nColorTableSize + 1;

        // Add entries for pixel values which apparently will not occur.
        for( iColor = nPCTSize; iColor < 256; iColor++ )
            anRemap[iColor] = 1;
    }

/* -------------------------------------------------------------------- */
/*      Boil out all duplicate entries.                                 */
/* -------------------------------------------------------------------- */
    int  i;

//.........这里部分代码省略.........

//.........这里部分代码省略.........

		//fetch raster bands
		bool zRasterProcessed = false;
		unsigned zInvalid = 0;
		double zMinMax[2] = {0, 0};

		for (int i=1; i<=rasterCount; ++i)
		{
			ccLog::Print( "Reading band #%i", i);
			GDALRasterBand* poBand = poDataset->GetRasterBand(i);

			GDALColorInterp colorInterp = poBand->GetColorInterpretation();
			GDALDataType bandType = poBand->GetRasterDataType();

			int nBlockXSize, nBlockYSize;
			poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );
			ccLog::Print( "Block=%dx%d Type=%s, ColorInterp=%s", nBlockXSize, nBlockYSize, GDALGetDataTypeName(poBand->GetRasterDataType()), GDALGetColorInterpretationName(colorInterp) );

			//fetching raster scan-line
			int nXSize = poBand->GetXSize();
			int nYSize = poBand->GetYSize();
			assert(nXSize == rasterX);
			assert(nYSize == rasterY);
			
			int bGotMin, bGotMax;
			double adfMinMax[2] = {0, 0};
			adfMinMax[0] = poBand->GetMinimum( &bGotMin );
			adfMinMax[1] = poBand->GetMaximum( &bGotMax );
			if (!bGotMin || !bGotMax )
				//DGM FIXME: if the file is corrupted (e.g. ASCII ArcGrid with missing rows) this method will enter in a infinite loop!
				GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);
			ccLog::Print( "Min=%.3fd, Max=%.3f", adfMinMax[0], adfMinMax[1] );

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

			if( poBand->GetOverviewCount() > 0 )
				printf( "Band has %d overviews", poBand->GetOverviewCount() );

			if (colorInterp == GCI_Undefined && !zRasterProcessed/*&& !colTable*/) //probably heights?
			{
				zRasterProcessed = true;
				zMinMax[0] = adfMinMax[0];
				zMinMax[1] = adfMinMax[1];

				double* scanline = (double*) CPLMalloc(sizeof(double)*nXSize);
				//double* scanline = new double[nXSize];
				memset(scanline,0,sizeof(double)*nXSize);

				for (int j=0; j<nYSize; ++j)
				{
					if (poBand->RasterIO( GF_Read, /*xOffset=*/0, /*yOffset=*/j, /*xSize=*/nXSize, /*ySize=*/1, /*buffer=*/scanline, /*bufferSizeX=*/nXSize, /*bufferSizeY=*/1, /*bufferType=*/GDT_Float64, /*x_offset=*/0, /*y_offset=*/0 ) != CE_None)
					{
						delete pc;
						CPLFree(scanline);
						GDALClose(poDataset);
						return CC_FERR_READING;
					}

					for (int k=0; k<nXSize; ++k)
					{
						double z = static_cast<double>(scanline[k]) + Pshift[2];
						unsigned pointIndex = static_cast<unsigned>(k + j * rasterX);
						if (pointIndex <= pc->size())
						{

//.........这里部分代码省略.........
                {
                    GDALRasterBand * band = dataset_.GetRasterBand(i + 1);
#ifdef MAPNIK_LOG
                    get_overview_meta(band);
#endif
                    GDALColorInterp color_interp = band->GetColorInterpretation();
                    switch (color_interp)
                    {
                    case GCI_RedBand:
                        red = band;
                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found red band";
                        break;
                    case GCI_GreenBand:
                        green = band;
                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found green band";
                        break;
                    case GCI_BlueBand:
                        blue = band;
                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found blue band";
                        break;
                    case GCI_AlphaBand:
                        alpha = band;
                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found alpha band";
                        break;
                    case GCI_GrayIndex:
                        grey = band;
                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found gray band";
                        break;
                    case GCI_PaletteIndex:
                    {
                        grey = band;
#ifdef MAPNIK_LOG
                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found gray band, and colortable...";
                        GDALColorTable *color_table = band->GetColorTable();

                        if (color_table)
                        {
                            int count = color_table->GetColorEntryCount();
                            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Color Table count=" << count;

                            for (int j = 0; j < count; j++)
                            {
                                const GDALColorEntry *ce = color_table->GetColorEntry (j);
                                if (! ce) continue;
                                MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Color entry RGB=" << ce->c1 << "," <<ce->c2 << "," << ce->c3;
                            }
                        }
#endif
                        break;
                    }
                    case GCI_Undefined:
#if GDAL_VERSION_NUM <= 1730
                        if (nbands_ == 4)
                        {
                            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming alpha band)";
                            alpha = band;
                        }
                        else
                        {
                            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)";
                            grey = band;
                        }
#else
                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found undefined band (assumming gray band)";
                        grey = band;
#endif

//.........这里部分代码省略.........
    // --------------------------------------------------------------------

    GDALRasterBand *poSrcBand;
    GDALRasterBand *poDstBand;
    double dfMin;
    double dfMax;
    double dfMean;
    double dfStdDev = -1;

    for( int i = 1; i <= poDstDS->nBands; i++)
    {
        delete poDstDS->GetRasterBand(i);
    }
    poDstDS->nBands = 0;

    if( poDstDS->hHeaderOne.DataTypeCode == Uncompressed24bit )
    {
        poDstDS->SetBand( 1, new IntergraphRGBBand( poDstDS, 1, 0, 3 ) );
        poDstDS->SetBand( 2, new IntergraphRGBBand( poDstDS, 2, 0, 2 ) );
        poDstDS->SetBand( 3, new IntergraphRGBBand( poDstDS, 3, 0, 1 ) );
        poDstDS->nBands = 3;
    }
    else
    {
        for( int i = 1; i <= poSrcDS->GetRasterCount(); i++ )
        {
            poSrcBand = poSrcDS->GetRasterBand(i);
            eType = poSrcDS->GetRasterBand(i)->GetRasterDataType();

            poDstBand = new IntergraphRasterBand( poDstDS, i, 0, eType );
            poDstDS->SetBand( i, poDstBand );

            poDstBand->SetCategoryNames( poSrcBand->GetCategoryNames() );
            poDstBand->SetColorTable( poSrcBand->GetColorTable() );
            poSrcBand->GetStatistics( false, true, &dfMin, &dfMax, &dfMean, &dfStdDev );
            poDstBand->SetStatistics( dfMin, dfMax, dfMean, dfStdDev );
        }
    }

    // --------------------------------------------------------------------
    // Copy image data
    // --------------------------------------------------------------------

    int nXSize = poDstDS->GetRasterXSize();
    int nYSize = poDstDS->GetRasterYSize();

    int nBlockXSize;
    int nBlockYSize;

    CPLErr eErr = CE_None;

    for( int iBand = 1; iBand <= poSrcDS->GetRasterCount(); iBand++ )
    {
        GDALRasterBand *poDstBand = poDstDS->GetRasterBand( iBand );
        GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( iBand );

        // ------------------------------------------------------------
        // Copy Untiled / Uncompressed
        // ------------------------------------------------------------

        int   iYOffset, iXOffset;
        void *pData;

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

        nBlockXSize = nXSize;

//.........这里部分代码省略.........
/* -------------------------------------------------------------------- */
/*      Geotransform.                                                   */
/* -------------------------------------------------------------------- */
    if( poAuxDS->GetGeoTransform( psPam->adfGeoTransform ) == CE_None )
        psPam->bHaveGeoTransform = TRUE;

/* -------------------------------------------------------------------- */
/*      GCPs                                                            */
/* -------------------------------------------------------------------- */
    if( poAuxDS->GetGCPCount() > 0 )
    {
        psPam->nGCPCount = poAuxDS->GetGCPCount();
        psPam->pasGCPList = GDALDuplicateGCPs( psPam->nGCPCount, 
                                               poAuxDS->GetGCPs() );
    }

/* -------------------------------------------------------------------- */
/*      Apply metadata. We likely ought to be merging this in rather    */
/*      than overwriting everything that was there.                     */
/* -------------------------------------------------------------------- */
    char **papszMD = poAuxDS->GetMetadata();
    if( CSLCount(papszMD) > 0 )
    {
        char **papszMerged = 
            CSLMerge( CSLDuplicate(GetMetadata()), papszMD );
        GDALPamDataset::SetMetadata( papszMerged );
        CSLDestroy( papszMerged );
    }

    papszMD = poAuxDS->GetMetadata("XFORMS");
    if( CSLCount(papszMD) > 0 )
    {
        char **papszMerged = 
            CSLMerge( CSLDuplicate(GetMetadata("XFORMS")), papszMD );
        GDALPamDataset::SetMetadata( papszMerged, "XFORMS" );
        CSLDestroy( papszMerged );
    }

/* ==================================================================== */
/*      Process bands.                                                  */
/* ==================================================================== */
    int iBand;

    for( iBand = 0; iBand < poAuxDS->GetRasterCount(); iBand++ )
    {
        if( iBand >= GetRasterCount() )
            break;

        GDALRasterBand *poAuxBand = poAuxDS->GetRasterBand( iBand+1 );
        GDALRasterBand *poBand = GetRasterBand( iBand+1 );

        papszMD = poAuxBand->GetMetadata();
        if( CSLCount(papszMD) > 0 )
        {
            char **papszMerged = 
                CSLMerge( CSLDuplicate(poBand->GetMetadata()), papszMD );
            poBand->SetMetadata( papszMerged );
            CSLDestroy( papszMerged );
        }

        if( poAuxBand->GetCategoryNames() != NULL )
            poBand->SetCategoryNames( poAuxBand->GetCategoryNames() );

        if( poAuxBand->GetColorTable() != NULL 
            && poBand->GetColorTable() == NULL )
            poBand->SetColorTable( poAuxBand->GetColorTable() );

        // histograms?
        double dfMin, dfMax;
        int nBuckets, *panHistogram=NULL;

        if( poAuxBand->GetDefaultHistogram( &dfMin, &dfMax, 
                                            &nBuckets, &panHistogram,
                                            FALSE, NULL, NULL ) == CE_None )
        {
            poBand->SetDefaultHistogram( dfMin, dfMax, nBuckets, 
                                         panHistogram );
            CPLFree( panHistogram );
        }

        // RAT 
        if( poAuxBand->GetDefaultRAT() != NULL )
            poBand->SetDefaultRAT( poAuxBand->GetDefaultRAT() );

        // NoData
        int bSuccess = FALSE;
        double dfNoDataValue = poAuxBand->GetNoDataValue( &bSuccess );
        if( bSuccess )
            poBand->SetNoDataValue( dfNoDataValue );
    }

    GDALClose( poAuxDS );
    
/* -------------------------------------------------------------------- */
/*      Mark PAM info as clean.                                         */
/* -------------------------------------------------------------------- */
    nPamFlags &= ~GPF_DIRTY;

    return CE_Failure;
}

CMapRaster::CMapRaster(const QString& fn, CCanvas * parent)
: IMap(eRaster, "",parent)
, x(0)
, y(0)
, zoomlevel(1)
, zoomfactor(1.0)
, rasterBandCount(0)
{
    filename = fn;
#ifdef WIN32
    dataset = (GDALDataset*)GDALOpen(filename.toLocal8Bit(),GA_ReadOnly);
#else
    dataset = (GDALDataset*)GDALOpen(filename.toUtf8(),GA_ReadOnly);
#endif
    if(dataset == 0)
    {
        QMessageBox::warning(0, tr("Error..."), tr("Failed to load file: %1").arg(filename));
        return;
    }

    rasterBandCount = dataset->GetRasterCount();
    if(rasterBandCount == 1)
    {
        GDALRasterBand * pBand;
        pBand = dataset->GetRasterBand(1);
        if(pBand == 0)
        {
            delete dataset; dataset = 0;
            QMessageBox::warning(0, tr("Error..."), tr("Failed to load file: %1").arg(filename));
            return;
        }

        if(pBand->GetColorInterpretation() !=  GCI_PaletteIndex && pBand->GetColorInterpretation() !=  GCI_GrayIndex)
        {
            delete dataset; dataset = 0;
            QMessageBox::warning(0, tr("Error..."), tr("File must be 8 bit palette or gray indexed."));
            return;
        }

        if(pBand->GetColorInterpretation() ==  GCI_PaletteIndex )
        {
            GDALColorTable * pct = pBand->GetColorTable();
            for(int i=0; i < pct->GetColorEntryCount(); ++i)
            {
                const GDALColorEntry& e = *pct->GetColorEntry(i);
                colortable << qRgba(e.c1, e.c2, e.c3, e.c4);
            }
        }
        else if(pBand->GetColorInterpretation() ==  GCI_GrayIndex )
        {
            for(int i=0; i < 256; ++i)
            {
                colortable << qRgba(i, i, i, 255);
            }
        }
        else
        {
            delete dataset; dataset = 0;
            QMessageBox::warning(0, tr("Error..."), tr("File must be 8 bit palette or gray indexed."));
            return;
        }

        int success = 0;
        double idx = pBand->GetNoDataValue(&success);

        if(success)
        {
            QColor tmp(colortable[idx]);
            tmp.setAlpha(0);
            colortable[idx] = tmp.rgba();
        }
    }
    maparea.setWidth(dataset->GetRasterXSize());
    maparea.setHeight(dataset->GetRasterYSize());

}

//.........这里部分代码省略.........
			  NULL, pProgressData ) )
	{
	    VSIFCloseL( fp );
	    VSIFree( pfData );
	    CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
	    return NULL;
	}
    }

    VSIFree( pfData );

    /* write out the min and max values */
    eErr = WriteHeader( fp, nXSize, nYSize,
			dfMinX, dfMaxX, dfMinY, dfMaxY, dfMinZ, dfMaxZ );

    if( eErr != CE_None )
    {
	VSIFCloseL( fp );
        return NULL;
    }

    VSIFCloseL( fp );

    GDALPamDataset *poDstDS = (GDALPamDataset *)GDALOpen( pszFilename,
							  GA_Update );
    if( poDstDS == NULL )
    {
	VSIUnlink( pszFilename );
	CPLError( CE_Failure, CPLE_FileIO,
		  "Unable to open copy of dataset.\n" );
	return NULL;
    }
    else if( dynamic_cast<GSBGDataset *>(poDstDS) == NULL )
    {
	VSIUnlink( pszFilename );
	delete poDstDS;
	CPLError( CE_Failure, CPLE_FileIO,
		  "Copy dataset not opened as Golden Surfer Binary Grid!?\n" );
	return NULL;
    }

    GDALRasterBand *poDstBand = poSrcDS->GetRasterBand(1);
    if( poDstBand == NULL )
    {
	VSIUnlink( pszFilename );
	delete poDstDS;
	CPLError( CE_Failure, CPLE_FileIO,
		  "Unable to open copy of raster band?\n" );
	return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Attempt to copy metadata.					*/
/* -------------------------------------------------------------------- */
    if( !bStrict )
	CPLPushErrorHandler( CPLQuietErrorHandler );

    /* non-zero transform 2 or 4 or negative 1 or 5  not supported natively */
    /*if( adfGeoTransform[2] != 0.0 || adfGeoTransform[4] != 0.0
	|| adfGeoTransform[1] < 0.0 || adfGeoTransform[5] < 0.0 )
	poDstDS->GDALPamDataset::SetGeoTransform( adfGeoTransform );*/

    const char *szProjectionRef = poSrcDS->GetProjectionRef();
    if( *szProjectionRef != '\0' )
	poDstDS->SetProjection( szProjectionRef );

    char **pszMetadata = poSrcDS->GetMetadata();
    if( pszMetadata != NULL )
	poDstDS->SetMetadata( pszMetadata );

    /* FIXME:  Should the dataset description be copied as well, or is it
     *         always the file name? */
    poDstBand->SetDescription( poSrcBand->GetDescription() );

    int bSuccess;
    double dfOffset = poSrcBand->GetOffset( &bSuccess );
    if( bSuccess && dfOffset != 0.0 )
	poDstBand->SetOffset( dfOffset );

    double dfScale = poSrcBand->GetScale( &bSuccess );
    if( bSuccess && dfScale != 1.0 )
	poDstBand->SetScale( dfScale );

    GDALColorInterp oColorInterp = poSrcBand->GetColorInterpretation();
    if( oColorInterp != GCI_Undefined )
        poDstBand->SetColorInterpretation( oColorInterp );

    char **pszCatNames = poSrcBand->GetCategoryNames();
    if( pszCatNames != NULL)
	poDstBand->SetCategoryNames( pszCatNames );

    GDALColorTable *poColorTable = poSrcBand->GetColorTable();
    if( poColorTable != NULL )
	poDstBand->SetColorTable( poColorTable );

    if( !bStrict )
	CPLPopErrorHandler();

    return poDstDS;
}

static GDALDataset *
XPMCreateCopy( const char * pszFilename, GDALDataset *poSrcDS, 
               int bStrict, char ** papszOptions, 
               GDALProgressFunc pfnProgress, void * pProgressData )

{
    int  nBands = poSrcDS->GetRasterCount();
    int  nXSize = poSrcDS->GetRasterXSize();
    int  nYSize = poSrcDS->GetRasterYSize();
    GDALColorTable *poCT;

/* -------------------------------------------------------------------- */
/*      Some some rudimentary checks                                    */
/* -------------------------------------------------------------------- */
    if( nBands != 1 )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "XPM driver only supports one band images.\n" );

        return NULL;
    }

    if( poSrcDS->GetRasterBand(1)->GetRasterDataType() != GDT_Byte 
        && bStrict )
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "XPM driver doesn't support data type %s. "
                  "Only eight bit bands supported.\n", 
                  GDALGetDataTypeName( 
                      poSrcDS->GetRasterBand(1)->GetRasterDataType()) );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      If there is no colortable on the source image, create a         */
/*      greyscale one with 64 levels of grey.                           */
/* -------------------------------------------------------------------- */
    GDALRasterBand	*poBand = poSrcDS->GetRasterBand(1);
    int                 i;
    GDALColorTable      oGreyTable;

    poCT = poBand->GetColorTable();
    if( poCT == NULL )
    {
        poCT = &oGreyTable;

        for( i = 0; i < 256; i++ )
        {
            GDALColorEntry sColor;

            sColor.c1 = (short) i;
            sColor.c2 = (short) i;
            sColor.c3 = (short) i;
            sColor.c4 = 255;

            poCT->SetColorEntry( i, &sColor );
        }
    }

/* -------------------------------------------------------------------- */
/*      Build list of active colors, and the mapping from pixels to     */
/*      our active colormap.                                            */
/* -------------------------------------------------------------------- */
    const char *pszColorCodes = " [email protected]#$%^&*()-+=[]|:;,.<>?/";

    int  anPixelMapping[256];
    GDALColorEntry asPixelColor[256];
    int  nActiveColors = MIN(poCT->GetColorEntryCount(),256);

    // Setup initial colortable and pixel value mapping. 
    memset( anPixelMapping+0, 0, sizeof(int) * 256 );
    for( i = 0; i < nActiveColors; i++ )
    {
        poCT->GetColorEntryAsRGB( i, asPixelColor + i );
        anPixelMapping[i] = i;
    }

/* ==================================================================== */
/*      Iterate merging colors until we are under our limit (about 85). */
/* ==================================================================== */
    while( nActiveColors > (int) strlen(pszColorCodes) )
    {
        int nClosestDistance = 768;
        int iClose1 = -1, iClose2 = -1;
        int iColor1, iColor2;

        // Find the closest pair of colors. 
        for( iColor1 = 0; iColor1 < nActiveColors; iColor1++ )
        {
            for( iColor2 = iColor1+1; iColor2 < nActiveColors; iColor2++ )
            {
                int	nDistance;

                if( asPixelColor[iColor1].c4 < 128 
                    && asPixelColor[iColor2].c4 < 128 )
                    nDistance = 0;
                else
                    nDistance = 
                        ABS(asPixelColor[iColor1].c1-asPixelColor[iColor2].c1)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            nBandBits = 128;
            nBandFormat = SAMPLEFORMAT_COMPLEXIEEEFP;
            break;

          default:
            CPLAssert( false );
            return CE_Failure;
        }

        if( hBand->GetMetadataItem( "NBITS", "IMAGE_STRUCTURE" ) )
        {
            nBandBits =
                atoi(hBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"));

            if( nBandBits == 1
                && STARTS_WITH_CI(pszResampling, "AVERAGE_BIT2") )
                nBandBits = 8;
        }

        if( iBand == 0 )
        {
            nBitsPerPixel = nBandBits;
            nSampleFormat = nBandFormat;
            nXSize = hBand->GetXSize();
            nYSize = hBand->GetYSize();
        }
        else if( nBitsPerPixel != nBandBits || nSampleFormat != nBandFormat )
        {
            CPLError( CE_Failure, CPLE_NotSupported,
                      "GTIFFBuildOverviews() doesn't support a mixture of band"
                      " data types." );
            return CE_Failure;
        }
        else if( hBand->GetColorTable() != nullptr )
        {
            CPLError( CE_Failure, CPLE_NotSupported,
                      "GTIFFBuildOverviews() doesn't support building"
                      " overviews of multiple colormapped bands." );
            return CE_Failure;
        }
        else if( hBand->GetXSize() != nXSize
                 || hBand->GetYSize() != nYSize )
        {
            CPLError( CE_Failure, CPLE_NotSupported,
                      "GTIFFBuildOverviews() doesn't support building"
                      " overviews of different sized bands." );
            return CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      Use specified compression method.                               */
/* -------------------------------------------------------------------- */
    const char *pszCompress = CPLGetConfigOption( "COMPRESS_OVERVIEW", nullptr );

    if( pszCompress != nullptr && pszCompress[0] != '\0' )
    {
        nCompression =
            GTIFFGetCompressionMethod(pszCompress, "COMPRESS_OVERVIEW");
        if( nCompression < 0 )
            return CE_Failure;
    }

    if( nCompression == COMPRESSION_JPEG && nBitsPerPixel > 8 )
    {
        if( nBitsPerPixel > 16 )

bool vtImageGeo::ReadTIF(const char *filename, bool progress_callback(int))
{
	// Use GDAL to read a TIF file (or any other format that GDAL is
	//  configured to read) into this OSG image.
	bool bRet = true;
	vtString message;

	setFileName(filename);

	g_GDALWrapper.RequestGDALFormats();

	GDALDataset *pDataset = NULL;
	GDALRasterBand *pBand;
	GDALRasterBand *pRed = NULL;
	GDALRasterBand *pGreen = NULL;
	GDALRasterBand *pBlue = NULL;
	GDALRasterBand *pAlpha = NULL;
	GDALColorTable *pTable;
	uchar *pScanline = NULL;
	uchar *pRedline = NULL;
	uchar *pGreenline = NULL;
	uchar *pBlueline = NULL;
	uchar *pAlphaline = NULL;

	CPLErr Err;
	bool bColorPalette = false;
	int iXSize, iYSize;
	int nxBlocks, nyBlocks;
	int xBlockSize, yBlockSize;

	try
	{
		pDataset = (GDALDataset *) GDALOpen(filename, GA_ReadOnly);
		if(pDataset == NULL )
			throw "Couldn't open that file.";

		// Get size
		iXSize = pDataset->GetRasterXSize();
		iYSize = pDataset->GetRasterYSize();

		// Try getting CRS
		vtProjection temp;
		bool bHaveProj = false;
		const char *pProjectionString = pDataset->GetProjectionRef();
		if (pProjectionString)
		{
			OGRErr err = temp.importFromWkt((char**)&pProjectionString);
			if (err == OGRERR_NONE)
			{
				m_proj = temp;
				bHaveProj = true;
			}
		}
		if (!bHaveProj)
		{
			// check for existence of .prj file
			bool bSuccess = temp.ReadProjFile(filename);
			if (bSuccess)
			{
				m_proj = temp;
				bHaveProj = true;
			}
		}

		// Try getting extents
		double affineTransform[6];
		if (pDataset->GetGeoTransform(affineTransform) == CE_None)
		{
			m_extents.left = affineTransform[0];
			m_extents.right = m_extents.left + affineTransform[1] * iXSize;
			m_extents.top = affineTransform[3];
			m_extents.bottom = m_extents.top + affineTransform[5] * iYSize;
		}

		// Raster count should be 3 for colour images (assume RGB)
		int iRasterCount = pDataset->GetRasterCount();

		if (iRasterCount != 1 && iRasterCount != 3 && iRasterCount != 4)
		{
			message.Format("Image has %d bands (not 1, 3, or 4).", iRasterCount);
			throw (const char *)message;
		}

		if (iRasterCount == 1)
		{
			pBand = pDataset->GetRasterBand(1);

			// Check the band's data type
			GDALDataType dtype = pBand->GetRasterDataType();
			if (dtype != GDT_Byte)
			{
				message.Format("Band is of type %s, but we support type Byte.", GDALGetDataTypeName(dtype));
				throw (const char *)message;
			}
			GDALColorInterp ci = pBand->GetColorInterpretation();

			if (ci == GCI_PaletteIndex)
			{
				if (NULL == (pTable = pBand->GetColorTable()))
					throw "Couldn't get color table.";
//.........这里部分代码省略.........

int readRaster(const char* pszFilename)
{
    // Set config options for GDAL (needs >= 2.0). Setting GTIFF_VIRTUAL_MEM_IO to "YES" can cause things bleed to
    // swap if enough RAM is not available. Use "IF_ENOUGH_RAM" for safer performance if unsure. NOTE that faster
    // mem io only works with *uncompressed* GeoTIFFs.
    // New in GDAL 2.0, from https://2015.foss4g-na.org/sites/default/files/slides/GDAL%202.0%20overview.pdf
    // GeoTIFF driver (with i7-4700 HQ (8 vCPUs)):with i7-4700 HQ (8 vCPUs)
    //   - Default: time ./testblockcache -ondisk: 7.5s
    //   - GTIFF_DIRECT_IO=YES: short circuit the block cache&libtiff for most RasterIO() operations (restricted to
    //     uncompressed stripped GeoTIFF). time ./testblockcache -ondisk: 2s
    //   - GTIFF_VIRTUAL_MEM_IO=YES: same as above, with tiled GeoTIFF as well. Uses memory-mapped file access. Linux
    //     only for now, 64bit recommended (could be extended to other platforms possible).
    //     time ./testblockcache -ondisk: 0.3s
    //
    CPLSetConfigOption("GTIFF_VIRTUAL_MEM_IO", "YES" );
    // Open the dataset
    GDALDataset *poDataset;
    GDALAllRegister();
    poDataset = (GDALDataset *) GDALOpen( pszFilename, GA_ReadOnly );

    if( poDataset == NULL )
    {
        std::cout << "Dataset " << pszFilename << " could not be opened." << std::endl;
        return -1;
    }
    else
    {
        GDALRasterBand  *poBand;
        int             nBlockXSize, nBlockYSize;
        int             bGotMin, bGotMax;
        double          adfMinMax[2];

        //  Get raster band and its size
        poBand = poDataset->GetRasterBand(1);
        poBand->GetBlockSize( &nBlockXSize, &nBlockYSize);

        std::cout << "Dataset: " << pszFilename << std::endl;

        std::cout << "Block=" << nBlockXSize << "x" << nBlockYSize << " Type=" <<
                GDALGetDataTypeName(poBand->GetRasterDataType()) << " ColorInterp=" <<
                GDALGetColorInterpretationName(poBand->GetColorInterpretation()) << std::endl;

        // Calculate some stats
        adfMinMax[0] = poBand->GetMinimum(&bGotMin);
        adfMinMax[1] = poBand->GetMaximum(&bGotMax);

        if(!(bGotMin && bGotMax)) {
            GDALComputeRasterMinMax((GDALRasterBandH) poBand, TRUE, adfMinMax);
        }
        std::cout << "Min=" << adfMinMax[0] << " Max=" << adfMinMax[1] << std::endl;

        if(poBand->GetOverviewCount() > 0) {
            std::cout << "Band has " << poBand->GetOverviewCount() << " overviews." << std::endl;
        }
        if( poBand->GetColorTable() != NULL ) {
            std::cout << "Band has a color table with " << poBand->GetColorTable()->GetColorEntryCount() <<
                    " entries." << std::endl;
        }

        // Get the actual data
        float *pafScanline;
        int   nXSize = poBand->GetXSize();
        pafScanline = (float *) CPLMalloc(sizeof(float) * nXSize);

        // RasterIO has a new argument psExtraArg in GDAL > 2.0. NOTE: that  GDALRasterBand::ReadBlock() probably has
        // better performance for reading the whole data at one go.
        #ifdef USE_GDAL_2
            GDALRasterIOExtraArg* arg = NULL;
            poBand->RasterIO(GF_Read, 0, 0, nXSize, 1, pafScanline, nXSize, 1, GDT_Float3GDALRasterBand::ReadBlock 2,
               0, 0, arg);
        #else
            poBand->RasterIO(GF_Read, 0, 0, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
        #endif

        // ... do something with the data ...

        // Free resources
        CPLFree(pafScanline);
        GDALClose((GDALDatasetH) poDataset);
        std::cout << std::endl;
    }
    return 0;
}

 virtual GDALColorTable* GetColorTable()
 {
     return poUnderlyingBand->GetColorTable();
 }