自学教程：C++ GDALGetRasterBand函数代码示例

    void object::test<6>()    {        // Index of test file being tested        const std::size_t fileIdx = 1;        std::string file(data_ + SEP);        file += grids_.at(fileIdx).file_;        GDALDatasetH ds = GDALOpen(file.c_str(), GA_ReadOnly);        ensure("Can't open dataset: " + file, nullptr != ds);        GDALRasterBandH band = GDALGetRasterBand(ds, grids_.at(fileIdx).band_);        ensure("Can't get raster band", nullptr != band);        const double noData = GDALGetRasterNoDataValue(band, nullptr);        ensure_equals("Grid NODATA value wrong or missing", noData, -99999);        ensure_equals("Data type is not GDT_Float32", GDALGetRasterDataType(band), GDT_Float32);        GDALClose(ds);    }

void generic_read(struct FI *f, char *filename){	if (filename_corresponds_to_tiffo(filename)) {#ifdef FANCY_TIFF		f->tiffo = true;		tiff_octaves_init0(f->t, filename, f->megabytes,f->max_octaves);		if (f->option_write) f->t->option_write = true;		f->w = f->t->i->w;		f->h = f->t->i->h;		f->pd = f->t->i->spp;		f->no = f->t->noctaves;#else		assert(false);#endif	} else if (!f->option_write && FORCE_GDAL()) {#ifdef FANCY_GDAL		f->gdal = true;		GDALAllRegister();		char buf[2*FILENAME_MAX];		snprintf(buf, 2*FILENAME_MAX, has_prefix(filename, "http://") ||				has_prefix(filename, "https://") ?				"/vsicurl/%s" : "%s", filename);		f->gdal_img = GDALOpen(buf, GA_ReadOnly);		fprintf(stderr, "gdal_dataset = %p/n", f->gdal_img);		f->pd = GDALGetRasterCount(f->gdal_img);		f->w = GDALGetRasterXSize(f->gdal_img);		f->h = GDALGetRasterYSize(f->gdal_img);		f->no = 1;		for (int i = 0; i < f->pd; i++)			f->gdal_band[i] = GDALGetRasterBand(f->gdal_img, i+1);#else		assert(false);#endif	} else {		f->x = iio_read_image_float_vec(filename, &f->w, &f->h, &f->pd);		f->no = build_pyramid(f, f->max_octaves);		snprintf(f->x_filename, FILENAME_MAX, "%s", filename);		f->x_changed = false;	}}

HRESULT CImage::SetPixel(int nRows, int nCols, BYTE* pPixel){	if(m_pGdalImage==NULL)	{		return S_FALSE;	}	if(m_bFlip)	{		nRows=m_nRows-nRows-1;	}	if(nRows<0||nRows>=m_nRows||nCols<0||nCols>=m_nCols)	{		return S_OK;	}	for(int i=0;i<m_nBandNum;i++)	{		int band=i+1;		GDALDataType datatype=m_datatype;		int nBPB=m_nBPB;				BYTE* pGray=new BYTE[nBPB];		memcpy(pGray,pPixel+i*nBPB,nBPB);		((GDALRasterBand*)GDALGetRasterBand(m_pGdalImage,band))->		RasterIO(GF_Write,				 nCols,nRows,				 1,1,				 pGray,				 1,1,				 datatype,				 0,0);		delete [] pGray; pGray=NULL;	}	return S_OK;}

void writeGeoTiffF(char * fileName, float * result, int nRow, int nCol, double xMin, double yMax, double cellSize){		GDALAllRegister();	OGRRegisterAll();	GDALDatasetH hDstDS;	GDALDriverH hDriver;	GDALRasterBandH hBand;	double adfGeoTransform[6];	char *papszOptions[] = {"COMPRESS=LZW",NULL};	const char *pszFormat="GTiff";	if(NULL == (hDriver = GDALGetDriverByName(pszFormat)))	{		printf("ERROR: hDriver is null cannot output using GDAL/n");		exit(1);	}		hDstDS = GDALCreate(hDriver, fileName, nCol, nRow, 1, GDT_Float32, papszOptions);	adfGeoTransform[0] = xMin;	adfGeoTransform[1] = cellSize;	adfGeoTransform[2] = 0;	adfGeoTransform[3] = yMax;	adfGeoTransform[4] = 0;	adfGeoTransform[5] = -cellSize;	GDALSetGeoTransform(hDstDS,adfGeoTransform);	hBand=GDALGetRasterBand(hDstDS,1);	GDALSetRasterNoDataValue(hBand,-1);	GDALRasterIO(hBand, GF_Write, 0, 0, nCol, nRow, result, nCol, nRow, GDT_Float32, 0, 0 );		GDALClose(hDstDS);	return;}

GDALDatasetH CPL_STDCALLGDALCreateWarpedVRT( 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;}

void CDialog3D::OnBnClickedBtnOpen(){	// TODO: 在此添加控件通知处理程序代码	CFileDialog fileDlg(TRUE);	if(fileDlg.DoModal()!=IDOK)		return;	CString strExt	   =fileDlg.GetFileExt();	CString strPathName=fileDlg.GetPathName();	if(strExt=="txt")	{		ifstream ifs(strPathName,ios_base::in);		char tmpchr[2048];		ifs.getline(tmpchr,2048);		do 		{			PNT3D tmpPnt;			ifs.getline(tmpchr,2048);			sscanf(tmpchr,"%f,%f,%f",&tmpPnt.DX,&tmpPnt.DY,&tmpPnt.DZ);			m_vec_PNT3Ds.push_back(tmpPnt);		} while (!ifs.eof());		ifs.close();	}	if(strExt=="BMP"||strExt=="bmp"||strExt=="JPG"||strExt=="jpg"||strExt=="TIF"||strExt=="tif")	{		GDALAllRegister();		GDALDatasetH hSrcDS=GDALOpen(strPathName,GA_ReadOnly);		double adfGeoTrans[6];		double scalex=1,scaley=1;		GDALGetGeoTransform(hSrcDS,adfGeoTrans);		int xsize=GDALGetRasterXSize(hSrcDS);		int ysize=GDALGetRasterYSize(hSrcDS);		int tmpxsize=xsize,tmpysize=ysize;		if(xsize>800)		{			tmpxsize=800;			scalex=xsize/tmpxsize;		}		if(ysize>600)		{			tmpysize=600;			scaley=ysize/tmpysize;		}		float *dataIn=new float[tmpxsize*tmpysize];		GDALRasterIO(GDALGetRasterBand(hSrcDS,1),GF_Read,0,0,xsize,ysize,dataIn,tmpxsize,tmpysize,GDT_Float32,0,0);		for(int i=0;i<tmpxsize;i++)		{			for (int j=0;j<tmpysize;j++)			{				PNT3D tmpPnt;				tmpPnt.DX=adfGeoTrans[0]+adfGeoTrans[1]*i*scalex+adfGeoTrans[2]*j*scaley;				tmpPnt.DY=adfGeoTrans[3]+adfGeoTrans[4]*i*scalex+adfGeoTrans[5]*j*scaley;				//tmpPnt.DX=i;				//tmpPnt.DY=j;				tmpPnt.DZ=dataIn[j*tmpxsize+i];				m_vec_PNT3Ds.push_back(tmpPnt);			}		}		delete[]dataIn;		GDALClose(hSrcDS);	}	if(!m_vec_PNT3Ds.empty())	{		m_min_DX=m_max_DX=m_vec_PNT3Ds[0].DX;		m_min_DY=m_max_DY=m_vec_PNT3Ds[0].DY;		m_min_DZ=m_max_DZ=m_vec_PNT3Ds[0].DZ;		for (int i=0;i<m_vec_PNT3Ds.size();i++)		{			m_max_DX=max(m_vec_PNT3Ds[i].DX,m_max_DX);			m_min_DX=min(m_vec_PNT3Ds[i].DX,m_min_DX);			m_max_DY=max(m_vec_PNT3Ds[i].DY,m_max_DY);			m_min_DY=min(m_vec_PNT3Ds[i].DY,m_min_DY);			m_max_DZ=max(m_vec_PNT3Ds[i].DZ,m_max_DZ);			m_min_DZ=min(m_vec_PNT3Ds[i].DZ,m_min_DZ);		}		AfxMessageBox("数据读取成功!/n");		InvalidateRect(NULL,FALSE);	}}

int main( int argc, char ** argv ){    GDALDatasetH hSrcDS;    int         iY, iX, nOutLevel=0, nXSize, nYSize, iArg, nFillDist=0;    void        *pStream;    GInt16      *panData;    const char  *pszFilename = NULL;    GDALRasterBandH hSrcBand;    double       adfGeoTransform[6];    int          bEnableTrim = FALSE;    GInt16       noDataValue = 0;    int          bHasNoData;/* -------------------------------------------------------------------- *//*      Identify arguments.                                             *//* -------------------------------------------------------------------- */    for( iArg = 1; iArg < argc; iArg++ )    {        if( EQUAL(argv[iArg],"-trim") )            bEnableTrim = TRUE;        else if( EQUAL(argv[iArg],"-fill") )            nFillDist = atoi(argv[++iArg]);        else if( EQUAL(argv[iArg],"-level") )            nOutLevel = atoi(argv[++iArg]);        else        {            if( pszFilename != NULL )                Usage();            pszFilename = argv[iArg];        }    }    if( pszFilename == NULL )        Usage();/* -------------------------------------------------------------------- *//*      Open input file.                                                *//* -------------------------------------------------------------------- */    GDALAllRegister();    hSrcDS = GDALOpen( pszFilename, GA_ReadOnly );    if( hSrcDS == NULL )        exit(1);    hSrcBand = GDALGetRasterBand( hSrcDS, 1 );    noDataValue = (GInt16)GDALGetRasterNoDataValue(hSrcBand, &bHasNoData);    nXSize = GDALGetRasterXSize( hSrcDS );    nYSize = GDALGetRasterYSize( hSrcDS );    GDALGetGeoTransform( hSrcDS, adfGeoTransform );/* -------------------------------------------------------------------- *//*      Create output stream.                                           *//* -------------------------------------------------------------------- */    pStream = DTEDCreatePtStream( ".", nOutLevel );    if( pStream == NULL )        exit( 1 );/* -------------------------------------------------------------------- *//*      Process all the profiles.                                       *//* -------------------------------------------------------------------- */    panData = (GInt16 *) malloc(sizeof(GInt16) * nXSize);    for( iY = 0; iY < nYSize; iY++ )    {        GDALRasterIO( hSrcBand, GF_Read, 0, iY, nXSize, 1,                      panData, nXSize, 1, GDT_Int16, 0, 0 );        if (bHasNoData)        {            for( iX = 0; iX < nXSize; iX++ )            {                if (panData[iX] == noDataValue)                    panData[iX] = DTED_NODATA_VALUE;            }        }        for( iX = 0; iX < nXSize; iX++ )        {            DTEDWritePt( pStream,                         adfGeoTransform[0]                         + adfGeoTransform[1] * (iX + 0.5)                         + adfGeoTransform[2] * (iY + 0.5),                         adfGeoTransform[3]                         + adfGeoTransform[4] * (iX + 0.5)                         + adfGeoTransform[5] * (iY + 0.5),                         panData[iX] );        }    }    free( panData );/* -------------------------------------------------------------------- *//*      Cleanup.                                                        *///.........这里部分代码省略.........

//.........这里部分代码省略.........        printf( "Subdatasets:/n" );        for( i = 0; papszMetadata[i] != NULL; i++ )        {            printf( "  %s/n", papszMetadata[i] );        }    }/* -------------------------------------------------------------------- *//*      Report corners.                                                 *//* -------------------------------------------------------------------- */    printf( "Corner Coordinates:/n" );    GDALInfoReportCorner( hDataset, "Upper Left",                           0.0, 0.0 );    GDALInfoReportCorner( hDataset, "Lower Left",                           0.0, GDALGetRasterYSize(hDataset));    GDALInfoReportCorner( hDataset, "Upper Right",                           GDALGetRasterXSize(hDataset), 0.0 );    GDALInfoReportCorner( hDataset, "Lower Right",                           GDALGetRasterXSize(hDataset),                           GDALGetRasterYSize(hDataset) );    GDALInfoReportCorner( hDataset, "Center",                           GDALGetRasterXSize(hDataset)/2.0,                           GDALGetRasterYSize(hDataset)/2.0 );/* ==================================================================== *//*      Loop over bands.                                                *//* ==================================================================== */    for( iBand = 0; iBand < GDALGetRasterCount( hDataset ); iBand++ )    {        double      dfMin, dfMax, adfCMinMax[2], dfNoData;        int         bGotMin, bGotMax, bGotNodata;        int         nBlockXSize, nBlockYSize;        hBand = GDALGetRasterBand( hDataset, iBand+1 );        GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );        printf( "Band %d Block=%dx%d Type=%d, ColorInterp=%d/n", iBand+1,                nBlockXSize, nBlockYSize,                GDALGetRasterDataType(hBand),                GDALGetRasterColorInterpretation(hBand) );        dfMin = GDALGetRasterMinimum( hBand, &bGotMin );        dfMax = GDALGetRasterMaximum( hBand, &bGotMax );        printf( "  Min=%.3f/%d, Max=%.3f/%d",  dfMin, bGotMin, dfMax, bGotMax);                if( bComputeMinMax )        {            GDALComputeRasterMinMax( hBand, TRUE, adfCMinMax );            printf( ", Computed Min/Max=%.3f,%.3f",                     adfCMinMax[0], adfCMinMax[1] );        }        printf( "/n" );        dfNoData = GDALGetRasterNoDataValue( hBand, &bGotNodata );        if( bGotNodata )        {            printf( "  NoData Value=%g/n", dfNoData );        }        if( GDALGetOverviewCount(hBand) > 0 )        {            int		iOverview;            printf( "  Overviews: " );            for( iOverview = 0;                  iOverview < GDALGetOverviewCount(hBand);                 iOverview++ )

int main( int argc, char ** argv ){    const char *pszSrcFilename = NULL;    int anReqOverviews[1000];    int nReqOverviewCount = 0;    bool bMasks = false;    GDALAllRegister();    argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );    if( argc < 1 )        exit( -argc );/* -------------------------------------------------------------------- *//*      Process arguments.                                              *//* -------------------------------------------------------------------- */    for( int iArg = 1; iArg < argc; iArg++ )    {        if( EQUAL(argv[iArg],"-masks") )        {            bMasks = true;        }        else if( pszSrcFilename == NULL )        {            pszSrcFilename = argv[iArg];        }        else if( atoi(argv[iArg]) > 0 || EQUAL(argv[iArg],"0") )        {            anReqOverviews[nReqOverviewCount++] = atoi(argv[iArg]);        }        else        {            Usage();        }    }    if( pszSrcFilename == NULL )        Usage();/* -------------------------------------------------------------------- *//*      Open the input file.                                            *//* -------------------------------------------------------------------- */    GDALDatasetH hSrcDS = GDALOpen( pszSrcFilename, GA_ReadOnly );    if( hSrcDS == NULL )        exit( 1 );/* ==================================================================== *//*      Process all bands.                                              *//* ==================================================================== */    int iBand;    int nBandCount = GDALGetRasterCount( hSrcDS );    for( iBand = 0; iBand < nBandCount; iBand++ )    {        GDALRasterBandH hBaseBand = GDALGetRasterBand( hSrcDS, iBand+1 );/* -------------------------------------------------------------------- *//*      Process all overviews.                                          *//* -------------------------------------------------------------------- */        int iOverview;        int nOverviewCount = GDALGetOverviewCount( hBaseBand );        for( iOverview = 0; iOverview < nOverviewCount; iOverview++ )        {            GDALRasterBandH hSrcOver = GDALGetOverview( hBaseBand, iOverview );            if (hSrcOver == NULL)            {                fprintf(stderr, "skipping overview %d as being null/n", iOverview);                continue;            }/* -------------------------------------------------------------------- *//*      Is this a requested overview?                                   *//* -------------------------------------------------------------------- */            if( nReqOverviewCount > 0 )            {                int i;                for( i = 0; i < nReqOverviewCount; i++ )                {                    if( anReqOverviews[i] == iOverview )                        break;                }                if( i == nReqOverviewCount )                    continue;            }/* -------------------------------------------------------------------- *//*      Create matching output file.                                    *//* -------------------------------------------------------------------- */            CPLString osFilename;            osFilename.Printf( "%s_%d_%d.tif",                               CPLGetBasename(pszSrcFilename),                               iBand+1, iOverview );            DumpBand( hSrcDS, hSrcOver, osFilename );//.........这里部分代码省略.........

static gint print_handler( void * cb_data, void * scanline_in ){    unsigned char *scanline = (unsigned char *) scanline_in;    static GDALDatasetH  working_ds = NULL;    static int next_scanline = 0;    static GDALRasterBandH red_band, green_band, blue_band;    gint   width;    if( cb_data == NULL && scanline_in == NULL )    {        next_scanline = 0;        working_ds = NULL;        return -1;    }    if( working_ds != cb_data )    {        working_ds = (GDALDatasetH) cb_data;        next_scanline = 0;        red_band = GDALGetRasterBand( working_ds, 1 );        if( GDALGetRasterCount( working_ds ) >= 3 )        {            green_band = GDALGetRasterBand( working_ds, 2 );            blue_band = GDALGetRasterBand( working_ds, 3 );        }        else        {            green_band = blue_band = NULL;        }        if( red_band == NULL )            return -1;    }    width = GDALGetRasterXSize( working_ds );    if( green_band == NULL )    {        GByte	*grey;        int     i, value;        grey = g_new( GByte, width );                for( i = 0; i < width; i++ )        {            value = *(scanline++);            value += *(scanline++);            value += *(scanline++);            value = (value+1) / 3;            grey[i] = value;        }                GDALRasterIO( red_band, GF_Write, 0, next_scanline, width, 1,                       grey, width, 1, GDT_Byte, 1, 0 );        g_free( grey );    }    else    {        GDALRasterIO( red_band, GF_Write, 0, next_scanline, width, 1,                       scanline+0, width, 1, GDT_Byte, 3, 0 );        GDALRasterIO( green_band, GF_Write, 0, next_scanline, width, 1,                       scanline+1, width, 1, GDT_Byte, 3, 0 );        GDALRasterIO( blue_band, GF_Write, 0, next_scanline, width, 1,                       scanline+2, width, 1, GDT_Byte, 3, 0 );    }        next_scanline++;    return 0;}

//.........这里部分代码省略.........            dfXMax = sEnvelope.MaxX;            bIsXExtentSet = TRUE;        }        if ( !bIsYExtentSet )        {            dfYMin = sEnvelope.MinY;            dfYMax = sEnvelope.MaxY;            bIsYExtentSet = TRUE;        }    }/* -------------------------------------------------------------------- *//*      Perform gridding.                                               *//* -------------------------------------------------------------------- */    const double    dfDeltaX = ( dfXMax - dfXMin ) / nXSize;    const double    dfDeltaY = ( dfYMax - dfYMin ) / nYSize;    if ( !bQuiet )    {        printf( "Grid data type is /"%s/"/n", GDALGetDataTypeName(eType) );        printf( "Grid size = (%lu %lu)./n",                (unsigned long)nXSize, (unsigned long)nYSize );        printf( "Corner coordinates = (%f %f)-(%f %f)./n",                dfXMin - dfDeltaX / 2, dfYMax + dfDeltaY / 2,                dfXMax + dfDeltaX / 2, dfYMin - dfDeltaY / 2 );        printf( "Grid cell size = (%f %f)./n", dfDeltaX, dfDeltaY );        printf( "Source point count = %lu./n", (unsigned long)adfX.size() );        PrintAlgorithmAndOptions( eAlgorithm, pOptions );        printf("/n");    }    GDALRasterBandH hBand = GDALGetRasterBand( hDstDS, nBand );    if (adfX.size() == 0)    {        // FIXME: Shoulda' set to nodata value instead        GDALFillRaster( hBand, 0.0 , 0.0 );        return CE_None;    }    GUInt32 nXOffset, nYOffset;    int     nBlockXSize, nBlockYSize;    int     nDataTypeSize = GDALGetDataTypeSize(eType) / 8;    // Try to grow the work buffer up to 16 MB if it is smaller    GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );    const GUInt32 nDesiredBufferSize = 16*1024*1024;    if( (GUInt32)nBlockXSize < nXSize && (GUInt32)nBlockYSize < nYSize &&        (GUInt32)nBlockXSize < nDesiredBufferSize / (nBlockYSize * nDataTypeSize) )    {        int nNewBlockXSize  = nDesiredBufferSize / (nBlockYSize * nDataTypeSize);        nBlockXSize = (nNewBlockXSize / nBlockXSize) * nBlockXSize;        if( (GUInt32)nBlockXSize > nXSize )            nBlockXSize = nXSize;    }    else if( (GUInt32)nBlockXSize == nXSize && (GUInt32)nBlockYSize < nYSize &&             (GUInt32)nBlockYSize < nDesiredBufferSize / (nXSize * nDataTypeSize) )    {        int nNewBlockYSize = nDesiredBufferSize / (nXSize * nDataTypeSize);        nBlockYSize = (nNewBlockYSize / nBlockYSize) * nBlockYSize;        if( (GUInt32)nBlockYSize > nYSize )            nBlockYSize = nYSize;    }    CPLDebug("GDAL_GRID", "Work buffer: %d * %d", nBlockXSize, nBlockYSize);

bool CUtils::isFloat32DataType(GDALDatasetH hDataset){	return ( GDALGetRasterDataType( GDALGetRasterBand(hDataset, 1) ) == GDT_Float32 );}

void CUtils::calculateByteGeoTIFFStatistics(GDALDatasetH hDataset, int userBandNumber, byte flNoDataValueAsBackground, byte NoDataValue){	fputs("/nCalculate statistics.../n", stderr);	GDALRasterBandH hBand =  GDALGetRasterBand(hDataset, 1);	int cols  = GDALGetRasterBandXSize(hBand);	int rows  = GDALGetRasterBandYSize(hBand);	int bands = GDALGetRasterCount(hDataset);		byte * pbuf = NULL;	pbuf = (byte *)CPLMalloc(sizeof(byte)*cols);		byte min = 0, max = 0, mean = 0;	double stddev = 0;    double summ = 0;    int count = 0;		for(int band=1; band<=bands; band++)	{		if(userBandNumber != -1) fprintf(stderr, "Band %d.../n", userBandNumber);		else fprintf(stderr, "Band %d.../n", band);				hBand =  GDALGetRasterBand(hDataset, band);				if(flNoDataValueAsBackground) NoDataValue = getFloatNoDataValueAsBackground(hBand);				min = max = mean = stddev = summ = 0;		count = 0;		bool flFirst = true;				int pr = CUtils::progress_ln_ex(stderr, 0, 0, START_PROGRESS);		for(int i=0; i<rows; i++)		{			GDALRasterIO(hBand, GF_Read, 0, i, cols, 1, pbuf, cols, 1, GDT_Byte, 0, 0 );			for(int j=0; j<cols; j++) if(pbuf[j]!=NoDataValue) 			{				if(flFirst)				{					mean = pbuf[j];					min = max = mean;					flFirst = false;				}				else				{					mean += pbuf[j];					if( min > pbuf[j] ) min = pbuf[j];					if( max < pbuf[j] ) max = pbuf[j];				}				count++;			}			pr = CUtils::progress_ln_ex(stderr, i, rows, pr);		}		CUtils::progress_ln_ex(stderr, 0, 0, END_PROGRESS);				double dmean = 0;		if(count > 0) dmean = mean / (double)count;				pr = CUtils::progress_ln_ex(stderr, 0, 0, START_PROGRESS);		for(int i=0; i<rows; i++)		{			GDALRasterIO(hBand, GF_Read, 0, i, cols, 1, pbuf, cols, 1, GDT_Byte, 0, 0 );			for(int j=0; j<cols; j++) if(pbuf[j]!=NoDataValue) summ += ((double)pbuf[j]-dmean)*((double)pbuf[j]-dmean);									pr = CUtils::progress_ln_ex(stderr, i, rows, pr);		}		CUtils::progress_ln_ex(stderr, 0, 0, END_PROGRESS);				summ = 0; stddev = 0;		if((count-1)>0)		{			summ /= (double)(count-1);			if(summ!=0) stddev = sqrt(summ);		}					GDALSetRasterStatistics(hBand, min, max, mean, stddev);		GDALSetRasterNoDataValue(hBand, NoDataValue);	}		CPLFree(pbuf);}

//.........这里部分代码省略.........                if (eErr != CE_None)                {                    break;                }            }            GDALDatasetH hMemDS = GDALCreate(hMemDriver, "MEM:::",                                              nReqXSize, nReqYSize, 0,                                               eDataType, NULL);            if (hMemDS == NULL)            {                eErr = CE_Failure;                break;            }                        int iBand;            for(iBand = 0; iBand < nBands; iBand ++)            {                char** papszOptions = NULL;                char szTmp[64];                memset(szTmp, 0, sizeof(szTmp));                CPLPrintPointer(szTmp,                                pabyMEMDSBuffer + iBand * nDataTypeSize *                                nReqXSize * nReqYSize, sizeof(szTmp));                papszOptions = CSLSetNameValue(papszOptions, "DATAPOINTER", szTmp);                GDALAddBand(hMemDS, eDataType, papszOptions);                CSLDestroy(papszOptions);            }            if( hPrevOvrMemDS != NULL )            {                for(iBand = 0; iBand < nBands; iBand ++)                {                    GDALRasterBandH hDstOvrBand = GDALGetRasterBand(hMemDS, iBand+1);                    eErr = GDALRegenerateOverviews( GDALGetRasterBand(hPrevOvrMemDS, iBand+1),                                                    1, &hDstOvrBand,                                                    pszResampling,                                                    NULL, NULL );                    if( eErr != CE_None )                        break;                }                GDALClose(hPrevOvrMemDS);            }            GDALDatasetH hOutDS = GDALCreateCopy(hTileDriver,                                        osTempFileName.c_str(), hMemDS, FALSE,                                        papszTileDriverOptions, NULL, NULL);            GDALClose(hMemDS);            if (hOutDS)                GDALClose(hOutDS);            else            {                eErr = CE_Failure;                break;            }/* -------------------------------------------------------------------- *//*      Insert new entry into raster table                              *//* -------------------------------------------------------------------- */            vsi_l_offset nDataLength;            GByte *pabyData = VSIGetMemFileBuffer( osTempFileName.c_str(),                                                   &nDataLength, FALSE);

void convert(const input_arguments& args){// determine the subdataset pattern name using the first available file  if(args.verbose)    std::cout << "/tlooking for the names of subdatasets (or variables)... " << std::flush;    std::vector<std::string> band_names = modis2scidb::extract_subdatasets_pattern_names(args.source_file_name);  if(args.verbose)    std::cout << "OK!" << std::endl;  // check if band numbers are in the valid range  if(args.verbose)    std::cout << "/tchecking the range for choosed band numbers... " << std::flush;    std::size_t num_bands = args.bands.size();    for(std::size_t i = 0; i != num_bands; ++i)    if(args.bands[i] >= band_names.size())      throw modis2scidb::invalid_arg_value() << modis2scidb::error_description("band number is invalid!");    if(args.verbose)    std::cout << "OK!" << std::endl;// let's buffer each subdataset  if(args.verbose)    std::cout << "/tbuffering data... " << std::flush;  std::vector<boost::shared_array<unsigned char> > data_buffers;  std::vector<unsigned char*> aux_data_buffers;  std::vector<std::size_t> band_datatype_size;    int64_t ncols = 0;  int64_t nrows = 0;  for(std::size_t i = 0; i != num_bands; ++i)  {    if(args.verbose)      std::cout << "/n/t/tband #" << args.bands[i] << "... " << std::flush;    boost::format subdataset(band_names[args.bands[i]]);    subdataset.bind_arg(1, args.source_file_name);        GDALDatasetH dataset = GDALOpen(subdataset.str().c_str(), GA_ReadOnly);        if(dataset == 0)    {      boost::format err_msg("could not open subdataset: '%1%', for input hdf file: '%2%'!");      throw modis2scidb::gdal_error() << modis2scidb::error_description((err_msg % subdataset.str() % args.source_file_name).str());    }        GDALRasterBandH band = GDALGetRasterBand(dataset, 1);        if(band == 0)    {      GDALClose(dataset);      boost::format err_msg("could not access band: %1%!");      throw modis2scidb::gdal_error() << modis2scidb::error_description((err_msg % args.bands[i]).str());    }        if(i == 0)    {      ncols = GDALGetRasterBandXSize(band);      nrows = GDALGetRasterBandYSize(band);    }    else    {      if((GDALGetRasterBandXSize(band) != ncols) ||         (GDALGetRasterBandYSize(band) != nrows))      {        GDALClose(dataset);        throw modis2scidb::gdal_error() << modis2scidb::error_description("selected bands must have the same dimension (rows and cols)!");      }    }        GDALDataType pixel_type = GDALGetRasterDataType(band);        std::size_t pixel_size = modis2scidb::num_bytes(pixel_type);        band_datatype_size.push_back(pixel_size);        boost::shared_array<unsigned char> buffer(new unsigned char[ncols * nrows * pixel_size]);        data_buffers.push_back(buffer);        aux_data_buffers.push_back(buffer.get());        CPLErr result = GDALRasterIO(band, GF_Read, 0, 0, static_cast<int>(ncols), static_cast<int>(nrows), buffer.get(), static_cast<int>(ncols), static_cast<int>(nrows), pixel_type, 0, 0);    if(result == CE_Failure)    {      GDALClose(dataset);      boost::format err_msg("could not read subdataset: '%1%', for input hdf file: '%2%'!");      throw modis2scidb::gdal_error() << modis2scidb::error_description((err_msg % subdataset.str() % args.source_file_name).str());    }    GDALClose(dataset);    if(args.verbose)//.........这里部分代码省略.........

/** * @param theBandNumber the number of the band for which you want a color table * @param theList a pointer the object that will hold the color table * @return true of a color table was able to be read, false otherwise */QList<QgsColorRampShader::ColorRampItem> QgsGdalProviderBase::colorTable( GDALDatasetH gdalDataset, int theBandNumber )const{  QgsDebugMsg( "entered." );  QList<QgsColorRampShader::ColorRampItem> ct;  //Invalid band number, segfault prevention  if ( 0 >= theBandNumber )  {    QgsDebugMsg( "Invalid parameter" );    return ct;  }  GDALRasterBandH myGdalBand = GDALGetRasterBand( gdalDataset, theBandNumber );  GDALColorTableH myGdalColorTable = GDALGetRasterColorTable( myGdalBand );  if ( myGdalColorTable )  {    QgsDebugMsg( "Color table found" );    // load category labels    char ** categoryNames = GDALGetRasterCategoryNames( myGdalBand );    QVector<QString> labels;    if ( categoryNames )    {      int i = 0;      while ( categoryNames[i] )      {        labels.append( QString( categoryNames[i] ) );        i++;      }    }    int myEntryCount = GDALGetColorEntryCount( myGdalColorTable );    GDALColorInterp myColorInterpretation =  GDALGetRasterColorInterpretation( myGdalBand );    QgsDebugMsg( "Color Interpretation: " + QString::number(( int )myColorInterpretation ) );    GDALPaletteInterp myPaletteInterpretation  = GDALGetPaletteInterpretation( myGdalColorTable );    QgsDebugMsg( "Palette Interpretation: " + QString::number(( int )myPaletteInterpretation ) );    const GDALColorEntry* myColorEntry = 0;    for ( int myIterator = 0; myIterator < myEntryCount; myIterator++ )    {      myColorEntry = GDALGetColorEntry( myGdalColorTable, myIterator );      if ( !myColorEntry )      {        continue;      }      else      {        QString label = labels.value( myIterator );        if ( label.isEmpty() )        {          label = QString::number( myIterator );        }        //Branch on the color interpretation type        if ( myColorInterpretation == GCI_GrayIndex )        {          QgsColorRampShader::ColorRampItem myColorRampItem;          myColorRampItem.value = ( double )myIterator;          myColorRampItem.label = label;          myColorRampItem.color = QColor::fromRgb( myColorEntry->c1, myColorEntry->c1, myColorEntry->c1, myColorEntry->c4 );          ct.append( myColorRampItem );        }        else if ( myColorInterpretation == GCI_PaletteIndex )        {          QgsColorRampShader::ColorRampItem myColorRampItem;          myColorRampItem.value = ( double )myIterator;          myColorRampItem.label = label;          //Branch on palette interpretation          if ( myPaletteInterpretation  == GPI_RGB )          {            myColorRampItem.color = QColor::fromRgb( myColorEntry->c1, myColorEntry->c2, myColorEntry->c3, myColorEntry->c4 );          }          else if ( myPaletteInterpretation  == GPI_CMYK )          {            myColorRampItem.color = QColor::fromCmyk( myColorEntry->c1, myColorEntry->c2, myColorEntry->c3, myColorEntry->c4 );          }          else if ( myPaletteInterpretation  == GPI_HLS )          {            myColorRampItem.color = QColor::fromHsv( myColorEntry->c1, myColorEntry->c3, myColorEntry->c2, myColorEntry->c4 );          }          else          {            myColorRampItem.color = QColor::fromRgb( myColorEntry->c1, myColorEntry->c1, myColorEntry->c1, myColorEntry->c4 );          }          ct.append( myColorRampItem );        }        else        {          QgsDebugMsg( "Color interpretation type not supported yet" );          return ct;        }      }    }  }//.........这里部分代码省略.........

CPLErr CPL_STDCALLGDALFillNodata( GDALRasterBandH hTargetBand,                GDALRasterBandH hMaskBand,                double dfMaxSearchDist,                CPL_UNUSED int bDeprecatedOption,                int nSmoothingIterations,                char **papszOptions,                GDALProgressFunc pfnProgress,                void * pProgressArg ){    VALIDATE_POINTER1( hTargetBand, "GDALFillNodata", CE_Failure );    const int nXSize = GDALGetRasterBandXSize(hTargetBand);    const int nYSize = GDALGetRasterBandYSize(hTargetBand);    if( dfMaxSearchDist == 0.0 )        dfMaxSearchDist = std::max(nXSize, nYSize) + 1;    const int nMaxSearchDist = static_cast<int>(floor(dfMaxSearchDist));    // Special "x" pixel values identifying pixels as special.    GDALDataType eType = GDT_UInt16;    GUInt32 nNoDataVal = 65535;    if( nXSize > 65533 || nYSize > 65533 )    {        eType = GDT_UInt32;        nNoDataVal = 4000002;    }    if( hMaskBand == nullptr )        hMaskBand = GDALGetMaskBand( hTargetBand );    // If there are smoothing iterations, reserve 10% of the progress for them.    const double dfProgressRatio = nSmoothingIterations > 0 ? 0.9 : 1.0;    const char* pszNoData = CSLFetchNameValue(papszOptions, "NODATA");    bool bHasNoData = false;    float fNoData = 0.0f;    if( pszNoData )    {        bHasNoData = true;        fNoData = static_cast<float>(CPLAtof(pszNoData));    }/* -------------------------------------------------------------------- *//*      Initialize progress counter.                                    *//* -------------------------------------------------------------------- */    if( pfnProgress == nullptr )        pfnProgress = GDALDummyProgress;    if( !pfnProgress( 0.0, "Filling...", pProgressArg ) )    {        CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );        return CE_Failure;    }/* -------------------------------------------------------------------- *//*      Determine format driver for temp work files.                    *//* -------------------------------------------------------------------- */    CPLString osTmpFileDriver = CSLFetchNameValueDef(            papszOptions, "TEMP_FILE_DRIVER", "GTiff");    GDALDriverH hDriver = GDALGetDriverByName(osTmpFileDriver.c_str());    if( hDriver == nullptr )    {        CPLError(CE_Failure, CPLE_AppDefined,                 "Given driver is not registered");        return CE_Failure;    }    if( GDALGetMetadataItem(hDriver, GDAL_DCAP_CREATE, nullptr) == nullptr )    {        CPLError(CE_Failure, CPLE_AppDefined,                 "Given driver is incapable of creating temp work files");        return CE_Failure;    }    char **papszWorkFileOptions = nullptr;    if( osTmpFileDriver == "GTiff" )    {        papszWorkFileOptions = CSLSetNameValue(                papszWorkFileOptions, "COMPRESS", "LZW");        papszWorkFileOptions = CSLSetNameValue(                papszWorkFileOptions, "BIGTIFF", "IF_SAFER");    }/* -------------------------------------------------------------------- *//*      Create a work file to hold the Y "last value" indices.          *//* -------------------------------------------------------------------- */    const CPLString osTmpFile = CPLGenerateTempFilename("");    const CPLString osYTmpFile = osTmpFile + "fill_y_work.tif";    GDALDatasetH hYDS =        GDALCreate( hDriver, osYTmpFile, nXSize, nYSize, 1,                    eType, papszWorkFileOptions );    if( hYDS == nullptr )    {//.........这里部分代码省略.........

static void DumpBand( GDALDatasetH hBaseDS, GDALRasterBandH hSrcOver,                      const char *pszName ){/* -------------------------------------------------------------------- *//*      Get base ds info.                                               *//* -------------------------------------------------------------------- */    double adfGeoTransform[6];    bool bHaveGT = GDALGetGeoTransform( hBaseDS, adfGeoTransform ) == CE_None;    int nOrigXSize = GDALGetRasterXSize( hBaseDS );    int nOrigYSize = GDALGetRasterYSize( hBaseDS );/* -------------------------------------------------------------------- *//*      Create matching output file.                                    *//* -------------------------------------------------------------------- */    int nXSize = GDALGetRasterBandXSize( hSrcOver );    int nYSize = GDALGetRasterBandYSize( hSrcOver );    GDALDataType eDT = GDALGetRasterDataType( hSrcOver );    GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );    GDALDatasetH hDstDS = GDALCreate( hDriver, pszName, nXSize, nYSize,                                      1, eDT, NULL );    if( hDstDS == NULL )        exit( 1 );/* -------------------------------------------------------------------- *//*      Apply corresponding georeferencing, scaled to size.             *//* -------------------------------------------------------------------- */    if( bHaveGT )    {        double adfOvGeoTransform[6];        memcpy( adfOvGeoTransform, adfGeoTransform,                sizeof(double) * 6 );        adfOvGeoTransform[1] *= (nOrigXSize / (double) nXSize);        adfOvGeoTransform[2] *= (nOrigXSize / (double) nXSize);        adfOvGeoTransform[4] *= (nOrigYSize / (double) nYSize);        adfOvGeoTransform[5] *= (nOrigYSize / (double) nYSize);        GDALSetGeoTransform( hDstDS, adfOvGeoTransform );        GDALSetProjection( hDstDS, GDALGetProjectionRef( hBaseDS ) );    }/* -------------------------------------------------------------------- *//*      Copy over all the image data.                                   *//* -------------------------------------------------------------------- */    void *pData = CPLMalloc(64 * nXSize);    for( int iLine = 0; iLine < nYSize; iLine++ )    {        GDALRasterIO( hSrcOver, GF_Read, 0, iLine, nXSize, 1,                      pData, nXSize, 1, eDT, 0, 0 );        GDALRasterIO( GDALGetRasterBand( hDstDS, 1 ), GF_Write,                      0, iLine, nXSize, 1,                      pData, nXSize, 1, eDT, 0, 0 );    }    CPLFree( pData );    GDALClose( hDstDS );}

bool QgsAlignRaster::createAndWarp( const Item& raster ){  GDALDriverH hDriver = GDALGetDriverByName( "GTiff" );  if ( !hDriver )  {    mErrorMessage = QString( "GDALGetDriverByName(GTiff) failed." );    return false;  }  // Open the source file.  GDALDatasetH hSrcDS = GDALOpen( raster.inputFilename.toLocal8Bit().constData(), GA_ReadOnly );  if ( !hSrcDS )  {    mErrorMessage = QObject::tr( "Unable to open input file: " ) + raster.inputFilename;    return false;  }  // Create output with same datatype as first input band.  int bandCount = GDALGetRasterCount( hSrcDS );  GDALDataType eDT = GDALGetRasterDataType( GDALGetRasterBand( hSrcDS, 1 ) );  // Create the output file.  GDALDatasetH hDstDS;  hDstDS = GDALCreate( hDriver, raster.outputFilename.toLocal8Bit().constData(), mXSize, mYSize,                       bandCount, eDT, NULL );  if ( !hDstDS )  {    GDALClose( hSrcDS );    mErrorMessage = QObject::tr( "Unable to create output file: " ) + raster.outputFilename;    return false;  }  // Write out the projection definition.  GDALSetProjection( hDstDS, mCrsWkt.toAscii().constData() );  GDALSetGeoTransform( hDstDS, ( double* )mGeoTransform );  // Copy the color table, if required.  GDALColorTableH hCT = GDALGetRasterColorTable( GDALGetRasterBand( hSrcDS, 1 ) );  if ( hCT != NULL )    GDALSetRasterColorTable( GDALGetRasterBand( hDstDS, 1 ), hCT );  // -----------------------------------------------------------------------  // Setup warp options.  GDALWarpOptions* psWarpOptions = GDALCreateWarpOptions();  psWarpOptions->hSrcDS = hSrcDS;  psWarpOptions->hDstDS = hDstDS;  psWarpOptions->nBandCount = GDALGetRasterCount( hSrcDS );  psWarpOptions->panSrcBands = ( int * ) CPLMalloc( sizeof( int ) * psWarpOptions->nBandCount );  psWarpOptions->panDstBands = ( int * ) CPLMalloc( sizeof( int ) * psWarpOptions->nBandCount );  for ( int i = 0; i < psWarpOptions->nBandCount; ++i )  {    psWarpOptions->panSrcBands[i] = i + 1;    psWarpOptions->panDstBands[i] = i + 1;  }  psWarpOptions->eResampleAlg = ( GDALResampleAlg ) raster.resampleMethod;  // our progress function  psWarpOptions->pfnProgress = _progress;  psWarpOptions->pProgressArg = this;  // Establish reprojection transformer.  psWarpOptions->pTransformerArg =    GDALCreateGenImgProjTransformer( hSrcDS, GDALGetProjectionRef( hSrcDS ),                                     hDstDS, GDALGetProjectionRef( hDstDS ),                                     FALSE, 0.0, 1 );  psWarpOptions->pfnTransformer = GDALGenImgProjTransform;  double rescaleArg[2];  if ( raster.rescaleValues )  {    rescaleArg[0] = raster.srcCellSizeInDestCRS; // source cell size    rescaleArg[1] = mCellSizeX * mCellSizeY;  // destination cell size    psWarpOptions->pfnPreWarpChunkProcessor = rescalePreWarpChunkProcessor;    psWarpOptions->pfnPostWarpChunkProcessor = rescalePostWarpChunkProcessor;    psWarpOptions->pPreWarpProcessorArg = rescaleArg;    psWarpOptions->pPostWarpProcessorArg = rescaleArg;    // force use of float32 data type as that is what our pre/post-processor uses    psWarpOptions->eWorkingDataType = GDT_Float32;  }  // Initialize and execute the warp operation.  GDALWarpOperation oOperation;  oOperation.Initialize( psWarpOptions );  oOperation.ChunkAndWarpImage( 0, 0, mXSize, mYSize );  GDALDestroyGenImgProjTransformer( psWarpOptions->pTransformerArg );  GDALDestroyWarpOptions( psWarpOptions );  GDALClose( hDstDS );  GDALClose( hSrcDS );  return true;}

//.........这里部分代码省略.........    source = parm.source->answer;    output = parm.output->answer;    flip = 0;    if (flag.h->answer)	flip |= FLIP_H;    if (flag.v->answer)	flip |= FLIP_V;    if (parm.title->answer) {	title = G_store(parm.title->answer);	G_strip(title);    }    else	title = NULL;    if (!input && !source)	G_fatal_error(_("One of options <%s> or <%s> must be given"),		      parm.input->key, parm.source->key);    if (input && source)	G_fatal_error(_("Option <%s> and <%s> are mutually exclusive"),		      parm.input->key, parm.source->key);        if (input && !G_is_absolute_path(input)) {	char path[GPATH_MAX];	getcwd(path, sizeof(path));	strcat(path, "/");	strcat(path, input);	input = G_store(path);    }    if (!input)	input = source;    hDS = GDALOpen(input, GA_ReadOnly);    if (hDS == NULL)	return 1;    setup_window(&cellhd, hDS, &flip);    check_projection(&cellhd, hDS, flag.o->answer);    Rast_set_window(&cellhd);    if (parm.band->answer)	min_band = max_band = atoi(parm.band->answer);    else	min_band = 1, max_band = GDALGetRasterCount(hDS);    G_verbose_message(_("Proceeding with import..."));    if (max_band > min_band) {	if (I_find_group(output) == 1)	    G_warning(_("Imagery group <%s> already exists and will be overwritten."), output);	I_init_group_ref(&reference);    }    for (band = min_band; band <= max_band; band++) {	char *output2, *title2 = NULL;	G_message(_("Reading band %d of %d..."),		  band, GDALGetRasterCount( hDS ));	hBand = GDALGetRasterBand(hDS, band);	if (!hBand)	    G_fatal_error(_("Selected band (%d) does not exist"), band);	if (max_band > min_band) {	    G_asprintf(&output2, "%s.%d", output, band);	    if (title)		G_asprintf(&title2, "%s (band %d)", title, band);	    G_debug(1, "Adding raster map <%s> to group <%s>", output2, output);	    I_add_file_to_group_ref(output2, G_mapset(), &reference);	}	else {	    output2 = G_store(output);	    if (title)		title2 = G_store(title);	}	query_band(hBand, output2, flag.r->answer, &cellhd, &info);	create_map(input, band, output2, &cellhd, &info, title, flip);	G_free(output2);	G_free(title2);    }    if (flag.e->answer)	update_default_window(&cellhd);    /* Create the imagery group if multiple bands are imported */    if (max_band > min_band) {    	I_put_group_ref(output, &reference);	I_put_group(output);	G_message(_("Imagery group <%s> created"), output);    }    exit(EXIT_SUCCESS);}

int main( int argc, const char* argv[] ){    GDALDriverH   hDriver;    double        adfGeoTransform[6];    GDALDatasetH  in_Dataset;    GDALDatasetH  mask_Dataset;    GDALDatasetH  out_Dataset;    GDALRasterBandH mask_band;    unsigned char   *out_scan_line, *data_scan_line;    int             nBlockXSize, nBlockYSize;    int             bGotMin, bGotMax;    int             bands;    int             xsize;    double          adfMinMax[2];    int             valid_data_pixels[10];    int             saturated_data_pixels[10];    int             y_index, x;    GDALRasterBandH  out_band;            if ( argc != 3 ) {        ussage(argv[0]);    }        GDALAllRegister();        /* Set cache to something reasonable.. - 1/2 gig*/    CPLSetConfigOption( "GDAL_CACHEMAX", "512" );    /* open datasets..*/    in_Dataset = GDAL_open_read( argv[1]);    out_Dataset = make_me_a_sandwitch(&in_Dataset, argv[2]);        /* Basic info on source dataset..*/    GDALGetBlockSize(GDALGetRasterBand( in_Dataset, 1 ) , &nBlockXSize, &nBlockYSize );        /* Loop though bands, checking for saturated pixels .... */    xsize = GDALGetRasterXSize( in_Dataset );    data_scan_line = (char *) CPLMalloc(sizeof(char)*xsize);    out_scan_line = (char *) CPLMalloc(sizeof(char)*xsize);        /* The output band... */    out_band =  GDALGetRasterBand( out_Dataset, 1);       /* wipe counters.. */     for (bands=1; bands <= GDALGetRasterCount( in_Dataset ); bands ++ ) {        valid_data_pixels[bands] = 0;        saturated_data_pixels[bands] = 0;    }        /* loop though the lines of the data, looking for no data and saturated pixels..*/    for (y_index = 0; y_index <GDALGetRasterYSize( in_Dataset ); y_index ++ ) {        for (bands=1; bands <= GDALGetRasterCount( in_Dataset ); bands ++ ) {            GDALRasterBandH data_band;            /* Read data..*/            data_band =  GDALGetRasterBand( in_Dataset, bands);            GDALRasterIO( data_band, GF_Read, 0, y_index, xsize , 1, data_scan_line, xsize , 1, GDT_Byte, 0, 0 );            /* If first band, then copy into output slice.. */            if (bands==1) {                unsigned char  data_value;                for(x=0; x < xsize; x++) {                    /*shift to make darker...*/                   out_scan_line[x] = data_scan_line[x] >> 1 ;                   if ( out_scan_line[x] ==0 && data_scan_line[x] != 0) {out_scan_line[x] = 1;}                }            }                        /* Loop though, looking for saturated pixels and no-data values.. */            for(x=0; x < xsize; x++) {                if (  data_scan_line[x] != 0 )  {                   valid_data_pixels[bands] += 1;                   if ( data_scan_line[x] == 255 ) {                    saturated_data_pixels[bands] += 1;                    out_scan_line[x] = 255;                   }                }            }        }        GDALRasterIO( out_band, GF_Write, 0, y_index, xsize , 1, out_scan_line, xsize , 1, GDT_Byte, 0, 0 );    }

int QgsZonalStatistics::calculateStatistics( QProgressDialog* p ){  if ( !mPolygonLayer || mPolygonLayer->geometryType() != QGis::Polygon )  {    return 1;  }  QgsVectorDataProvider* vectorProvider = mPolygonLayer->dataProvider();  if ( !vectorProvider )  {    return 2;  }  //open the raster layer and the raster band  GDALAllRegister();  GDALDatasetH inputDataset = GDALOpen( mRasterFilePath.toLocal8Bit().data(), GA_ReadOnly );  if ( inputDataset == NULL )  {    return 3;  }  if ( GDALGetRasterCount( inputDataset ) < ( mRasterBand - 1 ) )  {    GDALClose( inputDataset );    return 4;  }  GDALRasterBandH rasterBand = GDALGetRasterBand( inputDataset, mRasterBand );  if ( rasterBand == NULL )  {    GDALClose( inputDataset );    return 5;  }  mInputNodataValue = GDALGetRasterNoDataValue( rasterBand, NULL );  //get geometry info about raster layer  int nCellsX = GDALGetRasterXSize( inputDataset );  int nCellsY = GDALGetRasterYSize( inputDataset );  double geoTransform[6];  if ( GDALGetGeoTransform( inputDataset, geoTransform ) != CE_None )  {    GDALClose( inputDataset );    return 6;  }  double cellsizeX = geoTransform[1];  if ( cellsizeX < 0 )  {    cellsizeX = -cellsizeX;  }  double cellsizeY = geoTransform[5];  if ( cellsizeY < 0 )  {    cellsizeY = -cellsizeY;  }  QgsRectangle rasterBBox( geoTransform[0], geoTransform[3] - ( nCellsY * cellsizeY ), geoTransform[0] + ( nCellsX * cellsizeX ), geoTransform[3] );  //add the new count, sum, mean fields to the provider  QList<QgsField> newFieldList;  QgsField countField( mAttributePrefix + "count", QVariant::Double );  QgsField sumField( mAttributePrefix + "sum", QVariant::Double );  QgsField meanField( mAttributePrefix + "mean", QVariant::Double );  newFieldList.push_back( countField );  newFieldList.push_back( sumField );  newFieldList.push_back( meanField );  if ( !vectorProvider->addAttributes( newFieldList ) )  {    return 7;  }  //index of the new fields  int countIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "count" );  int sumIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "sum" );  int meanIndex = vectorProvider->fieldNameIndex( mAttributePrefix + "mean" );  if ( countIndex == -1 || sumIndex == -1 || meanIndex == -1 )  {    return 8;  }  //progress dialog  long featureCount = vectorProvider->featureCount();  if ( p )  {    p->setMaximum( featureCount );  }  //iterate over each polygon  vectorProvider->select( QgsAttributeList(), QgsRectangle(), true, false );  vectorProvider->rewind();  QgsFeature f;  double count = 0;  double sum = 0;  double mean = 0;  int featureCounter = 0;  while ( vectorProvider->nextFeature( f ) )  {    qWarning( "%d", featureCounter );    if ( p )//.........这里部分代码省略.........

//.........这里部分代码省略........./* -------------------------------------------------------------------- */    printf( "Corner Coordinates:/n" );    GDALInfoReportCorner( hDataset, hTransform, "Upper Left",                           0.0, 0.0 );    GDALInfoReportCorner( hDataset, hTransform, "Lower Left",                           0.0, GDALGetRasterYSize(hDataset));    GDALInfoReportCorner( hDataset, hTransform, "Upper Right",                           GDALGetRasterXSize(hDataset), 0.0 );    GDALInfoReportCorner( hDataset, hTransform, "Lower Right",                           GDALGetRasterXSize(hDataset),                           GDALGetRasterYSize(hDataset) );    GDALInfoReportCorner( hDataset, hTransform, "Center",                           GDALGetRasterXSize(hDataset)/2.0,                           GDALGetRasterYSize(hDataset)/2.0 );    if( hTransform != NULL )    {        OCTDestroyCoordinateTransformation( hTransform );        hTransform = NULL;    }    /* ==================================================================== *//*      Loop over bands.                                                *//* ==================================================================== */    for( iBand = 0; iBand < GDALGetRasterCount( hDataset ); iBand++ )    {        double      dfMin, dfMax, adfCMinMax[2], dfNoData;        int         bGotMin, bGotMax, bGotNodata, bSuccess;        int         nBlockXSize, nBlockYSize, nMaskFlags;        double      dfMean, dfStdDev;        GDALColorTableH	hTable;        CPLErr      eErr;        hBand = GDALGetRasterBand( hDataset, iBand+1 );        if( bSample )        {            float afSample[10000];            int   nCount;            nCount = GDALGetRandomRasterSample( hBand, 10000, afSample );            printf( "Got %d samples./n", nCount );        }                GDALGetBlockSize( hBand, &nBlockXSize, &nBlockYSize );        printf( "Band %d Block=%dx%d Type=%s, ColorInterp=%s/n", iBand+1,                nBlockXSize, nBlockYSize,                GDALGetDataTypeName(                    GDALGetRasterDataType(hBand)),                GDALGetColorInterpretationName(                    GDALGetRasterColorInterpretation(hBand)) );        if( GDALGetDescription( hBand ) != NULL             && strlen(GDALGetDescription( hBand )) > 0 )            printf( "  Description = %s/n", GDALGetDescription(hBand) );        dfMin = GDALGetRasterMinimum( hBand, &bGotMin );        dfMax = GDALGetRasterMaximum( hBand, &bGotMax );        if( bGotMin || bGotMax || bComputeMinMax )        {            printf( "  " );            if( bGotMin )                printf( "Min=%.3f ", dfMin );            if( bGotMax )                printf( "Max=%.3f ", dfMax );        

static intComputeEqualizationLUTs( GDALDatasetH hDataset, int nLUTBins,                         double **ppadfScaleMin, double **ppadfScaleMax,                          int ***ppapanLUTs,                          GDALProgressFunc pfnProgress ){    int iBand;    int nBandCount = GDALGetRasterCount(hDataset);    int nHistSize = 0;    int *panHistogram = NULL;    // For now we always compute min/max    *ppadfScaleMin = (double *) CPLCalloc(sizeof(double),nBandCount);    *ppadfScaleMax = (double *) CPLCalloc(sizeof(double),nBandCount);    *ppapanLUTs = (int **) CPLCalloc(sizeof(int *),nBandCount);/* ==================================================================== *//*      Process all bands.                                              *//* ==================================================================== */    for( iBand = 0; iBand < nBandCount; iBand++ )    {        GDALRasterBandH hBand = GDALGetRasterBand( hDataset, iBand+1 );        CPLErr eErr;/* -------------------------------------------------------------------- *//*      Get a reasonable histogram.                                     *//* -------------------------------------------------------------------- */        eErr =            GDALGetDefaultHistogram( hBand,                                      *ppadfScaleMin + iBand,                                     *ppadfScaleMax + iBand,                                     &nHistSize, &panHistogram,                                      TRUE, pfnProgress, NULL );        if( eErr != CE_None )            return FALSE;        panHistogram[0] = 0; // zero out extremes (nodata, etc)        panHistogram[nHistSize-1] = 0;/* -------------------------------------------------------------------- *//*      Total histogram count, and build cumulative histogram.          *//*      We take care to use big integers as there may be more than 4    *//*      Gigapixels.                                                     *//* -------------------------------------------------------------------- */        GIntBig *panCumHist = (GIntBig *) CPLCalloc(sizeof(GIntBig),nHistSize);        GIntBig nTotal = 0;        int iHist;        for( iHist = 0; iHist < nHistSize; iHist++ )        {            panCumHist[iHist] = nTotal + panHistogram[iHist]/2;            nTotal += panHistogram[iHist];        }        CPLFree( panHistogram );        if( nTotal == 0 )        {            CPLError( CE_Warning, CPLE_AppDefined,                       "Zero value entries in histogram, results will not be meaningful." );            nTotal = 1;        }/* -------------------------------------------------------------------- *//*      Now compute a LUT from the cumulative histogram.                *//* -------------------------------------------------------------------- */        int *panLUT = (int *) CPLCalloc(sizeof(int),nLUTBins);        int iLUT;        for( iLUT = 0; iLUT < nLUTBins; iLUT++ )        {            iHist = (iLUT * nHistSize) / nLUTBins;            int nValue = (int) ((panCumHist[iHist] * nLUTBins) / nTotal);            panLUT[iLUT] = MAX(0,MIN(nLUTBins-1,nValue));        }         (*ppapanLUTs)[iBand] = panLUT;    }    return TRUE;}

void doTif2Con(char fileName[], int headerFormat, int BandNumber, char separator, char SHPMaskFile[]){	bool flMask = false;		OGRDataSourceH poDS = NULL;	OGRSFDriverH poDriver = NULL;			if(strlen(SHPMaskFile)>0) 	{		flMask = true;		poDS = OGROpen(SHPMaskFile, FALSE, &poDriver);	}	if(flMask && (poDS == NULL)) 	{		fputs("/nError open mask file!!!/n/n", stderr);		return ;	}    GDALDatasetH  pDataset;    GDALRasterBandH pBand;				    pDataset = GDALOpen( fileName, GA_ReadOnly );    if(pDataset!=NULL)    {		int bands = 0;						if(0 == BandNumber) bands = GDALGetRasterCount( pDataset );		else bands = 1;				int cols = GDALGetRasterXSize(pDataset);		int rows = GDALGetRasterYSize(pDataset);						double adfGeoTransform[6];				float xOrigin 	  = 0;		float yOrigin 	  = 0;		float pixelWidth  = 0;		float pixelHeight = 0;				if( GDALGetGeoTransform( pDataset, adfGeoTransform ) == CE_None )				{			xOrigin 	= adfGeoTransform[0];			yOrigin 	= adfGeoTransform[3];			pixelWidth 	= adfGeoTransform[1];			pixelHeight = adfGeoTransform[5];		}				float *** pdata = NULL;		pdata = new float**[bands];		for(int i=0; i<bands; i++) pdata[i] = new float*[rows];		for(int i=0; i<bands; i++) for(int j=0; j<rows; j++) pdata[i][j] = new float[cols];		for(int i=0; i<bands; i++) for(int j=0; j<rows; j++) for(int k=0; k<cols; k++) pdata[i][j][k] = 0;				void *pbuf = NULL;		pBand = GDALGetRasterBand(pDataset, 1);		pbuf = CUtils::mallocData(pBand, pbuf, cols);				printHeader(headerFormat, BandNumber, bands, rows, cols, separator);				if(0 == BandNumber)		{			for(int i=1; i<=bands; i++)			{				pBand = GDALGetRasterBand(pDataset, i);				for(int j=0; j<rows; j++)				{										CUtils::getRasterLine(pBand, j, cols, pbuf);					for(int k=0; k<cols; k++) pdata[i-1][j][k] = CUtils::getDataAsFloat(pBand, pbuf, k);				}			}		}		else		{				pBand = GDALGetRasterBand(pDataset, BandNumber);				for(int j=0; j<rows; j++)				{					CUtils::getRasterLine(pBand, j, cols, pbuf);					for(int k=0; k<cols; k++) pdata[0][j][k] = CUtils::getDataAsFloat(pBand, pbuf, k);				}		}				CPLFree(pbuf);				GDALClose(pDataset);			printData((const float ***) pdata, headerFormat, bands, rows, cols,				  xOrigin, yOrigin, pixelWidth, pixelHeight, separator, poDS);				if(poDS != NULL) OGR_DS_Destroy(poDS);				for(int i=0; i<bands; i++) for(int j=0; j<rows; j++) delete [] pdata[i][j];		for(int i=0; i<bands; i++) delete [] pdata[i];		delete [] pdata;		pdata = NULL;		fputs("/nProcessing COMPLETE./n/n", stderr);	}	else	fputs("/nError open input image!!!/n/n", stderr);//.........这里部分代码省略.........

void MDAL::LoaderGdal::parseRasterBands( const MDAL::GdalDataset *cfGDALDataset ){  for ( unsigned int i = 1; i <= cfGDALDataset->mNBands; ++i ) // starts with 1 .... ehm....  {    // Get Band    GDALRasterBandH gdalBand = GDALGetRasterBand( cfGDALDataset->mHDataset, static_cast<int>( i ) );    if ( !gdalBand )    {      throw MDAL_Status::Err_InvalidData;    }    // Reference time    metadata_hash global_metadata = parseMetadata( cfGDALDataset->mHDataset );    parseGlobals( global_metadata );    // Get metadata    metadata_hash metadata = parseMetadata( gdalBand );    std::string band_name;    double time = std::numeric_limits<double>::min();    bool is_vector;    bool is_x;    if ( parseBandInfo( cfGDALDataset, metadata, band_name, &time, &is_vector, &is_x ) )    {      continue;    }    // Add to data structures    std::vector<GDALRasterBandH>::size_type data_count = is_vector ? 2 : 1;    std::vector<GDALRasterBandH>::size_type data_index = is_x ? 0 : 1;    if ( mBands.find( band_name ) == mBands.end() )    {      // this Face is not yet added at all      // => create new map      timestep_map qMap;      std::vector<GDALRasterBandH> raster_bands( data_count );      raster_bands[data_index] = gdalBand;      qMap[time] = raster_bands;      mBands[band_name] = qMap;    }    else    {      timestep_map::iterator timestep = mBands[band_name].find( time );      if ( timestep == mBands[band_name].end() )      {        // Face is there, but new timestep        // => create just new map entry        std::vector<GDALRasterBandH> raster_bands( data_count );        raster_bands[data_index] = gdalBand;        mBands[band_name][time] = raster_bands;      }      else      {        // Face is there, and timestep too, this must be other part        // of the existing vector        timestep->second[data_index] = gdalBand;      }    }  }}

//.........这里部分代码省略.........                sEnvelop.MaxY = sLayerEnvelop.MaxY;                if (hSRS == NULL)                    hSRS = OGR_L_GetSpatialRef(hLayer);                bFirstLayer = FALSE;            }            else            {                sEnvelop.MinX = MIN(sEnvelop.MinX, sLayerEnvelop.MinX);                sEnvelop.MinY = MIN(sEnvelop.MinY, sLayerEnvelop.MinY);                sEnvelop.MaxX = MAX(sEnvelop.MaxX, sLayerEnvelop.MaxX);                sEnvelop.MaxY = MAX(sEnvelop.MaxY, sLayerEnvelop.MaxY);            }        }        else        {            if (bFirstLayer)            {                if (hSRS == NULL)                    hSRS = OGR_L_GetSpatialRef(hLayer);                bFirstLayer = FALSE;            }        }    }    if (dfXRes == 0 && dfYRes == 0)    {        dfXRes = (sEnvelop.MaxX - sEnvelop.MinX) / nXSize;        dfYRes = (sEnvelop.MaxY - sEnvelop.MinY) / nYSize;    }    else if (bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)    {        sEnvelop.MinX = floor(sEnvelop.MinX / dfXRes) * dfXRes;        sEnvelop.MaxX = ceil(sEnvelop.MaxX / dfXRes) * dfXRes;        sEnvelop.MinY = floor(sEnvelop.MinY / dfYRes) * dfYRes;        sEnvelop.MaxY = ceil(sEnvelop.MaxY / dfYRes) * dfYRes;    }    double adfProjection[6];    adfProjection[0] = sEnvelop.MinX;    adfProjection[1] = dfXRes;    adfProjection[2] = 0;    adfProjection[3] = sEnvelop.MaxY;    adfProjection[4] = 0;    adfProjection[5] = -dfYRes;    if (nXSize == 0 && nYSize == 0)    {        nXSize = (int)(0.5 + (sEnvelop.MaxX - sEnvelop.MinX) / dfXRes);        nYSize = (int)(0.5 + (sEnvelop.MaxY - sEnvelop.MinY) / dfYRes);    }    hDstDS = GDALCreate(hDriver, pszDest, nXSize, nYSize,                        nBandCount, eOutputType, papszCreationOptions);    if (hDstDS == NULL)    {        CPLError(CE_Failure, CPLE_AppDefined, "Cannot create %s", pszDest);        return NULL;    }    GDALSetGeoTransform(hDstDS, adfProjection);    if (hSRS)        OSRExportToWkt(hSRS, &pszWKT);    if (pszWKT)        GDALSetProjection(hDstDS, pszWKT);    CPLFree(pszWKT);    int iBand;    /*if( nBandCount == 3 || nBandCount == 4 )    {        for(iBand = 0; iBand < nBandCount; iBand++)        {            GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);            GDALSetRasterColorInterpretation(hBand, (GDALColorInterp)(GCI_RedBand + iBand));        }    }*/    if (bNoDataSet)    {        for(iBand = 0; iBand < nBandCount; iBand++)        {            GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);            GDALSetRasterNoDataValue(hBand, dfNoData);        }    }    if (adfInitVals.size() != 0)    {        for(iBand = 0; iBand < MIN(nBandCount,(int)adfInitVals.size()); iBand++)        {            GDALRasterBandH hBand = GDALGetRasterBand(hDstDS, iBand + 1);            GDALFillRaster(hBand, adfInitVals[iBand], 0);        }    }    return hDstDS;}

GByte *gdal_to_rgba( GDALDatasetH hDS ){    int  nXSize, nYSize;    GByte *pabyRGBABuf = NULL;    /* validation of input parameters */    g_return_val_if_fail( hDS != NULL, NULL );/* -------------------------------------------------------------------- *//*      Allocate RGBA Raster buffer.                                    *//* -------------------------------------------------------------------- */    nXSize = GDALGetRasterXSize( hDS );    nYSize = GDALGetRasterYSize( hDS );    CPLDebug( "OpenEV", "creating buffer of (%d,%d)", nXSize, nYSize );    pabyRGBABuf = (GByte *) CPLMalloc( nXSize * nYSize * 4 );/* -------------------------------------------------------------------- *//*      Handle case where source is already presumed to be RGBA.        *//* -------------------------------------------------------------------- */    if( GDALGetRasterCount(hDS) == 4 )    {        GDALRasterIO( GDALGetRasterBand( hDS, 1 ), GF_Read,                       0, 0, nXSize, nYSize,                       pabyRGBABuf+0, nXSize, nYSize, GDT_Byte,                       4, nXSize * 4 );                              GDALRasterIO( GDALGetRasterBand( hDS, 2 ), GF_Read,                       0, 0, nXSize, nYSize,                       pabyRGBABuf+1, nXSize, nYSize, GDT_Byte, 4,                       nXSize * 4 );                              GDALRasterIO( GDALGetRasterBand( hDS, 3 ), GF_Read,                       0, 0, nXSize, nYSize,                       pabyRGBABuf+2, nXSize, nYSize, GDT_Byte, 4,                       nXSize * 4 );                              GDALRasterIO( GDALGetRasterBand( hDS, 4 ), GF_Read,                       0, 0, nXSize, nYSize,                       pabyRGBABuf+3, nXSize, nYSize, GDT_Byte, 4,                       nXSize * 4 );    }/* -------------------------------------------------------------------- *//*      Source is RGB.  Set Alpha to 255.                               *//* -------------------------------------------------------------------- */    else if( GDALGetRasterCount(hDS) == 3 )    {        memset( pabyRGBABuf, 255, 4 * nXSize * nYSize );                GDALRasterIO( GDALGetRasterBand( hDS, 1 ), GF_Read,                       0, 0, nXSize, nYSize,                       pabyRGBABuf+0, nXSize, nYSize, GDT_Byte,                       4, nXSize * 4 );                              GDALRasterIO( GDALGetRasterBand( hDS, 2 ), GF_Read,                       0, 0, nXSize, nYSize,                       pabyRGBABuf+1, nXSize, nYSize, GDT_Byte, 4,                       nXSize * 4 );                              GDALRasterIO( GDALGetRasterBand( hDS, 3 ), GF_Read,                       0, 0, nXSize, nYSize,                       pabyRGBABuf+2, nXSize, nYSize, GDT_Byte, 4,                       nXSize * 4 );    }/* -------------------------------------------------------------------- *//*      Source is pseudocolored.  Load and then convert to RGBA.        *//* -------------------------------------------------------------------- */    else if( GDALGetRasterCount(hDS) == 1              && GDALGetRasterColorTable( GDALGetRasterBand( hDS, 1 )) != NULL )    {        int        i;        GDALColorTableH hTable;        GByte      abyPCT[1024];        /* Load color table, and produce 256 entry table to RGBA. */        hTable = GDALGetRasterColorTable( GDALGetRasterBand( hDS, 1 ) );        for( i = 0; i < MIN(256,GDALGetColorEntryCount( hTable )); i++ )        {            GDALColorEntry sEntry;            GDALGetColorEntryAsRGB( hTable, i, &sEntry );            abyPCT[i*4+0] = sEntry.c1;            abyPCT[i*4+1] = sEntry.c2;            abyPCT[i*4+2] = sEntry.c3;            abyPCT[i*4+3] = sEntry.c4;        }        /* Fill in any missing colors with greyscale. */        for( i = GDALGetColorEntryCount( hTable ); i < 256; i++ )        {            abyPCT[i*4+0] = i;            abyPCT[i*4+1] = i;            abyPCT[i*4+2] = i;            abyPCT[i*4+3] = 255;        }        /* Read indexed raster */        GDALRasterIO( GDALGetRasterBand( hDS, 1 ), GF_Read, //.........这里部分代码省略.........

//.........这里部分代码省略.........    ctx->set_error(ctx,MAPCACHE_SOURCE_GDAL_ERROR,"gdal failed to create suggested warp output");    return;  }  GDALDestroyGenImgProjTransformer( hTransformArg );  double dfXRes = (tile->bbox[2] - tile->bbox[0]) / tile->sx;  double dfYRes = (tile->bbox[3] - tile->bbox[1]) / tile->sy;  adfDstGeoTransform[0] = tile->bbox[0];  adfDstGeoTransform[3] = tile->bbox[3];  adfDstGeoTransform[1] = dfXRes;  adfDstGeoTransform[5] = -dfYRes;  hDstDS = GDALCreate( hDriver, "tempd_gdal_image", tile->sx, tile->sy, 4, GDT_Byte, NULL );  /* -------------------------------------------------------------------- */  /*      Write out the projection definition.                            */  /* -------------------------------------------------------------------- */  GDALSetProjection( hDstDS, dstSRS );  GDALSetGeoTransform( hDstDS, adfDstGeoTransform );  char               **papszWarpOptions = NULL;  papszWarpOptions = CSLSetNameValue( papszWarpOptions, "INIT", "0" );  /* -------------------------------------------------------------------- */  /*      Create a transformation object from the source to               */  /*      destination coordinate system.                                  */  /* -------------------------------------------------------------------- */  GDALTransformerFunc pfnTransformer = NULL;  void               *hGenImgProjArg=NULL, *hApproxArg=NULL;  hTransformArg = hGenImgProjArg =                    GDALCreateGenImgProjTransformer( hDataset, srcSRS,                        hDstDS, dstSRS,                        TRUE, 1000.0, 0 );  if( hTransformArg == NULL )    exit( 1 );  pfnTransformer = GDALGenImgProjTransform;  hTransformArg = hApproxArg =                    GDALCreateApproxTransformer( GDALGenImgProjTransform,                        hGenImgProjArg, 0.125 );  pfnTransformer = GDALApproxTransform;  /* -------------------------------------------------------------------- */  /*      Now actually invoke the warper to do the work.                  */  /* -------------------------------------------------------------------- */  GDALSimpleImageWarp( hDataset, hDstDS, 0, NULL,                       pfnTransformer, hTransformArg,                       GDALDummyProgress, NULL, papszWarpOptions );  CSLDestroy( papszWarpOptions );  if( hApproxArg != NULL )    GDALDestroyApproxTransformer( hApproxArg );  if( hGenImgProjArg != NULL )    GDALDestroyGenImgProjTransformer( hGenImgProjArg );  if(GDALGetRasterCount(hDstDS) != 4) {    ctx->set_error(ctx,MAPCACHE_SOURCE_GDAL_ERROR,"gdal did not create a 4 band image");    return;  }  GDALRasterBandH *redband, *greenband, *blueband, *alphaband;  redband = GDALGetRasterBand(hDstDS,1);  greenband = GDALGetRasterBand(hDstDS,2);  blueband = GDALGetRasterBand(hDstDS,3);  alphaband = GDALGetRasterBand(hDstDS,4);  unsigned char *rasterdata = apr_palloc(ctx->pool,tile->sx*tile->sy*4);  data->buf = rasterdata;  data->avail = tile->sx*tile->sy*4;  data->size = tile->sx*tile->sy*4;  GDALRasterIO(redband,GF_Read,0,0,tile->sx,tile->sy,(void*)(rasterdata),tile->sx,tile->sy,GDT_Byte,4,4*tile->sx);  GDALRasterIO(greenband,GF_Read,0,0,tile->sx,tile->sy,(void*)(rasterdata+1),tile->sx,tile->sy,GDT_Byte,4,4*tile->sx);  GDALRasterIO(blueband,GF_Read,0,0,tile->sx,tile->sy,(void*)(rasterdata+2),tile->sx,tile->sy,GDT_Byte,4,4*tile->sx);  if(GDALGetRasterCount(hDataset)==4)    GDALRasterIO(alphaband,GF_Read,0,0,tile->sx,tile->sy,(void*)(rasterdata+3),tile->sx,tile->sy,GDT_Byte,4,4*tile->sx);  else {    unsigned char *alphaptr;    int i;    for(alphaptr = rasterdata+3, i=0; i<tile->sx*tile->sy; i++, alphaptr+=4) {      *alphaptr = 255;    }  }  tile->imdata = mapcache_image_create(ctx);  tile->imdata->w = tile->sx;  tile->imdata->h = tile->sy;  tile->imdata->stride = tile->sx * 4;  tile->imdata->data = rasterdata;  GDALClose( hDstDS );  GDALClose( hDataset);}