C++ GDALDriver::GetMetadata方法代码示例

void Raster<T>::OutputGTiff(const char* rasterName)
{
	const char *pszFormat = "GTiff";
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

	char **papszOptions = poDriver->GetMetadata();
	//papszOptions = CSLSetNameValue( papszOptions, "TILED", "YES" );
	//papszOptions = CSLSetNameValue( papszOptions, "COMPRESS", "PACKBITS" );
	GDALDataset *poDstDS = poDriver->Create(rasterName, m_nCols, m_nRows, 1, m_dType, papszOptions);
	
	//write the data to new file
	GDALRasterBand  *poDstBand= poDstDS->GetRasterBand(1);
	poDstBand->RasterIO(GF_Write, 0, 0,  m_nCols, m_nRows, m_data,  m_nCols, m_nRows, m_dType, 0, 0);
	poDstBand->SetNoDataValue(m_noDataValue);

	double geoTrans[6];
	geoTrans[0] = m_xMin;
	geoTrans[1] = m_dx;
	geoTrans[2] = 0;
	geoTrans[3] = m_yMax;
	geoTrans[4] = 0;
	geoTrans[5] = -m_dy;
	poDstDS->SetGeoTransform(geoTrans);

	poDstDS->SetProjection(m_proj.c_str());

	GDALClose(poDstDS);
}

void write_map(fs::path file_path, GDALDataType data_type, boost::shared_ptr<Map_Matrix<DataFormat> > data, std::string WKTprojection, GeoTransform transform, std::string driverName) throw(std::runtime_error)
{
	GDALAllRegister(); //This registers all availble raster file formats for use with this utility. How neat is that. We can input any GDAL supported rater file format.
    
    const char *pszFormat = driverName.c_str();
    GDALDriver * poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if (poDriver == NULL)
    {
        throw std::runtime_error("No driver for file tyle found");
    }
    
    char ** papszMetadata = poDriver->GetMetadata();
    if (!(CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE)))
    {
        throw std::runtime_error("Driver does not support raster creation");
    }
    
    char **papszOptions = NULL;
	papszOptions = CSLSetNameValue(papszOptions, "COMPRESS", "LZW");

    GDALDataset *poDstDS = poDriver->Create(file_path.string().c_str(), (int)data->NCols(), (int)data->NRows(), 1, data_type, papszOptions);
    
    double adfGeoTransform[6] = {1, 1, 1, 1, 1, 1};
    adfGeoTransform[0] = transform.x_origin;
    adfGeoTransform[1] = transform.pixel_width;
    adfGeoTransform[2] = transform.x_line_space;
    adfGeoTransform[3] = transform.y_origin;
    adfGeoTransform[4] = transform.pixel_height;
    adfGeoTransform[5] = transform.y_line_space;
    
    const char * psz_WKT = WKTprojection.c_str();
    poDstDS->SetGeoTransform(adfGeoTransform);             
    poDstDS->SetProjection(psz_WKT);
    
    DataFormat * pafScanline = new DataFormat[data->NCols() * data->NRows()];
    int pafIterator = 0;
	for (int i = 0; i < data->NRows(); i++)
    {
		for (int j = 0; j < data->NCols(); j++)
        {
            pafScanline[pafIterator] = data->Get(i, j);
            pafIterator++;
        }
    }
    
    GDALRasterBand * poBand = poDstDS->GetRasterBand(1);
    poBand->SetNoDataValue(data->NoDataValue());
    poBand->RasterIO(GF_Write, 0, 0, (int) data->NCols(), (int) data->NRows(), pafScanline, (int) data->NCols(), (int) data->NRows(), data_type, 0, 0);
    
    GDALClose( (GDALDatasetH) poDstDS);
}

int Raster<T>::OutputGeoTiff(const char* filename)
{
    /// Output GeoTiff all set as float datatype. by LJ.
    int n = m_nRows * m_nCols;
    float *data = new float[n];
    int index = 0;
    for (int i = 0; i < m_nRows; ++i)
    {
        for (int j = 0; j < m_nCols; ++j)
        {
            data[i*m_nCols+j] = float(m_data[i][j]);
        }
    }
    const char *pszFormat = "GTiff";
    GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    char **papszOptions = poDriver->GetMetadata();
    GDALDataset *poDstDS = poDriver->Create(filename, m_nCols, m_nRows, 1, GDT_Float32, papszOptions );

    /// Write the data to new file
    GDALRasterBand  *poDstBand= poDstDS->GetRasterBand(1);
    poDstBand->RasterIO(GF_Write, 0, 0,  m_nCols, m_nRows, data,  m_nCols, m_nRows, GDT_Float32, 0, 0);
    poDstBand->SetNoDataValue(m_noDataValue);

    double geoTrans[6];
    geoTrans[0] = m_xllCenter;
    geoTrans[1] = m_dx;
    geoTrans[2] = 0;
    geoTrans[3] = m_yllCenter + m_nRows*m_dx;
    geoTrans[4] = 0;
    geoTrans[5] = -m_dx;
    poDstDS->SetGeoTransform(geoTrans);
    poDstDS->SetProjection(m_srs.c_str());
    GDALClose(poDstDS);
    delete[] data;
    return 0;
}

void clsRasterData::outputGTiff(map<string,float> header,int nValidCells,float** position, float* value,string filename)
{
	float noDataValue = header["NODATA_VALUE"];

	int nCols = header["NCOLS"];
	int nRows = header["NROWS"];
	float xll = header["XLLCENTER"];
	float yll = header["YLLCENTER"];
	float dx = header["CELLSIZE"];

	int n = nRows * nCols;
	float *data = new float[n];

	int index = 0;
	for (int i = 0; i < nRows; ++i)
	{
		for (int j = 0; j < nCols; ++j)
		{
			if(index < nValidCells)
			{
				if(position[index][0] == i && position[index][1] == j) 
				{
					data[i*nCols+j] = value[index];	
					index++;
				}
				else 
					data[i*nCols+j] = noDataValue;
			}
			else 
				data[i*nCols+j] = noDataValue;				
		}
	}

	const char *pszFormat = "GTiff";
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

	char **papszOptions = poDriver->GetMetadata();
	GDALDataset *poDstDS = poDriver->Create(filename.c_str(), nCols, nRows, 1, GDT_Float32, papszOptions );
	
	//write the data to new file
	GDALRasterBand  *poDstBand= poDstDS->GetRasterBand(1);
	poDstBand->RasterIO(GF_Write, 0, 0,  nCols, nRows, data,  nCols, nRows, GDT_Float32, 0, 0);
	poDstBand->SetNoDataValue(noDataValue);

	double geoTrans[6];
	geoTrans[0] = xll;
	geoTrans[1] = dx;
	geoTrans[2] = 0;
	geoTrans[3] = yll + nRows*dx;
	geoTrans[4] = 0;
	geoTrans[5] = -dx;
	poDstDS->SetGeoTransform(geoTrans);

	OGRSpatialReference srs;
	srs.SetACEA(25, 47, 0, 105, 0, 0);
	srs.SetWellKnownGeogCS("WGS84");
	
	char *pSrsWkt = NULL;
	srs.exportToWkt(&pSrsWkt);
	poDstDS->SetProjection(pSrsWkt);
	CPLFree(pSrsWkt);

	GDALClose(poDstDS);

	delete[] data;
}

int meaningful_change(string cd_map_file, string result_file)
{
	GdalRasterApp cd_map;
	cd_map.open(cd_map_file.c_str());
	int iBandCount = cd_map.nBand();
	int iTileCountX = cd_map.getTileCountX();
	int iTileCountY = cd_map.getTileCountY();
	int iWidth = cd_map.width();
	int iHeight = cd_map.height();

	GDALDriver *poDriver;	//驱动，用于创建新的文件
	poDriver=GetGDALDriverManager()->GetDriverByName("GTIFF");
	char **papszMetadata = poDriver->GetMetadata();//获取格式类型
	GDALDataset *poDatasetNew;
	// 输出栅格
	poDatasetNew = poDriver->Create(result_file.c_str(), iWidth, iHeight, 1, GDT_Byte, papszMetadata);//根据文件路径文件名，图像宽，高，波段数，数据类型，文件类型，创建新的数据集				
	poDatasetNew->SetProjection(cd_map.getGetProjectionRef());
	poDatasetNew->SetGeoTransform(cd_map.getGeoTransform());//坐标赋值,与全色相同

	int nCount = 0;
	double norm_threshold = 20.0;
	int bgColor = 0;
	int fgColor = 255;
	for (int i = 0;i < iTileCountX;++i)
	{
		for (int j = 0;j < iTileCountY;++j)
		{
			GdalRasterApp::RasterBuf *pBuf = cd_map.getTileData(i, j, iBandCount);
			int bufWidth = pBuf->iBufWidth;
			int bufHeight = pBuf->iBufHeight;
			int bufBand = pBuf->iBandCount;
			int offsetX, offsetY;
			cd_map.getTileOffset(i, j, offsetX, offsetY);
			cv::Mat change_image(cv::Size(bufWidth, bufHeight), CV_8UC1, cv::Scalar(fgColor));

			//cv::Mat lbp_change(cv::Size(bufWidth, bufHeight), CV_8UC1);
			cv::Mat lbp_change = img_int2byte_band(pBuf, 0);
			//cv::Mat lbp_change = img_float2byte_band(pBuf, 0);
			//memcpy(lbp_change.data, pBuf->data, bufWidth*bufHeight);
			//cv::imwrite("E:\\minus.tif", lbp_change);
			int lbp_change_threshold = cvThresholdOtsu(lbp_change);

			int step = change_image.step;
			// change detection
			for (int row = 0; row < bufHeight; ++row) {
				for (int col = 0; col < bufWidth; ++col) {
					int lbp_change_ = lbp_change.data[row*lbp_change.step+col];
					// 判断是否变化很小
					if (lbp_change_ < lbp_change_threshold)
					{
						change_image.data[row*step+col] = bgColor;
						continue;
					}

					std::vector<double> change_vector(pFeatureBuf->iBandCount);
					for (int k = 0;k < pFeatureBuf->iBandCount;++k)
					{
						int pos = row * bufWidth * pFeatureBuf->iBandCount + col*pFeatureBuf->iBandCount + k;
						change_vector[k] = ((float*)pFeatureBuf->data)[pos];
					}

					//如果变化是否和样本相似
					for (int m = 0; m < (int)samples.size(); m++)
					{
						double angle_ = samples[m][0] - change_vector[0];
						double norm_ = samples[m][1] - change_vector[1];
						//double similarity = VectorSimilarity(samples[m], change_vector);
						//double distance = VectorDistance(samples[m], change_vector);
						//double angle = VectorAngle(samples[m], change_vector);
						//if (fabs(similarity) > similarityThreshold && distance < 15.0)
						//if (fabs(angle) < 10.0 && distance < 20.0)
						//if (fabs(similarity) < similarityThreshold)
						if (fabs(angle_) < similarityThreshold && fabs(norm_) < 10.0)
							//if (fabs(norm) < 10.0)
						{
							change_image.data[row*step+col] = bgColor;
							break;
						}
					}
				}
			}

			//RasterBuf2Opencv(pBuf, change_image);
			//cv::Mat change_image(bufHeight, bufWidth, CV_8U(bufBand));
			//int nBandDataSize = GDALGetDataTypeSize( pBuf->eDataType ) / 8;
			//memcpy(change_image.data, (GByte*)(pBuf->data), bufWidth*bufHeight*bufBand*nBandDataSize);
			//cvtColor( change_image, change_image, CV_BGR2GRAY );
			//int threshold = cvThresholdOtsu(change_image);
			//threshold = 60;
			//cv::threshold( change_image, change_image, threshold, 255, CV_THRESH_BINARY);
			//cv::Mat element = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(3, 3));
			//cv::morphologyEx(change_image, change_image, cv::MORPH_CLOSE, element);
			//element = cv::getStructuringElement(cv::MORPH_CROSS , cv::Size(7, 7));
			//cv::morphologyEx(change_image, change_image, cv::MORPH_OPEN, element);
			//element = cv::getStructuringElement(cv::MORPH_ELLIPSE, cv::Size(7, 7));
			//cv::morphologyEx(change_image, change_image, cv::MORPH_OPEN, element);
			//cv::GaussianBlur(change_image, change_image, cv::Size(5,5), 1.5);
			//cv::morphologyEx(change_image, change_image, cv::MORPH_CLOSE, element);

			//cv::imwrite("E:\\test0.tif", change_image);
//.........这里部分代码省略.........

static OGRErr CreateAndFillOutputDataset(OGRLayer* poSrcLayer,
                                         const char* pszDestDataSource,
                                         const char* pszFormat,
                                         const char* pszLayer,
                                         char** papszDSCO,
                                         char** papszLCO,
                                         int bQuiet)
{
    GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
    if( poDriver == NULL )
    {
        fprintf( stderr, "%s driver not available\n", pszFormat );
        return OGRERR_FAILURE;
    }

    if( !CSLTestBoolean(
                CSLFetchNameValueDef(poDriver->GetMetadata(), GDAL_DCAP_CREATE,
                                     "FALSE") ) )
    {
        fprintf( stderr,  "%s driver does not support data source creation.\n",
                pszFormat );
        return OGRERR_FAILURE;
    }

    GDALDataset* poODS = poDriver->Create( pszDestDataSource, 0, 0, 0,
                                         GDT_Unknown, papszDSCO );
    if( poODS == NULL )
    {
        fprintf( stderr,  "%s driver failed to create %s\n",
                pszFormat, pszDestDataSource );
        return OGRERR_FAILURE;
    }

    if(NULL == pszLayer)
        pszLayer = poSrcLayer->GetName();
    int nError;
    GetLayerAndOverwriteIfNecessary(poODS, pszLayer, TRUE, &nError);
    if(nError == TRUE)
    {
        return OGRERR_FAILURE;
    }

    // create layer
    OGRLayer * poLayer = poODS->CopyLayer(poSrcLayer, pszLayer, papszLCO);
    if (NULL == poLayer)
    {
        fprintf(stderr, "\nFAILURE: Can not copy path to %s\n",
                pszDestDataSource);
        GDALClose(poODS);

        return OGRERR_FAILURE;
    }

    if (bQuiet == FALSE)
    {
        printf("\nPath successfully copied and added to the network at %s\n",
            pszDestDataSource);
    }

    GDALClose(poODS);

    return OGRERR_NONE;
}

bool GDALUtilities::createRasterFile(QString& theFilename, 
									 QString& theFormat,  
									 GDALDataType theType, 
									 int theBands, 
									 int theRows, 
									 int theCols, 
									 void ** theData,
									 double * theGeoTransform, 
									 const QgsCoordinateReferenceSystem * theCrs,
									 double theNodataValue)
{
	if ( theBands <= 0 )
		return false;

	if ( theRows <= 0 )
		return false;

	if ( theCols <= 0 )
		return false;

	if ( !theData )
		return false;

/*	bool formatSupported = false;
	QMapIterator<QString, QString> i(mSupportedFormats);
	while (i.hasNext()) 
	{
		i.next();
		if( theFormat == i.key())
		{
			formatSupported = true;
			break;
		}
	}

	if ( !formatSupported )
		return false;

*/

	//GDALAllRegister();

	GDALDriver * driver;

	//set format
	char * format = new char[theFormat.size() + 1];
	strcpy( format, theFormat.toLocal8Bit().data() );

    driver = GetGDALDriverManager()->GetDriverByName(format);

	if( driver == NULL )
        return false;

	char ** metadata = driver->GetMetadata();
    if( !CSLFetchBoolean( metadata, GDAL_DCAP_CREATE, FALSE ) )
        return false;
		
	GDALDataset * dstDS;     

	//set options
	char ** options = NULL;
	options = CSLSetNameValue( options, "COMPRESS", "LZW" );

	//if it is a GeoTIFF format set correct compression options
	if ( !strcmp( format, "GTiff" ) )
	{
		if( theType == GDT_Byte )
		{
			options = CSLSetNameValue( options, "PREDICTOR", "1" );
		}
		else
		{
			if ( theType == GDT_UInt16 || theType == GDT_Int16  
				|| theType == GDT_UInt32 || theType == GDT_Int32 )
			{
				options = CSLSetNameValue( options, "PREDICTOR", "2" );
			} 
			else 
			{
				options = CSLSetNameValue( options, "PREDICTOR", "3" );
			}
		}
	}
	
	//set filename
	char * dstFilename = new char[theFilename.size() + 1];
	strcpy( dstFilename, theFilename.toLocal8Bit().data() );

	dstDS = driver->Create( dstFilename, theCols, theRows, theBands, theType, 
								options );
	delete dstFilename;
	delete [] options;

	//set geotransform
	dstDS->SetGeoTransform( theGeoTransform );

	//set CRS
	char * crsWkt = new char[theCrs->toWkt().size() + 1];
	strcpy( crsWkt, theCrs->toWkt().toLocal8Bit().data());
    dstDS->SetProjection( crsWkt );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                    fprintf(gridFiles[k], "\n");
            }
    }

#ifdef HAVE_GDAL
    GDALDataset **gdalFiles;
    char gdalFileName[1024];

    // open GDAL GeoTIFF files
    if(outputFormat == OUTPUT_FORMAT_GDAL_GTIFF || outputFormat == OUTPUT_FORMAT_ALL)
    {
        GDALAllRegister();

        if((gdalFiles = (GDALDataset **)malloc(sizeof(GDALDataset *) *  numTypes)) == NULL)
        {
            cerr << "File array allocation error" << endl;
            return -1;
        }

        for(i = 0; i < numTypes; i++)
        {
            if(outputType & type[i])
            {
                strncpy(gdalFileName, outputName.c_str(), sizeof(gdalFileName));
                strncat(gdalFileName, ext[i], strlen(ext[i]));
                strncat(gdalFileName, ".tif", strlen(".tif"));

                char **papszMetadata;
                const char *pszFormat = "GTIFF";
                GDALDriver* tpDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);

                if (tpDriver)
                {
                    papszMetadata = tpDriver->GetMetadata();
                    if (CSLFetchBoolean(papszMetadata, GDAL_DCAP_CREATE, FALSE))
                    {
                        char **papszOptions = NULL;
                        gdalFiles[i] = tpDriver->Create(gdalFileName, GRID_SIZE_X, GRID_SIZE_Y, 1, GDT_Float32, papszOptions);
                        if (gdalFiles[i] == NULL)
                        {
                            cerr << "File open error: " << gdalFileName << endl;
                            return -1;
                        } else {
                            if (adfGeoTransform)
                            {
                                gdalFiles[i]->SetGeoTransform(adfGeoTransform);
                            }
                            else
                            {
                                double defaultTransform [6] = { min_x - 0.5*GRID_DIST_X,                            // top left x
                                                                (double)GRID_DIST_X,                                // w-e pixel resolution
                                                                0.0,                                                // no rotation/shear
                                                                min_y - 0.5*GRID_DIST_Y + GRID_DIST_Y*GRID_SIZE_Y,  // top left y
                                                                0.0,                                                // no rotation/shear
                                                                -(double)GRID_DIST_Y };                             // n-x pixel resolution (negative value)
                                gdalFiles[i]->SetGeoTransform(defaultTransform);
                            }
                            if (wkt)
                                gdalFiles[i]->SetProjection(wkt);
                        }
                    }
                }
            } else {
                gdalFiles[i] = NULL;
            }
        }

//.........这里部分代码省略.........
    poDstDS = (GDALDataset*) OGROpen( pszOutputName, TRUE, NULL );

/* -------------------------------------------------------------------- */
/*      If that failed, find the driver so we can create the tile index.*/
/* -------------------------------------------------------------------- */
    if( poDstDS == NULL )
    {
        OGRSFDriverRegistrar     *poR = OGRSFDriverRegistrar::GetRegistrar();
        GDALDriver              *poDriver = NULL;
        int                      iDriver;

        for( iDriver = 0;
             iDriver < poR->GetDriverCount() && poDriver == NULL;
             iDriver++ )
        {
            if( EQUAL(poR->GetDriver(iDriver)->GetDescription(),pszFormat) )
            {
                poDriver = poR->GetDriver(iDriver);
            }
        }

        if( poDriver == NULL )
        {
            fprintf( stderr, "Unable to find driver `%s'.\n", pszFormat );
            fprintf( stderr, "The following drivers are available:\n" );
        
            for( iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
            {
                fprintf( stderr, "  -> `%s'\n", poR->GetDriver(iDriver)->GetDescription() );
            }
            exit( 1 );
        }

        if( !CSLTestBoolean( CSLFetchNameValueDef(poDriver->GetMetadata(), GDAL_DCAP_CREATE, "FALSE") ) )
        {
            fprintf( stderr, "%s driver does not support data source creation.\n",
                    pszFormat );
            exit( 1 );
        }

/* -------------------------------------------------------------------- */
/*      Now create it.                                                  */
/* -------------------------------------------------------------------- */
        
        poDstDS = poDriver->Create( pszOutputName, 0, 0, 0, GDT_Unknown, NULL );
        if( poDstDS == NULL )
        {
            fprintf( stderr, "%s driver failed to create %s\n", 
                    pszFormat, pszOutputName );
            exit( 1 );
        }

        if( poDstDS->GetLayerCount() == 0 )
        {
            OGRFieldDefn oLocation( pszTileIndexField, OFTString );
            
            oLocation.SetWidth( 200 );
            
            if( nFirstSourceDataset < nArgc && papszArgv[nFirstSourceDataset][0] == '-' )
            {
                nFirstSourceDataset++;
            }
            
            OGRSpatialReference* poSrcSpatialRef = NULL;
            
            /* Fetches the SRS of the first layer and use it when creating the tileindex layer */

//---------------------------------------------------------
bool CGDAL_Export::On_Execute(void)
{
	char					**pOptions	= NULL;
	int						x, y, n;
	double					*zLine;
	CSG_String				File_Name;
	CSG_Parameter_Grid_List	*pGrids;
	CSG_Grid				*pGrid;
	GDALDataType			gdal_Type;
	GDALDriver				*pDriver;
	GDALDataset				*pDataset;
	GDALRasterBand			*pBand;

	//-----------------------------------------------------
	pGrids		= Parameters("GRIDS")	->asGridList();
	File_Name	= Parameters("FILE")	->asString();

	//-----------------------------------------------------
	switch( Parameters("TYPE")->asInt() )
	{
	default:
	case 0:	gdal_Type	= g_GDAL_Driver.Get_GDAL_Type(pGrids);	break;	// match input data
	case 1:	gdal_Type	= GDT_Byte;		break;	// Eight bit unsigned integer
	case 2:	gdal_Type	= GDT_UInt16;	break;	// Sixteen bit unsigned integer
	case 3:	gdal_Type	= GDT_Int16;	break;	// Sixteen bit signed integer
	case 4:	gdal_Type	= GDT_UInt32;	break;	// Thirty two bit unsigned integer
	case 5:	gdal_Type	= GDT_Int32;	break;	// Thirty two bit signed integer
	case 6:	gdal_Type	= GDT_Float32;	break;	// Thirty two bit floating point
	case 7:	gdal_Type	= GDT_Float64;	break;	// Sixty four bit floating point
	}

	//-----------------------------------------------------
	if( (pDriver = g_GDAL_Driver.Get_Driver(SG_STR_SGTOMB(m_DriverNames[Parameters("FORMAT")->asInt()]))) == NULL )
	{
		Message_Add(_TL("Driver not found."));
	}
	else if( CSLFetchBoolean(pDriver->GetMetadata(), GDAL_DCAP_CREATE, false) == false )
	{
		Message_Add(_TL("Driver does not support file creation."));
	}
	else if( (pDataset = pDriver->Create(File_Name.b_str(), Get_NX(), Get_NY(), pGrids->Get_Count(), gdal_Type, pOptions)) == NULL )
	{
		Message_Add(_TL("Could not create dataset."));
	}
	else
	{
		g_GDAL_Driver.Set_Transform(pDataset, Get_System());

		if( pGrids->asGrid(0)->Get_Projection().Get_Type() != SG_PROJ_TYPE_CS_Undefined )
		{
			pDataset->SetProjection(SG_STR_SGTOMB(pGrids->asGrid(0)->Get_Projection().Get_WKT()));
		}

		zLine	= (double *)SG_Malloc(Get_NX() * sizeof(double));

		for(n=0; n<pGrids->Get_Count(); n++)
		{
			Process_Set_Text(CSG_String::Format(SG_T("%s %d"), _TL("Band"), n + 1));

			pGrid	= pGrids->asGrid(n);
			pBand	= pDataset->GetRasterBand(n + 1);

			for(y=0; y<Get_NY() && Set_Progress(y, Get_NY()); y++)
			{
				for(x=0; x<Get_NX(); x++)
				{
					zLine[x]	= pGrid->asDouble(x, Get_NY() - 1 - y);
				}

				pBand->RasterIO(GF_Write, 0, y, Get_NX(), 1, zLine, Get_NX(), 1, GDT_Float64, 0, 0);
			}
		}

		//-------------------------------------------------
		SG_Free(zLine);

		GDALClose(pDataset);

		return( true );
	}

	//-----------------------------------------------------
	return( false );
}

  int Delaunay(maps*& conf,maps*& inputs,maps*& outputs){
#ifdef DEBUG
    fprintf(stderr,"\nService internal print\nStarting\n");
#endif
    maps* cursor=inputs;
    OGRGeometryH geometry,res;
    int bufferDistance;
    map* tmpm=NULL;
    OGRRegisterAll();

    std::vector<Point> points;
    if(int res=parseInput(conf,inputs,&points,"/vsimem/tmp")!=SERVICE_SUCCEEDED)
      return res;

    DelaunayT T;
    T.insert(points.begin(), points.end());

    /* -------------------------------------------------------------------- */
    /*      Try opening the output datasource as an existing, writable      */
    /* -------------------------------------------------------------------- */
#if GDAL_VERSION_MAJOR >= 2
    GDALDataset *poODS;
    GDALDriverManager* poR=GetGDALDriverManager();
    GDALDriver          *poDriver = NULL;
#else
    OGRDataSource       *poODS;    
    OGRSFDriverRegistrar *poR = OGRSFDriverRegistrar::GetRegistrar();
    OGRSFDriver          *poDriver = NULL;
#endif
    int                  iDriver;
    map *tmpMap=getMapFromMaps(outputs,"Result","mimeType");
    const char *oDriver;
    oDriver="GeoJSON";
    if(tmpMap!=NULL){
      if(strcmp(tmpMap->value,"text/xml")==0){
	oDriver="GML";
      }
    }
    
    for( iDriver = 0;
	 iDriver < poR->GetDriverCount() && poDriver == NULL;
	 iDriver++ )
      {
#ifdef DEBUG
#if GDAL_VERSION_MAJOR >= 2
	fprintf(stderr,"D:%s\n",poR->GetDriver(iDriver)->GetDescription());
#else
	fprintf(stderr,"D:%s\n",poR->GetDriver(iDriver)->GetName());
#endif
#endif
	if( EQUAL(
#if GDAL_VERSION_MAJOR >= 2
		  poR->GetDriver(iDriver)->GetDescription()
#else
		  poR->GetDriver(iDriver)->GetName()
#endif
		  ,
		  oDriver) )
	  {
	    poDriver = poR->GetDriver(iDriver);
	  }
      }

    if( poDriver == NULL )
      {
	char emessage[8192];
	sprintf( emessage, "Unable to find driver `%s'.\n", oDriver );
	sprintf( emessage,  "%sThe following drivers are available:\n",emessage );
        
	for( iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
	  {
#if GDAL_VERSION_MAJOR >= 2
	    sprintf( emessage,  "%s  -> `%s'\n", emessage, poR->GetDriver(iDriver)->GetDescription() );
#else
	    sprintf( emessage,  "%s  -> `%s'\n", emessage, poR->GetDriver(iDriver)->GetName() );
#endif
	  }

	setMapInMaps(conf,"lenv","message",emessage);
	return SERVICE_FAILED;

      }

#if GDAL_VERSION_MAJOR >=2
    if( !CPLTestBool( CSLFetchNameValueDef(poDriver->GetMetadata(), GDAL_DCAP_CREATE, "FALSE") ) )
#else
    if( !poDriver->TestCapability( ODrCCreateDataSource ) )
#endif
      {
	char emessage[1024];
	sprintf( emessage,  "%s driver does not support data source creation.\n",
		 "json" );
	setMapInMaps(conf,"lenv","message",emessage);
	return SERVICE_FAILED;
      }

    /* -------------------------------------------------------------------- */
    /*      Create the output data source.                                  */
    /* -------------------------------------------------------------------- */
    char *pszDestDataSource=(char*)malloc(100);
//.........这里部分代码省略.........

MAIN_START(argc, argv)
{
    /* Check strict compilation and runtime library version as we use C++ API */
    if (! GDAL_CHECK_VERSION(argv[0]))
        exit(1);

    EarlySetConfigOptions(argc, argv);

/* -------------------------------------------------------------------- */
/*      Generic arg processing.                                         */
/* -------------------------------------------------------------------- */
    GDALAllRegister();
    argc = GDALGeneralCmdLineProcessor(argc, &argv, 0);
    if( argc < 1 )
        exit( -argc );

    for( int i = 0; i < argc; i++ )
    {
        if( EQUAL(argv[i], "--utility_version") )
        {
            printf("%s was compiled against GDAL %s and "
                   "is running against GDAL %s\n",
                   argv[0], GDAL_RELEASE_NAME, GDALVersionInfo("RELEASE_NAME"));
            CSLDestroy(argv);
            return 0;
        }
        else if( EQUAL(argv[i], "--help") )
        {
            Usage();
        }
    }

    GDALRasterizeOptionsForBinary* psOptionsForBinary =
        GDALRasterizeOptionsForBinaryNew();
    // coverity[tainted_data]
    GDALRasterizeOptions *psOptions =
        GDALRasterizeOptionsNew(argv + 1, psOptionsForBinary);
    CSLDestroy(argv);

    if( psOptions == nullptr )
    {
        Usage();
    }

    if( !(psOptionsForBinary->bQuiet) )
    {
        GDALRasterizeOptionsSetProgress(psOptions, GDALTermProgress, nullptr);
    }

    if( psOptionsForBinary->pszSource == nullptr )
        Usage("No input file specified.");

    if( psOptionsForBinary->pszDest == nullptr )
        Usage("No output file specified.");

/* -------------------------------------------------------------------- */
/*      Open input file.                                                */
/* -------------------------------------------------------------------- */
    GDALDatasetH hInDS = GDALOpenEx(
        psOptionsForBinary->pszSource, GDAL_OF_VECTOR | GDAL_OF_VERBOSE_ERROR,
        nullptr, nullptr, nullptr);

    if( hInDS == nullptr )
        exit(1);

/* -------------------------------------------------------------------- */
/*      Open output file if it exists.                                  */
/* -------------------------------------------------------------------- */
    GDALDatasetH hDstDS = nullptr;
    if( !(psOptionsForBinary->bCreateOutput) )
    {
        CPLPushErrorHandler(CPLQuietErrorHandler);
        hDstDS = GDALOpenEx(
            psOptionsForBinary->pszDest,
            GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR | GDAL_OF_UPDATE,
            nullptr, nullptr, nullptr );
        CPLPopErrorHandler();
    }

    if( psOptionsForBinary->pszFormat != nullptr &&
        (psOptionsForBinary->bCreateOutput || hDstDS == nullptr) )
    {
        GDALDriverManager *poDM = GetGDALDriverManager();
        GDALDriver *poDriver =
            poDM->GetDriverByName(psOptionsForBinary->pszFormat);
        char** papszDriverMD = (poDriver) ? poDriver->GetMetadata(): nullptr;
        if( poDriver == nullptr ||
            !CPLTestBool(CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_RASTER,
                                              "FALSE")) ||
            !CPLTestBool(CSLFetchNameValueDef(papszDriverMD, GDAL_DCAP_CREATE,
                                              "FALSE")) )
        {
            fprintf(stderr,
                    "Output driver `%s' not recognised or does not support "
                    "direct output file creation.\n",
                    psOptionsForBinary->pszFormat);
            fprintf(stderr,
                    "The following format drivers are configured and "
                    "support direct output:\n" );

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

void ccRasterizeTool::generateRaster() const
{
#ifdef CC_GDAL_SUPPORT

	if (!m_cloud || !m_grid.isValid())
		return;

	GDALAllRegister();
	ccLog::PrintDebug("(GDAL drivers: %i)", GetGDALDriverManager()->GetDriverCount());

	const char *pszFormat = "GTiff";
	GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
	if (!poDriver)
	{
		ccLog::Error("[GDAL] Driver %s is not supported", pszFormat);
		return;
	}

	char** papszMetadata = poDriver->GetMetadata();
	if( !CSLFetchBoolean( papszMetadata, GDAL_DCAP_CREATE, FALSE ) )
	{
		ccLog::Error("[GDAL] Driver %s doesn't support Create() method", pszFormat);
		return;
	}

	//which (and how many) bands shall we create?
	bool heightBand = true; //height by default
	bool densityBand = false;
	bool allSFBands = false;
	int sfBandIndex = -1; //scalar field index
	int totalBands = 0;

	bool interpolateSF = (getTypeOfSFInterpolation() != INVALID_PROJECTION_TYPE);
	ccPointCloud* pc = m_cloud->isA(CC_TYPES::POINT_CLOUD) ? static_cast<ccPointCloud*>(m_cloud) : 0;

	bool hasSF =  interpolateSF && pc && !m_grid.scalarFields.empty();
	
	RasterExportOptionsDlg reoDlg;
	reoDlg.dimensionsLabel->setText(QString("%1 x %2").arg(m_grid.width).arg(m_grid.height));
	reoDlg.exportHeightsCheckBox->setChecked(heightBand);
	reoDlg.exportDensityCheckBox->setChecked(densityBand);
	reoDlg.exportDisplayedSFCheckBox->setEnabled(hasSF);
	reoDlg.exportAllSFCheckBox->setEnabled(hasSF);
	reoDlg.exportAllSFCheckBox->setChecked(allSFBands);

	if (!reoDlg.exec())
		return;

	//we ask the output filename AFTER displaying the export parameters ;)
	QString outputFilename;
	{
		QSettings settings;
		settings.beginGroup(ccPS::HeightGridGeneration());
		QString imageSavePath = settings.value("savePathImage",QApplication::applicationDirPath()).toString();
		outputFilename = QFileDialog::getSaveFileName(0,"Save height grid raster",imageSavePath+QString("/raster.tif"),"geotiff (*.tif)");

		if (outputFilename.isNull())
			return;

		//save current export path to persistent settings
		settings.setValue("savePathImage",QFileInfo(outputFilename).absolutePath());
	}

	heightBand = reoDlg.exportHeightsCheckBox->isChecked();
	densityBand = reoDlg.exportDensityCheckBox->isChecked();
	if (hasSF)
	{
		assert(pc);
		allSFBands = reoDlg.exportAllSFCheckBox->isChecked() && hasSF;
		if (!allSFBands && reoDlg.exportDisplayedSFCheckBox->isChecked())
		{
			sfBandIndex = pc->getCurrentDisplayedScalarFieldIndex();
			if (sfBandIndex < 0)
				ccLog::Warning("[Rasterize] Cloud has no active (displayed) SF!");
		}
	}

	totalBands = heightBand ? 1 : 0;
	if (densityBand)
	{
		++totalBands;
	}
	if (allSFBands)
	{
		assert(hasSF);
		for (size_t i=0; i<m_grid.scalarFields.size(); ++i)
			if (m_grid.scalarFields[i])
				++totalBands;
	}
	else if (sfBandIndex >= 0)
	{
		++totalBands;
	}
	
	if (totalBands == 0)
	{
		ccLog::Warning("[Rasterize] Warning, can't output a raster with no band! (check export parameters)");
		return;
	}

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

static void Usage(const char* pszAdditionalMsg, int bShort)

{
    OGRSFDriverRegistrar        *poR = OGRSFDriverRegistrar::GetRegistrar();


    printf( "Usage: ogr2ogr [--help-general] [-skipfailures] [-append] [-update]\n"
            "               [-select field_list] [-where restricted_where|@filename]\n"
            "               [-progress] [-sql <sql statement>|@filename] [-dialect dialect]\n"
            "               [-preserve_fid] [-fid FID]\n"
            "               [-spat xmin ymin xmax ymax] [-spat_srs srs_def] [-geomfield field]\n"
            "               [-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]\n"
            "               [-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]\n"
            "               dst_datasource_name src_datasource_name\n"
            "               [-lco NAME=VALUE] [-nln name] \n"
            "               [-nlt type|PROMOTE_TO_MULTI|CONVERT_TO_LINEAR|CONVERT_TO_CURVE]\n"
            "               [-dim 2|3|layer_dim] [layer [layer ...]]\n"
            "\n"
            "Advanced options :\n"
            "               [-gt n] [-ds_transaction]\n"
            "               [[-oo NAME=VALUE] ...] [[-doo NAME=VALUE] ...]\n"
            "               [-clipsrc [xmin ymin xmax ymax]|WKT|datasource|spat_extent]\n"
            "               [-clipsrcsql sql_statement] [-clipsrclayer layer]\n"
            "               [-clipsrcwhere expression]\n"
            "               [-clipdst [xmin ymin xmax ymax]|WKT|datasource]\n"
            "               [-clipdstsql sql_statement] [-clipdstlayer layer]\n"
            "               [-clipdstwhere expression]\n"
            "               [-wrapdateline][-datelineoffset val]\n"
            "               [[-simplify tolerance] | [-segmentize max_dist]]\n"
            "               [-addfields] [-unsetFid]\n"
            "               [-relaxedFieldNameMatch] [-forceNullable] [-unsetDefault]\n"
            "               [-fieldTypeToString All|(type1[,type2]*)] [-unsetFieldWidth]\n"
            "               [-mapFieldType srctype|All=dsttype[,srctype2=dsttype2]*]\n"
            "               [-fieldmap identity | index1[,index2]*]\n"
            "               [-splitlistfields] [-maxsubfields val]\n"
            "               [-explodecollections] [-zfield field_name]\n"
            "               [-gcp pixel line easting northing [elevation]]* [-order n | -tps]\n"
            "               [-nomd] [-mo \"META-TAG=VALUE\"]* [-noNativeData]\n");

    if (bShort)
    {
        printf( "\nNote: ogr2ogr --long-usage for full help.\n");
        if( pszAdditionalMsg )
            fprintf(stderr, "\nFAILURE: %s\n", pszAdditionalMsg);
        exit( 1 );
    }

    printf("\n -f format_name: output file format name, possible values are:\n");

    std::vector<CPLString> aoSetDrivers;
    for( int iDriver = 0; iDriver < poR->GetDriverCount(); iDriver++ )
    {
        GDALDriver *poDriver = poR->GetDriver(iDriver);

        if( CPLTestBool( CSLFetchNameValueDef(poDriver->GetMetadata(), GDAL_DCAP_CREATE, "FALSE") ) )
            aoSetDrivers.push_back( poDriver->GetDescription() );
    }
    std::sort (aoSetDrivers.begin(), aoSetDrivers.end(), StringCISortFunction);
    for( size_t i = 0; i < aoSetDrivers.size(); i++ )
    {
        printf( "     -f \"%s\"\n", aoSetDrivers[i].c_str() );
    }

    printf( " -append: Append to existing layer instead of creating new if it exists\n"
            " -overwrite: delete the output layer and recreate it empty\n"
            " -update: Open existing output datasource in update mode\n"
            " -progress: Display progress on terminal. Only works if input layers have the \n"
            "                                          \"fast feature count\" capability\n"
            " -select field_list: Comma-delimited list of fields from input layer to\n"
            "                     copy to the new layer (defaults to all)\n"
            " -where restricted_where: Attribute query (like SQL WHERE)\n"
            " -wrapdateline: split geometries crossing the dateline meridian\n"
            "                (long. = +/- 180deg)\n"
            " -datelineoffset: offset from dateline in degrees\n"
            "                (default long. = +/- 10deg,\n"
            "                geometries within 170deg to -170deg will be split)\n"
            " -sql statement: Execute given SQL statement and save result.\n"
            " -dialect value: select a dialect, usually OGRSQL to avoid native sql.\n"
            " -skipfailures: skip features or layers that fail to convert\n"
            " -gt n: group n features per transaction (default 20000). n can be set to unlimited\n"
            " -spat xmin ymin xmax ymax: spatial query extents\n"
            " -simplify tolerance: distance tolerance for simplification.\n"
            " -segmentize max_dist: maximum distance between 2 nodes.\n"
            "                       Used to create intermediate points\n"
            " -dsco NAME=VALUE: Dataset creation option (format specific)\n"
            " -lco  NAME=VALUE: Layer creation option (format specific)\n"
            " -oo   NAME=VALUE: Input dataset open option (format specific)\n"
            " -doo  NAME=VALUE: Destination dataset open option (format specific)\n"
            " -nln name: Assign an alternate name to the new layer\n"
            " -nlt type: Force a geometry type for new layer.  One of NONE, GEOMETRY,\n"
            "      POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION, MULTIPOINT,\n"
            "      MULTIPOLYGON, or MULTILINESTRING, or PROMOTE_TO_MULTI or CONVERT_TO_LINEAR.  Add \"25D\" for 3D layers.\n"
            "      Default is type of source layer.\n"
            " -dim dimension: Force the coordinate dimension to the specified value.\n"
            " -fieldTypeToString type1,...: Converts fields of specified types to\n"
            "      fields of type string in the new layer. Valid types are : Integer,\n"
            "      Integer64, Real, String, Date, Time, DateTime, Binary, IntegerList, Integer64List, RealList,\n"
            "      StringList. Special value All will convert all fields to strings.\n"
            " -fieldmap index1,index2,...: Specifies the list of field indexes to be\n"
            "      copied from the source to the destination. The (n)th value specified\n"
//.........这里部分代码省略.........