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

void write_map(fs::path file_path, GDALDataType data_type, boost::shared_ptr<Map_Matrix<DataFormat> > data, std::string WKTprojection, GeoTransform transform, std::string driverName) throw(std::runtime_error)
{
	GDALAllRegister(); //This registers all availble raster file formats for use with this utility. How neat is that. We can input any GDAL supported rater file format.
    
    const char *pszFormat = driverName.c_str();
    GDALDriver * poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if (poDriver == NULL)
    {
        throw std::runtime_error("No driver for file tyle found");
    }
    
    char ** papszMetadata = poDriver->GetMetadata();
    if (!(CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE)))
    {
        throw std::runtime_error("Driver does not support raster creation");
    }
    
    char **papszOptions = NULL;
	papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "LZW");

    GDALDataset *poDstDS = poDriver->Create(file_path.string().c_str(), (int)data->NCols(), (int)data->NRows(), 1, data_type, papszOptions);
    
    double adfGeoTransform[6] = {1, 1, 1, 1, 1, 1};
    adfGeoTransform[0] = transform.x_origin;
    adfGeoTransform[1] = transform.pixel_width;
    adfGeoTransform[2] = transform.x_line_space;
    adfGeoTransform[3] = transform.y_origin;
    adfGeoTransform[4] = transform.pixel_height;
    adfGeoTransform[5] = transform.y_line_space;
    
    const char * psz_WKT = WKTprojection.c_str();
    poDstDS->SetGeoTransform(adfGeoTransform);             
    poDstDS->SetProjection(psz_WKT);
    
    DataFormat * pafScanline = new DataFormat[data->NCols() * data->NRows()];
    int pafIterator = 0;
	for (int i = 0; i < data->NRows(); i++)
    {
		for (int j = 0; j < data->NCols(); j++)
        {
            pafScanline[pafIterator] = data->Get(i, j);
            pafIterator++;
        }
    }
    
    GDALRasterBand * poBand = poDstDS->GetRasterBand(1);
    poBand->SetNoDataValue(data->NoDataValue());
    poBand->RasterIO(GF_Write, 0, 0, (int) data->NCols(), (int) data->NRows(), pafScanline, (int) data->NCols(), (int) data->NRows(), data_type, 0, 0);
    
    GDALClose( (GDALDatasetH) poDstDS);
}

SEXP
RGDAL_SetNoDataValue(SEXP sxpRasterBand, SEXP NoDataValue) {
  CPLErr err;

  GDALRasterBand *pRasterBand = getGDALRasterPtr(sxpRasterBand);

  err = pRasterBand->SetNoDataValue(NUMERIC_POINTER(NoDataValue)[0]);

  if (err == CE_Failure)
	warning("setting of missing value not supported by this driver");

  return(sxpRasterBand);

}

  void saveGDAL(const std::string &filename, const std::string &template_name, int xoffset, int yoffset){
    GDALDataset *fintempl = (GDALDataset*)GDALOpen(template_name.c_str(), GA_ReadOnly);
    assert(fintempl!=NULL); //TODO: Error handle

    GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
    assert(poDriver!=NULL); //TODO: Error handle
    GDALDataset *fout    = poDriver->Create(filename.c_str(), viewWidth(), viewHeight(), 1, myGDALType(), NULL);
    assert(fout!=NULL);     //TODO: Error handle

    GDALRasterBand *oband = fout->GetRasterBand(1);
    oband->SetNoDataValue(no_data);

    //The geotransform maps each grid cell to a point in an affine-transformed
    //projection of the actual terrain. The geostransform is specified as follows:
    //    Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
    //    Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
    //In case of north up images, the GT(2) and GT(4) coefficients are zero, and
    //the GT(1) is pixel width, and GT(5) is pixel height. The (GT(0),GT(3))
    //position is the top left corner of the top left pixel of the raster.
    double geotrans[6];
    fintempl->GetGeoTransform(geotrans);

    //We shift the top-left pixel of hte image eastward to the appropriate
    //coordinate
    geotrans[0] += xoffset*geotrans[1];

    //We shift the top-left pixel of the image southward to the appropriate
    //coordinate
    geotrans[3] += yoffset*geotrans[5];

    #ifdef DEBUG
      std::cerr<<"Filename: "<<std::setw(20)<<filename<<" Xoffset: "<<std::setw(6)<<xoffset<<" Yoffset: "<<std::setw(6)<<yoffset<<" Geotrans0: "<<std::setw(10)<<std::setprecision(10)<<std::fixed<<geotrans[0]<<" Geotrans3: "<<std::setw(10)<<std::setprecision(10)<<std::fixed<<geotrans[3]<< std::endl;
    #endif

    fout->SetGeoTransform(geotrans);

    const char* projection_string=fintempl->GetProjectionRef();
    fout->SetProjection(projection_string);

    GDALClose(fintempl);

    for(int y=0;y<view_height;y++)
      oband->RasterIO(GF_Write, 0, y, viewWidth(), 1, data[y].data(), viewWidth(), 1, myGDALType(), 0, 0);

    GDALClose(fout);
  }

//.........这里部分代码省略.........
	if ( theCols <= 0 )
		return false;

	if ( !theData )
		return false;

/*	bool formatSupported = false;
	QMapIterator<QString, QString> i(mSupportedFormats);
	while (i.hasNext()) 
	{
		i.next();
		if( theFormat == i.key())
		{
			formatSupported = true;
			break;
		}
	}

	if ( !formatSupported )
		return false;

*/

	//GDALAllRegister();

	GDALDriver * driver;

	//set format
	char * format = new char[theFormat.size() + 1];
	strcpy( format, theFormat.toLocal8Bit().data() );

    driver = GetGDALDriverManager()->GetDriverByName(format);

	if( driver == NULL )
        return false;

	char ** metadata = driver->GetMetadata();
    if( !CSLFetchBoolean( metadata, GDAL_DCAP_CREATE, FALSE ) )
        return false;
		
	GDALDataset * dstDS;     

	//set options
	char ** options = NULL;
	options = CSLSetNameValue( options, "COMPRESS", "LZW" );

	//if it is a GeoTIFF format set correct compression options
	if ( !strcmp( format, "GTiff" ) )
	{
		if( theType == GDT_Byte )
		{
			options = CSLSetNameValue( options, "PREDICTOR", "1" );
		}
		else
		{
			if ( theType == GDT_UInt16 || theType == GDT_Int16  
				|| theType == GDT_UInt32 || theType == GDT_Int32 )
			{
				options = CSLSetNameValue( options, "PREDICTOR", "2" );
			} 
			else 
			{
				options = CSLSetNameValue( options, "PREDICTOR", "3" );
			}
		}
	}
	
	//set filename
	char * dstFilename = new char[theFilename.size() + 1];
	strcpy( dstFilename, theFilename.toLocal8Bit().data() );

	dstDS = driver->Create( dstFilename, theCols, theRows, theBands, theType, 
								options );
	delete dstFilename;
	delete [] options;

	//set geotransform
	dstDS->SetGeoTransform( theGeoTransform );

	//set CRS
	char * crsWkt = new char[theCrs->toWkt().size() + 1];
	strcpy( crsWkt, theCrs->toWkt().toLocal8Bit().data());
    dstDS->SetProjection( crsWkt );
    delete crsWkt;

	GDALRasterBand * band;

	for( int i=1; i <= theBands; i++ )
	{
		band = dstDS->GetRasterBand( i );
		band->SetNoDataValue( theNodataValue );
		band->RasterIO( GF_Write, 0, 0, theCols, theRows, theData[ i-1 ], 
						theCols, theRows, theType, 0, 0);
	}

	GDALClose( (GDALDatasetH) dstDS );

	return true;

}

//.........这里部分代码省略.........
        {
            poODS = static_cast<GDALDataset *>(
                GDALOpen( osOvrFilename, GA_Update ) );
            if( poODS == NULL )
                eErr = CE_Failure;
        }
    }

    GDALDestroyScaledProgress( pScaledProgress );

/* -------------------------------------------------------------------- */
/*      Refresh old overviews that were listed.                         */
/* -------------------------------------------------------------------- */
    GDALRasterBand **papoOverviewBands = static_cast<GDALRasterBand **>(
        CPLCalloc(sizeof(void*), nOverviews) );

    for( int iBand = 0; iBand < nBands && eErr == CE_None; iBand++ )
    {
        poBand = poDS->GetRasterBand( panBandList[iBand] );

        nNewOverviews = 0;
        for( int i = 0; i < nOverviews && poBand != NULL; i++ )
        {
            for( int j = 0; j < poBand->GetOverviewCount(); j++ )
            {
                GDALRasterBand * poOverview = poBand->GetOverview( j );
                if( poOverview == NULL )
                    continue;

                int bHasNoData = FALSE;
                double noDataValue = poBand->GetNoDataValue(&bHasNoData);

                if( bHasNoData )
                  poOverview->SetNoDataValue(noDataValue);

                const int nOvFactor =
                    GDALComputeOvFactor(poOverview->GetXSize(),
                                        poBand->GetXSize(),
                                        poOverview->GetYSize(),
                                        poBand->GetYSize());

                if( nOvFactor == - panOverviewList[i]
                    || (panOverviewList[i] < 0 &&
                        nOvFactor == GDALOvLevelAdjust2( -panOverviewList[i],
                                                       poBand->GetXSize(),
                                                       poBand->GetYSize() )) )
                {
                    papoOverviewBands[nNewOverviews++] = poOverview;
                    break;
                }
            }
        }

        if( nNewOverviews > 0 )
        {
            const double dfOffset = dfAreaNewOverviews / dfAreaRefreshedOverviews;
            const double dfScale = 1.0 - dfOffset;
            pScaledProgress = GDALCreateScaledProgress(
                    dfOffset + dfScale * iBand / nBands,
                    dfOffset + dfScale * (iBand+1) / nBands,
                    pfnProgress, pProgressData );
            eErr = GDALRegenerateOverviews( (GDALRasterBandH) poBand,
                                            nNewOverviews,
                                            (GDALRasterBandH*)papoOverviewBands,
                                            pszResampling,
                                            GDALScaledProgress, pScaledProgress );

//.........这里部分代码省略.........
                                gdalFiles[i]->SetGeoTransform(defaultTransform);
                            }
                            if (wkt)
                                gdalFiles[i]->SetProjection(wkt);
                        }
                    }
                }
            } else {
                gdalFiles[i] = NULL;
            }
        }
    } else {
      gdalFiles = NULL;
    }

    if (gdalFiles != NULL)
    {
        for (i = 0; i < numTypes; i++)
        {
            if (gdalFiles[i] != NULL)
            {
                float *poRasterData = new float[GRID_SIZE_X*GRID_SIZE_Y];
                for (int j = 0; j < GRID_SIZE_X*GRID_SIZE_Y; j++)
                {
                    poRasterData[j] = 0;
                }

                for(j = GRID_SIZE_Y - 1; j >= 0; j--)
                {
                    for(k = 0; k < GRID_SIZE_X; k++)
                    {
                        int index = (GRID_SIZE_Y - 1 - j) * GRID_SIZE_X + k;

                        if(interp[k][j].empty == 0 &&
                                interp[k][j].filled == 0)
                        {
                            poRasterData[index] = -9999.f;
                        } else {
                            switch (i)
                            {
                                case 0:
                                    poRasterData[index] = interp[k][j].Zmin;
                                    break;

                                case 1:
                                    poRasterData[index] = interp[k][j].Zmax;
                                    break;

                                case 2:
                                    poRasterData[index] = interp[k][j].Zmean;
                                    break;

                                case 3:
                                    poRasterData[index] = interp[k][j].Zidw;
                                    break;

                                case 4:
                                    poRasterData[index] = interp[k][j].count;
                                    break;

                                case 5:
                                    poRasterData[index] = interp[k][j].Zstd;
                                    break;
                            }
                        }
                    }
                }
                GDALRasterBand *tBand = gdalFiles[i]->GetRasterBand(1);
                tBand->SetNoDataValue(-9999.f);

                if (GRID_SIZE_X > 0 && GRID_SIZE_Y > 0)
                    tBand->RasterIO(GF_Write, 0, 0, GRID_SIZE_X, GRID_SIZE_Y, poRasterData, GRID_SIZE_X, GRID_SIZE_Y, GDT_Float32, 0, 0);
                GDALClose((GDALDatasetH) gdalFiles[i]);
                delete [] poRasterData;
            }
        }
    }
#endif // HAVE_GDAL

    // close files
    if(gridFiles != NULL)
    {
        for(i = 0; i < numTypes; i++)
        {
            if(gridFiles[i] != NULL)
                fclose(gridFiles[i]);
        }
    }

    if(arcFiles != NULL)
    {
        for(i = 0; i < numTypes; i++)
        {
            if(arcFiles[i] != NULL)
                fclose(arcFiles[i]);
        }
    }

    return 0;
}

// Slot called when the menu item is triggered
// If you created more menu items / toolbar buttons in initiGui, you should
// create a separate handler for each action - this single run() method will
// not be enough
void Heatmap::run()
{
  HeatmapGui d( mQGisIface->mainWindow(), QgisGui::ModalDialogFlags, &mSessionSettings );

  if ( d.exec() == QDialog::Accepted )
  {
    // everything runs here

    // Get the required data from the dialog
    QgsRectangle myBBox = d.bbox();
    int columns = d.columns();
    int rows = d.rows();
    double cellsize = d.cellSizeX(); // or d.cellSizeY();  both have the same value
    mDecay = d.decayRatio();
    int kernelShape = d.kernelShape();

    // Start working on the input vector
    QgsVectorLayer* inputLayer = d.inputVectorLayer();

    // Getting the rasterdataset in place
    GDALAllRegister();

    GDALDataset *emptyDataset;
    GDALDriver *myDriver;

    myDriver = GetGDALDriverManager()->GetDriverByName( d.outputFormat().toUtf8() );
    if ( myDriver == NULL )
    {
      QMessageBox::information( 0, tr( "GDAL driver error" ), tr( "Cannot open the driver for the specified format" ) );
      return;
    }

    double geoTransform[6] = { myBBox.xMinimum(), cellsize, 0, myBBox.yMinimum(), 0, cellsize };
    emptyDataset = myDriver->Create( d.outputFilename().toUtf8(), columns, rows, 1, GDT_Float32, NULL );
    emptyDataset->SetGeoTransform( geoTransform );
    // Set the projection on the raster destination to match the input layer
    emptyDataset->SetProjection( inputLayer->crs().toWkt().toLocal8Bit().data() );

    GDALRasterBand *poBand;
    poBand = emptyDataset->GetRasterBand( 1 );
    poBand->SetNoDataValue( NO_DATA );

    float* line = ( float * ) CPLMalloc( sizeof( float ) * columns );
    for ( int i = 0; i < columns ; i++ )
    {
      line[i] = NO_DATA;
    }
    // Write the empty raster
    for ( int i = 0; i < rows ; i++ )
    {
      poBand->RasterIO( GF_Write, 0, i, columns, 1, line, columns, 1, GDT_Float32, 0, 0 );
    }

    CPLFree( line );
    //close the dataset
    GDALClose(( GDALDatasetH ) emptyDataset );

    // open the raster in GA_Update mode
    GDALDataset *heatmapDS;
    heatmapDS = ( GDALDataset * ) GDALOpen( d.outputFilename().toUtf8(), GA_Update );
    if ( !heatmapDS )
    {
      QMessageBox::information( 0, tr( "Raster update error" ), tr( "Could not open the created raster for updating. The heatmap was not generated." ) );
      return;
    }
    poBand = heatmapDS->GetRasterBand( 1 );

    QgsAttributeList myAttrList;
    int rField = 0;
    int wField = 0;

    // Handle different radius options
    double radius;
    double radiusToMapUnits = 1;
    int myBuffer = 0;
    if ( d.variableRadius() )
    {
      rField = d.radiusField();
      myAttrList.append( rField );
      QgsDebugMsg( QString( "Radius Field index received: %1" ).arg( rField ) );

      // If not using map units, then calculate a conversion factor to convert the radii to map units
      if ( d.radiusUnit() == HeatmapGui::Meters )
      {
        radiusToMapUnits = mapUnitsOf( 1, inputLayer->crs() );
      }
    }
    else
    {
      radius = d.radius(); // radius returned by d.radius() is already in map units
      myBuffer = bufferSize( radius, cellsize );
    }

    if ( d.weighted() )
    {
      wField = d.weightField();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    CPLErr eErr = CE_None;

    if(  ((bSourceIsPixelInterleaved && bSourceIsJPEG2000) ||
          (nCompression != COMPRESSION_NONE)) &&
         nPlanarConfig == PLANARCONFIG_CONTIG &&
         !GDALDataTypeIsComplex(papoBandList[0]->GetRasterDataType()) &&
         papoBandList[0]->GetColorTable() == nullptr &&
         (STARTS_WITH_CI(pszResampling, "NEAR") ||
          EQUAL(pszResampling, "AVERAGE") ||
          EQUAL(pszResampling, "GAUSS") ||
          EQUAL(pszResampling, "CUBIC") ||
          EQUAL(pszResampling, "CUBICSPLINE") ||
          EQUAL(pszResampling, "LANCZOS") ||
          EQUAL(pszResampling, "BILINEAR")) )
    {
        // In the case of pixel interleaved compressed overviews, we want to
        // generate the overviews for all the bands block by block, and not
        // band after band, in order to write the block once and not loose
        // space in the TIFF file.
        GDALRasterBand ***papapoOverviewBands =
            static_cast<GDALRasterBand ***>(
                CPLCalloc(sizeof(void *), nBands) );
        for( int iBand = 0; iBand < nBands && eErr == CE_None; iBand++ )
        {
            GDALRasterBand *poSrcBand = papoBandList[iBand];
            GDALRasterBand *poDstBand = hODS->GetRasterBand( iBand + 1 );
            papapoOverviewBands[iBand] =
                static_cast<GDALRasterBand **>(
                    CPLCalloc(sizeof(void *), nOverviews) );

            int bHasNoData = FALSE;
            const double noDataValue = poSrcBand->GetNoDataValue(&bHasNoData);
            if( bHasNoData )
                poDstBand->SetNoDataValue(noDataValue);

            for( int i = 0; i < nOverviews && eErr == CE_None; i++ )
            {
                for( int j = -1; j < poDstBand->GetOverviewCount() &&
                                 eErr == CE_None; j++ )
                {
                    GDALRasterBand * poOverview =
                            (j < 0 ) ? poDstBand : poDstBand->GetOverview( j );
                    if( poOverview == nullptr )
                    {
                        eErr = CE_Failure;
                        continue;
                    }

                    const int nOvFactor =
                        GDALComputeOvFactor(poOverview->GetXSize(),
                                            poSrcBand->GetXSize(),
                                            poOverview->GetYSize(),
                                            poSrcBand->GetYSize());

                    if( nOvFactor == panOverviewListSorted[i]
                        || nOvFactor == GDALOvLevelAdjust2(
                                            panOverviewListSorted[i],
                                            poSrcBand->GetXSize(),
                                            poSrcBand->GetYSize() ) )
                    {
                        papapoOverviewBands[iBand][i] = poOverview;
                        if( bHasNoData )
                            poOverview->SetNoDataValue(noDataValue);
                        break;
                    }
                }

//.........这里部分代码省略.........
        {
            poODS = (GDALDataset *) GDALOpen( osOvrFilename, GA_Update );
            if( poODS == NULL )
                eErr = CE_Failure;
        }
    }

/* -------------------------------------------------------------------- */
/*      Refresh old overviews that were listed.                         */
/* -------------------------------------------------------------------- */
    GDALRasterBand **papoOverviewBands;

    papoOverviewBands = (GDALRasterBand **) 
        CPLCalloc(sizeof(void*),nOverviews);

    for( int iBand = 0; iBand < nBands && eErr == CE_None; iBand++ )
    {
        poBand = poDS->GetRasterBand( panBandList[iBand] );

        nNewOverviews = 0;
        for( i = 0; i < nOverviews && poBand != NULL; i++ )
        {
            int   j;
            
            for( j = 0; j < poBand->GetOverviewCount(); j++ )
            {
                int    nOvFactor;
                GDALRasterBand * poOverview = poBand->GetOverview( j );

                int bHasNoData;
                double noDataValue = poBand->GetNoDataValue(&bHasNoData);

                if (bHasNoData)
                  poOverview->SetNoDataValue(noDataValue);

                nOvFactor = (int) 
                    (0.5 + poBand->GetXSize() / (double) poOverview->GetXSize());

                if( nOvFactor == - panOverviewList[i] 
                    || nOvFactor == GDALOvLevelAdjust( -panOverviewList[i], 
                                                       poBand->GetXSize() ) )
                {
                    papoOverviewBands[nNewOverviews++] = poOverview;
                    break;
                }
            }
        }

        if( nNewOverviews > 0 )
        {
            eErr = GDALRegenerateOverviews( (GDALRasterBandH) poBand, 
                                            nNewOverviews, 
                                            (GDALRasterBandH*)papoOverviewBands,
                                            pszResampling, 
                                            pfnProgress, pProgressData );
        }
    }

/* -------------------------------------------------------------------- */
/*      Cleanup                                                         */
/* -------------------------------------------------------------------- */
    CPLFree( papoOverviewBands );
    CPLFree( panNewOverviewList );
    CPLFree( pahBands );

/* -------------------------------------------------------------------- */

int main(int argc, char* argv[])
{
  if (argc < 2)
  {
    cout << "void-filing-color <infile> <outfile>" << endl;;
    exit(1);
  }

  const char* InFilename = argv[1];
  const char* OutFilename = argv[2];

  GDALAllRegister();

  // Open dataset and get raster band
  GDALDataset* poDataset = (GDALDataset*) GDALOpen(InFilename, GA_ReadOnly);
  if(poDataset == NULL)
  {
    cout << "Couldn't open dataset " << InFilename << endl;
  }

  GDALRasterBand *poInBandr;
  GDALRasterBand *poInBandg;
  GDALRasterBand *poInBandb;
  poInBandr = poDataset->GetRasterBand(1);
  poInBandg = poDataset->GetRasterBand(2);
  poInBandb = poDataset->GetRasterBand(3);
  double adfGeoTransform[6];
  poDataset->GetGeoTransform(adfGeoTransform);

  // Get variables from input dataset
  const int nXSize = poInBandr->GetXSize();
  const int nYSize = poInBandr->GetYSize();

  // Create the output dataset and copy over relevant metadata
  const char*  Format = "GTiff";
  GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(Format);
  char** Options = NULL;

  GDALDataset* poDS = poDriver->Create(OutFilename,nXSize,nYSize,3,GDT_Byte,Options);
  poDS->SetGeoTransform(adfGeoTransform);
  poDS->SetProjection(poDataset->GetProjectionRef());
  GDALRasterBand* poBandr = poDS->GetRasterBand(1);
  GDALRasterBand* poBandg = poDS->GetRasterBand(2);
  GDALRasterBand* poBandb = poDS->GetRasterBand(3);
  poBandr->SetNoDataValue(0);
  poBandg->SetNoDataValue(0);
  poBandb->SetNoDataValue(0);

  GDALAllRegister();

  cout << "Read image." << endl;
  CImg<unsigned char>* ReadPixels = new CImg<unsigned char>(nXSize, nYSize, 1, 3, 0);
  for (int i = 0; i < nYSize; i++) {
    cout << "\r" << i << "/" << nYSize;
    for (int j = 0; j < nXSize; j++) {
      unsigned char InPixel;
      poInBandr->RasterIO(GF_Read, j, i, 1, 1, &InPixel, 1, 1, GDT_Byte, 0, 0);
      (*ReadPixels)(j, i, 0, 0) = InPixel;
      poInBandg->RasterIO(GF_Read, j, i, 1, 1, &InPixel, 1, 1, GDT_Byte, 0, 0);
      (*ReadPixels)(j, i, 0, 1) = InPixel;
      poInBandb->RasterIO(GF_Read, j, i, 1, 1, &InPixel, 1, 1, GDT_Byte, 0, 0);
      (*ReadPixels)(j, i, 0, 2) = InPixel;
    }
  }
  cout << endl;

  cout << "Void filling." << endl;
  CImg<unsigned char>* InPixels = new CImg<unsigned char>(nXSize, nYSize, 1, 3, 0);
  for (int i = 0; i < nYSize; i++) {
    cout << "\r" << i << "/" << nYSize;
    
    for (int j = 0; j < nXSize; j++) {
    
      if (isEmpty(ReadPixels, j, i) && isEmpty(InPixels, j, i)) {
        int is = i;
        int ie = i;
        int js = j;
        int je = j;
        const int maxsize = 2;
        js = max(0, j-maxsize);
        je = min(nXSize-1, j+maxsize);
        is = max(0, i-maxsize);
        ie = min(nYSize-1, i+maxsize);
//        cout << endl << js << ", " << je << ", " << is << ", " << ie;

        float fact = 0;
        float sumr = 0;
        float sumg = 0;
        float sumb = 0;
        for (int ia = is ; ia <= ie ; ia++) {
          for (int ja = js ; ja <= je ; ja++) {
//            cout << endl << ia << ", " << ja;
            if (!isEmpty(ReadPixels, ja, ia)) {
              int ik = ia - i;
              int jk = ja - j;
              float length = ik*ik+jk*jk;
              float coef = 1/(length*length);
              sumr += ((*ReadPixels)(ja, ia, 0, 0))*coef;
              sumg += ((*ReadPixels)(ja, ia, 0, 1))*coef;
              sumb += ((*ReadPixels)(ja, ia, 0, 2))*coef;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        delete poDS;
        return NULL;
    }

    poDS->eAccess = poOpenInfo->eAccess;

/* -------------------------------------------------------------------- */
/*      Compute the line offsets.                                        */
/* -------------------------------------------------------------------- */

    const long int nItemSize = GDALGetDataTypeSize(eDataType)/8;
    const long int nPixelOffset = nItemSize;
    const long int nLineOffset = nPixelOffset * nCols + atoi(poDS->GetKeyword("NBB")) ;
    const long int nBandOffset = nLineOffset * nRows;

    nSkipBytes = atoi(poDS->GetKeyword("LBLSIZE"));

/* -------------------------------------------------------------------- */
/*      Create band information objects.                                */
/* -------------------------------------------------------------------- */
    for( int i = 0; i < nBands; i++ )
    {
        GDALRasterBand	*poBand
            = new RawRasterBand( poDS, i+1, poDS->fpImage, nSkipBytes + nBandOffset * i,
                                 nPixelOffset, nLineOffset, eDataType,
#ifdef CPL_LSB
                                   chByteOrder == 'I' || chByteOrder == 'L',
#else
                                   chByteOrder == 'M',
#endif
                                   TRUE );

        poDS->SetBand( i+1, poBand );
        poBand->SetNoDataValue( dfNoData );
        if (bIsDTM) {
            poBand->SetScale( (double) CPLAtof(poDS->GetKeyword( "DTM.DTM_SCALING_FACTOR") ) );
            poBand->SetOffset( (double) CPLAtof(poDS->GetKeyword( "DTM.DTM_OFFSET") ) );
            const char* pszMin = poDS->GetKeyword( "DTM.DTM_MINIMUM_DN", NULL );
            const char* pszMax = poDS->GetKeyword( "DTM.DTM_MAXIMUM_DN", NULL );
            if (pszMin != NULL && pszMax != NULL )
                poBand->SetStatistics(CPLAtofM(pszMin),CPLAtofM(pszMax),0,0);
            const char* pszNoData = poDS->GetKeyword( "DTM.DTM_MISSING_DN", NULL );
            if (pszNoData != NULL )
                poBand->SetNoDataValue( CPLAtofM(pszNoData) );
        } else if (EQUAL( poDS->GetKeyword( "BLTYPE"), "M94_HRSC" )) {
            float scale=CPLAtof(poDS->GetKeyword("DLRTO8.REFLECTANCE_SCALING_FACTOR","-1."));
            if (scale < 0.) {
                scale = CPLAtof(poDS->GetKeyword( "HRCAL.REFLECTANCE_SCALING_FACTOR","1."));
            }
            poBand->SetScale( scale );
            float offset=CPLAtof(poDS->GetKeyword("DLRTO8.REFLECTANCE_OFFSET","-1."));
            if (offset < 0.) {
                offset = CPLAtof(poDS->GetKeyword( "HRCAL.REFLECTANCE_OFFSET","0."));
            }
            poBand->SetOffset( offset );
        }
        const char* pszMin = poDS->GetKeyword( "STATISTICS.MINIMUM", NULL );
        const char* pszMax = poDS->GetKeyword( "STATISTICS.MAXIMUM", NULL );
        const char* pszMean = poDS->GetKeyword( "STATISTICS.MEAN", NULL );
        const char* pszStdDev = poDS->GetKeyword( "STATISTICS.STANDARD_DEVIATION", NULL );
        if (pszMin != NULL && pszMax != NULL && pszMean != NULL && pszStdDev != NULL )
                poBand->SetStatistics(CPLAtofM(pszMin),CPLAtofM(pszMax),CPLAtofM(pszMean),CPLAtofM(pszStdDev));
    }


/* -------------------------------------------------------------------- */

/*
 * Algorithm taken from: http://www.quantdec.com/SYSEN597/GTKAV/section9/map_algebra.htm
 *
*/
int Raster::ReSample(const char * pOutputRaster, double fNewCellSize,
                     double fNewLeft, double fNewTop, int nNewRows, int nNewCols)
{
    if (fNewCellSize <= 0)
        return CELL_SIZE_ERROR;

    if (fNewLeft <=0)
        return LEFT_ERROR;

    if (fNewTop <=0)
        return TOP_ERROR;

    if (nNewRows <=0)
        return ROWS_ERROR;

    if (nNewCols <=0)
        return COLS_ERROR;


    /*************************************************************************************************
    * Open the original dataset and retrieve its basic properties
    */
    GDALDataset * pDSOld = (GDALDataset*) GDALOpen(m_sFilePath, GA_ReadOnly);
    if (pDSOld  == NULL)
        return INPUT_FILE_ERROR;

    GDALRasterBand * pRBInput = pDSOld->GetRasterBand(1);

    /*************************************************************************************************
     * Create the new dataset. Determine the driver from the output file extension.
     * Enforce LZW compression for TIFs. The predictor 3 is used for floating point prediction.
     * Not using this value defaults the LZW to prediction to 1 which causes striping.
     */
    char **papszOptions = NULL;
    GDALDriver * pDR = NULL;
    const char * pSuffix = ExtractFileExt(pOutputRaster);
    if (pSuffix == NULL)
        return OUTPUT_FILE_EXT_ERROR;
    else
    {
        if (strcmp(pSuffix, ".tif") == 0)
        {
            pDR = GetGDALDriverManager()->GetDriverByName("GTiff");
            papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "LZW");
            //papszOptions = CSLSetNameValue(papszOptions, "PREDICTOR", "3");
        }
        else if (strcmp(pSuffix, ".img") == 0)
            pDR = GetGDALDriverManager()->GetDriverByName("HFA");
        else
            return OUTPUT_UNHANDLED_DRIVER;
    }

    GDALDataset * pDSOutput = pDR->Create(pOutputRaster, nNewCols, nNewRows, 1,  *GetGDALDataType(), papszOptions);
    CSLDestroy( papszOptions );
    if (pDSOutput == NULL)
        return OUTPUT_FILE_ERROR;

    GDALRasterBand * pRBOutput = pDSOutput->GetRasterBand(1);

    if (HasNoDataValue())
    {
        CPLErr er = pRBOutput->SetNoDataValue(this->GetNoDataValue());
        if (er == CE_Failure || er == CE_Fatal)
            return OUTPUT_NO_DATA_ERROR;
    }

    double newTransform[6];
    newTransform[0] = fNewLeft;
    newTransform[1] = fNewCellSize;
    newTransform[2] = 0;
    newTransform[3] = fNewTop;
    newTransform[4] = 0;
    newTransform[5] = -1 * fNewCellSize;
    pDSOutput->SetGeoTransform(newTransform);
    pDSOutput->SetProjection(GetProjectionRef());

    ReSampleRaster(pRBInput, pRBOutput, fNewCellSize, fNewLeft, fNewTop, nNewRows, nNewCols);

    CalculateStats(pDSOutput->GetRasterBand(1));

    GDALClose(pDSOld);
    GDALClose(pDSOutput);

    return PROCESS_OK;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    return CE_Failure;
}

//.........这里部分代码省略.........
									0,			//0
									shiftY,		//top left y
									0,			//0
									stepY		//n-s pixel resolution (can be negative)
	};

	poDstDS->SetGeoTransform( adfGeoTransform );

	//OGRSpatialReference oSRS;
	//oSRS.SetUTM( 11, TRUE );
	//oSRS.SetWellKnownGeogCS( "NAD27" );
	//char *pszSRS_WKT = NULL;
	//oSRS.exportToWkt( &pszSRS_WKT );
	//poDstDS->SetProjection( pszSRS_WKT );
	//CPLFree( pszSRS_WKT );

	double* scanline = (double*) CPLMalloc(sizeof(double)*m_grid.width);
	int currentBand = 0;

	//exort height band?
	if (heightBand)
	{
		GDALRasterBand* poBand = poDstDS->GetRasterBand(++currentBand);
		assert(poBand);
		poBand->SetColorInterpretation(GCI_Undefined);

		EmptyCellFillOption fillEmptyCellsStrategy = getFillEmptyCellsStrategy(fillEmptyCellsComboBox);

		double emptyCellHeight = 0;
		switch (fillEmptyCellsStrategy)
		{
		case LEAVE_EMPTY:
			emptyCellHeight = m_grid.minHeight-1.0;
			poBand->SetNoDataValue(emptyCellHeight); //should be transparent!
			break;
		case FILL_MINIMUM_HEIGHT:
			emptyCellHeight = m_grid.minHeight;
			break;
		case FILL_MAXIMUM_HEIGHT:
			emptyCellHeight = m_grid.maxHeight;
			break;
		case FILL_CUSTOM_HEIGHT:
			emptyCellHeight = getCustomHeightForEmptyCells();
			break;
		case FILL_AVERAGE_HEIGHT:
			emptyCellHeight = m_grid.meanHeight;
			break;
		default:
			assert(false);
		}

		for (unsigned j=0; j<m_grid.height; ++j)
		{
			const RasterCell* aCell = m_grid.data[j];
			for (unsigned i=0; i<m_grid.width; ++i,++aCell)
			{
				scanline[i] = aCell->h == aCell->h ? aCell->h : emptyCellHeight;
			}

			if (poBand->RasterIO( GF_Write, 0, static_cast<int>(j), static_cast<int>(m_grid.width), 1, scanline, static_cast<int>(m_grid.width), 1, GDT_Float64, 0, 0 ) != CE_None)
			{
				ccLog::Error("[GDAL] An error occurred while writing the height band!");
				if (scanline)
					CPLFree(scanline);
				GDALClose( (GDALDatasetH) poDstDS );
				return;

//.........这里部分代码省略.........
                      nJpegQuality );
        GTIFFSetJpegQuality((GDALDatasetH)hODS, nJpegQuality);
    }

/* -------------------------------------------------------------------- */
/*      Loop writing overview data.                                     */
/* -------------------------------------------------------------------- */

    if (nCompression != COMPRESSION_NONE &&
        nPlanarConfig == PLANARCONFIG_CONTIG &&
        GDALDataTypeIsComplex(papoBandList[0]->GetRasterDataType()) == FALSE &&
        papoBandList[0]->GetColorTable() == NULL &&
        (EQUALN(pszResampling, "NEAR", 4) || EQUAL(pszResampling, "AVERAGE") ||
         EQUAL(pszResampling, "GAUSS") || EQUAL(pszResampling, "CUBIC") ||
         EQUAL(pszResampling, "CUBICSPLINE") || EQUAL(pszResampling, "LANCZOS") ||
         EQUAL(pszResampling, "BILINEAR")))
    {
        /* In the case of pixel interleaved compressed overviews, we want to generate */
        /* the overviews for all the bands block by block, and not band after band, */
        /* in order to write the block once and not loose space in the TIFF file */
        GDALRasterBand ***papapoOverviewBands;

        papapoOverviewBands = (GDALRasterBand ***) CPLCalloc(sizeof(void*),nBands);
        for( iBand = 0; iBand < nBands && eErr == CE_None; iBand++ )
        {
            GDALRasterBand    *hSrcBand = papoBandList[iBand];
            GDALRasterBand    *hDstBand = hODS->GetRasterBand( iBand+1 );
            papapoOverviewBands[iBand] = (GDALRasterBand **) CPLCalloc(sizeof(void*),nOverviews);
            papapoOverviewBands[iBand][0] = hDstBand;

            int bHasNoData;
            double noDataValue = hSrcBand->GetNoDataValue(&bHasNoData);
            if (bHasNoData)
                hDstBand->SetNoDataValue(noDataValue);

            for( int i = 0; i < nOverviews-1 && eErr == CE_None; i++ )
            {
                papapoOverviewBands[iBand][i+1] = hDstBand->GetOverview(i);
                if (papapoOverviewBands[iBand][i+1] == NULL)
                    eErr = CE_Failure;
                else
                {
                    if (bHasNoData)
                        papapoOverviewBands[iBand][i+1]->SetNoDataValue(noDataValue);
                }
            }
        }

        if (eErr == CE_None)
            eErr = GDALRegenerateOverviewsMultiBand(nBands, papoBandList,
                                            nOverviews, papapoOverviewBands,
                                            pszResampling, pfnProgress, pProgressData );

        for( iBand = 0; iBand < nBands; iBand++ )
        {
            CPLFree(papapoOverviewBands[iBand]);
        }
        CPLFree(papapoOverviewBands);
    }
    else
    {
        GDALRasterBand   **papoOverviews;

        papoOverviews = (GDALRasterBand **) CPLCalloc(sizeof(void*),128);

        for( iBand = 0; iBand < nBands && eErr == CE_None; iBand++ )