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

int main(int argc, char* argv[]){    CPLJoinableThread* hThread;    printf("main thread %p/n", (void*)CPLGetPID());    argc = GDALGeneralCmdLineProcessor( argc, &argv, 0 );    CPLSetConfigOption("GDAL_CACHEMAX", "0");    CPLSetConfigOption("GDAL_DEBUG_BLOCK_CACHE", "ON");    MyDataset* poDS = new MyDataset();    char buf1[] = { 1 } ;    CPL_IGNORE_RET_VAL(GDALRasterIO(GDALGetRasterBand(poDS, 1), GF_Write, 0, 0, 1, 1, buf1, 1, 1, GDT_Byte, 0, 0));    hThread = CPLCreateJoinableThread(thread_func, NULL);    CPLSleep(0.3);    CPL_IGNORE_RET_VAL(GDALRasterIO(GDALGetRasterBand(poDS, 1), GF_Write, 1, 0, 1, 1, buf1, 1, 1, GDT_Byte, 0, 0));    GDALFlushCacheBlock();    CPLJoinThread(hThread);    delete poDS;    GDALDestroyDriverManager();    CSLDestroy( argv );    return 0;}

//辐射校正处理=====================================================================================================================================================================================//绝对辐射校正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;}

void CDialog3D::OnBnClickedBtnColor(){	// TODO: 在此添加控件通知处理程序代码	CFileDialog fileDlg(TRUE);	if(fileDlg.DoModal()!=IDOK)		return;	CString strExt	   =fileDlg.GetFileExt();	CString strPathName=fileDlg.GetPathName();	if(strExt=="BMP"||strExt=="bmp"||strExt=="JPG"||strExt=="jpg"||strExt=="TIF"||strExt=="tif"||strExt=="png"||strExt=="PNG")	{		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 *dataIn1=new float[tmpxsize*tmpysize];		float *dataIn2=new float[tmpxsize*tmpysize];		float *dataIn3=new float[tmpxsize*tmpysize];		GDALRasterIO(GDALGetRasterBand(hSrcDS,1),GF_Read,0,0,xsize,ysize,dataIn1,tmpxsize,tmpysize,GDT_Float32,0,0);		GDALRasterIO(GDALGetRasterBand(hSrcDS,2),GF_Read,0,0,xsize,ysize,dataIn2,tmpxsize,tmpysize,GDT_Float32,0,0);		GDALRasterIO(GDALGetRasterBand(hSrcDS,3),GF_Read,0,0,xsize,ysize,dataIn3,tmpxsize,tmpysize,GDT_Float32,0,0);				for(int i=0;i<tmpxsize;i++)		{			for (int j=0;j<tmpysize;j++)			{				PNT3D tmpPnt;				tmpPnt.DX=dataIn1[j*tmpxsize+i]/255;				tmpPnt.DY=dataIn2[j*tmpxsize+i]/255;				tmpPnt.DZ=dataIn3[j*tmpxsize+i]/255;				m_vec_colors.push_back(tmpPnt);			}		}		delete[]dataIn1;		delete[]dataIn2;		delete[]dataIn3;		GDALClose(hSrcDS);	}}

IplImage  * GetDataFromImg(const char *pszImg,bool bIsColor){	GDALAllRegister();	GDALDatasetH hDs = GDALOpen(pszImg,GA_ReadOnly);	if (hDs == NULL)		return NULL;	int nXSize = GDALGetRasterXSize(hDs);	int nYSize = GDALGetRasterYSize(hDs);	GDALRasterBandH hBand = GDALGetRasterBand(hDs,1);	unsigned char *pData = NULL;	if (!bIsColor)	{		pData = new unsigned char [ROW*COL];		GDALRasterIO(hBand, GF_Read, STAT_ROW, START_COL, ROW, COL, pData, ROW, COL, GDT_Byte, NULL, NULL);	}	else	{		pData = new unsigned char [ROW*COL*3];		int pBandMap[3]={1,2,3};		GDALDatasetRasterIO(hDs,GF_Read, STAT_ROW, START_COL, ROW, COL, pData, ROW, COL, GDT_Byte,3,pBandMap,3,3*COL,1);	}		GDALClose(hDs);	cv::Mat mat(ROW,COL,CV_8UC3,pData);	IplImage  limg =(IplImage)mat;	IplImage * img = (IplImage *)calloc(1,sizeof(IplImage) ) ;	memcpy(img,&limg,sizeof(IplImage));	return img;}

void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, QgsGeometry* poly, int pixelOffsetX,    int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, double& sum, double& count ){  double cellCenterX, cellCenterY;  QgsPoint currentCellCenter;  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );  cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;  count = 0;  sum = 0;  for ( int i = 0; i < nCellsY; ++i )  {    GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 );    cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;    for ( int j = 0; j < nCellsX; ++j )    {      currentCellCenter = QgsPoint( cellCenterX, cellCenterY );      if ( poly->contains( &currentCellCenter ) )      {        if ( scanLine[j] != mInputNodataValue ) //don't consider nodata values        {          sum += scanLine[j];          ++count;        }      }      cellCenterX += cellSizeX;    }    cellCenterY -= cellSizeY;  }  CPLFree( scanLine );}

/* Function to read the density data from the GDAL file into the array rho. * Returns 1 if there was a problem, 0 otherwise */int readpop(GDALRasterBandH hBand, double **rho, int xsize, int ysize){	int ix,iy;	double mean;	double sum;	double *rasterData, *rasterDataPtr, *rhoPtr;		// Read the raster band into rho.	rasterData = malloc(xsize*ysize*sizeof(double));	GDALRasterIO(hBand, GF_Read, 0, 0, xsize, ysize, rasterData, xsize, ysize, GDT_Float64, 0, 0);			// Compute the mean value.	sum = 0.0;	rasterDataPtr = rasterData;	for (iy = 0; iy < ysize; iy++) {		for (ix = 0; ix < xsize; ix++) {			rho[ix][iy] = *rasterDataPtr;			sum += *rasterDataPtr;			rasterDataPtr++;		}	}	mean = sum / ((double)xsize*(double)ysize);		// Add an bias to all raster values.	rhoPtr = *rho;	for (iy = 0; iy < ysize; iy++) {		for (ix = 0; ix < xsize; ix++) {			*rhoPtr += OFFSET*mean;			rhoPtr++;		}	}		return 0;}

float* get_tif_data(char* tif_file, int* tif_width, int* tif_height) {	GDALDatasetH hDataset;	hDataset = GDALOpen(tif_file, GA_ReadOnly);	if (hDataset == NULL ) {		fprintf(stderr, "ERROR: Failed to open the file %s/n", tif_file);		return NULL ;	}	GDALRasterBandH hBand;	hBand = GDALGetRasterBand(hDataset, 1);	int width1 = GDALGetRasterXSize(hDataset);	int height1 = GDALGetRasterYSize(hDataset);	double nodata1 = GDALGetRasterNoDataValue(hBand, NULL );	float* data = (float *) CPLMalloc(sizeof(float) * width1 * height1);	if (!data) {		fprintf(stderr, "ERROR: Failed to allocate data of size %d /n", width1 * height1);		return NULL;	}	GDALRasterIO(hBand, GF_Read, 0, 0, width1, height1, data, width1, height1,			GDT_Float32, 0, 0);	*tif_width = width1;	*tif_height = height1;	return data;}

void ColorizationFilter::filter(PointView& view){    int32_t pixel(0);    int32_t line(0);    std::array<double, 2> pix = { {0.0, 0.0} };    for (PointId idx = 0; idx < view.size(); ++idx)    {        double x = view.getFieldAs<double>(Dimension::Id::X, idx);        double y = view.getFieldAs<double>(Dimension::Id::Y, idx);        if (!getPixelAndLinePosition(x, y, m_inverse_transform, pixel,                line, m_ds))            continue;        for (auto bi = m_bands.begin(); bi != m_bands.end(); ++bi)        {            BandInfo& b = *bi;            GDALRasterBandH hBand = GDALGetRasterBand(m_ds, b.m_band);            if (hBand == NULL)            {                std::ostringstream oss;                oss << "Unable to get band " << b.m_band <<                    " from data source!";                throw pdal_error(oss.str());            }            if (GDALRasterIO(hBand, GF_Read, pixel, line, 1, 1,                &pix[0], 1, 1, GDT_CFloat64, 0, 0) == CE_None)                view.setField(b.m_dim, idx, pix[0] * b.m_scale);        }    }}

void MDAL::DriverGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<MemoryDataset> tos, bool is_vector, bool is_x ){  assert( raster_band );  double nodata =  GDALGetRasterNoDataValue( raster_band, nullptr );  unsigned int mXSize = meshGDALDataset()->mXSize;  unsigned int mYSize = meshGDALDataset()->mYSize;  double *values = tos->values();  for ( unsigned int y = 0; y < mYSize; ++y )  {    // buffering per-line    CPLErr err = GDALRasterIO(                   raster_band,                   GF_Read,                   0, //nXOff                   static_cast<int>( y ), //nYOff                   static_cast<int>( mXSize ), //nXSize                   1, //nYSize                   mPafScanline, //pData                   static_cast<int>( mXSize ), //nBufXSize                   1, //nBufYSize                   GDT_Float64, //eBufType                   0, //nPixelSpace                   0 //nLineSpace                 );    if ( err != CE_None )    {      throw MDAL_Status::Err_InvalidData;    }    for ( unsigned int x = 0; x < mXSize; ++x )    {      unsigned int idx = x + mXSize * y;      double val = mPafScanline[x];      if ( !MDAL::equals( val, nodata ) )      {        // values is prepolulated with NODATA values, so store only legal values        if ( is_vector )        {          if ( is_x )          {            values[2 * idx] = val;          }          else          {            values[2 * idx + 1] = val;          }        }        else        {          values[idx] = val;        }      }    }  }}

void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeometry& poly, int pixelOffsetX,    int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats ){  double cellCenterX, cellCenterY;  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );  cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;  stats.reset();  GEOSGeometry* polyGeos = poly.exportToGeos();  if ( !polyGeos )  {    return;  }  GEOSContextHandle_t geosctxt = QgsGeometry::getGEOSHandler();  const GEOSPreparedGeometry* polyGeosPrepared = GEOSPrepare_r( geosctxt, polyGeos );  if ( !polyGeosPrepared )  {    GEOSGeom_destroy_r( geosctxt, polyGeos );    return;  }  GEOSCoordSequence* cellCenterCoords = nullptr;  GEOSGeometry* currentCellCenter = nullptr;  for ( int i = 0; i < nCellsY; ++i )  {    if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )         != CPLE_None )    {      continue;    }    cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;    for ( int j = 0; j < nCellsX; ++j )    {      if ( validPixel( scanLine[j] ) )      {        GEOSGeom_destroy_r( geosctxt, currentCellCenter );        cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );        GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX );        GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY );        currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );        if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )        {          stats.addValue( scanLine[j] );        }      }      cellCenterX += cellSizeX;    }    cellCenterY -= cellSizeY;  }  GEOSGeom_destroy_r( geosctxt, currentCellCenter );  CPLFree( scanLine );  GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );  GEOSGeom_destroy_r( geosctxt, polyGeos );}

void QgsZonalStatistics::statisticsFromPreciseIntersection( void* band, const QgsGeometry* poly, int pixelOffsetX,    int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats ){  stats.reset();  double currentY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;  float* pixelData = ( float * ) CPLMalloc( sizeof( float ) );  QgsGeometry* pixelRectGeometry = nullptr;  double hCellSizeX = cellSizeX / 2.0;  double hCellSizeY = cellSizeY / 2.0;  double pixelArea = cellSizeX * cellSizeY;  double weight = 0;  for ( int row = 0; row < nCellsY; ++row )  {    double currentX = rasterBBox.xMinimum() + cellSizeX / 2.0 + pixelOffsetX * cellSizeX;    for ( int col = 0; col < nCellsX; ++col )    {      if ( GDALRasterIO( band, GF_Read, pixelOffsetX + col, pixelOffsetY + row, nCellsX, 1, pixelData, 1, 1, GDT_Float32, 0, 0 ) != CE_None )        QgsDebugMsg( "Raster IO Error" );      if ( !validPixel( *pixelData ) )        continue;      pixelRectGeometry = QgsGeometry::fromRect( QgsRectangle( currentX - hCellSizeX, currentY - hCellSizeY, currentX + hCellSizeX, currentY + hCellSizeY ) );      if ( pixelRectGeometry )      {        //intersection        QgsGeometry *intersectGeometry = pixelRectGeometry->intersection( poly );        if ( intersectGeometry )        {          double intersectionArea = intersectGeometry->area();          if ( intersectionArea >= 0.0 )          {            weight = intersectionArea / pixelArea;            stats.addValue( *pixelData, weight );          }          delete intersectGeometry;        }        delete pixelRectGeometry;        pixelRectGeometry = nullptr;      }      currentX += cellSizeX;    }    currentY -= cellSizeY;  }  CPLFree( pixelData );}

boost::uint32_t Colorization::readBufferImpl(PointBuffer& data){    const boost::uint32_t numRead = getPrevIterator().read(data);#ifdef PDAL_HAVE_GDAL    boost::int32_t pixel(0);    boost::int32_t line(0);    double x(0.0);    double y(0.0);    bool fetched(false);    boost::array<double, 2> pix;    pix.assign(0.0);    for (boost::uint32_t pointIndex=0; pointIndex<numRead; pointIndex++)    {        x = getScaledValue(data, *m_dimX, pointIndex);        y = getScaledValue(data, *m_dimY, pointIndex);        fetched = getPixelAndLinePosition(x, y, m_inverse_transform, pixel, line, m_ds);        if (!fetched)            continue;        for (std::vector<boost::int32_t>::size_type i = 0;                i < m_bands.size(); i++)        {            GDALRasterBandH hBand = GDALGetRasterBand(m_ds, m_bands[i]);            if (hBand == NULL)            {                std::ostringstream oss;                oss << "Unable to get band " << m_bands[i] << " from data source!";                throw pdal_error(oss.str());            }            if (GDALRasterIO(hBand, GF_Read, pixel, line, 1, 1,                             &pix[0], 1, 1, GDT_CFloat64, 0, 0) == CE_None)            {                double output = pix[0];                output = output * m_scales[i];                setScaledValue(data, output, *m_dimensions[i], pointIndex);            }        }    }#endif    return numRead;}

double QgsAlignRaster::RasterInfo::identify( double mx, double my ){  GDALRasterBandH hBand = GDALGetRasterBand( mDataset, 1 );  // must not be rotated in order for this to work  int px = int(( mx - mGeoTransform[0] ) / mGeoTransform[1] );  int py = int(( my - mGeoTransform[3] ) / mGeoTransform[5] );  float* pafScanline = ( float * ) CPLMalloc( sizeof( float ) );  CPLErr err = GDALRasterIO( hBand, GF_Read, px, py, 1, 1,                             pafScanline, 1, 1, GDT_Float32, 0, 0 );  double value = err == CE_None ? pafScanline[0] : NAN;  CPLFree( pafScanline );  return value;}

static intvips__gdal_generate( VipsRegion *out, 	void *_seq, void *_read, void *unused, gboolean *stop ){	Read *read = _read;	VipsRect *r = &out->valid;	int n = r->width * r->height;	unsigned char *buf = (unsigned char *) VIPS_REGION_ADDR( out, r->left, r->top );	VIPS_DEBUG_MSG( "vips__gdal_generate: %dx%d @ %dx%d/n",		r->width, r->height, r->left, r->top );	/* We're inside a cache, so requests should always be	 * tile_width by tile_height pixels and on a tile boundary.	 */	g_assert( (r->left % read->tile_width) == 0 );	g_assert( (r->top % read->tile_height) == 0 );	g_assert( r->width <= read->tile_width );	g_assert( r->height <= read->tile_height );	g_assert( read->data_type == GDT_Byte );	g_assert( GDALGetRasterCount );  for (int channel = 0; channel < 3; ++channel) {		GDALRasterBandH hBand;		hBand = GDALGetRasterBand( read->hDataset, channel + 1 );		unsigned char* gdal_data;		gdal_data = (unsigned char*) CPLMalloc(sizeof(unsigned char) * r->width * r->height);		GDALRasterIO( hBand, GF_Read, r->left, r->top, r->width, r->height, 		              gdal_data, r->width, r->height, GDT_Byte, 		              0, 0 );		for (int i = 0; i < r->width * r->height; ++i) {			buf[i * 3 + channel] = gdal_data[i];		}	  CPLFree (gdal_data);	}	return( 0 );}

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

// API: get a sample from an image, at the requested octavefloat fancy_image_getsample_oct(struct fancy_image *fi,		int octave, int i,int j, int l){	struct FI *f = (void*)fi;	if (octave < 0 || octave >= f->no)		return NAN;	if (l < 0) l = 0;	if (l >= f->pd) l = f->pd - 1;	if (f->tiffo) {#ifdef FANCY_TIFF		uint8_t *p_pixel = tiff_octaves_getpixel(f->t, octave, i, j);		if (!p_pixel) return NAN;		uint8_t *p_sample = p_pixel + (l * f->t->i->bps) / 8;		return convert_sample_to_float(f->t->i, p_sample);#else		assert(false);#endif	} else if (f->gdal) {#ifdef FANCY_GDAL		if (octave != 0) return NAN;		static float *roi = NULL;		if (!roi) roi = CPLMalloc(1*1*sizeof*roi);		GDALRasterBandH img = f->gdal_band[l];		int r = GDALRasterIO(img, GF_Read, i,j,1, 1, roi,1,1,				GDT_Float32, 0,0);		return roi[0*0+0];#else		assert(false);#endif	} else {		float *x = f->pyr_x[octave];		int    w = f->pyr_w[octave];		int    h = f->pyr_h[octave];		if (i < 0 || j < 0 || i >= w || j >= h)			return NAN;		int  idx = (j * w + i) * f->pd + l;		return x[idx];	}	return NAN;}

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

CPLErr QgsGdalProviderBase::gdalRasterIO( GDALRasterBandH hBand, GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize, void * pData, int nBufXSize, int nBufYSize, GDALDataType eBufType, int nPixelSpace, int nLineSpace ){  // See http://hub.qgis.org/issues/8356 and http://trac.osgeo.org/gdal/ticket/5170#if GDAL_VERSION_MAJOR == 1 && ( (GDAL_VERSION_MINOR == 9 && GDAL_VERSION_REV <= 2) || (GDAL_VERSION_MINOR == 10 && GDAL_VERSION_REV <= 0) )  char* pszOldVal = CPLStrdup( CPLGetConfigOption( "VSI_CACHE", "FALSE" ) );  CPLSetThreadLocalConfigOption( "VSI_CACHE", "FALSE" );  QgsDebugMsg( "Disabled VSI_CACHE" );#endif  CPLErr err = GDALRasterIO( hBand, eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize, eBufType, nPixelSpace, nLineSpace );#if GDAL_VERSION_MAJOR == 1 && ( (GDAL_VERSION_MINOR == 9 && GDAL_VERSION_REV <= 2) || (GDAL_VERSION_MINOR == 10 && GDAL_VERSION_REV <= 0) )  CPLSetThreadLocalConfigOption( "VSI_CACHE", pszOldVal );  CPLFree( pszOldVal );  QgsDebugMsg( "Reset VSI_CACHE" );#endif  return err;}

static CPLErr GPMaskImageData( GDALRasterBandH hMaskBand, GByte* pabyMaskLine, int iY, int nXSize,                  float *pafImageLine ){    CPLErr eErr;    eErr = GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,                          pabyMaskLine, nXSize, 1, GDT_Byte, 0, 0 );    if( eErr == CE_None )    {        int i;        for( i = 0; i < nXSize; i++ )        {            if( pabyMaskLine[i] == 0 )                pafImageLine[i] = GP_NODATA_MARKER;        }    }    return eErr;}

int GDALFillBandNoData(GDALDataset *poDS, int nBand, int nSearchPixels){    (void)nBand;    (void)nSearchPixels;    if(poDS == NULL)    {        fprintf(stderr, "Invalid GDAL Dataset Handle, cannot fill no data/n");        return -1;    }    int nPixels, nLines;    nPixels = poDS->GetRasterXSize();    nLines = poDS->GetRasterYSize();    GDALRasterBand *poBand;    poBand = poDS->GetRasterBand(1);    GDALFillNodata(poBand, NULL, 100, 0, 0, NULL, NULL, NULL);    double dfNoData = poBand->GetNoDataValue(NULL);    double *padfScanline;    padfScanline = (double *) CPLMalloc(sizeof(double)*nPixels);    int nNoDataCount = 0;    for(int i = 0;i < nLines;i++)    {        GDALRasterIO(poBand, GF_Read, 0, i, nPixels, 1,                     padfScanline, nPixels, 1, GDT_Float64, 0, 0);        for(int j = 0; j < nPixels;j++)        {            if(CPLIsEqual(padfScanline[j], dfNoData))                nNoDataCount++;        }    }    CPLFree(padfScanline);    return nNoDataCount;}

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

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

void MDAL::LoaderGdal::addDataToOutput( GDALRasterBandH raster_band, std::shared_ptr<Dataset> tos, bool is_vector, bool is_x ){  assert( raster_band );  double nodata =  GDALGetRasterNoDataValue( raster_band, nullptr );  unsigned int mXSize = meshGDALDataset()->mXSize;  unsigned int mYSize = meshGDALDataset()->mYSize;  for ( unsigned int y = 0; y < mYSize; ++y )  {    // buffering per-line    CPLErr err = GDALRasterIO(                   raster_band,                   GF_Read,                   0, //nXOff                   static_cast<int>( y ), //nYOff                   static_cast<int>( mXSize ), //nXSize                   1, //nYSize                   mPafScanline, //pData                   static_cast<int>( mXSize ), //nBufXSize                   1, //nBufYSize                   GDT_Float64, //eBufType                   0, //nPixelSpace                   0 //nLineSpace                 );    if ( err != CE_None )    {      throw MDAL_Status::Err_InvalidData;    }    for ( unsigned int x = 0; x < mXSize; ++x )    {      unsigned int idx = x + mXSize * y;      double val = mPafScanline[x];      bool noData = false;      if ( MDAL::equals( val, nodata ) )      {        // store all nodata value as this hardcoded number        val = MDAL_NODATA;        noData = true;      }      if ( is_vector )      {        if ( is_x )        {          tos->values[idx].x = val;          tos->values[idx].noData = noData;        }        else        {          tos->values[idx].y = val;          tos->values[idx].noData = noData;        }      }      else      {        tos->values[idx].x = val;        tos->values[idx].noData = noData;      }    }  }}

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

CPLErr CPL_STDCALLGDALFPolygonize( GDALRasterBandH hSrcBand,                GDALRasterBandH hMaskBand,                OGRLayerH hOutLayer, int iPixValField,                 char **papszOptions,                GDALProgressFunc pfnProgress,                 void * pProgressArg ){#ifndef OGR_ENABLED    CPLError(CE_Failure, CPLE_NotSupported, "GDALFPolygonize() unimplemented in a non OGR build");    return CE_Failure;#else    VALIDATE_POINTER1( hSrcBand, "GDALFPolygonize", CE_Failure );    VALIDATE_POINTER1( hOutLayer, "GDALFPolygonize", CE_Failure );    if( pfnProgress == NULL )        pfnProgress = GDALDummyProgress;    int nConnectedness = CSLFetchNameValue( papszOptions, "8CONNECTED" ) ? 8 : 4;/* -------------------------------------------------------------------- *//*      Confirm our output layer will support feature creation.         *//* -------------------------------------------------------------------- */    if( !OGR_L_TestCapability( hOutLayer, OLCSequentialWrite ) )    {        CPLError( CE_Failure, CPLE_AppDefined,                  "Output feature layer does not appear to support creation/n"                  "of features in GDALFPolygonize()." );        return CE_Failure;    }/* -------------------------------------------------------------------- *//*      Allocate working buffers.                                       *//* -------------------------------------------------------------------- */    CPLErr eErr = CE_None;    int nXSize = GDALGetRasterBandXSize( hSrcBand );    int nYSize = GDALGetRasterBandYSize( hSrcBand );    float *pafLastLineVal = (float *) VSIMalloc2(sizeof(float),nXSize + 2);    float *pafThisLineVal = (float *) VSIMalloc2(sizeof(float),nXSize + 2);    GInt32 *panLastLineId =  (GInt32 *) VSIMalloc2(sizeof(GInt32),nXSize + 2);    GInt32 *panThisLineId =  (GInt32 *) VSIMalloc2(sizeof(GInt32),nXSize + 2);    GByte *pabyMaskLine = (hMaskBand != NULL) ? (GByte *) VSIMalloc(nXSize) : NULL;    if (pafLastLineVal == NULL || pafThisLineVal == NULL ||        panLastLineId == NULL || panThisLineId == NULL ||        (hMaskBand != NULL && pabyMaskLine == NULL))    {        CPLError(CE_Failure, CPLE_OutOfMemory,                 "Could not allocate enough memory for temporary buffers");        CPLFree( panThisLineId );        CPLFree( panLastLineId );        CPLFree( pafThisLineVal );        CPLFree( pafLastLineVal );        CPLFree( pabyMaskLine );        return CE_Failure;    }/* -------------------------------------------------------------------- *//*      Get the geotransform, if there is one, so we can convert the    *//*      vectors into georeferenced coordinates.                         *//* -------------------------------------------------------------------- */    GDALDatasetH hSrcDS = GDALGetBandDataset( hSrcBand );    double adfGeoTransform[6] = { 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };    if( hSrcDS )        GDALGetGeoTransform( hSrcDS, adfGeoTransform );/* -------------------------------------------------------------------- *//*      The first pass over the raster is only used to build up the     *//*      polygon id map so we will know in advance what polygons are     *//*      what on the second pass.                                        *//* -------------------------------------------------------------------- */    int iY;    GDALRasterFPolygonEnumerator oFirstEnum(nConnectedness);    for( iY = 0; eErr == CE_None && iY < nYSize; iY++ )    {        eErr = GDALRasterIO(             hSrcBand,            GF_Read, 0, iY, nXSize, 1,             pafThisLineVal, nXSize, 1, GDT_Float32, 0, 0 );                if( eErr == CE_None && hMaskBand != NULL )            eErr = GPMaskImageData( hMaskBand, pabyMaskLine, iY, nXSize,                    pafThisLineVal );        if( iY == 0 )            oFirstEnum.ProcessLine(                 NULL, pafThisLineVal, NULL, panThisLineId, nXSize );        else            oFirstEnum.ProcessLine(                pafLastLineVal, pafThisLineVal,                panLastLineId,  panThisLineId,                 nXSize );        // swap lines        float * pafTmp = pafLastLineVal;        pafLastLineVal = pafThisLineVal;        pafThisLineVal = pafTmp;//.........这里部分代码省略.........

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