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

int
makeGeotiff (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;
}