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

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

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

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

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

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

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

SEXP
RGDAL_GetNoDataValue(SEXP sxpRasterBand) {

  GDALRasterBand *pRasterBand = getGDALRasterPtr(sxpRasterBand);

  int hasNoDataValue;

  double noDataValue = pRasterBand->GetNoDataValue(&hasNoDataValue);

  return(hasNoDataValue ? ScalarReal(noDataValue) : R_NilValue);

}

feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
{
    CPLErr raster_io_error = CE_None;

    if (band_ > 0)
    {
        unsigned raster_xsize = dataset_.GetRasterXSize();
        unsigned raster_ysize = dataset_.GetRasterYSize();

        double gt[6];
        dataset_.GetGeoTransform(gt);

        double det = gt[1] * gt[5] - gt[2] * gt[4];
        // subtract half a pixel width & height because gdal coord reference
        // is the top-left corner of a pixel, not the center.
        double X = pt.x - gt[0] - gt[1]/2;
        double Y = pt.y - gt[3] - gt[5]/2;
        double det1 = gt[1]*Y + gt[4]*X;
        double det2 = gt[2]*Y + gt[5]*X;
        unsigned x = static_cast<unsigned>(det2/det);
        unsigned y = static_cast<unsigned>(det1/det);

        if (x < raster_xsize && y < raster_ysize)
        {
            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: pt.x=" << pt.x << " pt.y=" << pt.y;
            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: x=" << x << " y=" << y;

            GDALRasterBand* band = dataset_.GetRasterBand(band_);
            int raster_has_nodata;
            double nodata = band->GetNoDataValue(&raster_has_nodata);
            double value;
            raster_io_error = band->RasterIO(GF_Read, x, y, 1, 1, &value, 1, 1, GDT_Float64, 0, 0);
            if (raster_io_error == CE_Failure) {
                throw datasource_exception(CPLGetLastErrorMsg());
            }
            if (! raster_has_nodata || value != nodata)
            {
                // construct feature
                feature_ptr feature = feature_factory::create(ctx_,1);
                feature->set_geometry(mapnik::geometry::point<double>(pt.x,pt.y));
                feature->put_new("value",value);
                if (raster_has_nodata)
                {
                    feature->put_new("nodata",nodata);
                }
                return feature;
            }
        }
    }
    return feature_ptr();
}

void getGDALHeader(const std::string &filename, int &height, int &width, T &no_data, double *geotrans){
  GDALAllRegister();
  GDALDataset *fin = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
  assert(fin!=NULL);

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

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

  fin->GetGeoTransform(geotrans);

  GDALClose(fin);
}

feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
{
    if (band_ > 0)
    {
        unsigned raster_xsize = dataset_.GetRasterXSize();
        unsigned raster_ysize = dataset_.GetRasterYSize();

        double gt[6];
        dataset_.GetGeoTransform(gt);

        double det = gt[1] * gt[5] - gt[2] * gt[4];
        // subtract half a pixel width & height because gdal coord reference
        // is the top-left corner of a pixel, not the center.
        double X = pt.x - gt[0] - gt[1]/2;
        double Y = pt.y - gt[3] - gt[5]/2;
        double det1 = gt[1]*Y + gt[4]*X;
        double det2 = gt[2]*Y + gt[5]*X;
        unsigned x = static_cast<unsigned>(det2/det);
        unsigned y = static_cast<unsigned>(det1/det);

        if (x < raster_xsize && y < raster_ysize)
        {
            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: pt.x=" << pt.x << " pt.y=" << pt.y;
            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: x=" << x << " y=" << y;

            GDALRasterBand* band = dataset_.GetRasterBand(band_);
            int raster_has_nodata;
            double nodata = band->GetNoDataValue(&raster_has_nodata);
            double value;
            band->RasterIO(GF_Read, x, y, 1, 1, &value, 1, 1, GDT_Float64, 0, 0);
            if (! raster_has_nodata || value != nodata)
            {
                // construct feature
                feature_ptr feature = feature_factory::create(ctx_,1);
                std::unique_ptr<geometry_type> point = std::make_unique<geometry_type>(mapnik::geometry_type::types::Point);
                point->move_to(pt.x, pt.y);
                feature->add_geometry(point.release());
                feature->put_new("value",value);
                if (raster_has_nodata)
                {
                    feature->put_new("nodata",nodata);
                }
                return feature;
            }
        }
    }
    return feature_ptr();
}

feature_ptr gdal_featureset::get_feature_at_point(mapnik::coord2d const& pt)
{
    if (band_ > 0)
    {
        unsigned raster_xsize = dataset_.GetRasterXSize();
        unsigned raster_ysize = dataset_.GetRasterYSize();

        double gt[6];
        dataset_.GetGeoTransform(gt);

        double det = gt[1] * gt[5] - gt[2] * gt[4];
        // subtract half a pixel width & height because gdal coord reference
        // is the top-left corner of a pixel, not the center.
        double X = pt.x - gt[0] - gt[1]/2;
        double Y = pt.y - gt[3] - gt[5]/2;
        double det1 = gt[1]*Y + gt[4]*X;
        double det2 = gt[2]*Y + gt[5]*X;
        unsigned x = det2/det, y = det1/det;

        if (x < raster_xsize && y < raster_ysize)
        {
#ifdef MAPNIK_DEBUG
            std::clog << boost::format("GDAL Plugin: pt.x=%f pt.y=%f") % pt.x % pt.y << std::endl;
            std::clog << boost::format("GDAL Plugin: x=%f y=%f") % x % y << std::endl;
#endif
            GDALRasterBand* band = dataset_.GetRasterBand(band_);
            int hasNoData;
            double nodata = band->GetNoDataValue(&hasNoData);
            double value;
            band->RasterIO(GF_Read, x, y, 1, 1, &value, 1, 1, GDT_Float64, 0, 0);

            if (! hasNoData || value != nodata)
            {
                // construct feature
                feature_ptr feature(new Feature(1));
                geometry_type * point = new geometry_type(mapnik::Point);
                point->move_to(pt.x, pt.y);
                feature->add_geometry(point);
                (*feature)["value"] = value;
                return feature;
            }
        }
    }
    return feature_ptr();
}

void getGDALHeader(
  const   std::string &filename,
  int32_t &height,
  int32_t &width,
  T       &no_data,
  double  geotransform[6]
){
  GDALAllRegister();
  GDALDataset *fin = (GDALDataset*)GDALOpen(filename.c_str(), GA_ReadOnly);
  if(fin==NULL)
    throw std::runtime_error("Could not get GDAL header: file '" + filename + "'' did not open!");

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

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

  fin->GetGeoTransform(geotransform);

  GDALClose(fin);
}

Dataset::ElevationDataPtr Dataset::convertToElevationData() const
{
  if ( 0x0 == _data )
    return 0x0;

  GDALRasterBand* band ( _data->GetRasterBand ( 1 ) );
  if ( 0x0 == band )
    return 0x0;

  // Get the width and height.
  const int width ( _data->GetRasterXSize() );
  const int height ( _data->GetRasterYSize() );

  // Read the values.
  const int size ( width * height );
  std::vector<IElevationData::ValueType> bytes ( size, 0 );
  if ( CE_None == band->RasterIO ( GF_Read, 0, 0, width, height, &bytes[0], width, height, GDT_Float32, 0, 0 ) )
  {
    Minerva::Core::ElevationData::RefPtr elevationData ( new Minerva::Core::ElevationData ( width, height ) );

    // The no data value.
    const double noDataValue ( band->GetNoDataValue() );
    elevationData->noDataValue ( static_cast<IElevationData::ValueType> ( noDataValue ) );

    for ( int row = 0; row < height; ++row )
    {
      for ( int column = 0; column < width; ++column )
      {
        const unsigned int index ( column + ( row * width ) );
        elevationData->value ( column, ( height - row - 1 ), bytes.at ( index ) ); 
      }
    }

    return ElevationDataPtr ( elevationData );
  }

  return 0x0;

}

int GDALFillBandNoData(GDALDataset *poDS, int nBand, int nSearchPixels)
{
    (void)nBand;
    (void)nSearchPixels;
    if(poDS == NULL)
    {
        fprintf(stderr, "Invalid GDAL Dataset Handle, cannot fill no data\n");
        return -1;
    }
    int nPixels, nLines;
    nPixels = poDS->GetRasterXSize();
    nLines = poDS->GetRasterYSize();

    GDALRasterBand *poBand;

    poBand = poDS->GetRasterBand(1);

    GDALFillNodata(poBand, NULL, 100, 0, 0, NULL, NULL, NULL);

    double dfNoData = poBand->GetNoDataValue(NULL);

    double *padfScanline;
    padfScanline = (double *) CPLMalloc(sizeof(double)*nPixels);
    int nNoDataCount = 0;
    for(int i = 0;i < nLines;i++)
    {
        GDALRasterIO(poBand, GF_Read, 0, i, nPixels, 1,
                     padfScanline, nPixels, 1, GDT_Float64, 0, 0);
        for(int j = 0; j < nPixels;j++)
        {
            if(CPLIsEqual(padfScanline[j], dfNoData))
                nNoDataCount++;
        }
    }

    CPLFree(padfScanline);

    return nNoDataCount;
}

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

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

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

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

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

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

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

    view_xoff = xOffset;
    view_yoff = yOffset;

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

    GDALClose(fin);
  }

/** Check for no data values in a band for an image
 * @param poDS a pointer to a valid GDALDataset
 * @param band an integer representation of which band in the image
 * @return true if the band contains any no data values
 * @warning May be cpu intensive on large images
 */
bool GDALHasNoData( GDALDataset *poDS, int band )
{
    bool hasNDV = false;
    //check if poDS is NULL #lm
    int ncols = poDS->GetRasterXSize();
    int nrows = poDS->GetRasterYSize();

    double nDV;

    GDALRasterBand *poBand = poDS->GetRasterBand( band );
    if( poBand == NULL )
	return false;

    int pbSuccess = 0;
    nDV = poBand->GetNoDataValue( &pbSuccess );
    if( pbSuccess == false )
	nDV = -9999.0;

    double *padfScanline;
    padfScanline = new double[ncols];
    for( int i = 0;i < nrows;i++ ) {
	poBand->RasterIO( GF_Read, 0, i, ncols, 1, padfScanline, ncols, 1,
			  GDT_Float64, 0, 0 );
	for( int j = 0;j < ncols;j++ ) {
	    if( CPLIsEqual( (float)padfScanline[j], (float)nDV ) )
            {
		hasNDV = true;
                goto done;
            }
	}
    }
done:
    delete[] padfScanline;

    return hasNDV;
}

//.........这里部分代码省略.........
        {
            im_width = width;
            im_height = height;
        }

        if (im_width > 0 && im_height > 0)
        {
            mapnik::raster_ptr raster = boost::make_shared<mapnik::raster>(intersect, im_width, im_height);
            feature->set_raster(raster);
            mapnik::image_data_32 & image = raster->data_;
            image.set(0xffffffff);

            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Image Size=(" << im_width << "," << im_height << ")";
            MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Reading band=" << band_;

            if (band_ > 0) // we are querying a single band
            {
                if (band_ > nbands_)
                {
                    throw datasource_exception((boost::format("GDAL Plugin: '%d' is an invalid band, dataset only has '%d' bands\n") % band_ % nbands_).str());
                }

                float* imageData = (float*)image.getBytes();
                GDALRasterBand * band = dataset_.GetRasterBand(band_);
                int hasNoData(0);
                double nodata(0);
                if (nodata_value_)
                {
                    hasNoData = 1;
                    nodata = *nodata_value_;
                }
                else
                {
                    nodata = band->GetNoDataValue(&hasNoData);
                }
                band->RasterIO(GF_Read, x_off, y_off, width, height,
                               imageData, image.width(), image.height(),
                               GDT_Float32, 0, 0);

                if (hasNoData)
                {
                    feature->put("NODATA",nodata);
                }
            }
            else // working with all bands
            {
                for (int i = 0; i < nbands_; ++i)
                {
                    GDALRasterBand * band = dataset_.GetRasterBand(i + 1);

#ifdef MAPNIK_LOG
                    get_overview_meta(band);
#endif

                    GDALColorInterp color_interp = band->GetColorInterpretation();
                    switch (color_interp)
                    {
                    case GCI_RedBand:
                        red = band;

                        MAPNIK_LOG_DEBUG(gdal) << "gdal_featureset: Found red band";

                        break;
                    case GCI_GreenBand:
                        green = band;

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

/* -------------------------------------------------------------------- */
/*      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, NULL );

        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.                                         */
/* -------------------------------------------------------------------- */
    GInt16      anProfData[3601];
    int         dfNodataCount=0;
    GByte       iPartialCell;

    for( int iProfile = 0; iProfile < psDTED->nXSize; iProfile++ )
    {
        for( int iY = 0; iY < psDTED->nYSize; iY++ )
        {
            anProfData[iY] = panData[iProfile + iY * psDTED->nXSize];
            if ( bSrcBandHasNoData && anProfData[iY] == srcBandNoData)
            {
                anProfData[iY] = DTED_NODATA_VALUE;
                dfNodataCount++;
            }
            else if ( anProfData[iY] == DTED_NODATA_VALUE )
                dfNodataCount++;
        }
        DTEDWriteProfile( psDTED, iProfile, anProfData );

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

//.........这里部分代码省略.........
        int nJpegQuality = atoi(CPLGetConfigOption("JPEG_QUALITY_OVERVIEW","75"));
        TIFFSetField( hTIFF, TIFFTAG_JPEGQUALITY, 
                      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);

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