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

gint gv_view_area_print_to_file(GvViewArea *view, int width, int height, const char * filename, const char * format, int is_rgb){    GDALDriverH   driver;    GDALDatasetH  dataset;    gint          errcode;    driver = GDALGetDriverByName( format );    if( driver == NULL )        return -1;    if( is_rgb )        dataset = GDALCreate( driver, filename, width, height, 3, GDT_Byte,                               NULL );    else        dataset = GDALCreate( driver, filename, width, height, 1, GDT_Byte,                               NULL );    if( dataset == NULL )        return -1;    errcode = gv_view_area_render_to_func( view, width, height,                                            print_handler, dataset );    GDALClose( dataset );    print_handler( NULL, NULL );        return errcode;}

/* Makes a copy of a dataset, and opens it for writing.. */GDALDatasetHmake_me_a_sandwitch (GDALDatasetH * in_dataset, const char *copy_file_name){  char **papszOptions = NULL;  const char *pszFormat = "GTiff";  double adfGeoTransform[6];  GDALDriverH hDriver;  GDALDatasetH out_gdalfile;  hDriver = GDALGetDriverByName (pszFormat);  papszOptions = CSLSetNameValue (papszOptions, "TILED", "YES");  papszOptions = CSLSetNameValue (papszOptions, "COMPRESS", "DEFLATE");  /*Perhaps controversal - default to bigtiff... */  /*papszOptions = CSLSetNameValue( papszOptions, "BIGTIFF", "YES" ); */  /*return GDALCreateCopy( hDriver, copy_file_name, *in_dataset, FALSE, papszOptions, NULL, NULL ); */  out_gdalfile = GDALCreate (hDriver, copy_file_name,			     GDALGetRasterXSize (*in_dataset),			     GDALGetRasterYSize (*in_dataset),			     GDALGetRasterCount (*in_dataset),			     GDT_Byte, papszOptions);  /* Set geotransform */  GDALGetGeoTransform (*in_dataset, adfGeoTransform);  GDALSetGeoTransform (out_gdalfile, adfGeoTransform);  /* Set projection */  GDALSetProjection (out_gdalfile, GDALGetProjectionRef (*in_dataset));  return out_gdalfile;}

int sat_save_ch(s_sat * sat, const char * fname, unsigned ch){	try;		int ret = 0;		unsigned height = sat->height, width = sat->width;		GDALDatasetH ds = NULL;		throw_null((ds = GDALCreate(drv_tiff, fname, width, height, 1, GDT_Byte, NULL)));		throw((GDALSetProjection(ds, sat->proj_ref) == CE_Failure));		throw((GDALSetGeoTransform(ds, sat->gt_coef) == CE_Failure));		throw(GDALDatasetRasterIO(ds, GF_Write, 0, 0, width, height, sat->pixel[ch], width, height, GDT_Byte, 1, NULL, 0, 0, 0) == CE_Failure);	catch;		ret = -1;	finally;		if(ds != NULL)			GDALClose(ds);	return ret;}

void write_image_file(const char * filename, int nx, int ny, float * image){    int i;    int nBands = 1;    printf("writing image file %s (%d,%d)/n", filename, nx, ny);    unsigned char * buf = (unsigned char*) malloc(nx*ny*sizeof(unsigned char));    for (i = 0; i < nx*ny; i++) {        buf[i] = (unsigned char) image[i];    }    GDALAllRegister();    GDALDriverH driver = GDALGetDriverByName("GTiff");    if (driver == NULL) {        exit(1);    }    GDALDatasetH dataset = GDALCreate(driver, filename, nx, ny, nBands, GDT_Byte, NULL);    if (dataset == NULL) {        fprintf(stderr, "write_image_file: failed to open file %s/n", filename);    }    GDALDatasetRasterIO(dataset, GF_Write, 0, 0, nx, ny, buf, nx, ny, GDT_Byte,                        nBands, NULL, 1, nx, nx*ny);    GDALClose(dataset);    free(buf);}

//void CUtils::createNewByteGeoTIFF(const char* fileName, int bands, int rows, int cols, double adfGeoTransform[6], char szProjection[512], byte fillData, byte noDataValue)void CUtils::createNewByteGeoTIFF(const char* fileName, int bands, int rows, int cols, double adfGeoTransform[6], const char * szProjection, byte fillData, byte noDataValue){	char **papszOptions = NULL;    GDALDriverH hDriver;    GDALRasterBandH hDataset;//	GDALRasterBandH hBand;	if( (hDriver = GDALGetDriverByName("GTiff")) != NULL)	{		fprintf(stderr, "Create image %s.../n", fileName);        if( (hDataset = GDALCreate( hDriver, fileName, cols, rows, bands, GDT_Byte, papszOptions )) !=NULL )		{			GDALSetGeoTransform(hDataset, adfGeoTransform );			GDALSetProjection(hDataset, szProjection );			/*			int pr = CUtils::progress_ln_ex(stderr, 0, 0, START_PROGRESS);			for(int band = 1; band<=bands; band++)			{				if( (hBand = GDALGetRasterBand(hDataset, band)) != NULL )				{					byte *pline = (byte *)CPLMalloc(sizeof(byte)*cols);					for(int i=0; i<cols; i++) pline[i] = fillData;					for(int i=0; i<rows; i++) GDALRasterIO(hBand, GF_Write, 0, i, cols, 1, pline, cols, 1, GDT_Byte, 0, 0 );					CPLFree(pline);					GDALSetRasterNoDataValue(hBand, noDataValue);				}				pr = CUtils::progress_ln_ex(stderr, band-1, bands, pr);			}			CUtils::progress_ln_ex(stderr, 0, 0, END_PROGRESS);			*/			GDALClose(hDataset);		}	}}

//辐射校正处理=====================================================================================================================================================================================//绝对辐射校正long QPDLevel1Process::Level1Proc_RadiationAbsolute(const char* pathImg, const char* pathImgRad, const char* pathAbsRegFile){	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");	//中文路径	GDALAllRegister();	long lError = 0;	unsigned short *imgBuffer = NULL;	//影像数据	float* parametersA = NULL, *parametersB = NULL, *parametersAux = NULL;//校正系数	GDALDatasetH m_dataset = GDALOpen(pathImg, GA_ReadOnly);	int xsize = GDALGetRasterXSize(m_dataset);	int ysize = GDALGetRasterYSize(m_dataset);	int bands = GDALGetRasterCount(m_dataset);	char **papszOptions = NULL;	papszOptions = CSLSetNameValue(papszOptions, "INTERLEAVE", "BAND");	GDALDatasetH m_datasetdst = GDALCreate(GDALGetDriverByName("GTiff"), pathImgRad, xsize, ysize, bands, GDT_UInt16, papszOptions);	//int nSamples, nLines, nLevels;	//LevelProc_GetParameterInfo(pathImgRad, nSamples, nLines, nLevels);	try	{		parametersA = new float[bands];		parametersB = new float[bands];		imgBuffer = new unsigned short[xsize*ysize];	}	catch (bad_alloc)	{		printf("allocate memory error/n");		exit(-1);	}	Level1Proc_AbsoluteParameters(pathAbsRegFile, parametersA, parametersB);	for (int i = 0; i < bands; i++)	{		GDALRasterIO(GDALGetRasterBand(m_dataset, i + 1), GF_Read, 0, 0, xsize, ysize, imgBuffer, xsize, ysize, GDT_UInt16, 0, 0);		for (int j = 0; j < ysize; j++)		{			for (int k = 0; k < xsize; k++)			{				//扩大100倍精度为0.01				imgBuffer[k] = (unsigned short)((imgBuffer[j*xsize + k] * parametersA[i] + parametersB[i]) * 100);			}		}		GDALRasterIO(GDALGetRasterBand(m_datasetdst, i + 1), GF_Write, 0, 0, xsize, ysize, imgBuffer, xsize, ysize, GDT_UInt16, 0, 0);	}	delete[]parametersA; parametersA = NULL;	delete[]parametersB; parametersB = NULL;	delete[]imgBuffer;	 imgBuffer = NULL;	GDALClose(m_dataset);	GDALClose(m_datasetdst);	return lError;}

GDALDataset* geGdalVSI::VsiGdalCreateWrap(GDALDriverH hdriver,                                          std::string* const vsifile,                                          int nx, int ny, int nbands,                                          GDALDataType bandtype,                                          char **papszoptions) {  GDALDatasetH hvsi_ds;  khMutex &mutex = GetMutex();  khLockGuard lock(mutex);  *vsifile = UniqueVSIFilename();  hvsi_ds = GDALCreate(hdriver, (*vsifile).c_str(), nx, ny, nbands,                       bandtype, papszoptions);  return static_cast<GDALDataset*>(hvsi_ds);}

GDALDatasetH QgsRelief::openOutputFile( GDALDatasetH inputDataset, GDALDriverH outputDriver ){  if ( !inputDataset )  {    return nullptr;  }  int xSize = GDALGetRasterXSize( inputDataset );  int ySize = GDALGetRasterYSize( inputDataset );  //open output file  char **papszOptions = nullptr;  //use PACKBITS compression for tiffs by default  papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );  //create three band raster (reg, green, blue)  GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), xSize, ySize, 3, GDT_Byte, papszOptions );  if ( !outputDataset )  {    return outputDataset;  }  //get geotransform from inputDataset  double geotransform[6];  if ( GDALGetGeoTransform( inputDataset, geotransform ) != CE_None )  {    GDALClose( outputDataset );    return nullptr;  }  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;}

GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver ){  //open output file  char **papszOptions = NULL;  GDALDatasetH outputDataset = GDALCreate( outputDriver, mOutputFile.toLocal8Bit().data(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );  if ( outputDataset == NULL )  {    return outputDataset;  }  //assign georef information  double geotransform[6];  outputGeoTransform( geotransform );  GDALSetGeoTransform( outputDataset, geotransform );  return outputDataset;}

GDALDatasetH QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver ){  //open output file  char **papszOptions = nullptr;  GDALDatasetH outputDataset = GDALCreate( outputDriver, TO8F( mOutputFile ), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions );  if ( !outputDataset )  {    return outputDataset;  }  //assign georef information  double geotransform[6];  outputGeoTransform( geotransform );  GDALSetGeoTransform( outputDataset, geotransform );  return outputDataset;}

gdal::dataset_unique_ptr QgsRasterCalculator::openOutputFile( GDALDriverH outputDriver ){  //open output file  char **papszOptions = nullptr;  gdal::dataset_unique_ptr outputDataset( GDALCreate( outputDriver, mOutputFile.toUtf8().constData(), mNumOutputColumns, mNumOutputRows, 1, GDT_Float32, papszOptions ) );  if ( !outputDataset )  {    return nullptr;  }  //assign georef information  double geotransform[6];  outputGeoTransform( geotransform );  GDALSetGeoTransform( outputDataset.get(), geotransform );  return outputDataset;}

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

static dErr JakoGDALDatasetCreateMem(OGRSpatialReferenceH ref,const double geo[6],dInt n,dInt nlines,GDALDataType dtype,GDALDatasetH *dset,void *bandmem){  char *wkt;  GDALDriverH memdriver;  CPLErr cplerr;  OGRErr oerr;  dErr err;  dFunctionBegin;  oerr = OSRExportToWkt(ref,&wkt);dOGRCHK(oerr);  memdriver = GDALGetDriverByName("MEM");  *dset = GDALCreate(memdriver,"MEM:::",n,nlines,0,dtype,NULL);  cplerr = GDALSetProjection(*dset,wkt);dCPLCHK(cplerr);  cplerr = GDALSetGeoTransform(*dset,(double*)geo);dCPLCHK(cplerr); /* const-incorrect interface */  OGRFree(wkt);  if (bandmem) {err = JakoGDALMemAddBand(*dset,GDT_Float64,&bandmem);dCHK(err);}  dFunctionReturn(0);}

s_sat * sat_rasterize_copy(s_sat * c_sat, OGRGeometryH geometry){	try;		int band = 1;		unsigned height, width, height_width;		double burn_value = 255;		GDALDatasetH ds = NULL;		s_sat * sat = NULL;		throw_null((sat = sat_init(1)));		height = sat->height = c_sat->height;		width = sat->width = c_sat->width;		height_width = height * width;		sat->proj_ref = strdup(c_sat->proj_ref);		sat->ch_num = 0;		memcpy(sat->gt_coef, c_sat->gt_coef, sizeof(double) * 6);		throw_null((sat->pixel[0] = sfire_alloc(sizeof(uint8_t), 1, height_width)));		sat->ch_num = 1;		throw_null((ds = GDALCreate(drv_r_mem, "", width, height, 1, GDT_Byte, NULL)));		throw((GDALSetProjection(ds, sat->proj_ref) == CE_Failure));		throw((GDALSetGeoTransform(ds, sat->gt_coef) == CE_Failure));		throw((GDALRasterizeGeometries(ds, 1, & band, 1, & geometry, NULL, NULL, & burn_value, NULL, NULL, NULL) == CE_Failure));		throw(GDALDatasetRasterIO(ds, GF_Read, 0, 0, width, height, sat->pixel[0], width, height, GDT_Byte, 1, NULL, 0, 0, 0) == CE_Failure);	catch;		sat_destroy(sat);		sat = NULL;	finally;		if(ds != NULL)			GDALClose(ds);	return sat;}

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;}

intmakeGeotiff (struct deminfo *d0, char *outpath,int nodata){  GDALAllRegister ();  GDALDataType band_type = GDT_Float32;  int bands = 1;  int dsn_xsize = (d0->highx - d0->lowx + 1);  int dsn_ysize = (d0->highy - d0->lowy + 1);  char **papszCreateOptions = NULL;  papszCreateOptions = CSLSetNameValue (papszCreateOptions, "PROFILE", "GeoTIFF");  //papszCreateOptions = CSLSetNameValue( papszCreateOptions, "TFW", "YES" );  //papszCreateOptions = CSLSetNameValue (papszCreateOptions, "INTERLEAVE", "PIXEL");  //papszCreateOptions = CSLSetNameValue (papszCreateOptions, "TILED", "YES");  //papszCreateOptions = CSLSetNameValue (papszCreateOptions, "COMPRESS", "LZW");  GDALDriverH hDriver = GDALGetDriverByName ("GTiff");  GDALDatasetH hDsnDS = GDALCreate (hDriver, outpath, dsn_xsize, dsn_ysize, bands, band_type, papszCreateOptions);  double dsnGeoTransform[6];  dsnGeoTransform[0] = d0->W;  dsnGeoTransform[1] = (d0->E - d0->W) / dsn_xsize;  dsnGeoTransform[2] = 0;  dsnGeoTransform[3] = d0->N;  dsnGeoTransform[4] = 0;  dsnGeoTransform[5] = -1.0 * (d0->N - d0->S) / dsn_ysize;  GDALSetGeoTransform (hDsnDS, dsnGeoTransform);  char pszSRS_WKT[1024] = "GEOGCS[/"JGD2000/", DATUM[/"Japanese Geodetic Datum 2000/", SPHEROID[/"GRS 1980/", 6378137.0, 298.257222101, AUTHORITY[/"EPSG/",/"7019/"]], TOWGS84[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],AUTHORITY[/"EPSG/",/"6612/"]], PRIMEM[/"Greenwich/", 0.0, AUTHORITY[/"EPSG/",/"8901/"]], UNIT[/"degree/", 0.017453292519943295], AXIS[/"Geodetic longitude/", EAST], AXIS[/"Geodetic latitude/", NORTH], AUTHORITY[/"EPSG/",/"4612/"]]";  GDALSetProjection (hDsnDS, pszSRS_WKT);  GDALRasterBandH t_band = GDALGetRasterBand (hDsnDS, 1);  if(nodata==1){   GDALSetRasterNoDataValue (t_band, -9999);  }  GDALRasterIO (t_band, GF_Write, 0, 0, dsn_xsize, dsn_ysize, d0->alti, dsn_xsize, dsn_ysize, band_type, 0, 0);  CSLDestroy (papszCreateOptions);  GDALClose (hDsnDS);  return 0;}

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

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 );}

//.........这里部分代码省略.........    if( hOutDS != NULL && papszCreationOptions != NULL)    {        CPLError(CE_Warning, CPLE_AppDefined,                  "Warning: creation options are ignored when writing to an existing file.");    }/* -------------------------------------------------------------------- *//*      Do we need to create output file?                               *//* -------------------------------------------------------------------- */    if( hOutDS == NULL )    {        GDALDriverH hDriver = GDALGetDriverByName( pszDriverName );        if (hDriver == NULL)            exit(1);        if (!bQuiet && !bFormatExplicitelySet)            CheckExtensionConsistency(pszOutFile, pszDriverName);        if (bSetAlpha)        {            /***** fixme there should be a way to preserve alpha band data not in the collar *****/            if (nBands == 4)                nBands --;            else                nDstBands ++;        }        if (bSetMask)        {            if (nBands == 4)                nDstBands = nBands = 3;        }        hOutDS = GDALCreate( hDriver, pszOutFile,                              nXSize, nYSize, nDstBands, GDT_Byte,                              papszCreationOptions );        if( hOutDS == NULL )            exit( 1 );        double adfGeoTransform[6];        if( GDALGetGeoTransform( hInDS, adfGeoTransform ) == CE_None )        {            GDALSetGeoTransform( hOutDS, adfGeoTransform );            GDALSetProjection( hOutDS, GDALGetProjectionRef( hInDS ) );        }    }    else    {        if (bSetAlpha)        {            if (nBands != 4 &&                (nBands < 2 ||                 GDALGetRasterColorInterpretation(GDALGetRasterBand(hOutDS, nBands)) != GCI_AlphaBand))            {                CPLError(CE_Failure, CPLE_AppDefined,                        "Last band is not an alpha band.");                exit(1);            }            nBands --;        }        if (bSetMask)        {            if (nBands == 4)

/** * /private /memberof mapcache_source_gdal * /sa mapcache_source::render_metatile() */void _mapcache_source_gdal_render_metatile(mapcache_context *ctx, mapcache_metatile *tile){  mapcache_source_gdal *gdal = (mapcache_source_gdal*)tile->tile.tileset->source;  char *srcSRS = "", *dstSRS;  mapcache_buffer *data = mapcache_buffer_create(0,ctx->pool);  GC_CHECK_ERROR(ctx);  GDALDatasetH  hDataset;  GDALAllRegister();  OGRSpatialReferenceH hSRS;  CPLErrorReset();  hSRS = OSRNewSpatialReference( NULL );  if( OSRSetFromUserInput( hSRS, tile->tile.grid->srs ) == OGRERR_NONE )    OSRExportToWkt( hSRS, &dstSRS );  else {    ctx->set_error(ctx,MAPCACHE_SOURCE_GDAL_ERROR,"failed to parse gdal srs %s",tile->tile.grid->srs);    return;  }  OSRDestroySpatialReference( hSRS );  hDataset = GDALOpen( gdal->datastr, GA_ReadOnly );  if( hDataset == NULL ) {    ctx->set_error(ctx,MAPCACHE_SOURCE_GDAL_ERROR,"GDAL failed to open %s",gdal->datastr);    return;  }  /* -------------------------------------------------------------------- */  /*      Check that there's at least one raster band                     */  /* -------------------------------------------------------------------- */  if ( GDALGetRasterCount(hDataset) == 0 ) {    ctx->set_error(ctx,MAPCACHE_SOURCE_GDAL_ERROR,"raster %s has no bands",gdal->datastr);    return;  }  if( GDALGetProjectionRef( hDataset ) != NULL      && strlen(GDALGetProjectionRef( hDataset )) > 0 )    srcSRS = apr_pstrdup(ctx->pool,GDALGetProjectionRef( hDataset ));  else if( GDALGetGCPProjection( hDataset ) != NULL           && strlen(GDALGetGCPProjection(hDataset)) > 0           && GDALGetGCPCount( hDataset ) > 1 )    srcSRS = apr_pstrdup(ctx->pool,GDALGetGCPProjection( hDataset ));  GDALDriverH hDriver = GDALGetDriverByName( "MEM" );  GDALDatasetH hDstDS;  /* -------------------------------------------------------------------- */  /*      Create a transformation object from the source to               */  /*      destination coordinate system.                                  */  /* -------------------------------------------------------------------- */  void *hTransformArg =    GDALCreateGenImgProjTransformer( hDataset, srcSRS,                                     NULL, dstSRS,                                     TRUE, 1000.0, 0 );  if( hTransformArg == NULL ) {    ctx->set_error(ctx,MAPCACHE_SOURCE_GDAL_ERROR,"gdal failed to create SRS transformation object");    return;  }  /* -------------------------------------------------------------------- */  /*      Get approximate output definition.                              */  /* -------------------------------------------------------------------- */  int nPixels, nLines;  double adfDstGeoTransform[6];  if( GDALSuggestedWarpOutput( hDataset,                               GDALGenImgProjTransform, hTransformArg,                               adfDstGeoTransform, &nPixels, &nLines )      != CE_None ) {    ctx->set_error(ctx,MAPCACHE_SOURCE_GDAL_ERROR,"gdal failed to create suggested warp output");    return;  }  GDALDestroyGenImgProjTransformer( hTransformArg );  double dfXRes = (tile->bbox[2] - tile->bbox[0]) / tile->sx;  double dfYRes = (tile->bbox[3] - tile->bbox[1]) / tile->sy;  adfDstGeoTransform[0] = tile->bbox[0];  adfDstGeoTransform[3] = tile->bbox[3];  adfDstGeoTransform[1] = dfXRes;  adfDstGeoTransform[5] = -dfYRes;  hDstDS = GDALCreate( hDriver, "tempd_gdal_image", tile->sx, tile->sy, 4, GDT_Byte, NULL );  /* -------------------------------------------------------------------- */  /*      Write out the projection definition.                            */  /* -------------------------------------------------------------------- */  GDALSetProjection( hDstDS, dstSRS );  GDALSetGeoTransform( hDstDS, adfDstGeoTransform );  char               **papszWarpOptions = NULL;  papszWarpOptions = CSLSetNameValue( papszWarpOptions, "INIT", "0" );  /* -------------------------------------------------------------------- */  /*      Create a transformation object from the source to               *///.........这里部分代码省略.........

staticGDALDatasetH CreateOutputDataset(std::vector<OGRLayerH> ahLayers,                                 OGRSpatialReferenceH hSRS,                                 int bGotBounds, OGREnvelope sEnvelop,                                 GDALDriverH hDriver, const char* pszDest,                                 int nXSize, int nYSize, double dfXRes, double dfYRes,                                 int bTargetAlignedPixels,                                 int nBandCount, GDALDataType eOutputType,                                 char** papszCreationOptions, std::vector<double> adfInitVals,                                 int bNoDataSet, double dfNoData){    int bFirstLayer = TRUE;    char* pszWKT = NULL;    GDALDatasetH hDstDS = NULL;    unsigned int i;    for( i = 0; i < ahLayers.size(); i++ )    {        OGRLayerH hLayer = ahLayers[i];        if (!bGotBounds)        {            OGREnvelope sLayerEnvelop;            if (OGR_L_GetExtent(hLayer, &sLayerEnvelop, TRUE) != OGRERR_NONE)            {                CPLError(CE_Failure, CPLE_AppDefined, "Cannot get layer extent");                return NULL;            }            /* Voluntarily increase the extent by a half-pixel size to avoid */            /* missing points on the border */            if (!bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)            {                sLayerEnvelop.MinX -= dfXRes / 2;                sLayerEnvelop.MaxX += dfXRes / 2;                sLayerEnvelop.MinY -= dfYRes / 2;                sLayerEnvelop.MaxY += dfYRes / 2;            }            if (bFirstLayer)            {                sEnvelop.MinX = sLayerEnvelop.MinX;                sEnvelop.MinY = sLayerEnvelop.MinY;                sEnvelop.MaxX = sLayerEnvelop.MaxX;                sEnvelop.MaxY = sLayerEnvelop.MaxY;                if (hSRS == NULL)                    hSRS = OGR_L_GetSpatialRef(hLayer);                bFirstLayer = FALSE;            }            else            {                sEnvelop.MinX = MIN(sEnvelop.MinX, sLayerEnvelop.MinX);                sEnvelop.MinY = MIN(sEnvelop.MinY, sLayerEnvelop.MinY);                sEnvelop.MaxX = MAX(sEnvelop.MaxX, sLayerEnvelop.MaxX);                sEnvelop.MaxY = MAX(sEnvelop.MaxY, sLayerEnvelop.MaxY);            }        }        else        {            if (bFirstLayer)            {                if (hSRS == NULL)                    hSRS = OGR_L_GetSpatialRef(hLayer);                bFirstLayer = FALSE;            }        }    }    if (dfXRes == 0 && dfYRes == 0)    {        dfXRes = (sEnvelop.MaxX - sEnvelop.MinX) / nXSize;        dfYRes = (sEnvelop.MaxY - sEnvelop.MinY) / nYSize;    }    else if (bTargetAlignedPixels && dfXRes != 0 && dfYRes != 0)    {        sEnvelop.MinX = floor(sEnvelop.MinX / dfXRes) * dfXRes;        sEnvelop.MaxX = ceil(sEnvelop.MaxX / dfXRes) * dfXRes;        sEnvelop.MinY = floor(sEnvelop.MinY / dfYRes) * dfYRes;        sEnvelop.MaxY = ceil(sEnvelop.MaxY / dfYRes) * dfYRes;    }    double adfProjection[6];    adfProjection[0] = sEnvelop.MinX;    adfProjection[1] = dfXRes;    adfProjection[2] = 0;    adfProjection[3] = sEnvelop.MaxY;    adfProjection[4] = 0;    adfProjection[5] = -dfYRes;    if (nXSize == 0 && nYSize == 0)    {        nXSize = (int)(0.5 + (sEnvelop.MaxX - sEnvelop.MinX) / dfXRes);        nYSize = (int)(0.5 + (sEnvelop.MaxY - sEnvelop.MinY) / dfYRes);    }    hDstDS = GDALCreate(hDriver, pszDest, nXSize, nYSize,//.........这里部分代码省略.........

//.........这里部分代码省略........./* -------------------------------------------------------------------- *//*      Iterate over blocks to add data into raster and metadata tables *//* -------------------------------------------------------------------- */    OGR_DS_ExecuteSQL(hDS, "BEGIN", NULL, NULL);        CPLErr eErr = CE_None;    int nBlockXOff, nBlockYOff;    for(nBlockYOff=0;eErr == CE_None && nBlockYOff<nYBlocks;nBlockYOff++)    {        for(nBlockXOff=0;eErr == CE_None && nBlockXOff<nXBlocks;nBlockXOff++)        {/* -------------------------------------------------------------------- *//*      Create in-memory tile                                           *//* -------------------------------------------------------------------- */            int nReqXSize = nBlockXSize, nReqYSize = nBlockYSize;            if ((nBlockXOff+1) * nBlockXSize > nOvrXSize)                nReqXSize = nOvrXSize - nBlockXOff * nBlockXSize;            if ((nBlockYOff+1) * nBlockYSize > nOvrYSize)                nReqYSize = nOvrYSize - nBlockYOff * nBlockYSize;                        eErr = RasterIO(GF_Read,                            nBlockXOff * nBlockXSize * nOvrFactor,                            nBlockYOff * nBlockYSize * nOvrFactor,                            nReqXSize * nOvrFactor, nReqYSize * nOvrFactor,                            pabyMEMDSBuffer, nReqXSize, nReqYSize,                            eDataType, nBands, NULL,                            0, 0, 0);            if (eErr != CE_None)            {                break;            }                        GDALDatasetH hMemDS = GDALCreate(hMemDriver, "MEM:::",                                              nReqXSize, nReqYSize, 0,                                               eDataType, NULL);            if (hMemDS == NULL)            {                eErr = CE_Failure;                break;            }                        int iBand;            for(iBand = 0; iBand < nBands; iBand ++)            {                char** papszOptions = NULL;                char szTmp[64];                memset(szTmp, 0, sizeof(szTmp));                CPLPrintPointer(szTmp,                                pabyMEMDSBuffer + iBand * nDataTypeSize *                                nReqXSize * nReqYSize, sizeof(szTmp));                papszOptions = CSLSetNameValue(papszOptions, "DATAPOINTER", szTmp);                GDALAddBand(hMemDS, eDataType, papszOptions);                CSLDestroy(papszOptions);            }                        GDALDatasetH hOutDS = GDALCreateCopy(hTileDriver,                                        osTempFileName.c_str(), hMemDS, FALSE,                                        papszTileDriverOptions, NULL, NULL);            GDALClose(hMemDS);            if (hOutDS)                GDALClose(hOutDS);            else            {                eErr = CE_Failure;

void Dust::MakeGrid(WindNinjaInputs &input, AsciiGrid<double> &grid){    /*------------------------------------------*/    /* Open grid as a GDAL dataset              */    /*------------------------------------------*/    int nXSize = grid.get_nCols();    int nYSize = grid.get_nRows();        GDALDriverH hDriver = GDALGetDriverByName( "MEM" );            GDALDatasetH hMemDS = GDALCreate(hDriver, "", nXSize, nYSize, 1, GDT_Float64, NULL);        double *padfScanline;    padfScanline = new double[nXSize];        double adfGeoTransform[6];    adfGeoTransform[0] = grid.get_xllCorner();    adfGeoTransform[1] = grid.get_cellSize();    adfGeoTransform[2] = 0;    adfGeoTransform[3] = grid.get_yllCorner()+(grid.get_nRows()*grid.get_cellSize());    adfGeoTransform[4] = 0;    adfGeoTransform[5] = -grid.get_cellSize();            char* pszDstWKT = (char*)grid.prjString.c_str();    GDALSetProjection(hMemDS, pszDstWKT);    GDALSetGeoTransform(hMemDS, adfGeoTransform);            GDALRasterBandH hBand = GDALGetRasterBand( hMemDS, 1 );            GDALSetRasterNoDataValue(hBand, -9999.0);            for(int i=nYSize-1; i>=0; i--)    {        for(int j=0; j<nXSize; j++)        {              padfScanline[j] = grid.get_cellValue(nYSize-1-i, j);        }        GDALRasterIO(hBand, GF_Write, 0, i, nXSize, 1, padfScanline, nXSize,                     1, GDT_Float64, 0, 0);    }        /*------------------------------------------*/    /* Get the geometry info                    */    /*------------------------------------------*/        OGRDataSourceH hOGRDS = 0;    hOGRDS = OGROpen(input.dustFilename.c_str(), FALSE, 0);    if(hOGRDS == NULL)    {        throw std::runtime_error("Could not open the fire perimeter file '" +              input.dustFilename + "' for reading.");    }    OGRLayer *poLayer;    OGRFeature *poFeature;    OGRGeometry *poGeo;        poLayer = (OGRLayer*)OGR_DS_GetLayer(hOGRDS, 0);    poLayer->ResetReading();    poFeature = poLayer->GetNextFeature();    poGeo = poFeature->GetGeometryRef();    OGRGeometryH hPolygon = (OGRGeometryH) poGeo;    /* -------------------------------------------------------------------- */    /*  Check for same CRS in fire perimeter and DEM files                  */    /* -------------------------------------------------------------------- */    char *pszSrcWKT;    OGRSpatialReference *poSrcSRS, oDstSRS;    poSrcSRS = poLayer->GetSpatialRef(); //shapefile CRS    poSrcSRS->exportToWkt( &pszSrcWKT );    //printf("CRS of DEM is:/n %s/n", pszDstWKT);    //printf("WKT CRS of .shp is:/n %s/n", pszSrcWKT);        oDstSRS.importFromWkt( &pszDstWKT );        char *pszDstProj4, *pszSrcProj4;    oDstSRS.exportToProj4( &pszDstProj4 );    poSrcSRS->exportToProj4( &pszSrcProj4 );        //printf("proj4 of .shp is:/n %s/n", pszSrcProj4);    //printf("proj4 of dem is:/n %s/n", pszDstProj4);        /* -------------------------------------------------------------------- */    /*  If the CRSs are not equal, convert shapefile CRS to DEM CRS         */    /* -------------------------------------------------------------------- */    GDALTransformerFunc pfnTransformer = NULL;    if( !EQUAL( pszSrcProj4, pszDstProj4 ) ){ //tranform shp CRS to DEM CRS        poGeo->transformTo(&oDstSRS);    }         /* -------------------------------------------------------------------- */    /*  Rasterize the shapefile                                             */    /* -------------------------------------------------------------------- */        int nTargetBand = 1;    double BurnValue = 1.0;    CPLErr eErr;//.........这里部分代码省略.........

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

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

CPLErrGDALWarpCutlineMasker( void *pMaskFuncArg,                       CPL_UNUSED int nBandCount,                       CPL_UNUSED GDALDataType eType,                       int nXOff, int nYOff, int nXSize, int nYSize,                       GByte ** /*ppImageData */,                       int bMaskIsFloat, void *pValidityMask ){    GDALWarpOptions *psWO = (GDALWarpOptions *) pMaskFuncArg;    float *pafMask = (float *) pValidityMask;    CPLErr eErr;    GDALDriverH hMemDriver;    if( nXSize < 1 || nYSize < 1 )        return CE_None;/* -------------------------------------------------------------------- *//*      Do some minimal checking.                                       *//* -------------------------------------------------------------------- */    if( !bMaskIsFloat )    {        CPLAssert( FALSE );        return CE_Failure;    }    if( psWO == NULL || psWO->hCutline == NULL )    {        CPLAssert( FALSE );        return CE_Failure;    }    hMemDriver = GDALGetDriverByName("MEM");    if (hMemDriver == NULL)    {        CPLError(CE_Failure, CPLE_AppDefined, "GDALWarpCutlineMasker needs MEM driver");        return CE_Failure;    }/* -------------------------------------------------------------------- *//*      Check the polygon.                                              *//* -------------------------------------------------------------------- */    OGRGeometryH hPolygon = (OGRGeometryH) psWO->hCutline;    OGREnvelope  sEnvelope;    if( wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbPolygon        && wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbMultiPolygon )    {        CPLAssert( FALSE );        return CE_Failure;    }    OGR_G_GetEnvelope( hPolygon, &sEnvelope );    if( sEnvelope.MaxX + psWO->dfCutlineBlendDist < nXOff        || sEnvelope.MinX - psWO->dfCutlineBlendDist > nXOff + nXSize        || sEnvelope.MaxY + psWO->dfCutlineBlendDist < nYOff        || sEnvelope.MinY - psWO->dfCutlineBlendDist > nYOff + nYSize )    {        // We are far from the blend line - everything is masked to zero.        // It would be nice to realize no work is required for this whole        // chunk!        memset( pafMask, 0, sizeof(float) * nXSize * nYSize );        return CE_None;    }/* -------------------------------------------------------------------- *//*      Create a byte buffer into which we can burn the                 *//*      mask polygon and wrap it up as a memory dataset.                *//* -------------------------------------------------------------------- */    GByte *pabyPolyMask = (GByte *) CPLCalloc( nXSize, nYSize );    GDALDatasetH hMemDS;    double adfGeoTransform[6] = { 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };    char szDataPointer[100];    char *apszOptions[] = { szDataPointer, NULL };    memset( szDataPointer, 0, sizeof(szDataPointer) );    sprintf( szDataPointer, "DATAPOINTER=" );    CPLPrintPointer( szDataPointer+strlen(szDataPointer),                     pabyPolyMask,                      sizeof(szDataPointer) - strlen(szDataPointer) );    hMemDS = GDALCreate( hMemDriver, "warp_temp",                          nXSize, nYSize, 0, GDT_Byte, NULL );    GDALAddBand( hMemDS, GDT_Byte, apszOptions );    GDALSetGeoTransform( hMemDS, adfGeoTransform );/* -------------------------------------------------------------------- *//*      Burn the polygon into the mask with 1.0 values.                 *//* -------------------------------------------------------------------- */    int nTargetBand = 1;    double dfBurnValue = 255.0;    int    anXYOff[2];    char   **papszRasterizeOptions = NULL;        if( CSLFetchBoolean( psWO->papszWarpOptions, "CUTLINE_ALL_TOUCHED", FALSE ))        papszRasterizeOptions =             CSLSetNameValue( papszRasterizeOptions, "ALL_TOUCHED", "TRUE" );//.........这里部分代码省略.........

//.........这里部分代码省略.........    int nTotalBlocks = nXBlocks * nYBlocks;    char** papszTileDriverOptions = RasterliteGetTileDriverOptions(papszOptions);        OGR_DS_ExecuteSQL(hDS, "BEGIN", NULL, NULL);        CPLErr eErr = CE_None;    int nBlockXOff, nBlockYOff;    for(nBlockYOff=0;eErr == CE_None && nBlockYOff<nYBlocks;nBlockYOff++)    {        for(nBlockXOff=0;eErr == CE_None && nBlockXOff<nXBlocks;nBlockXOff++)        {/* -------------------------------------------------------------------- *//*      Create in-memory tile                                           *//* -------------------------------------------------------------------- */            int nReqXSize = nBlockXSize, nReqYSize = nBlockYSize;            if ((nBlockXOff+1) * nBlockXSize > nXSize)                nReqXSize = nXSize - nBlockXOff * nBlockXSize;            if ((nBlockYOff+1) * nBlockYSize > nYSize)                nReqYSize = nYSize - nBlockYOff * nBlockYSize;            eErr = poSrcDS->RasterIO(GF_Read,                                     nBlockXOff * nBlockXSize,                                     nBlockYOff * nBlockYSize,                                     nReqXSize, nReqYSize,                                     pabyMEMDSBuffer, nReqXSize, nReqYSize,                                     eDataType, nBands, NULL,                                     0, 0, 0);            if (eErr != CE_None)            {                break;            }                        GDALDatasetH hMemDS = GDALCreate(hMemDriver, "MEM:::",                                              nReqXSize, nReqYSize, 0,                                               eDataType, NULL);            if (hMemDS == NULL)            {                eErr = CE_Failure;                break;            }                        int iBand;            for(iBand = 0; iBand < nBands; iBand ++)            {                char** papszMEMDSOptions = NULL;                char szTmp[64];                memset(szTmp, 0, sizeof(szTmp));                CPLPrintPointer(szTmp,                                pabyMEMDSBuffer + iBand * nDataTypeSize *                                nReqXSize * nReqYSize, sizeof(szTmp));                papszMEMDSOptions = CSLSetNameValue(papszMEMDSOptions, "DATAPOINTER", szTmp);                GDALAddBand(hMemDS, eDataType, papszMEMDSOptions);                CSLDestroy(papszMEMDSOptions);            }                        GDALDatasetH hOutDS = GDALCreateCopy(hTileDriver,                                        osTempFileName.c_str(), hMemDS, FALSE,                                        papszTileDriverOptions, NULL, NULL);            GDALClose(hMemDS);            if (hOutDS)                GDALClose(hOutDS);            else            {                eErr = CE_Failure;

//.........这里部分代码省略........./* -------------------------------------------------------------------- *//*      Open input datasource.                                          *//* -------------------------------------------------------------------- */    OGRDataSourceH hSrcDS;    hSrcDS = OGROpen( pszSource, FALSE, NULL );    if( hSrcDS == NULL )    {        fprintf( stderr, "Unable to open input datasource /"%s/"./n",                 pszSource );        fprintf( stderr, "%s/n", CPLGetLastErrorMsg() );        exit( 3 );    }/* -------------------------------------------------------------------- *//*      Create target raster file.                                      *//* -------------------------------------------------------------------- */    GDALDatasetH    hDstDS;    int             nLayerCount = CSLCount(papszLayers);    int             nBands = nLayerCount;    if ( pszSQL )        nBands++;    // FIXME    if ( nXSize == 0 )        nXSize = 256;    if ( nYSize == 0 )        nYSize = 256;    if (!bQuiet && !bFormatExplicitelySet)        CheckExtensionConsistency(pszDest, pszFormat);    hDstDS = GDALCreate( hDriver, pszDest, nXSize, nYSize, nBands,                         eOutputType, papszCreateOptions );    if ( hDstDS == NULL )    {        fprintf( stderr, "Unable to create target dataset /"%s/"./n",                 pszDest );        fprintf( stderr, "%s/n", CPLGetLastErrorMsg() );        exit( 3 );    }/* -------------------------------------------------------------------- *//*      If algorithm was not specified assigh default one.              *//* -------------------------------------------------------------------- */    if ( !pOptions )        ParseAlgorithmAndOptions( szAlgNameInvDist, &eAlgorithm, &pOptions );/* -------------------------------------------------------------------- *//*      Process SQL request.                                            *//* -------------------------------------------------------------------- */    if( pszSQL != NULL )    {        OGRLayerH hLayer;        hLayer = OGR_DS_ExecuteSQL( hSrcDS, pszSQL,                                    (OGRGeometryH)poSpatialFilter, NULL );         if( hLayer != NULL )        {            // Custom layer will be rasterized in the first band.            ProcessLayer( hLayer, hDstDS, poSpatialFilter, nXSize, nYSize, 1,                          bIsXExtentSet, bIsYExtentSet,                          dfXMin, dfXMax, dfYMin, dfYMax, pszBurnAttribute,                          dfIncreaseBurnValue, dfMultiplyBurnValue, eOutputType, eAlgorithm, pOptions,                          bQuiet, pfnProgress );

//.........这里部分代码省略.........            if( aoDrivers.empty() )            {                CPLError( CE_Failure, CPLE_AppDefined,                        "Cannot guess driver for %s", index_filename);                exit( 10 );            }            else            {                if( aoDrivers.size() > 1 )                {                    CPLError( CE_Warning, CPLE_AppDefined,                            "Several drivers matching %s extension. Using %s",                            CPLGetExtension(index_filename), aoDrivers[0].c_str() );                }                osFormat = aoDrivers[0];            }        }        else        {            osFormat = pszDriverName;        }        if( !EQUAL(osFormat, "ESRI Shapefile") )            nMaxFieldSize = 0;        GDALDriverH hDriver = GDALGetDriverByName( osFormat.c_str() );        if( hDriver == nullptr )        {            printf( "%s driver not available./n", osFormat.c_str() );            exit( 1 );        }        hTileIndexDS =             GDALCreate( hDriver, index_filename, 0, 0, 0, GDT_Unknown, nullptr );    }    if( hTileIndexDS != nullptr && hLayer == nullptr )    {        OGRSpatialReferenceH hSpatialRef = nullptr;        char* pszLayerName = nullptr;        if( pszIndexLayerName == nullptr )        {            VSIStatBuf sStat;            if( EQUAL(osFormat, "ESRI Shapefile") ||                VSIStat(index_filename, &sStat) == 0 )            {                pszLayerName = CPLStrdup(CPLGetBasename(index_filename));            }            else            {                printf( "-lyr_name must be specified./n" );                exit( 1 );            }        }        else        {            pszLayerName = CPLStrdup(pszIndexLayerName);        }        /* get spatial reference for output file from target SRS (if set) */        /* or from first input file */        if( bSetTargetSRS )        {            hSpatialRef = OSRClone( hTargetSRS );        }        else

//.........这里部分代码省略.........        for(nBlockXOff=0;eErr == CE_None && nBlockXOff<nXBlocks;nBlockXOff++)        {            GDALDatasetH hPrevOvrMemDS = NULL;/* -------------------------------------------------------------------- *//*      Create in-memory tile                                           *//* -------------------------------------------------------------------- */            int nReqXSize = nBlockXSize, nReqYSize = nBlockYSize;            if ((nBlockXOff+1) * nBlockXSize > nOvrXSize)                nReqXSize = nOvrXSize - nBlockXOff * nBlockXSize;            if ((nBlockYOff+1) * nBlockYSize > nOvrYSize)                nReqYSize = nOvrYSize - nBlockYOff * nBlockYSize;            if( pabyPrevOvrMEMDSBuffer != NULL )            {                int nPrevOvrReqXSize =                    (int)(nReqXSize * dfRatioPrevOvr + 0.5);                int nPrevOvrReqYSize =                    (int)(nReqYSize * dfRatioPrevOvr + 0.5);                eErr = RasterIO(GF_Read,                                nBlockXOff * nBlockXSize * nOvrFactor,                                nBlockYOff * nBlockYSize * nOvrFactor,                                nReqXSize * nOvrFactor, nReqYSize * nOvrFactor,                                pabyPrevOvrMEMDSBuffer, nPrevOvrReqXSize, nPrevOvrReqYSize,                                eDataType, nBands, NULL,                                0, 0, 0, NULL);                if (eErr != CE_None)                {                    break;                }                hPrevOvrMemDS = GDALCreate(hMemDriver, "MEM:::",                                           nPrevOvrReqXSize, nPrevOvrReqYSize, 0,                                           eDataType, NULL);                if (hPrevOvrMemDS == NULL)                {                    eErr = CE_Failure;                    break;                }                int iBand;                for(iBand = 0; iBand < nBands; iBand ++)                {                    char** papszOptions = NULL;                    char szTmp[64];                    memset(szTmp, 0, sizeof(szTmp));                    CPLPrintPointer(szTmp,                                    pabyPrevOvrMEMDSBuffer + iBand * nDataTypeSize *                                    nPrevOvrReqXSize * nPrevOvrReqYSize, sizeof(szTmp));                    papszOptions = CSLSetNameValue(papszOptions, "DATAPOINTER", szTmp);                    GDALAddBand(hPrevOvrMemDS, eDataType, papszOptions);                    CSLDestroy(papszOptions);                }            }            else            {                eErr = RasterIO(GF_Read,                                nBlockXOff * nBlockXSize * nOvrFactor,                                nBlockYOff * nBlockYSize * nOvrFactor,                                nReqXSize * nOvrFactor, nReqYSize * nOvrFactor,                                pabyMEMDSBuffer, nReqXSize, nReqYSize,                                eDataType, nBands, NULL,                                0, 0, 0, NULL);