C++ GDALRasterBand类代码示例

CPLErr VRTSourcedRasterBand::ComputeRasterMinMax( int bApproxOK, double* adfMinMax )
{
    double  dfMin = 0.0;
    double  dfMax = 0.0;

/* -------------------------------------------------------------------- */
/*      Does the driver already know the min/max?                       */
/* -------------------------------------------------------------------- */
    if( bApproxOK )
    {
        int          bSuccessMin, bSuccessMax;

        dfMin = GetMinimum( &bSuccessMin );
        dfMax = GetMaximum( &bSuccessMax );

        if( bSuccessMin && bSuccessMax )
        {
            adfMinMax[0] = dfMin;
            adfMinMax[1] = dfMax;
            return CE_None;
        }
    }
    
/* -------------------------------------------------------------------- */
/*      If we have overview bands, use them for min/max.                */
/* -------------------------------------------------------------------- */
    if ( bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews() )
    {
        GDALRasterBand *poBand;

        poBand = GetRasterSampleOverview( GDALSTAT_APPROX_NUMSAMPLES );

        if ( poBand != this )
            return poBand->ComputeRasterMinMax( FALSE, adfMinMax );
    }
    
/* -------------------------------------------------------------------- */
/*      Try with source bands.                                          */
/* -------------------------------------------------------------------- */
    if ( bAntiRecursionFlag )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "VRTSourcedRasterBand::ComputeRasterMinMax() called recursively on the same band. "
                  "It looks like the VRT is referencing itself." );
        return CE_Failure;
    }
    bAntiRecursionFlag = TRUE;

    adfMinMax[0] = 0.0;
    adfMinMax[1] = 0.0;
    for( int iSource = 0; iSource < nSources; iSource++ )
    {
        double adfSourceMinMax[2];
        CPLErr eErr = papoSources[iSource]->ComputeRasterMinMax(GetXSize(), GetYSize(), bApproxOK, adfSourceMinMax);
        if (eErr != CE_None)
        {
            eErr = GDALRasterBand::ComputeRasterMinMax(bApproxOK, adfMinMax);
            bAntiRecursionFlag = FALSE;
            return eErr;
        }

        if (iSource == 0 || adfSourceMinMax[0] < adfMinMax[0])
            adfMinMax[0] = adfSourceMinMax[0];
        if (iSource == 0 || adfSourceMinMax[1] > adfMinMax[1])
            adfMinMax[1] = adfSourceMinMax[1];
    }

    bAntiRecursionFlag = FALSE;

    return CE_None;
}

/** Private method to calculate statistics for each band. Populates rasterStatsMap. */
void ImageWriter::calculateStats(RasterBandStats * theRasterBandStats,GDALDataset * gdalDataset)
{
        std::cout << "Calculating statistics..." << std::endl;
	GDALRasterBand  *myGdalBand = gdalDataset->GetRasterBand( 1 );
	QString myColorInterpretation = GDALGetColorInterpretationName(myGdalBand->GetColorInterpretation());
	theRasterBandStats->bandName=myColorInterpretation;
	theRasterBandStats->bandNo=1;
	// get the dimensions of the raster
	int myColsInt = myGdalBand->GetXSize();
	int myRowsInt = myGdalBand->GetYSize();

	theRasterBandStats->elementCountInt=myColsInt*myRowsInt;
	theRasterBandStats->noDataDouble=myGdalBand->GetNoDataValue();
	//allocate a buffer to hold one row of ints
	int myAllocationSizeInt = sizeof(uint)*myColsInt;
	uint * myScanlineAllocInt = (uint*) CPLMalloc(myAllocationSizeInt);
	bool myFirstIterationFlag = true;
	//unfortunately we need to make two passes through the data to calculate stddev
	for (int myCurrentRowInt=0; myCurrentRowInt < myRowsInt;myCurrentRowInt++)
	{
		CPLErr myResult = myGdalBand->RasterIO(
				GF_Read, 0, myCurrentRowInt, myColsInt, 1, myScanlineAllocInt, myColsInt, 1, GDT_UInt32, 0, 0 );
		for (int myCurrentColInt=0; myCurrentColInt < myColsInt; myCurrentColInt++)
		{
			//get the nth element from the current row
			double myDouble=myScanlineAllocInt[myCurrentColInt];
			//only use this element if we have a non null element
			if (myDouble != theRasterBandStats->noDataDouble )
			{
				if (myFirstIterationFlag)
				{
					//this is the first iteration so initialise vars
					myFirstIterationFlag=false;
					theRasterBandStats->minValDouble=myDouble;
					theRasterBandStats->maxValDouble=myDouble;
				} //end of true part for first iteration check
				else
				{
					//this is done for all subsequent iterations
					if (myDouble < theRasterBandStats->minValDouble)
					{
						theRasterBandStats->minValDouble=myDouble;
					}
					if (myDouble > theRasterBandStats->maxValDouble)
					{
					//	printf ("Maxval updated to %f\n",myDouble);
						theRasterBandStats->maxValDouble=myDouble;
					}
					//only increment the running total if it is not a nodata value
					if (myDouble != theRasterBandStats->noDataDouble)
					{
						theRasterBandStats->sumDouble += myDouble;
						++theRasterBandStats->elementCountInt;
					}
				} //end of false part for first iteration check
			} //end of nodata chec
		} //end of column wise loop
	} //end of row wise loop
	//
	//end of first pass through data now calculate the range
	theRasterBandStats->rangeDouble = theRasterBandStats->maxValDouble-theRasterBandStats->minValDouble;
	//calculate the mean
	theRasterBandStats->meanDouble = theRasterBandStats->sumDouble / theRasterBandStats->elementCountInt;
	//for the second pass we will get the sum of the squares / mean
	for (int myCurrentRowInt=0; myCurrentRowInt < myRowsInt;myCurrentRowInt++)
	{
		CPLErr myResult = myGdalBand->RasterIO(
				GF_Read, 0, myCurrentRowInt, myColsInt, 1, myScanlineAllocInt, myColsInt, 1, GDT_UInt32, 0, 0 );
		for (int myCurrentColInt=0; myCurrentColInt < myColsInt; myCurrentColInt++)
		{
			//get the nth element from the current row
			double myDouble=myScanlineAllocInt[myCurrentColInt];
			theRasterBandStats->sumSqrDevDouble += static_cast<double>(pow(myDouble - theRasterBandStats->meanDouble,2));
		} //end of column wise loop
	} //end of row wise loop
	//divide result by sample size - 1 and get square root to get stdev
	theRasterBandStats->stdDevDouble = static_cast<double>(sqrt(theRasterBandStats->sumSqrDevDouble /
				(theRasterBandStats->elementCountInt - 1)));
	CPLFree(myScanlineAllocInt);
	//printf("CalculateStats::\n");
	//std::cout << "Band Name   : " << theRasterBandStats->bandName << std::endl;
	//printf("Band No   : %i\n",theRasterBandStats->bandNo);
	//printf("Band min  : %f\n",theRasterBandStats->minValDouble);
	//printf("Band max  : %f\n",theRasterBandStats->maxValDouble);
	//printf("Band range: %f\n",theRasterBandStats->rangeDouble);
	//printf("Band mean : %f\n",theRasterBandStats->meanDouble);
	//printf("Band sum : %f\n",theRasterBandStats->sumDouble);
	return ;
}

GDALDataset *GS7BGDataset::CreateCopy( const char *pszFilename,
                                       GDALDataset *poSrcDS,
                                       int bStrict,
                                       CPL_UNUSED char **papszOptions,
                                       GDALProgressFunc pfnProgress,
                                       void *pProgressData )
{
    if( pfnProgress == NULL )
        pfnProgress = GDALDummyProgress;

    int nBands = poSrcDS->GetRasterCount();
    if (nBands == 0)
    {
        CPLError( CE_Failure, CPLE_NotSupported,
                  "Driver does not support source dataset with zero band.\n");
        return NULL;
    }
    else if (nBands > 1)
    {
        if( bStrict )
        {
            CPLError( CE_Failure, CPLE_NotSupported,
                "Unable to create copy, "
                "format only supports one raster band.\n" );
            return NULL;
        }
        else
            CPLError( CE_Warning, CPLE_NotSupported,
                "Format only supports one "
                "raster band, first band will be copied.\n" );
    }

    GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( 1 );

    if( !pfnProgress( 0.0, NULL, pProgressData ) )
    {
        CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated\n" );
        return NULL;
    }

    VSILFILE    *fp = VSIFOpenL( pszFilename, "w+b" );

    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed,
                  "Attempt to create file '%s' failed.\n",
                  pszFilename );
        return NULL;
    }

    GInt32  nXSize = poSrcBand->GetXSize();
    GInt32  nYSize = poSrcBand->GetYSize();
    double  adfGeoTransform[6];

    poSrcDS->GetGeoTransform( adfGeoTransform );

    double dfMinX = adfGeoTransform[0] + adfGeoTransform[1] / 2;
    double dfMaxX = adfGeoTransform[1] * (nXSize - 0.5) + adfGeoTransform[0];
    double dfMinY = adfGeoTransform[5] * (nYSize - 0.5) + adfGeoTransform[3];
    double dfMaxY = adfGeoTransform[3] + adfGeoTransform[5] / 2;
    CPLErr eErr = WriteHeader( fp, nXSize, nYSize,
                   dfMinX, dfMaxX, dfMinY, dfMaxY, 0.0, 0.0 );

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

/* -------------------------------------------------------------------- */
/*      Copy band data.                                                 */
/* -------------------------------------------------------------------- */
    double *pfData = (double *)VSI_MALLOC2_VERBOSE( nXSize, sizeof( double ) );
    if( pfData == NULL )
    {
        VSIFCloseL( fp );
        return NULL;
    }

    int     bSrcHasNDValue;
    double   dfSrcNoDataValue = poSrcBand->GetNoDataValue( &bSrcHasNDValue );
    double  dfMinZ = DBL_MAX;
    double  dfMaxZ = -DBL_MAX;
    for( GInt32 iRow = nYSize - 1; iRow >= 0; iRow-- )
    {
        eErr = poSrcBand->RasterIO( GF_Read, 0, iRow,
                    nXSize, 1, pfData,
                    nXSize, 1, GDT_Float64, 0, 0, NULL );

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

        for( int iCol=0; iCol<nXSize; iCol++ )
        {
            if( bSrcHasNDValue && pfData[iCol] == dfSrcNoDataValue )
            {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if( CSLFetchNameValue(papszOptions,"FILLUNDEFVAL") != NULL ) {
        fillundefval = atoi(CSLFetchNameValue(papszOptions,"FILLUNDEFVAL"));
        if( (fillundefval < -32768) || (fillundefval > 32767) ) {
            CPLError( CE_Failure, CPLE_IllegalArg,
                      "FILLUNDEFVAL=%s is not a legal value in the range -32768, 32767.",
                      CSLFetchNameValue(papszOptions,"FILLUNDEFVAL") );
            return NULL;
        }
    }

    int endian = LITTLEENDIAN;
    if( CSLFetchNameValue(papszOptions,"BIGENDIAN") != NULL
        && !EQUAL(CSLFetchNameValue(papszOptions,"BIGENDIAN"),"NO") )
        endian = BIGENDIAN;

// --------------------------------------------------------------------
//      Create the dataset.
// --------------------------------------------------------------------

    // Create a BLX context
    blxcontext_t *ctx = blx_create_context();

    // Setup BLX parameters
    ctx->cell_rows = nYSize / ctx->cell_ysize;
    ctx->cell_cols = nXSize / ctx->cell_xsize;
    ctx->zscale = zscale;
    ctx->fillundef = fillundef;
    ctx->fillundefval = fillundefval;
    ctx->endian = endian;

    if(blxopen(ctx, pszFilename, "wb")) {
        CPLError( CE_Failure, CPLE_OpenFailed,
                  "Unable to create blx file %s.\n",
                  pszFilename );
        blx_free_context(ctx);
        return NULL;
    }

// --------------------------------------------------------------------
//      Loop over image, copying image data.
// --------------------------------------------------------------------

    GInt16 *pabyTile
        = (GInt16 *) VSI_MALLOC_VERBOSE( sizeof(GInt16)*ctx->cell_xsize*ctx->cell_ysize );
    if (pabyTile == NULL)
    {
        blxclose(ctx);
        blx_free_context(ctx);
        return NULL;
    }

    CPLErr eErr=CE_None;
    if( !pfnProgress( 0.0, NULL, pProgressData ) )
        eErr = CE_Failure;

    for(int i=0; (i < ctx->cell_rows) && (eErr == CE_None); i++)
        for(int j=0; j < ctx->cell_cols; j++) {
            blxdata *celldata;
            GDALRasterBand * poBand = poSrcDS->GetRasterBand( 1 );
            eErr = poBand->RasterIO( GF_Read, j*ctx->cell_xsize, i*ctx->cell_ysize,
                                     ctx->cell_xsize, ctx->cell_ysize,
                                     pabyTile, ctx->cell_xsize, ctx->cell_ysize, GDT_Int16,
                                     0, 0, NULL );
            if(eErr >= CE_Failure)
                 break;
            celldata = pabyTile;
            if (blx_writecell(ctx, celldata, i, j) != 0)
            {
                eErr = CE_Failure;
                break;
            }

            if ( ! pfnProgress( 1.0 * (i * ctx->cell_cols + j) / (ctx->cell_rows * ctx->cell_cols), NULL, pProgressData ))
            {
                eErr = CE_Failure;
                break;
            }
    }

    pfnProgress( 1.0, NULL, pProgressData );

    CPLFree( pabyTile );

    double adfGeoTransform[6];
    if( poSrcDS->GetGeoTransform( adfGeoTransform ) == CE_None )
    {
        ctx->lon = adfGeoTransform[0];
        ctx->lat = adfGeoTransform[3];
        ctx->pixelsize_lon = adfGeoTransform[1];
        ctx->pixelsize_lat = adfGeoTransform[5];
    }

    blxclose(ctx);
    blx_free_context(ctx);

    if (eErr == CE_None)
        return reinterpret_cast<GDALDataset *>( GDALOpen( pszFilename, GA_ReadOnly ) );

    return NULL;
}

GDALDataset *
GIFDataset::CreateCopy( 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();
    int	 bInterlace = FALSE;

/* -------------------------------------------------------------------- */
/*      Check for interlaced option.                                    */
/* -------------------------------------------------------------------- */
    bInterlace = CSLFetchBoolean(papszOptions, "INTERLACING", FALSE);

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

        return NULL;
    }

    if (nXSize > 65535 || nYSize > 65535)
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "GIF driver only supports datasets up to 65535x65535 size.\n" );

        return NULL;
    }

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

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Open the output file.                                           */
/* -------------------------------------------------------------------- */
    GifFileType *hGifFile;
    VSILFILE *fp;

    fp = VSIFOpenL( pszFilename, "wb" );
    if( fp == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "Failed to create %s:\n%s", 
                  pszFilename, VSIStrerror( errno ) );
        return NULL;
    }

#if defined(GIFLIB_MAJOR) && GIFLIB_MAJOR >= 5
    int nError;
    hGifFile = EGifOpen( fp, VSIGIFWriteFunc, &nError );
#else
    hGifFile = EGifOpen( fp, VSIGIFWriteFunc );
#endif
    if( hGifFile == NULL )
    {
        VSIFCloseL( fp );
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "EGifOpenFilename(%s) failed.  Does file already exist?",
                  pszFilename );

        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Prepare colortable.                                             */
/* -------------------------------------------------------------------- */
    GDALRasterBand	*poBand = poSrcDS->GetRasterBand(1);
    ColorMapObject	*psGifCT;
    int			iColor;

    if( poBand->GetColorTable() == NULL )
    {
        psGifCT = GifMakeMapObject( 256, NULL );
        for( iColor = 0; iColor < 256; iColor++ )
        {
            psGifCT->Colors[iColor].Red = (GifByteType) iColor;
            psGifCT->Colors[iColor].Green = (GifByteType) iColor;
            psGifCT->Colors[iColor].Blue = (GifByteType) iColor;
        }
    }
    else
    {
        GDALColorTable	*poCT = poBand->GetColorTable();
        int nFullCount = 1;

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

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;

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

CPLErr GDALMRFRasterBand::IWriteBlock(int xblk, int yblk, void *buffer)

{
    GInt32 cstride = img.pagesize.c;
    ILSize req(xblk, yblk, 0, m_band/cstride, m_l);
    GUIntBig infooffset = IdxOffset(req, img);

    CPLDebug("MRF_IB", "IWriteBlock %d,%d,0,%d, level  %d, stride %d\n", xblk, yblk, 
	m_band, m_l, cstride);

    if (1 == cstride) {     // Separate bands, we can write it as is
	// Empty page skip

	int success;
	double val = GetNoDataValue(&success);
	if (!success) val = 0.0;
	if (isAllVal(eDataType, buffer, img.pageSizeBytes, val))
	    return poDS->WriteTile(NULL, infooffset, 0);

	// Use the pbuffer to hold the compressed page before writing it
	poDS->tile = ILSize(); // Mark it corrupt

	buf_mgr src;
        src.buffer = (char *)buffer;
        src.size = static_cast<size_t>(img.pageSizeBytes);
	buf_mgr dst = {(char *)poDS->GetPBuffer(), poDS->GetPBufferSize()};

	// Swab the source before encoding if we need to 
	if (is_Endianess_Dependent(img.dt, img.comp) && (img.nbo != NET_ORDER)) 
	    swab_buff(src, img);

	// Compress functions need to return the compressed size in
	// the bytes in buffer field
	Compress(dst, src);
	void *usebuff = dst.buffer;
	if (deflatep) {
	    usebuff = DeflateBlock(dst, poDS->pbsize - dst.size, deflate_flags);
	    if (!usebuff) {
		CPLError(CE_Failure,CPLE_AppDefined, "MRF: Deflate error");
		return CE_Failure;
	    }
	}
	return poDS->WriteTile(usebuff, infooffset , dst.size);
    }

    // Multiple bands per page, use a temporary to assemble the page
    // Temporary is large because we use it to hold both the uncompressed and the compressed
    poDS->tile=req; poDS->bdirty=0;

    // Keep track of what bands are empty
    GUIntBig empties=0;

    void *tbuffer = VSIMalloc(img.pageSizeBytes + poDS->pbsize);

    if (!tbuffer) {
	CPLError(CE_Failure,CPLE_AppDefined, "MRF: Can't allocate write buffer");
	return CE_Failure;
    }
	
    // Get the other bands from the block cache
    for (int iBand=0; iBand < poDS->nBands; iBand++ )
    {
	const char *pabyThisImage=NULL;
	GDALRasterBlock *poBlock=NULL;

	if (iBand == m_band)
	{
	    pabyThisImage = (char *) buffer;
	    poDS->bdirty |= bandbit();
	} else {
	    GDALRasterBand *band = poDS->GetRasterBand(iBand +1);
	    // Pick the right overview
	    if (m_l) band = band->GetOverview(m_l -1);
	    poBlock = ((GDALMRFRasterBand *)band)
		->TryGetLockedBlockRef(xblk, yblk);
	    if (NULL==poBlock) continue;
	    // This is where the image data is for this band

	    pabyThisImage = (char*) poBlock->GetDataRef();
	    poDS->bdirty |= bandbit(iBand);
	}

	// Keep track of empty bands, but encode them anyhow, in case some are not empty
	int success;
	double val = GetNoDataValue(&success);
	if (!success) val = 0.0;
        if (isAllVal(eDataType, (char *)pabyThisImage, blockSizeBytes(), val))
	    empties |= bandbit(iBand);

	// Copy the data into the dataset buffer here
	// Just the right mix of templates and macros make this real tidy
#define CpySO(T) cpy_stride_out<T> (((T *)tbuffer)+iBand, pabyThisImage,\
		blockSizeBytes()/sizeof(T), cstride)


	// Build the page in tbuffer
	switch (GDALGetDataTypeSize(eDataType)/8)
	{
	    case 1: CpySO(GByte); break;
	    case 2: CpySO(GInt16); break;
//.........这里部分代码省略.........

GDALDataset *IntergraphDataset::CreateCopy( const char *pszFilename,
        GDALDataset *poSrcDS,
        int bStrict,
        char **papszOptions,
        GDALProgressFunc pfnProgress,
        void *pProgressData )
{
    (void) bStrict;

    int nBands = poSrcDS->GetRasterCount();
    if (nBands == 0)
    {
        CPLError( CE_Failure, CPLE_NotSupported,
                  "Intergraph driver does not support source dataset with zero band.\n");
        return NULL;
    }

    if( !pfnProgress( 0.0, NULL, pProgressData ) )
    {
        return NULL;
    }

    // --------------------------------------------------------------------
    // Query GDAL Data Type
    // --------------------------------------------------------------------

    GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();

    // --------------------------------------------------------------------
    // Copy metadata
    // --------------------------------------------------------------------

    char **papszCreateOptions = CSLDuplicate( papszOptions );
    const char  *pszValue;

    pszValue = CSLFetchNameValue(papszCreateOptions, "RESOLUTION");
    if( pszValue == NULL )
    {
        const char *value = poSrcDS->GetMetadataItem("RESOLUTION");
        if (value)
        {
            papszCreateOptions = CSLSetNameValue( papszCreateOptions, "RESOLUTION",
                                                  value );
        }
    }

    // --------------------------------------------------------------------
    // Create IntergraphDataset
    // --------------------------------------------------------------------

    IntergraphDataset *poDstDS;

    poDstDS = (IntergraphDataset*) IntergraphDataset::Create( pszFilename,
              poSrcDS->GetRasterXSize(),
              poSrcDS->GetRasterYSize(),
              poSrcDS->GetRasterCount(),
              eType,
              papszCreateOptions );

    CSLDestroy( papszCreateOptions );

    if( poDstDS == NULL )
    {
        return NULL;
    }

    // --------------------------------------------------------------------
    // Copy Transformation Matrix to the dataset
    // --------------------------------------------------------------------

    double adfGeoTransform[6];

    poDstDS->SetProjection( poSrcDS->GetProjectionRef() );
    poSrcDS->GetGeoTransform( adfGeoTransform );
    poDstDS->SetGeoTransform( adfGeoTransform );

    // --------------------------------------------------------------------
    // Copy information to the raster band
    // --------------------------------------------------------------------

    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;
    }
//.........这里部分代码省略.........

static std::vector<GDALFeaturePoint> *
GatherFeaturePoints( GDALDataset* poDataset, int* panBands,
                     int nOctaveStart, int nOctaveEnd, double dfThreshold )
{
    if( poDataset == nullptr )
    {
        CPLError(CE_Failure, CPLE_AppDefined, "GDALDataset isn't specified");
        return nullptr;
    }

    if( panBands == nullptr )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Raster bands are not specified");
        return nullptr;
    }

    if( nOctaveStart <= 0 || nOctaveEnd < 0 ||
        nOctaveStart > nOctaveEnd )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Octave numbers are invalid");
        return nullptr;
    }

    if( dfThreshold < 0 )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Threshold have to be greater than zero");
        return nullptr;
    }

    GDALRasterBand *poRstRedBand = poDataset->GetRasterBand(panBands[0]);
    GDALRasterBand *poRstGreenBand = poDataset->GetRasterBand(panBands[1]);
    GDALRasterBand *poRstBlueBand = poDataset->GetRasterBand(panBands[2]);

    const int nWidth = poRstRedBand->GetXSize();
    const int nHeight = poRstRedBand->GetYSize();

    if( nWidth == 0 || nHeight == 0 )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Must have non-zero width and height.");
        return nullptr;
    }

    // Allocate memory for grayscale image.
    double **padfImg = new double*[nHeight];
    for( int i = 0; ; )
    {
        padfImg[i] = new double[nWidth];
        for( int j = 0; j < nWidth; ++j )
          padfImg[i][j] = 0.0;
        ++i;
        if( i == nHeight )
            break;
    }

    // Create grayscale image.
    GDALSimpleSURF::ConvertRGBToLuminosity(
        poRstRedBand, poRstGreenBand, poRstBlueBand, nWidth, nHeight,
        padfImg, nHeight, nWidth);

    // Prepare integral image.
    GDALIntegralImage *poImg = new GDALIntegralImage();
    poImg->Initialize(const_cast<const double**>(padfImg), nHeight, nWidth);

    // Get feature points.
    GDALSimpleSURF *poSurf = new GDALSimpleSURF(nOctaveStart, nOctaveEnd);

    std::vector<GDALFeaturePoint> *poCollection =
        poSurf->ExtractFeaturePoints(poImg, dfThreshold);

    // Clean up.
    delete poImg;
    delete poSurf;

    for( int i = 0; i < nHeight; ++i )
        delete[] padfImg[i];

    delete[] padfImg;

    return poCollection;
}

//.........这里部分代码省略.........
    GDALDataset *srcDS, *wrpDS;
    std::string temp;
    std::string srcWkt;

    std::vector<std::string> varList = getVariableList();

    /*
     * Set the initial values in the warped dataset to no data
     */
    GDALWarpOptions* psWarpOptions;

    for( unsigned int i = 0;i < varList.size();i++ ) {

        temp = "NETCDF:" + input.forecastFilename + ":" + varList[i];

        srcDS = (GDALDataset*)GDALOpenShared( temp.c_str(), GA_ReadOnly );
        if( srcDS == NULL ) {
            CPLDebug( "ncepNdfdInitialization::setSurfaceGrids()",
                    "Bad forecast file" );
        }

        srcWkt = srcDS->GetProjectionRef();

        if( srcWkt.empty() ) {
            CPLDebug( "ncepNdfdInitialization::setSurfaceGrids()",
                    "Bad forecast file" );
            //throw
        }

        /*
         * Grab the first band to get the nodata value for the variable,
         * assume all bands have the same ndv
         */
        GDALRasterBand *poBand = srcDS->GetRasterBand( 1 );
        int pbSuccess;
        double dfNoData = poBand->GetNoDataValue( &pbSuccess );

        psWarpOptions = GDALCreateWarpOptions();

        int nBandCount = srcDS->GetRasterCount();

        psWarpOptions->nBandCount = nBandCount;

        psWarpOptions->padfDstNoDataReal =
            (double*) CPLMalloc( sizeof( double ) * nBandCount );
        psWarpOptions->padfDstNoDataImag =
            (double*) CPLMalloc( sizeof( double ) * nBandCount );

        for( int b = 0;b < srcDS->GetRasterCount();b++ ) {
            psWarpOptions->padfDstNoDataReal[b] = dfNoData;
            psWarpOptions->padfDstNoDataImag[b] = dfNoData;
        }

        if( pbSuccess == false )
            dfNoData = -9999.0;

        psWarpOptions->papszWarpOptions =
            CSLSetNameValue( psWarpOptions->papszWarpOptions,
                            "INIT_DEST", "NO_DATA" );

        wrpDS = (GDALDataset*) GDALAutoCreateWarpedVRT( srcDS, srcWkt.c_str(),
                                                        dstWkt.c_str(),
                                                        GRA_NearestNeighbour,
                                                        1.0, psWarpOptions );

        if( varList[i] == "Temperature_height_above_ground" ) {

//.........这里部分代码省略.........
    {
        CPLError( CE_Warning, CPLE_AppDefined, 
               "The horizontal source size is not conformant with the one "
               "expected by DTED Level %d at this latitude (%d pixels found instead of %d).", nLevel,
                poSrcDS->GetRasterXSize(), expectedXSize);
    }

/* -------------------------------------------------------------------- */
/*      Create the output dted file.                                    */
/* -------------------------------------------------------------------- */
    const char *pszError;

    pszError = DTEDCreate( pszFilename, nLevel, nLLOriginLat, nLLOriginLong );

    if( pszError != NULL )
    {
        CPLError( CE_Failure, CPLE_AppDefined, 
                  "%s", pszError );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Open the DTED file so we can output the data to it.             */
/* -------------------------------------------------------------------- */
    DTEDInfo *psDTED;

    psDTED = DTEDOpen( pszFilename, "rb+", FALSE );
    if( psDTED == NULL )
        return NULL;

/* -------------------------------------------------------------------- */
/*      Read all the data in a single buffer.                           */
/* -------------------------------------------------------------------- */
    GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( 1 );
    GInt16      *panData;

    panData = (GInt16 *) 
        VSIMalloc(sizeof(GInt16) * psDTED->nXSize * psDTED->nYSize);
    if (panData == NULL)
    {
        CPLError( CE_Failure, CPLE_OutOfMemory, "Out of memory");
        DTEDClose(psDTED);
        return NULL;
    }

    for( int iY = 0; iY < psDTED->nYSize; iY++ )
    {
        poSrcBand->RasterIO( GF_Read, 0, iY, psDTED->nXSize, 1, 
                            (void *) (panData + iY * psDTED->nXSize), psDTED->nXSize, 1, 
                            GDT_Int16, 0, 0 );

        if( pfnProgress && !pfnProgress(0.5 * (iY+1) / (double) psDTED->nYSize, NULL, pProgressData ) )
        {
            CPLError( CE_Failure, CPLE_UserInterrupt, 
                        "User terminated CreateCopy()" );
            DTEDClose( psDTED );
            CPLFree( panData );
            return NULL;
        }
    }

    int bSrcBandHasNoData;
    double srcBandNoData = poSrcBand->GetNoDataValue(&bSrcBandHasNoData);

/* -------------------------------------------------------------------- */
/*      Write all the profiles.                                         */

/**
* Checks the downloaded data to see if it is all valid.
*/
void genericSurfInitialization::checkForValidData()
{
    //just make up a "dummy" timezone for use here
    boost::local_time::time_zone_ptr zone(new boost::local_time::posix_time_zone("MST-07"));

    //get time list
    std::vector<boost::local_time::local_date_time> timeList( getTimeList(zone) );

    boost::posix_time::ptime pt_low(boost::gregorian::date(1900,boost::gregorian::Jan,1), boost::posix_time::hours(12));
    boost::posix_time::ptime pt_high(boost::gregorian::date(2100,boost::gregorian::Jan,1), boost::posix_time::hours(12));
    boost::local_time::local_date_time low_time(pt_low, zone);
    boost::local_time::local_date_time high_time(pt_high, zone);

    //check times
    for(unsigned int i = 0; i < timeList.size(); i++)
    {
        if(timeList[i].is_special())    //if time is any special value (not_a_date_time, infinity, etc.)
            throw badForecastFile("Bad time in forecast file.");
        if(timeList[i] < low_time || timeList[i] > high_time)
            throw badForecastFile("Bad time in forecast file.");
    }

    // open ds variable by variable
    GDALDataset *srcDS;
    std::string temp;
    std::string srcWkt;
    int nBands = 0;
    bool noDataValueExists;
    bool noDataIsNan;

    std::vector<std::string> varList = getVariableList();

    //Acquire a lock to protect the non-thread safe netCDF library
#ifdef _OPENMP
    omp_guard netCDF_guard(netCDF_lock);
#endif

    for( unsigned int i = 0;i < varList.size();i++ ) {

        temp = "NETCDF:" + wxModelFileName + ":" + varList[i];

        srcDS = (GDALDataset*)GDALOpen( temp.c_str(), GA_ReadOnly );
        if( srcDS == NULL )
            throw badForecastFile("Cannot open forecast file.");

        srcWkt = srcDS->GetProjectionRef();

        if( srcWkt.empty() )
            throw badForecastFile("Forecast file doesn't have projection information.");

        //Get total bands (time steps)
        nBands = srcDS->GetRasterCount();
        int nXSize, nYSize;
        GDALRasterBand *poBand;
        int pbSuccess;
        double dfNoData;
        double *padfScanline;

        nXSize = srcDS->GetRasterXSize();
        nYSize = srcDS->GetRasterYSize();

        //loop over all bands for this variable (bands are time steps)
        for(int j = 1; j <= nBands; j++)
        {
            poBand = srcDS->GetRasterBand( j );

            pbSuccess = 0;
            dfNoData = poBand->GetNoDataValue( &pbSuccess );
            if( pbSuccess == false )
                noDataValueExists = false;
            else
            {
                noDataValueExists = true;
                noDataIsNan = CPLIsNan(dfNoData);
            }

            //set the data
            padfScanline = new double[nXSize*nYSize];
            poBand->RasterIO(GF_Read, 0, 0, nXSize, nYSize, padfScanline, nXSize, nYSize,
                    GDT_Float64, 0, 0);
            for(int k = 0;k < nXSize*nYSize; k++)
            {
                //Check if value is no data (if no data value was defined in file)
                if(noDataValueExists)
                {
                    if(noDataIsNan)
                    {
                        if(CPLIsNan(padfScanline[k]))
                            throw badForecastFile("Forecast file contains no_data values.");
                    }else
                    {
                        if(padfScanline[k] == dfNoData)
                            throw badForecastFile("Forecast file contains no_data values.");
                    }
                }

                if( varList[i] == "Temperature_height_above_ground" )   //units are Kelvin
//.........这里部分代码省略.........

CPLErr VRTSourcedRasterBand::GetHistogram( double dfMin, double dfMax, 
                                     int nBuckets, int *panHistogram, 
                                     int bIncludeOutOfRange, int bApproxOK,
                                     GDALProgressFunc pfnProgress, 
                                     void *pProgressData )

{
    if( nSources != 1 )
        return GDALRasterBand::GetHistogram( dfMin, dfMax,
                                             nBuckets, panHistogram,
                                             bIncludeOutOfRange, bApproxOK,
                                             pfnProgress, pProgressData );

    if( pfnProgress == NULL )
        pfnProgress = GDALDummyProgress;

/* -------------------------------------------------------------------- */
/*      If we have overviews, use them for the histogram.               */
/* -------------------------------------------------------------------- */
    if( bApproxOK && GetOverviewCount() > 0 && !HasArbitraryOverviews() )
    {
        // FIXME: should we use the most reduced overview here or use some
        // minimum number of samples like GDALRasterBand::ComputeStatistics()
        // does?
        GDALRasterBand *poBestOverview = GetRasterSampleOverview( 0 );
        
        if( poBestOverview != this )
        {
            return poBestOverview->GetHistogram( dfMin, dfMax, nBuckets,
                                                 panHistogram,
                                                 bIncludeOutOfRange, bApproxOK,
                                                 pfnProgress, pProgressData );
        }
    }

/* -------------------------------------------------------------------- */
/*      Try with source bands.                                          */
/* -------------------------------------------------------------------- */
    if ( bAntiRecursionFlag )
    {
        CPLError( CE_Failure, CPLE_AppDefined,
                  "VRTSourcedRasterBand::GetHistogram() called recursively on the same band. "
                  "It looks like the VRT is referencing itself." );
        return CE_Failure;
    }
    bAntiRecursionFlag = TRUE;

    CPLErr eErr = papoSources[0]->GetHistogram(GetXSize(), GetYSize(), dfMin, dfMax, nBuckets,
                                               panHistogram,
                                               bIncludeOutOfRange, bApproxOK,
                                               pfnProgress, pProgressData);
    if (eErr != CE_None)
    {
        eErr = GDALRasterBand::GetHistogram( dfMin, dfMax,
                                                  nBuckets, panHistogram,
                                                  bIncludeOutOfRange, bApproxOK,
                                                  pfnProgress, pProgressData );
        bAntiRecursionFlag = FALSE;
        return eErr;
    }

    bAntiRecursionFlag = FALSE;

    return CE_None;
}

void output_geotiff ( rgb_image & out )
{
    int i, r, c;
    int val;
    char s[2];
    OGRSpatialReference ref;
    GDALDataset  *df;
    char *wkt = NULL;
    GDALRasterBand *bd;
    double        trans[6];
    GDALDriver   *gt;
    char         **options = NULL;
    int          ov[] = { 2, 4, 8, 16, 32 };
    int          nov;
    int          n;
    int          bands[] = { 1, 2, 3 };
    char         file[1000];
    int          block, ir, rows;

    options = CSLSetNameValue ( options, "TILED", "NO" );
    options = CSLSetNameValue ( options, "BLOCKXSIZE", "256" );
    options = CSLSetNameValue ( options, "BLOCKYSIZE", "256" );
    options = CSLSetNameValue ( options, "COMPRESS", "LZW" );

    gt = GetGDALDriverManager()->GetDriverByName("GTiff");
    if ( !gt ) {
        fprintf(stderr,"Could not get GTiff driver\n");
        exit(1);
    }

    strcpy ( file, p.value("output_file").c_str() );
    df = gt->Create( file, out.cols, out.rows, 3, GDT_Byte, options );
    if( df == NULL ) {
        fprintf(stderr,"Could not create %s\n", file );
        exit(1);
    }

    trans[0] = p.dvalue("easting_left");
    trans[1] = p.dvalue("output_cell_size");
    trans[2] = 0.0;
    trans[3] = p.dvalue("northing_top");
    trans[4] = 0.0;
    trans[5] = -p.dvalue("output_cell_size");
    df->SetGeoTransform ( trans );
    ref.SetUTM ( p.ivalue("utm_zone") );
    ref.SetWellKnownGeogCS ( "NAD27" );
    ref.exportToWkt ( &wkt );
    df->SetProjection(wkt);
    CPLFree ( wkt );

    for ( ir = 0; ir < out.rows; ir += bs ) {
        rows = out.rows - ir;
        if ( rows > bs ) rows = bs;
        //printf("Writer waiting for %d\n",ir );
        pe.wait_for("data",ir);
        for ( i = 0; i < 3; i++ ) {
            bd = df->GetRasterBand(i+1);
            bd->RasterIO( GF_Write, 0, ir, out.cols, rows,
                          out.img[i].data+ir*out.cols,
                          out.cols, rows, GDT_Byte, 0, 0 );
        }
    }

    delete df;

    df = (GDALDataset *)GDALOpen ( file, GA_Update );
    if( df == NULL ) {
        fprintf(stderr,"Could not open for update %s\n", file );
        exit(1);
    }
    nov = p.ivalue("overviews");
    if ( nov > 5 ) nov = 5;
    if ( nov > 0 ) {
        df->BuildOverviews("NEAREST", nov, ov, 3, bands, NULL, NULL );
    }

}

GDALDataset *SAGADataset::CreateCopy( const char *pszFilename,
				      GDALDataset *poSrcDS,
				      int bStrict, char **papszOptions,
				      GDALProgressFunc pfnProgress,
				      void *pProgressData )
{
    if( pfnProgress == NULL )
        pfnProgress = GDALDummyProgress;

    int nBands = poSrcDS->GetRasterCount();
    if (nBands == 0)
    {
        CPLError( CE_Failure, CPLE_NotSupported, 
                  "SAGA driver does not support source dataset with zero band.\n");
        return NULL;
    }
    else if (nBands > 1)
    {
        if( bStrict )
        {
            CPLError( CE_Failure, CPLE_NotSupported,
                      "Unable to create copy, SAGA Binary Grid "
                      "format only supports one raster band.\n" );
            return NULL;
        }
        else
            CPLError( CE_Warning, CPLE_NotSupported,
                      "SAGA Binary Grid format only supports one "
                      "raster band, first band will be copied.\n" );
    }

    GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand( 1 );
    
    char** papszCreateOptions = NULL;
    papszCreateOptions = CSLSetNameValue(papszCreateOptions, "FILL_NODATA", "NO");

    int bHasNoDataValue = FALSE;
    double dfNoDataValue = poSrcBand->GetNoDataValue(&bHasNoDataValue);
    if (bHasNoDataValue)
        papszCreateOptions = CSLSetNameValue(papszCreateOptions, "NODATA_VALUE",
                                             CPLSPrintf("%.16g", dfNoDataValue));
    
    GDALDataset* poDstDS =
        Create(pszFilename, poSrcBand->GetXSize(), poSrcBand->GetYSize(),
               1, poSrcBand->GetRasterDataType(), papszCreateOptions);
    CSLDestroy(papszCreateOptions);
    
    if (poDstDS == NULL)
        return NULL;

    /* -------------------------------------------------------------------- */
    /*      Copy band data.	                                                */
    /* -------------------------------------------------------------------- */

    CPLErr	eErr;
    
    eErr = GDALDatasetCopyWholeRaster( (GDALDatasetH) poSrcDS, 
                                       (GDALDatasetH) poDstDS,
                                       NULL,
                                       pfnProgress, pProgressData );
                                       
    if (eErr == CE_Failure)
    {
        delete poDstDS;
        return NULL;
    }

    double  adfGeoTransform[6];

    poSrcDS->GetGeoTransform( adfGeoTransform );
    poDstDS->SetGeoTransform( adfGeoTransform );
    
    poDstDS->SetProjection( poSrcDS->GetProjectionRef() );

    return poDstDS;
}