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

QgsAlignRaster::RasterInfo::RasterInfo( const QString& layerpath ){  mDataset = GDALOpen( layerpath.toLocal8Bit().constData(), GA_ReadOnly );  if ( !mDataset )    return;  mXSize = GDALGetRasterXSize( mDataset );  mYSize = GDALGetRasterYSize( mDataset );  GDALGetGeoTransform( mDataset, mGeoTransform );  // TODO: may be null or empty string  mCrsWkt = QString::fromAscii( GDALGetProjectionRef( mDataset ) );  mBandCnt = GDALGetBandNumber( mDataset );}

void MDAL::GdalDataset::parseParameters(){  mNBands = static_cast<unsigned int>( GDALGetRasterCount( mHDataset ) );  if ( mNBands == 0 ) throw MDAL_Status::Err_InvalidData;  GDALGetGeoTransform( mHDataset, mGT ); // in case of error it returns Identid  mXSize = static_cast<unsigned int>( GDALGetRasterXSize( mHDataset ) ); //raster width in pixels  if ( mXSize == 0 ) throw MDAL_Status::Err_InvalidData;  mYSize = static_cast<unsigned int>( GDALGetRasterYSize( mHDataset ) ); //raster height in pixels  if ( mYSize == 0 ) throw MDAL_Status::Err_InvalidData;  mNPoints = mXSize * mYSize;  mNVolumes = ( mXSize - 1 ) * ( mYSize - 1 );}

static void WorkerFunc( void * ){    GDALDatasetH hDS;    int iIter, iOpenIter;    for( iOpenIter = 0; iOpenIter < nOpenIterations; iOpenIter++ )    {        if( bLockOnOpen )            CPLAcquireMutex( pGlobalMutex, 100.0 );        hDS = GDALOpen( pszFilename, GA_ReadOnly );        if( bLockOnOpen )            CPLReleaseMutex( pGlobalMutex );        for( iIter = 0; iIter < nIterations && hDS != NULL; iIter++ )        {            int nMyChecksum;                    nMyChecksum = GDALChecksumImage( GDALGetRasterBand( hDS, 1 ),                                              0, 0,                                              GDALGetRasterXSize( hDS ),                                              GDALGetRasterYSize( hDS ) );            if( nMyChecksum != nChecksum )            {                printf( "Checksum ERROR in worker thread!/n" );                break;            }        }        if( hDS )        {            if( bLockOnOpen )                CPLAcquireMutex( pGlobalMutex, 100.0 );            GDALClose( hDS );            if( bLockOnOpen )                CPLReleaseMutex( pGlobalMutex );        }    }    CPLAcquireMutex( pGlobalMutex, 100.0 );    nPendingThreads--;    CPLReleaseMutex( pGlobalMutex );}

//duplicated from QgsNineCellFilter. Todo: make common base classGDALDatasetH QgsRelief::openInputFile( int &nCellsX, int &nCellsY ){  GDALDatasetH inputDataset = GDALOpen( mInputFile.toUtf8().constData(), GA_ReadOnly );  if ( inputDataset )  {    nCellsX = GDALGetRasterXSize( inputDataset );    nCellsY = GDALGetRasterYSize( inputDataset );    //we need at least one band    if ( GDALGetRasterCount( inputDataset ) < 1 )    {      GDALClose( inputDataset );      return nullptr;    }  }  return inputDataset;}

GDALDatasetH QgsNineCellFilter::openInputFile( int& nCellsX, int& nCellsY ){  GDALDatasetH inputDataset = GDALOpen( mInputFile.toLocal8Bit().data(), GA_ReadOnly );  if ( inputDataset != NULL )  {    nCellsX = GDALGetRasterXSize( inputDataset );    nCellsY = GDALGetRasterYSize( inputDataset );    //we need at least one band    if ( GDALGetRasterCount( inputDataset ) < 1 )    {      GDALClose( inputDataset );      return NULL;    }  }  return inputDataset;}

//duplicated from QgsNineCellFilter. Todo: make common base classGDALDatasetH QgsRelief::openInputFile( int& nCellsX, int& nCellsY ){  GDALDatasetH inputDataset = GDALOpen( TO8F( mInputFile ), GA_ReadOnly );  if ( inputDataset != NULL )  {    nCellsX = GDALGetRasterXSize( inputDataset );    nCellsY = GDALGetRasterYSize( inputDataset );    //we need at least one band    if ( GDALGetRasterCount( inputDataset ) < 1 )    {      GDALClose( inputDataset );      return NULL;    }  }  return inputDataset;}

void TestQgsGdalUtils::testCreateSingleBandMemoryDataset(){  gdal::dataset_unique_ptr ds1 = QgsGdalUtils::createSingleBandMemoryDataset( GDT_Float32, QgsRectangle( 1, 1, 21, 11 ), 40, 20, QgsCoordinateReferenceSystem( "EPSG:4326" ) );  QVERIFY( ds1 );  QCOMPARE( GDALGetRasterCount( ds1.get() ), 1 );  QCOMPARE( GDALGetRasterXSize( ds1.get() ), 40 );  QCOMPARE( GDALGetRasterYSize( ds1.get() ), 20 );  QCOMPARE( GDALGetProjectionRef( ds1.get() ), EPSG_4326_WKT );  double geoTransform[6];  double geoTransformExpected[] = { 1, 0.5, 0, 11, 0, -0.5 };  QCOMPARE( GDALGetGeoTransform( ds1.get(), geoTransform ), CE_None );  QVERIFY( memcmp( geoTransform, geoTransformExpected, sizeof( double ) * 6 ) == 0 );  QCOMPARE( GDALGetRasterDataType( GDALGetRasterBand( ds1.get(), 1 ) ), GDT_Float32 );}

/* Makes a copy of a dataset, and opens it for writing.. */GDALDatasetH make_me_a_sandwitch(GDALDatasetH *in_dataset, char *filename){    char **papszOptions = NULL;    const char *pszFormat = "GTiff";    GDALDriverH hDriver;    GDALDatasetH out_gdalfile;    hDriver = GDALGetDriverByName( pszFormat );    papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );    papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "DEFLATE" );        /*Create copy..*/    /*return GDALCreateCopy( hDriver, filename, *in_dataset, FALSE, papszOptions, NULL, NULL );*/    return GDALCreate(hDriver, filename,        GDALGetRasterXSize( *in_dataset ),        GDALGetRasterYSize( *in_dataset ),        1,        GDT_Byte, papszOptions );}

GDALDatasetH QgsNineCellFilter::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver ){  if ( inputDataset == NULL )  {    return NULL;  }  int xSize = GDALGetRasterXSize( inputDataset );  int ySize = GDALGetRasterYSize( inputDataset );;  //open output file  char **papszOptions = NULL;  GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), xSize, ySize, 1, GDT_Float32, papszOptions );  if ( outputDataset == NULL )  {    return outputDataset;  }  //get geotransform from inputDataset  double geotransform[6];  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )  {    GDALClose( outputDataset );    return NULL;  }  GDALSetGeoTransform( outputDataset, geotransform );  //make sure mCellSizeX and mCellSizeY are always > 0  mCellSizeX = geotransform[1];  if ( mCellSizeX < 0 )  {    mCellSizeX = -mCellSizeX;  }  mCellSizeY = geotransform[5];  if ( mCellSizeY < 0 )  {    mCellSizeY = -mCellSizeY;  }  const char* projection = GDALGetProjectionRef( inputDataset );  GDALSetProjection( outputDataset, projection );  return outputDataset;}

static dErr JakoFileDataView(GDALDatasetH filedata,const char *name,PetscViewer viewer){  dErr err;  CPLErr cplerr;  double geo[6],data[8*12];  int nx=8,ny=12,snx,sny;  GDALRasterBandH band;  dFunctionBegin;  cplerr = GDALGetGeoTransform(filedata,geo);  err = dRealTableView(2,3,geo,PETSC_VIEWER_STDOUT_WORLD,"%s:geo",name);dCHK(err);  snx = GDALGetRasterXSize(filedata);  sny = GDALGetRasterYSize(filedata);  err = PetscViewerASCIIPrintf(viewer,"%s: nx=%d ny=%d/n",name,snx,sny);dCHK(err);  band = GDALGetRasterBand(filedata,1);  cplerr = GDALRasterIO(band,GF_Read,snx/2,sny/2,nx,ny,data,nx,ny,GDT_Float64,0,0);dCPLCHK(cplerr);  err = dRealTableView(ny,nx,data,PETSC_VIEWER_STDOUT_WORLD,name);dCHK(err);  dFunctionReturn(0);}

static dErr JakoGDALMemAddBand(GDALDatasetH dset,GDALDataType dtype,void *memory){  char buf[256] = {0},**bandoptions = NULL;  int bytes,nx,ny;  CPLErr cplerr;  dErr err;  dFunctionBegin;  bytes = GDALGetDataTypeSize(dtype);  nx = GDALGetRasterXSize(dset);  ny = GDALGetRasterYSize(dset);  err = dMalloc(nx*ny*bytes,(void**)memory);dCHK(err);  // This is where the API moves from merely cumbersome to outright demeaning, like some twisted hazing ritual.  CPLPrintPointer(buf,*(void**)memory,sizeof(buf));  bandoptions = CSLSetNameValue(bandoptions,"DATAPOINTER",buf);  cplerr = GDALAddBand(dset,dtype,bandoptions);dCPLCHK(cplerr);  CSLDestroy(bandoptions);  dFunctionReturn(0);}

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

// Determines the pixel/line position given an x/y.// No reprojection is done at this time.bool Colorization::getPixelAndLinePosition(double x,        double y,        boost::array<double, 6> const& inverse,        boost::int32_t& pixel,        boost::int32_t& line,        void* ds){#ifdef PDAL_HAVE_GDAL    pixel = (boost::int32_t) std::floor(                inverse[0]                + inverse[1] * x                + inverse[2] * y);    line = (boost::int32_t) std::floor(               inverse[3]               + inverse[4] * x               + inverse[5] * y);        int xs = GDALGetRasterXSize(ds);    int ys = GDALGetRasterYSize(ds);        if (!xs || !ys)    {        throw pdal_error("Unable to get X or Y size from raster!");    }    if (pixel < 0 ||         line < 0 ||         pixel >= xs ||         line  >= ys       )    {        // The x, y is not coincident with this raster        return false;    }#endif    return true;}

int msUVRASTERLayerGetExtent(layerObj *layer, rectObj *extent){  char szPath[MS_MAXPATHLEN];  mapObj *map = layer->map;  double adfGeoTransform[6];  int nXSize, nYSize;  GDALDatasetH hDS;  shapefileObj *tileshpfile;  int tilelayerindex = -1;  CPLErr eErr = CE_Failure;  char *decrypted_path;  if( (!layer->data || strlen(layer->data) == 0)      && layer->tileindex == NULL) {    /* should we be issuing a specific error about not supporting       extents for tileindexed raster layers? */    return MS_FAILURE;  }  if( map == NULL )    return MS_FAILURE;  /* If the layer use a tileindex, return the extent of the tileindex shapefile/referenced layer */  if (layer->tileindex) {    tilelayerindex = msGetLayerIndex(map, layer->tileindex);    if(tilelayerindex != -1) /* does the tileindex reference another layer */      return msLayerGetExtent(GET_LAYER(map, tilelayerindex), extent);    else {      tileshpfile = (shapefileObj *) malloc(sizeof(shapefileObj));      MS_CHECK_ALLOC(tileshpfile, sizeof(shapefileObj), MS_FAILURE);      if(msShapefileOpen(tileshpfile, "rb", msBuildPath3(szPath, map->mappath, map->shapepath, layer->tileindex), MS_TRUE) == -1)        if(msShapefileOpen(tileshpfile, "rb", msBuildPath(szPath, map->mappath, layer->tileindex), MS_TRUE) == -1)          return MS_FAILURE;      *extent = tileshpfile->bounds;      msShapefileClose(tileshpfile);      free(tileshpfile);      return MS_SUCCESS;    }  }  msTryBuildPath3(szPath, map->mappath, map->shapepath, layer->data);  decrypted_path = msDecryptStringTokens( map, szPath );  msAcquireLock( TLOCK_GDAL );  if( decrypted_path ) {    hDS = GDALOpen(decrypted_path, GA_ReadOnly );    msFree( decrypted_path );  } else    hDS = NULL;  if( hDS != NULL ) {    nXSize = GDALGetRasterXSize( hDS );    nYSize = GDALGetRasterYSize( hDS );    eErr = GDALGetGeoTransform( hDS, adfGeoTransform );    GDALClose( hDS );  }  msReleaseLock( TLOCK_GDAL );  if( hDS == NULL || eErr != CE_None ) {    return MS_FAILURE;  }  /* If this appears to be an ungeoreferenced raster than flip it for     mapservers purposes. */  if( adfGeoTransform[5] == 1.0 && adfGeoTransform[3] == 0.0 ) {    adfGeoTransform[5] = -1.0;    adfGeoTransform[3] = nYSize;  }  extent->minx = adfGeoTransform[0];  extent->maxy = adfGeoTransform[3];  extent->maxx = adfGeoTransform[0] + nXSize * adfGeoTransform[1];  extent->miny = adfGeoTransform[3] + nYSize * adfGeoTransform[5];  return MS_SUCCESS;}

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

//.........这里部分代码省略.........            if( papanLUTs )            {                int iLUT;                for( iLUT = 0; iLUT < nLUTBins; iLUT++ )                    fprintf( fpConfig, "%d ", papanLUTs[iBand][iLUT] );            }            fprintf( fpConfig, "/n" );        }        if( pszConfigFile )            fclose( fpConfig );        exit( 0 );    }    if (padfScaleMin == NULL || padfScaleMax == NULL)    {        fprintf( stderr, "-equalize or -config filename command line options must be specified./n");        exit(1);    }/* ==================================================================== *//*      Create a virtual dataset.                                       *//* ==================================================================== */    VRTDataset *poVDS;    EnhanceCBInfo *pasEInfo = (EnhanceCBInfo *)         CPLCalloc(nBandCount, sizeof(EnhanceCBInfo));        /* -------------------------------------------------------------------- *//*      Make a virtual clone.                                           *//* -------------------------------------------------------------------- */    poVDS = new VRTDataset( GDALGetRasterXSize(hDataset),                            GDALGetRasterYSize(hDataset) );    if( GDALGetGCPCount(hDataset) == 0 )    {        const char *pszProjection;        double adfGeoTransform[6];        pszProjection = GDALGetProjectionRef( hDataset );        if( pszProjection != NULL && strlen(pszProjection) > 0 )            poVDS->SetProjection( pszProjection );        if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )            poVDS->SetGeoTransform( adfGeoTransform );    }    else    {        poVDS->SetGCPs( GDALGetGCPCount(hDataset),                         GDALGetGCPs(hDataset),                        GDALGetGCPProjection( hDataset ) );    }        poVDS->SetMetadata( ((GDALDataset*)hDataset)->GetMetadata() );    for( iBand = 0; iBand < nBandCount; iBand++ )    {        VRTSourcedRasterBand   *poVRTBand;        GDALRasterBand  *poSrcBand;        GDALDataType    eBandType;        poSrcBand = ((GDALDataset *) hDataset)->GetRasterBand(iBand+1);/* -------------------------------------------------------------------- */

int main( int argc, char ** argv ){    GDALDatasetH	hDataset;    GDALRasterBandH	hBand;    int			i, iBand;    double		adfGeoTransform[6];    GDALDriverH		hDriver;    char		**papszMetadata;    int                 bComputeMinMax = FALSE;    if( !GDALBridgeInitialize( "..", stderr ) )    {        fprintf( stderr, "Unable to intiailize GDAL bridge./n" );        exit( 10 );    }    if( argc > 1 && strcmp(argv[1],"-mm") == 0 )    {        bComputeMinMax = TRUE;        argv++;    }    GDALAllRegister();    hDataset = GDALOpen( argv[1], GA_ReadOnly );        if( hDataset == NULL )    {        fprintf( stderr,                 "GDALOpen failed - %d/n%s/n",                 CPLGetLastErrorNo(), CPLGetLastErrorMsg() );        exit( 1 );    }    /* -------------------------------------------------------------------- *//*      Report general info.                                            *//* -------------------------------------------------------------------- */    hDriver = GDALGetDatasetDriver( hDataset );    printf( "Driver: %s/%s/n",            GDALGetDriverShortName( hDriver ),            GDALGetDriverLongName( hDriver ) );    printf( "Size is %d, %d/n",            GDALGetRasterXSize( hDataset ),             GDALGetRasterYSize( hDataset ) );/* -------------------------------------------------------------------- *//*      Report projection.                                              *//* -------------------------------------------------------------------- */    if( GDALGetProjectionRef( hDataset ) != NULL )    {        OGRSpatialReferenceH  hSRS;        char		      *pszProjection;        pszProjection = (char *) GDALGetProjectionRef( hDataset );        hSRS = OSRNewSpatialReference(NULL);        if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )        {            char	*pszPrettyWkt = NULL;            OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );            printf( "Coordinate System is:/n%s/n", pszPrettyWkt );        }        else            printf( "Coordinate System is `%s'/n",                    GDALGetProjectionRef( hDataset ) );        OSRDestroySpatialReference( hSRS );    }/* -------------------------------------------------------------------- *//*      Report Geotransform.                                            *//* -------------------------------------------------------------------- */    if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )    {        printf( "Origin = (%.6f,%.6f)/n",                adfGeoTransform[0], adfGeoTransform[3] );        printf( "Pixel Size = (%.6f,%.6f)/n",                adfGeoTransform[1], adfGeoTransform[5] );    }/* -------------------------------------------------------------------- *//*      Report GCPs.                                                    *//* -------------------------------------------------------------------- */    if( GDALGetGCPCount( hDataset ) > 0 )    {        printf( "GCP Projection = %s/n", GDALGetGCPProjection(hDataset) );        for( i = 0; i < GDALGetGCPCount(hDataset); i++ )        {            const GDAL_GCP	*psGCP;                        psGCP = GDALGetGCPs( hDataset ) + i;            printf( "GCP[%3d]: Id=%s, Info=%s/n"                    "          (%g,%g) -> (%g,%g,%g)/n",                     i, psGCP->pszId, psGCP->pszInfo,                     psGCP->dfGCPPixel, psGCP->dfGCPLine, //.........这里部分代码省略.........

//.........这里部分代码省略.........        CSLDestroy( argv );            GDALDumpOpenDatasets( stderr );        GDALDestroyDriverManager();        CPLDumpSharedList( NULL );        exit( 1 );    }    /* -------------------------------------------------------------------- *//*      Report general info.                                            *//* -------------------------------------------------------------------- */    hDriver = GDALGetDatasetDriver( hDataset );    printf( "Driver: %s/%s/n",            GDALGetDriverShortName( hDriver ),            GDALGetDriverLongName( hDriver ) );    papszFileList = GDALGetFileList( hDataset );    if( CSLCount(papszFileList) == 0 )    {        printf( "Files: none associated/n" );    }    else    {        printf( "Files: %s/n", papszFileList[0] );        for( i = 1; papszFileList[i] != NULL; i++ )            printf( "       %s/n", papszFileList[i] );    }    CSLDestroy( papszFileList );    printf( "Size is %d, %d/n",            GDALGetRasterXSize( hDataset ),             GDALGetRasterYSize( hDataset ) );/* -------------------------------------------------------------------- *//*      Report projection.                                              *//* -------------------------------------------------------------------- */    if( GDALGetProjectionRef( hDataset ) != NULL )    {        OGRSpatialReferenceH  hSRS;        char		      *pszProjection;        pszProjection = (char *) GDALGetProjectionRef( hDataset );        hSRS = OSRNewSpatialReference(NULL);        if( OSRImportFromWkt( hSRS, &pszProjection ) == CE_None )        {            char	*pszPrettyWkt = NULL;            OSRExportToPrettyWkt( hSRS, &pszPrettyWkt, FALSE );            printf( "Coordinate System is:/n%s/n", pszPrettyWkt );            CPLFree( pszPrettyWkt );        }        else            printf( "Coordinate System is `%s'/n",                    GDALGetProjectionRef( hDataset ) );        OSRDestroySpatialReference( hSRS );    }/* -------------------------------------------------------------------- *//*      Report Geotransform.                                            *//* -------------------------------------------------------------------- */    if( GDALGetGeoTransform( hDataset, adfGeoTransform ) == CE_None )

GvSymbolObj *gv_symbol_manager_get_symbol(GvSymbolManager *manager, const char *symbol_name){    gchar   *pszOpenEVHome = NULL;    gchar   *pszSymbolsDir = NULL;    gchar   *pszAbsolutePath = NULL;    gchar   *pszPathSeparator = NULL;    GDALDatasetH hDataset;    GvSymbolObj *poSymbol;    CPLXMLNode  *xml_shape = NULL;    GByte *rgba_buffer;/* -------------------------------------------------------------------- *//*      Lookup the symbol in the hash table, and return it if found.    *//* -------------------------------------------------------------------- */    poSymbol = g_hash_table_lookup( manager->symbol_cache, symbol_name );    if( poSymbol != NULL )        return poSymbol;/* -------------------------------------------------------------------- *//*      We didn't already have it, so try to find a file to load it     *//*      from.                                                           *//* -------------------------------------------------------------------- */#ifndef WIN32    pszPathSeparator = "/";#else    pszPathSeparator = "//";#endif /* WIN32 */    /* validate inputs */    g_return_val_if_fail( manager != NULL, 0 );    g_return_val_if_fail( symbol_name != NULL, 0 );    /* get an absolute path */    if ( !g_path_is_absolute( symbol_name ) )    {        /* check configuration option first */        pszSymbolsDir = g_strdup( gv_manager_get_preference( gv_get_manager(),                                                             "symbols_dir" ) );        /* if not configured check $OPENEV_HOME */        if ( !pszSymbolsDir )        {            pszOpenEVHome = g_getenv( "OPENEV_HOME" );            if( pszOpenEVHome == NULL )                pszOpenEVHome = g_getenv( "OPENEVHOME" );            if ( pszOpenEVHome )                pszSymbolsDir = g_strjoin( pszPathSeparator, pszOpenEVHome,                                           "symbols", NULL );        }                /* get current directory as last resort */        if ( !pszSymbolsDir )            pszSymbolsDir = g_get_current_dir();        pszAbsolutePath = g_strjoin( pszPathSeparator, pszSymbolsDir,                                      symbol_name, NULL );        g_free( pszSymbolsDir );    }    else        pszAbsolutePath = g_strdup( symbol_name );/* -------------------------------------------------------------------- *//*      pszAbsolutePath contains a newly allocated string that is       *//*      suitable for using as a key in the hash table.  If a texture    *//*      is found in the hash table then this string needs to be         *//*      freed.  If one isn't found then the string is used in the       *//*      hash table and should be freed when the associated hash         *//*      table entry is released                                         *//* -------------------------------------------------------------------- */    CPLDebug( "OpenEV",               "gv_symbol_manager_get_symbol(%s) ... need to load.",                   pszAbsolutePath );        /*      * validate path by opening with GDAL and looking for an error      * Disable CPL error handler to supress error reporting     */        CPLErrorReset();    CPLPushErrorHandler( CPLQuietErrorHandler );    hDataset = GDALOpen( pszAbsolutePath, GA_ReadOnly );    CPLPopErrorHandler();        if ( hDataset )    {        rgba_buffer = gdal_to_rgba( hDataset );                if ( rgba_buffer )        {            gv_symbol_manager_inject_raster_symbol(                 manager, symbol_name,                GDALGetRasterXSize( hDataset ),                GDALGetRasterYSize( hDataset ),                rgba_buffer );            CPLFree( rgba_buffer );        }                GDALClose( hDataset );    }//.........这里部分代码省略.........

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

void generateTexture(string fname, GLuint& tex, int bandnum){   if(bandnum <= 0 )   {     bandnum = 1;   }   GDALDataset *poDataset;   GDALAllRegister();   poDataset= (GDALDataset*) GDALOpen(fname.c_str(),GA_ReadOnly);   if(poDataset == NULL)   {      cout << "OUCH!" << endl;      //exit(0);      return;   }   cout << "Data size: " << GDALGetRasterXSize(poDataset) << " " << GDALGetRasterYSize(poDataset) << endl;   GDALRasterBand  *poBand;   int             nBlockXSize, nBlockYSize;   int             bGotMin, bGotMax;   double          adfMinMax[2];  int bands = poDataset->GetRasterCount();   bandnum = bandnum % bands + 1;   if(bandnum > bands)   {      bandnum = 1;   }   poBand = poDataset->GetRasterBand( bandnum );   poBand->GetBlockSize( &nBlockXSize, &nBlockYSize );   printf( "Block=%dx%d Type=%s, ColorInterp=%s/n",            nBlockXSize, nBlockYSize,            GDALGetDataTypeName(poBand->GetRasterDataType()),            GDALGetColorInterpretationName(            poBand->GetColorInterpretation()) );   float max = adfMinMax[0] = poBand->GetMinimum( &bGotMin );   float min = adfMinMax[1] = poBand->GetMaximum( &bGotMax );   if( ! (bGotMin && bGotMax) )      GDALComputeRasterMinMax((GDALRasterBandH)poBand, TRUE, adfMinMax);   int width = poBand->GetXSize();    int height = poBand->GetYSize();   float *pafScanline;   std::cout << "Before allocation" << adfMinMax[0] << " " << adfMinMax[1] << endl;   min = adfMinMax[0];   max = adfMinMax[1];   int dsize = 256;   pafScanline = (float *) CPLMalloc(sizeof(float)*512*512);   vector<vector<float>> out = vector<vector<float>>(height,vector<float> (width,0));   //vector<vector<unsigned char>> texs = vector<vector<unsigned char>>(height,vector<unsigned char> (width,0));   unsigned char texs[512*512];   poBand->RasterIO(GF_Read,0,0,width,height,pafScanline,512,512,GDT_Float32,0,0);   float no = poBand->GetNoDataValue();   cout << "After allocation" << endl;   for(int i = 0; i < 512; i++)   {    for(int j = 0; j < 512; j++)    {      //cout << i << j << endl << pafS;      if(pafScanline[i*width+j] != no)      {        // set tex val        texs[i*512+j] = (unsigned char)(255*((pafScanline[i*512+j] - min)/(max-min)));        //if((int)texs[i*width] < 0)        //cout << (int)texs[i*512 +j] << " " << pafScanline[i*512+j] << " " << no << " " << fname << " " << min << " " << max << endl;      }      else      {        // Set zero val        texs[i*512+j] = 0;        //cout << (int)texs[i*512 +j] << fname << endl;      }      //texs[i*512+j] = 255;      //ut[i][j] = pafScanline[i*width+j];    }   }   CPLFree(pafScanline);   //exit(0);   // Create a texture   glGenTextures(1,&tex);   glBindTexture(GL_TEXTURE_2D,tex);   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);   glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);   glTexImage2D(GL_TEXTURE_2D, 0, GL_RED, 512,512, 0, GL_RED, GL_UNSIGNED_BYTE,texs);   GDALClose( (GDALDatasetH) poDataset);   return;} 

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

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

static void msTransformToGeospatialPDF(imageObj *img, mapObj *map, cairo_renderer *r){  /* We need a GDAL 1.10 PDF driver at runtime, but as far as the C API is concerned, GDAL 1.9 is */  /* largely sufficient. */#if defined(USE_GDAL) && defined(GDAL_VERSION_NUM) && GDAL_VERSION_NUM >= 1900  GDALDatasetH hDS = NULL;  const char* pszGEO_ENCODING = NULL;  GDALDriverH hPDFDriver = NULL;  const char* pszVirtualIO = NULL;  int bVirtualIO = FALSE;  char* pszTmpFilename = NULL;  VSILFILE* fp = NULL;  if (map == NULL)    return;  pszGEO_ENCODING = msGetOutputFormatOption(img->format, "GEO_ENCODING", NULL);  if (pszGEO_ENCODING == NULL)    return;  msGDALInitialize();  hPDFDriver = GDALGetDriverByName("PDF");  if (hPDFDriver == NULL)    return;  /* When compiled against libpoppler, the PDF driver is VirtualIO capable */  /* but not, when it is compiled against libpodofo. */  pszVirtualIO = GDALGetMetadataItem( hPDFDriver, GDAL_DCAP_VIRTUALIO, NULL );  if (pszVirtualIO)    bVirtualIO = CSLTestBoolean(pszVirtualIO);  if (bVirtualIO)    pszTmpFilename = msTmpFile(map, NULL, "/vsimem/mscairopdf/", "pdf");  else    pszTmpFilename = msTmpFile(map, map->mappath, NULL, "pdf");  /* Copy content of outputStream buffer into file */  fp = VSIFOpenL(pszTmpFilename, "wb");  if (fp == NULL) {    msFree(pszTmpFilename);    return;  }  VSIFWriteL(r->outputStream->data, 1, r->outputStream->size, fp);  VSIFCloseL(fp);  fp = NULL;  hDS = GDALOpen(pszTmpFilename, GA_Update);  if ( hDS != NULL ) {    char* pszWKT = msProjectionObj2OGCWKT( &(map->projection) );    if( pszWKT != NULL ) {      double adfGeoTransform[6];      int i;      /* Add user-specified options */      for( i = 0; i < img->format->numformatoptions; i++ ) {        const char* pszOption = img->format->formatoptions[i];        if( strncasecmp(pszOption,"METADATA_ITEM:",14) == 0 ) {          char* pszKey = NULL;          const char* pszValue = CPLParseNameValue(pszOption + 14,                                 &pszKey);          if( pszKey != NULL ) {            GDALSetMetadataItem(hDS, pszKey, pszValue, NULL);            CPLFree(pszKey);          }        }      }      /* We need to rescale the geotransform because GDAL will not necessary */      /* open the PDF with the DPI that was used to generate it */      memcpy(adfGeoTransform, map->gt.geotransform, 6 * sizeof(double));      adfGeoTransform[1] = adfGeoTransform[1] * map->width / GDALGetRasterXSize(hDS);      adfGeoTransform[5] = adfGeoTransform[5] * map->height / GDALGetRasterYSize(hDS);      GDALSetGeoTransform(hDS, adfGeoTransform);      GDALSetProjection(hDS, pszWKT);      msFree( pszWKT );      pszWKT = NULL;      CPLSetThreadLocalConfigOption("GDAL_PDF_GEO_ENCODING", pszGEO_ENCODING);      GDALClose(hDS);      hDS = NULL;      CPLSetThreadLocalConfigOption("GDAL_PDF_GEO_ENCODING", NULL);      /* We need to replace the buffer with the content of the GDAL file */      fp = VSIFOpenL(pszTmpFilename, "rb");      if( fp != NULL ) {        int nFileSize;        VSIFSeekL(fp, 0, SEEK_END);        nFileSize = (int)VSIFTellL(fp);        msBufferResize(r->outputStream, nFileSize);        VSIFSeekL(fp, 0, SEEK_SET);        r->outputStream->size = VSIFReadL(r->outputStream->data, 1, nFileSize, fp);        VSIFCloseL(fp);//.........这里部分代码省略.........

//.........这里部分代码省略.........    nfdr.SetVersion( 78 );    double one_hour, ten_hour, hundred_hour, thousand_hour, woody, herb;    /* What we actually get */    double nfdr_erc;    int nfdr_bi, nfdr_sc;;    double dfDummy;    int nDummy;    GDALDatasetH hFuelDS;    hFuelDS = GDALOpen( "/home/kyle/src/omfrr/trunk/data/"                        "small_fuel40.tif", GA_ReadOnly );    if( hFuelDS == NULL )    {        fprintf( stderr, "Cannot open fuels data for query./n" );        exit( 1 );    }    GDALDatasetH hSlopeDS;    hSlopeDS = GDALOpen( "/home/kyle/src/omfrr/trunk/data/"                         "small_slope.tif", GA_ReadOnly );    if( hSlopeDS == NULL )    {        fprintf( stderr, "Cannot open slope data for query./n" );        exit( 1 );    }    GDALDatasetH hNarrMaskDS;    hNarrMaskDS = GDALOpen( "/home/kyle/src/omfrr/trunk/data/"                            "narr_mask_byte.tif", GA_ReadOnly );    if( hNarrMaskDS == NULL )    {        fprintf( stderr, "Cannot open NARR mask for query./n" );        exit( 1 );    }    int nMaskXSize = GDALGetRasterXSize( hNarrMaskDS );    int nMaskYSize = GDALGetRasterYSize( hNarrMaskDS );    const char *pszAlbersWkt =        "PROJCS[/"USA_Contiguous_Albers_Equal_Area_Conic_USGS_version/","        "GEOGCS[/"GCS_North_American_1983/",DATUM[/"D_North_American_1983/","        "SPHEROID[/"GRS_1980/",6378137.0,298.257222101]],"        "PRIMEM[/"Greenwich/",0.0],UNIT[/"Degree/",0.0174532925199433]],"        "PROJECTION[/"Albers/"],PARAMETER[/"False_Easting/",0.0],"        "PARAMETER[/"False_Northing/",0.0],"        "PARAMETER[/"Central_Meridian/",-96.0],"        "PARAMETER[/"Standard_Parallel_1/",29.5],"        "PARAMETER[/"Standard_Parallel_2/",45.5],"        "PARAMETER[/"Latitude_Of_Origin/",23.0],UNIT[/"Meter/",1.0]]";    /* Assume same srs */    //const char *pszAlbersWkt = NULL;    char out_buf[2048];    int day, year, cause, time;    double x, y, erc;    int fuel, mask, slope, slope_class;    ifstream frisk_stream( "/home/kyle/src/omfrr/trunk/data/"                           "narr_32km_frisk.flt", ios::binary );    FRisk oFRisk;    int skip = 0;    int done = 0;    int wind_speed;    long offset;    int nPixel;    int nLine;    if( bProgress )    {

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

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

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

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