C++ OGRLayer::GetFeatureCount方法代码示例

void FeatureCursorCDBV::Read_Data_Directly(void)
{
	FeatureList preProcessList;

	OGR_SCOPED_LOCK;

	OGR_L_ResetReading(_layerHandle);
	OGRFeatureH feat_handle;

	OGRLayer * thislayer = (OGRLayer *)_layerHandle;
	int totalCount = thislayer->GetFeatureCount();
	int fcount = -1;
	while ((feat_handle = OGR_L_GetNextFeature(_layerHandle)) != NULL)
	{
		++fcount;
		if (feat_handle)
		{
			osg::ref_ptr<Feature> f = OgrUtils::createFeature(feat_handle, _profile.get());
			if (f.valid() && !_source->isBlacklisted(f->getFID()))
			{
				if (isGeometryValid(f->getGeometry()))
				{
					_queue.push(f);

					if (_filters.size() > 0)
					{
						preProcessList.push_back(f.release());
					}
				}
				else
				{
					OE_DEBUG << LC << "Skipping feature with invalid geometry: " << f->getGeoJSON() << std::endl;
				}
			}
			OGR_F_Destroy(feat_handle);
		}
	}

	// preprocess the features using the filter list:
	if (preProcessList.size() > 0)
	{
		FilterContext cx;
		cx.setProfile(_profile.get());

		for (FeatureFilterList::const_iterator i = _filters.begin(); i != _filters.end(); ++i)
		{
			FeatureFilter* filter = i->get();
			cx = filter->push(preProcessList, cx);
		}
	}

}

void MSN_Helper::LoadSamples(Config &config, Matrix &mat)
{
	if(config.bShapefile)
	{
		OGRRegisterAll();
		string pLayerName = StringGetFileName(config.sSamples);
		OGRDataSource* poDS = OGRSFDriverRegistrar::Open(config.sSamples.c_str(),FALSE);
		OGRLayer* poLayer = poDS->GetLayerByName(pLayerName.c_str());
	
		config.nSamples = poLayer->GetFeatureCount();
		mat.Resize(config.nSamples, 4);

		OGRPoint *poPoint;
		OGRFeature * pFeature =poLayer->GetNextFeature();
		for(unsigned long i=1; i<=config.nSamples; i++)
		{
			//样本取值字段名为value，double型
			poPoint = (OGRPoint *)pFeature->GetGeometryRef();
			mat[i][1] = poPoint->getX();
			mat[i][2] = poPoint->getY();
			mat[i][3] = pFeature->GetFieldAsInteger("stratum");
			mat[i][4] = pFeature->GetFieldAsDouble("value");
			pFeature = poLayer->GetNextFeature();
		}

		OGRDataSource::DestroyDataSource( poDS );
	}
	else
	{
		double x, y, v, stratum;
		string sline;
		ifstream infile(config.sSamples.c_str());
		infile >> config.nSamples;

		mat.Resize(config.nSamples, 4);
		for(unsigned long i=1; i<=config.nSamples; i++)
		{
			infile >> x >> y >> stratum >> v;
			mat[i][1] = x;
			mat[i][2] = y;
			mat[i][3] = stratum;
			mat[i][4] = v;
		}
		infile.close();
	}
}

bool ShpReader::Open(std::wstring fullPath, StudyControllerPtr studyController,
                     VectorMapControllerPtr vectorMapController, ProgressDlgPtr progressDlg) 
{
	// Registers all format drivers built into GDAL/OGR.
	OGRRegisterAll();

	OGRDataSource *OGRDataset;

	// Open vector file path
	std::string tempStr( fullPath.begin(), fullPath.end() );
	OGRDataset = OGRSFDriverRegistrar::Open( tempStr.c_str(), FALSE );

	// Return if no vector files are found
	if( OGRDataset == NULL )
	{
		Log::Inst().Warning("(Warning) Failed to open file at " + tempStr + ".");
		return false;
	}
	if ( App::Inst().GetLayerTreeController()->GetNumMapLayers() >0 )

		MapControllerPtr mapController= App::Inst().GetLayerTreeController()->GetLayerTreeModel()->GetStudy(0)->GetMapLayer(0)->GetMapController();

	// It appears that shapefiles (*.SHP) only support up to one layer per file
	// This will need to be further investigated for other vector filetypes (e.g., KML)
	// For now just grab layer at position 0
	
	OGRLayer *poLayer = OGRDataset->GetLayer( 0 );
	
	// Determine the XY boundaries for the entire vector dataset
	OGREnvelope psEnvelope;

	VectorMapModelPtr vectorMapModel = vectorMapController->GetVectorMapModel();
	poLayer->GetExtent( &psEnvelope );
	vectorMapModel->SetVectorBoundary(psEnvelope);

	if(!SetupVectorProjection(OGRDataset,studyController,poLayer,vectorMapController ))
	{
		OGRDataset->DestroyDataSource(OGRDataset);
		return false;
	}
	
	if(progressDlg)
	{
			if(!progressDlg->Update(0, _T("Reading Vector Map Information...")))
			{
				OGRDataset->DestroyDataSource(OGRDataset);
				return false;
			}
	}	
	GLdouble minX, minY, maxX, maxY;
	minX = minY = std::numeric_limits<float>::max();
	maxX = maxY = -std::numeric_limits<float>::max();

	// Retrieve features from the dataset
	OGRFeature *poFeature;
	poLayer->ResetReading();
	int numFeatures = poLayer->GetFeatureCount();
	int count=0;
    //Log::Inst().Write("Loading shapefile with the following meta data:");

	while ( ( poFeature = poLayer->GetNextFeature() ) != NULL )
	{
		/////////////////////////////////////////////////
		// PHASE 1: Retrieve METADATA from the dataset //
		/////////////////////////////////////////////////
		OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
		int iField;

		for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ )
		{
			OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField );

			//if( poFieldDefn->GetType() == OFTInteger )
			//    printf( "%d,", poFeature->GetFieldAsInteger( iField ) );
			//else if( poFieldDefn->GetType() == OFTReal )
			//    printf( "%.3f,", poFeature->GetFieldAsDouble(iField) );
			//else if( poFieldDefn->GetType() == OFTString )
			//    printf( "%s,", poFeature->GetFieldAsString(iField) );
			//else
			//    printf( "%s,", poFeature->GetFieldAsString(iField) );

			//ofs << poFeature->GetFieldAsString(iField) << ",";

			std::string metaData = poFeature->GetFieldAsString(iField);
			// do something with the meta data...
			//Log::Inst().Write(metaData);
		}
		count++;
		if(progressDlg)
		{
			if(!progressDlg->Update(int(50 + (float(count)/numFeatures)*50)))
				return false;
		}

		///////////////////////////////////////////////////
		// PHASE 2: Retrieve GEOMETRIES from the dataset //
		///////////////////////////////////////////////////
		OGRGeometry *poGeometry;
		poGeometry = poFeature->GetGeometryRef();

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

//.........这里部分代码省略.........
            CPLError(CE_Failure, CPLE_AppDefined,
                    "Cannot open file %s", pszFilename);
            return FALSE;
        }

        char** argv = GetArgv(bExplicitFeatures, bWaypoints, bRoutes,
                              bTracks, pszGPSBabelDriverName, "-");

        VSILFILE* tmpfp = VSIFOpenL(osTmpFileName.c_str(), "wb");

        CPLPushErrorHandler(CPLQuietErrorHandler);
        bRet = (CPLSpawn(argv, fp, tmpfp, TRUE) == 0);
        CPLPopErrorHandler();

        CSLDestroy(argv);
        argv = NULL;

        CPLErr nLastErrorType = CPLGetLastErrorType();
        CPLErrorNum nLastErrorNo = CPLGetLastErrorNo();
        CPLString osLastErrorMsg = CPLGetLastErrorMsg();

        VSIFCloseL(tmpfp);
        tmpfp = NULL;

        VSIFCloseL(fp);
        fp = NULL;

        if (!bRet)
        {
            if ( strstr(osLastErrorMsg.c_str(),
                        "This format cannot be used in piped commands") == NULL)
            {
                CPLError( nLastErrorType, nLastErrorNo, "%s",
                          osLastErrorMsg.c_str());
            }
            else
            {
                VSIStatBuf sStatBuf;
                if (VSIStat(pszFilename, &sStatBuf) != 0)
                {
                    CPLError( CE_Failure, CPLE_NotSupported,
                              "Driver %s only supports real (non virtual) "
                              "files",
                              pszGPSBabelDriverName );
                    return FALSE;
                }

                /* Try without piping in */
                argv = GetArgv(bExplicitFeatures, bWaypoints, bRoutes,
                              bTracks, pszGPSBabelDriverName, pszFilename);
                tmpfp = VSIFOpenL(osTmpFileName.c_str(), "wb");
                bRet = (CPLSpawn(argv, NULL, tmpfp, TRUE) == 0);
                VSIFCloseL(tmpfp);
                tmpfp = NULL;

                CSLDestroy(argv);
                argv = NULL;
            }
        }
    }


    if (bRet)
    {
        poGPXDS = static_cast<GDALDataset *>(
            GDALOpenEx( osTmpFileName.c_str(),
                        GDAL_OF_VECTOR, NULL, NULL, NULL ) );
        if (poGPXDS)
        {
            if (bWaypoints)
            {
                OGRLayer* poLayer = poGPXDS->GetLayerByName("waypoints");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
            }

            if (bRoutes)
            {
                OGRLayer* poLayer = poGPXDS->GetLayerByName("routes");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
                poLayer = poGPXDS->GetLayerByName("route_points");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
            }

            if (bTracks)
            {
                OGRLayer* poLayer = poGPXDS->GetLayerByName("tracks");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
                poLayer = poGPXDS->GetLayerByName("track_points");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
            }
        }
    }

    return nLayers > 0;
}

bool ShortNavigation::createObstaclesTables()
{
    //    OGRLayer* layer;
    status = true;
    if ( !status==false ) {
        qDebug()<<"bool ShortNavigation::createObstaclesTables()";

        /// POINT_OFFSET is the area around the obstacle,
        /// so is it's value is 5, it will draw 10 m^2

        //float POINT_OFFSET = settingsManager->getPOffset();
        float POINT_OFFSET = 5; //test values, above comment from old code
        //float NOT_ALLOWED_DEPTH = settingsManager->getMinDepth();
        float NOT_ALLOWED_DEPTH = 2; //test values, above comment from old code

        /// signsound with depth from 0 to this value are
        /// taken in consideration for the obstacles layer

        //float SIGNSOUND_THRESHOLD = settingsManager->getSignsound();
        float SIGNSOUND_THRESHOLD = 1; //test values, above comment from old code



        // ##################################################
        // list of polygon layers considered obstacles
        QList<QString> polygon_layers;
        polygon_layers.append("generarea_r");
        // list of point layers to be considered as obstacles
        QList<QString> point_layers;
        point_layers.append("rock_p");
        point_layers.append("wreck_p");
        point_layers.append("navaid_p");
        point_layers.append("signsound_p");
        // line layers to be considered as obstacles
        QList<QString> line_layers;
        line_layers.append("depthcont_l");
        //####################################################

        for (int i = 0 ; i < this->chartObjects.size(); i++) {
            //        qDebug() << "chartObjects size" << chartObjects.size();

            if(this->chartObjects.at(i)->getTableName() == "generarea_r") {
                for (int j = 0 ; j < this->chartObjects.at(i)->getCoordinateGeometry().size();j++) {
                    this->polyObstacles.append(this->chartObjects.at(i)->getCoordinateGeometry().at(j));
                }
            }
            if(this->chartObjects.at(i)->getTableName() == "depthcont_l") {
                for (int j = 0 ; j < this->chartObjects.at(i)->getCoordinateGeometry().size();j++) { //TODO: parse line from polygon
                    //this->lineObstacles.append(this->chartObjects.at(i)->getCoordinateGeometry().at(j)); //Not used in Shortnavigation -> can be removed!
                }
            }
            if(this->chartObjects.at(i)->getTableName() == "rock_p") {
                for (int j = 0 ; j < this->chartObjects.at(i)->getCoordinateGeometry().size();j++) { //need to check how vectors are written from DB. is getfeaturecount = polygonvector.size ?
                    for(int k = 0; k < this->chartObjects.at(i)->getCoordinateGeometry().at(j).size();k++) {
                        this->pointObstacles.append(this->chartObjects.at(i)->getCoordinateGeometry().at(j).at(k));
                    }
                }
            }
            if(this->chartObjects.at(i)->getTableName() == "wreck_p") {
                for (int j = 0 ; j < this->chartObjects.at(i)->getCoordinateGeometry().size();j++) {
                    for(int k = 0; k < this->chartObjects.at(i)->getCoordinateGeometry().at(j).size();k++) {
                        this->pointObstacles.append(this->chartObjects.at(i)->getCoordinateGeometry().at(j).at(k));
                    }
                }
            }
            if(this->chartObjects.at(i)->getTableName() == "navaid_p") {
                for (int j = 0 ; j < this->chartObjects.at(i)->getCoordinateGeometry().size();j++) {
                    for(int k = 0; k < this->chartObjects.at(i)->getCoordinateGeometry().at(j).size();k++) {
                        this->pointObstacles.append(this->chartObjects.at(i)->getCoordinateGeometry().at(j).at(k));

                    }
                }
            }
            if(this->chartObjects.at(i)->getTableName() == "signsound_p") {

                QVector<QPolygonF> qpoint;
                QList<QPointF> listqpoint;
                OGRLayer* layer;
                OGRFeatureDefn *poFDefn;
                OGRFeature *poFeature;

                qpoint = this->chartObjects.at(i)->getCoordinateGeometry();
                for(int n = 0; n < qpoint.size(); n++){
                    listqpoint = qpoint.value(n).toList();
                }

                layer = this->chartObjects.at(i)->getFeatureData();
                layer->ResetReading();

                QString /*sql("SELECT * FROM *"); //FIX THIS! //*/ sql("SELECT * FROM ( SELECT DISTINCT Intersects( wkb_geometry, ");
                for(int m= 0; m < layer->GetFeatureCount();m++ ) {
                    poFDefn = layer->GetLayerDefn();
                    poFeature = layer->GetNextFeature();
                    for(int j = 0; j < poFDefn->GetFieldCount(); j++ ){
                        OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( j );
                        QString str = poFieldDefn->GetNameRef();
                        if(str.contains("depth") == true) {
                            if (poFeature->GetFieldAsDouble(j) < SIGNSOUND_THRESHOLD) { //parse out depths above signsound threshold since they are not obstacles
                                pointObstacles.append(listqpoint.at(m));
                            }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    {
        VSILFILE* fp = VSIFOpenL(pszFilename, "rb");
        if (fp == NULL)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                    "Cannot open file %s", pszFilename);
            return FALSE;
        }

        char** argv = GetArgv(bExplicitFeatures, bWaypoints, bRoutes,
                              bTracks, pszGPSBabelDriverName, "-");

        VSILFILE* tmpfp = VSIFOpenL(osTmpFileName.c_str(), "wb");

        CPLPushErrorHandler(CPLQuietErrorHandler);
        bRet = (CPLSpawn(argv, fp, tmpfp, TRUE) == 0);
        CPLPopErrorHandler();

        CSLDestroy(argv);
        argv = NULL;

        CPLErr nLastErrorType = CPLGetLastErrorType();
        int nLastErrorNo = CPLGetLastErrorNo();
        CPLString osLastErrorMsg = CPLGetLastErrorMsg();

        VSIFCloseL(tmpfp);
        tmpfp = NULL;

        VSIFCloseL(fp);
        fp = NULL;

        if (!bRet)
        {
            if (strstr(osLastErrorMsg.c_str(), "This format cannot be used in piped commands") == NULL)
            {
                CPLError(nLastErrorType, nLastErrorNo, "%s", osLastErrorMsg.c_str());
            }
            else
            {
                VSIStatBuf sStatBuf;
                if (VSIStat(pszFilename, &sStatBuf) != 0)
                {
                    CPLError(CE_Failure, CPLE_NotSupported,
                            "Driver %s only supports real (non virtual) files", pszGPSBabelDriverName);
                    return FALSE;
                }

                /* Try without piping in */
                argv = GetArgv(bExplicitFeatures, bWaypoints, bRoutes,
                              bTracks, pszGPSBabelDriverName, pszFilename);
                tmpfp = VSIFOpenL(osTmpFileName.c_str(), "wb");
                bRet = (CPLSpawn(argv, NULL, tmpfp, TRUE) == 0);
                VSIFCloseL(tmpfp);
                tmpfp = NULL;

                CSLDestroy(argv);
                argv = NULL;
            }
        }
    }


    if (bRet)
    {
        poGPXDS = OGRSFDriverRegistrar::Open(osTmpFileName.c_str());
        if (poGPXDS)
        {
            OGRLayer* poLayer;

            if (bWaypoints)
            {
                poLayer = poGPXDS->GetLayerByName("waypoints");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
            }

            if (bRoutes)
            {
                poLayer = poGPXDS->GetLayerByName("routes");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
                poLayer = poGPXDS->GetLayerByName("route_points");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
            }

            if (bTracks)
            {
                poLayer = poGPXDS->GetLayerByName("tracks");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
                poLayer = poGPXDS->GetLayerByName("track_points");
                if (poLayer != NULL && poLayer->GetFeatureCount() != 0)
                    apoLayers[nLayers++] = poLayer;
            }
        }
    }

    return nLayers > 0;
}

//.........这里部分代码省略.........
        if (layer_by_name)
        {
            s << "cannot find layer by name '" << *layer_by_name;
        }
        else if (layer_by_index)
        {
            s << "cannot find layer by index number '" << *layer_by_index;
        }
        else if (layer_by_sql)
        {
            s << "cannot find layer by sql query '" << *layer_by_sql;
        }

        s << "' in dataset '" << dataset_name_ << "'";

        throw datasource_exception(s.str());
    }

    // work with real OGR layer
    OGRLayer* layer = layer_.layer();

    // initialize envelope
    boost::optional<std::string> ext = params.get<std::string>("extent");
    if (ext && !ext->empty())
    {
        extent_.from_string(*ext);
    }
    else
    {
        OGREnvelope envelope;
        OGRErr e = layer->GetExtent(&envelope);
        if (e == OGRERR_FAILURE)
        {
            if (layer->GetFeatureCount() == 0)
            {
                MAPNIK_LOG_ERROR(ogr) << "could not determine extent, layer '" << layer->GetLayerDefn()->GetName() << "' appears to have no features";
            }
            else
            {
                std::ostringstream s;
                s << "OGR Plugin: Cannot determine extent for layer '" << layer->GetLayerDefn()->GetName() << "'. Please provide a manual extent string (minx,miny,maxx,maxy).";
                throw datasource_exception(s.str());
            }
        }
        extent_.init(envelope.MinX, envelope.MinY, envelope.MaxX, envelope.MaxY);
    }

    // scan for index file
    // TODO - layer names don't match dataset name, so this will break for
    // any layer types of ogr than shapefiles, etc
    // fix here and in ogrindex
    size_t breakpoint = dataset_name_.find_last_of(".");
    if (breakpoint == std::string::npos)
    {
        breakpoint = dataset_name_.length();
    }
    index_name_ = dataset_name_.substr(0, breakpoint) + ".ogrindex";

#if defined (_WINDOWS)
    std::ifstream index_file(mapnik::utf8_to_utf16(index_name_), std::ios::in | std::ios::binary);
#else
    std::ifstream index_file(index_name_.c_str(), std::ios::in | std::ios::binary);
#endif

    if (index_file)
    {

//extern "C" {
SEXP ogrFIDs(SEXP filename, SEXP layer) {
    SEXP fids, nf, ii;
    int /*layerNum,*/i;
    int nFeatures, pc=0;
    OGRLayer *poLayer;
    OGRFeature *poFeature;
#ifdef GDALV2
    GDALDataset *poDS;
#else
    OGRDataSource *poDS;
    OGRSFDriver *poDriver;
#endif

    installErrorHandler();
#ifdef GDALV2
    poDS=(GDALDataset*) GDALOpenEx(CHAR(STRING_ELT(filename, 0)), GDAL_OF_VECTOR, NULL, NULL, NULL);
#else
    poDS=OGRSFDriverRegistrar::Open(CHAR(STRING_ELT(filename, 0)),
                                    FALSE, &poDriver);
#endif
    uninstallErrorHandlerAndTriggerError();

    if(poDS==NULL) {
        error("Cannot open file");
    }

    installErrorHandler();
    poLayer = poDS->GetLayerByName(CHAR(STRING_ELT(layer, 0)));
    uninstallErrorHandlerAndTriggerError();

    if(poLayer == NULL) {
        error("Cannot open layer");
    }
    installErrorHandler();
#ifdef GDALV2
    GIntBig nFIDs64 = poLayer->GetFeatureCount();
    nFeatures = (nFIDs64 > INT_MAX) ? INT_MAX :
                (nFIDs64 < INT_MIN) ? INT_MIN : (int) nFIDs64;
    if ((GIntBig) nFeatures != nFIDs64) {
        uninstallErrorHandlerAndTriggerError();
        error("ogrFIDs: feature count overflow");
    }
#else
    nFeatures=poLayer->GetFeatureCount();
#endif
    uninstallErrorHandlerAndTriggerError();

    if (nFeatures == -1) {
        i=0;
        installErrorHandler();
        while( ((poFeature = poLayer->GetNextFeature()) != NULL) && i <= INT_MAX) {
            i++;
            OGRFeature::DestroyFeature( poFeature );
//    delete poFeature;
        }
        uninstallErrorHandlerAndTriggerError();
        installErrorHandler();
        poLayer->ResetReading();
        uninstallErrorHandlerAndTriggerError();
        if (i == INT_MAX) {
            error("ogrFIDs: feature count overflow");
        } else {
            nFeatures = i;
        }
    }

    PROTECT(fids=allocVector(INTSXP,nFeatures));
    pc++;
    PROTECT(nf = NEW_INTEGER(1));
    pc++;
    INTEGER_POINTER(nf)[0] = nFeatures;
    PROTECT(ii = NEW_INTEGER(1));
    pc++;

    installErrorHandler();
    poLayer->ResetReading();
    uninstallErrorHandlerAndTriggerError();


    i=0;
    installErrorHandler();
    while( (poFeature = poLayer->GetNextFeature()) != NULL ) {
        INTEGER(fids)[i]= (int) poFeature->GetFID();
        i++;
        OGRFeature::DestroyFeature( poFeature );
//    delete poFeature;
    }
    uninstallErrorHandlerAndTriggerError();

    INTEGER_POINTER(ii)[0] = i;
    setAttrib(fids, install("nf"), nf);
    setAttrib(fids, install("i"), ii);

    installErrorHandler();
#ifdef GDALV2
    GDALClose( poDS );
#else
    OGRDataSource::DestroyDataSource( poDS );
#endif
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            return CE_Failure;
        }
    }

/* ==================================================================== */
/*      Read the specified layers transfoming and rasterizing           */
/*      geometries.                                                     */
/* ==================================================================== */
    CPLErr      eErr = CE_None;
    int         iLayer;
    const char  *pszBurnAttribute =
        CSLFetchNameValue( papszOptions, "ATTRIBUTE" );

    pfnProgress( 0.0, NULL, pProgressArg );

    for( iLayer = 0; iLayer < nLayerCount; iLayer++ )
    {
        int         iBurnField = -1;
        double      *padfBurnValues = NULL;
        OGRLayer    *poLayer = (OGRLayer *) pahLayers[iLayer];

        if ( !poLayer )
        {
            CPLError( CE_Warning, CPLE_AppDefined, 
                      "Layer element number %d is NULL, skipping.\n", iLayer );
            continue;
        }

/* -------------------------------------------------------------------- */
/*      If the layer does not contain any features just skip it.        */
/*      Do not force the feature count, so if driver doesn't know       */
/*      exact number of features, go down the normal way.               */
/* -------------------------------------------------------------------- */
        if ( poLayer->GetFeatureCount(FALSE) == 0 )
            continue;

        if ( pszBurnAttribute )
        {
            iBurnField =
                poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute );
            if ( iBurnField == -1 )
            {
                CPLError( CE_Warning, CPLE_AppDefined, 
                          "Failed to find field %s on layer %s, skipping.\n",
                          pszBurnAttribute, 
                          poLayer->GetLayerDefn()->GetName() );
                continue;
            }
        }
        else
            padfBurnValues = padfLayerBurnValues + iLayer * nBandCount;

/* -------------------------------------------------------------------- */
/*      If we have no transformer, create the one from input file       */
/*      projection. Note that each layer can be georefernced            */
/*      separately.                                                     */
/* -------------------------------------------------------------------- */
        int bNeedToFreeTransformer = FALSE;

        if( pfnTransformer == NULL )
        {
            char    *pszProjection = NULL;
            bNeedToFreeTransformer = TRUE;

            OGRSpatialReference *poSRS = poLayer->GetSpatialRef();
            if ( !poSRS )

//.........这里部分代码省略.........
        std::cerr << "Error: geometry must be Polygon or MultiPolygon" << std::endl;
        return 1;
      }
    }
    
    //-- reading from WKT stored in first line of a text file
    else if (strcmp(argv[argNum], "-f") == 0) {
      unsigned int bufferSize = 100000000;
      char *inputWKT = (char *)malloc(bufferSize*sizeof(char));
      if (argNum + 1 <= argc - 1 && argv[argNum+1][0] != '-') {
        std::ifstream infile(argv[argNum+1], std::ifstream::in);
        infile.getline(inputWKT, bufferSize);
        ++argNum;
      } else {
        std::cerr << "Error: Missing input file name." << std::endl;
        return 1;
      }
      OGRGeometryFactory::createFromWkt(&inputWKT, NULL, &geometry);
      if (geometry == NULL) {
        std::cout << "Error: WKT is not valid" << std::endl;
        return 1;
      }
    }
    
    //-- reading from a ogr dataset (most supported: shp, geojson, gml, etc)
    else if (strcmp(argv[argNum], "--ogr") == 0) {
      OGRRegisterAll();
      OGRDataSource *dataSource = OGRSFDriverRegistrar::Open(argv[argNum+1], false);
      ++argNum;
      if (dataSource == NULL) {
        std::cerr << "Error: Could not open file." << std::endl;
        return false;
      }
      OGRLayer *dataLayer = dataSource->GetLayer(0); //-- get first layer
      dataLayer->ResetReading();
      //-- Reads all features in this layer
      OGRFeature *feature;
      //-- get first feature
      if (dataLayer->GetFeatureCount(true) > 1)
        std::cout << "Reading only the first feature in the shapefile." << std::endl;
      feature = dataLayer->GetNextFeature();
      if (feature->GetGeometryRef()->getGeometryType() == wkbPolygon) {
        geometry = static_cast<OGRPolygon *>(feature->GetGeometryRef());feature->GetGeometryRef();
      }
      else {
        std::cout << "First feature ain't a POLYGON." << std::endl;
        return(0);
      }
      
    }
    
    else {
      usage();
      return(0);
    }
  }
  
  if (!timeResults) startTime = 0;
  
  PolygonRepair prepair;
  
  //-- compute robustness
  if (computeRobustness == true)
    std::cout << "Robustness of input polygon: " << sqrt(prepair.computeRobustness(geometry)) <<std::endl;
  
  
  OGRMultiPolygon *outPolygons;
  if (pointSet) {
    outPolygons = prepair.repairPointSet(geometry, startTime);
  } else {
    outPolygons = prepair.repairOddEven(geometry, startTime);
  }
  
  if (minArea > 0) {
    prepair.removeSmallPolygons(outPolygons, minArea);
  }
  
  //-- output well known text
  if (shpOut) {
    prepair.saveToShp(outPolygons, "out.shp");
  }
  else {
    char *outputWKT;
    outPolygons->exportToWkt(&outputWKT);
    std::cout << outputWKT << std::endl;
  }
  
  
  //-- compute robustness
  if (computeRobustness == true)
    std::cout << "Robustness of output polygon: " << sqrt(prepair.computeRobustness()) <<std::endl;
  
  //-- time results
  if (timeResults) {
    time_t totalTime = time(NULL)-startTime;
    std::cout << "Done! Process finished in " << totalTime/60 << " minutes " << totalTime%60 << " seconds." << std::endl;
  }
  
  return 0;
}

CPLErr GDALRasterizeLayersBuf( void *pData, int nBufXSize, int nBufYSize,
                               GDALDataType eBufType,
                               int nPixelSpace, int nLineSpace,
                               int nLayerCount, OGRLayerH *pahLayers,
                               const char *pszDstProjection,
                               double *padfDstGeoTransform,
                               GDALTransformerFunc pfnTransformer, 
                               void *pTransformArg, double dfBurnValue,
                               char **papszOptions,
                               GDALProgressFunc pfnProgress,
                               void *pProgressArg )

{
/* -------------------------------------------------------------------- */
/*      If pixel and line spaceing are defaulted assign reasonable      */
/*      value assuming a packed buffer.                                 */
/* -------------------------------------------------------------------- */
    if( nPixelSpace == 0 )
        nPixelSpace = GDALGetDataTypeSize( eBufType ) / 8;
    
    if( nLineSpace == 0 )
        nLineSpace = nPixelSpace * nBufXSize;

    if( pfnProgress == NULL )
        pfnProgress = GDALDummyProgress;

/* -------------------------------------------------------------------- */
/*      Do some rudimentary arg checking.                               */
/* -------------------------------------------------------------------- */
    if( nLayerCount == 0 )
        return CE_None;

    int bAllTouched = CSLFetchBoolean( papszOptions, "ALL_TOUCHED", FALSE );
    const char *pszOpt = CSLFetchNameValue( papszOptions, "BURN_VALUE_FROM" );
    GDALBurnValueSrc eBurnValueSource = GBV_UserBurnValue;
    if( pszOpt )
    {
        if( EQUAL(pszOpt,"Z"))
            eBurnValueSource = GBV_Z;
        /*else if( EQUAL(pszOpt,"M"))
            eBurnValueSource = GBV_M;*/
    }

/* ==================================================================== */
/*      Read thes pecified layers transfoming and rasterizing           */
/*      geometries.                                                     */
/* ==================================================================== */
    CPLErr      eErr = CE_None;
    int         iLayer;
    const char  *pszBurnAttribute =
        CSLFetchNameValue( papszOptions, "ATTRIBUTE" );

    pfnProgress( 0.0, NULL, pProgressArg );

    for( iLayer = 0; iLayer < nLayerCount; iLayer++ )
    {
        int         iBurnField = -1;
        OGRLayer    *poLayer = (OGRLayer *) pahLayers[iLayer];

        if ( !poLayer )
        {
            CPLError( CE_Warning, CPLE_AppDefined, 
                      "Layer element number %d is NULL, skipping.\n", iLayer );
            continue;
        }

/* -------------------------------------------------------------------- */
/*      If the layer does not contain any features just skip it.        */
/*      Do not force the feature count, so if driver doesn't know       */
/*      exact number of features, go down the normal way.               */
/* -------------------------------------------------------------------- */
        if ( poLayer->GetFeatureCount(FALSE) == 0 )
            continue;

        if ( pszBurnAttribute )
        {
            iBurnField =
                poLayer->GetLayerDefn()->GetFieldIndex( pszBurnAttribute );
            if ( iBurnField == -1 )
            {
                CPLError( CE_Warning, CPLE_AppDefined, 
                          "Failed to find field %s on layer %s, skipping.\n",
                          pszBurnAttribute, 
                          poLayer->GetLayerDefn()->GetName() );
                continue;
            }
        }

/* -------------------------------------------------------------------- */
/*      If we have no transformer, create the one from input file       */
/*      projection. Note that each layer can be georefernced            */
/*      separately.                                                     */
/* -------------------------------------------------------------------- */
        int bNeedToFreeTransformer = FALSE;

        if( pfnTransformer == NULL )
        {
            char    *pszProjection = NULL;
            bNeedToFreeTransformer = TRUE;

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

int main()
{
    GDALAllRegister();  //register all the format drivers
    cout << "GDAL All Registed!" << endl;

    GDALDataset *poDS;  //Data source
    poDS = (GDALDataset*) GDALOpenEx("./beijing/road/Nbeijing_point.shp", GDAL_OF_VECTOR, NULL, NULL, NULL);
    if(poDS == NULL)
    {
        cout << "Open shp file failed!" << endl;
        exit(1);
    }
    cout << "Data source open success!" << endl;

    int layerNum = poDS->GetLayerCount();   //a dataset may have many layers
    cout << "Layer number:" << layerNum << endl;

    OGRLayer *poLayer;
    poLayer = poDS->GetLayerByName("Nbeijing_point");
    
    //feature is a geometry and a set of attributes
    OGRFeature *poFeature;
    poLayer->ResetReading();    //Start at the beginning of the layer
    cout << "Feature number:" << poLayer->GetFeatureCount() << endl;

    stringstream ss;
    map<string, pair<double, double> > mip;
    map<string, pair<double, double> >::iterator imip;
    
	ofstream ofile("beijingNode");

    string id;
	int	crossFlag;
	string stmp, stmp2;
	vector<string> vs;
	vector<string>::iterator ivs;
    while((poFeature = poLayer->GetNextFeature()) != NULL)
    {
        //Defn contains the definition of all the fields
        OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
//        for(iField = 0; iField < poFDefn->GetFieldCount(); iField++)
 //       {
//        OGRFieldDefn * poFieldDefn = poFDefn->GetFieldDefn(iField);
        id = poFeature->GetFieldAsString(1);

		ofile << id;

		crossFlag = poFeature->GetFieldAsInteger(4);
		if(crossFlag == 0)
		{
			ofile << "\t" << crossFlag << "\t" << 0 << "\t" << 0 << "\t" << 0;
		}
		else if(crossFlag == 1)
		{	
			ofile << "\t" << crossFlag << "\t" << 0 << "\t" << poFeature->GetFieldAsInteger(7) << "\t" << 0;
		}
		else if(crossFlag == 2)
		{
			stmp = poFeature->GetFieldAsString(6);
			vs = split(stmp, "|");
			ofile << "\t" << crossFlag << "\t" << vs.size();
			for(ivs = vs.begin(); ivs != vs.end(); ivs++)
				ofile << "\t" << *ivs;
			vs.clear();
			ofile << "\t" << poFeature->GetFieldAsInteger(7) << "\t" << 0;
		}
		else if(crossFlag == 3)
		{
			stmp = poFeature->GetFieldAsString(6);
			vs = split(stmp, "|");
			ofile << "\t" << crossFlag << "\t" << vs.size();
			for(ivs = vs.begin(); ivs != vs.end(); ivs++)
				ofile << "\t" << *ivs;
			vs.clear();
			
			ofile << "\t" << poFeature->GetFieldAsInteger(7);
			stmp = poFeature->GetFieldAsString(8);
			stmp2 = poFeature->GetFieldAsString(9);
			if(stmp2 != "")
				stmp += "|" + stmp2;
			vs = split(stmp, "|");
			ofile << "\t" << vs.size();
			for(ivs = vs.begin(); ivs != vs.end(); ivs++)
				ofile << "\t" << *ivs;
			vs.clear();
		}
	
		ofile << "\t" << poFeature->GetFieldAsString(11);

		stmp = poFeature->GetFieldAsString(12);
		vs = split(stmp, "|");
		ofile << "\t" << vs.size();
		for(ivs = vs.begin(); ivs != vs.end(); ivs++)
			ofile << "\t" << *ivs;
		vs.clear();
/*            if(poFieldDefn->GetType() == OFTInteger)
                cout << poFeature->GetFieldAsInteger(iField) << ", ";
            else if(poFieldDefn->GetType() == OFTInteger64)
                cout << poFeature->GetFieldAsInteger64(iField) << ", ";
            else if(poFieldDefn->GetType() == OFTReal)
//.........这里部分代码省略.........

// extern "C" {
SEXP ogrInfo(SEXP ogrsourcename, SEXP Layer) {
    // return FIDs, nFields, fieldInfo

    SEXP ans, vec1, vec2, vec3,/*mat,*/drv, dvec;
    SEXP itemlist, itemnames, itemwidth, itemtype, itemTypeNames;
    SEXP itemlistmaxcount;
#ifdef GDALV2
    SEXP dFIDs;
#endif
    /*SEXP geotype;*/

    int nFIDs, nFields, iField, *nCount, pc=0;

#ifdef GDALV2
    GDALDriver *poDriver;
    GDALDataset *poDS;
#else
    OGRDataSource *poDS;
    OGRSFDriver *poDriver;
#endif
    OGRLayer *poLayer;
    OGRFeature *poFeature;
    OGRFeatureDefn *poDefn;
    /*  OGRGeometry *poGeom;*/

    installErrorHandler();
#ifdef GDALV2
    poDS=(GDALDataset*) GDALOpenEx(CHAR(STRING_ELT(ogrsourcename, 0)), GDAL_OF_VECTOR, NULL, NULL, NULL);
    if(poDS==NULL) {
        uninstallErrorHandlerAndTriggerError();
        error("Cannot open data source");
    }
    poDriver = poDS->GetDriver();
#else
    poDS=OGRSFDriverRegistrar::Open(CHAR(STRING_ELT(ogrsourcename, 0)),
                                    FALSE, &poDriver);
#endif
    uninstallErrorHandlerAndTriggerError();

    if(poDS==NULL) {
        installErrorHandler();
#ifdef GDALV2
        GDALClose( poDS );
#else
        OGRDataSource::DestroyDataSource( poDS );
#endif
        uninstallErrorHandlerAndTriggerError();
//    delete poDS;
        error("Cannot open file");
    }

    installErrorHandler();
    poLayer = poDS->GetLayerByName(CHAR(STRING_ELT(Layer, 0)));
    uninstallErrorHandlerAndTriggerError();

    if(poLayer == NULL) {
        installErrorHandler();
#ifdef GDALV2
        GDALClose( poDS );
#else
        OGRDataSource::DestroyDataSource( poDS );
#endif
        uninstallErrorHandlerAndTriggerError();
//    delete poDS;
        error("Cannot open layer");
    }

    // allocate a list for return values
    PROTECT(ans=allocVector(VECSXP,6));
    pc++;

    PROTECT(drv=allocVector(STRSXP,1));
    pc++;
    installErrorHandler();
#ifdef GDALV2
    SET_STRING_ELT(drv, 0, mkChar(poDriver->GetDescription()));
#else
    SET_STRING_ELT(drv, 0, mkChar(poDriver->GetName()));
#endif
    uninstallErrorHandlerAndTriggerError();
    SET_VECTOR_ELT(ans,3,drv);

    PROTECT(vec1=allocVector(INTSXP,1));
    pc++;
    installErrorHandler();
#ifdef GDALV2
    GIntBig nFIDs64 = poLayer->GetFeatureCount();
    nFIDs = (nFIDs64 > INT_MAX) ? INT_MAX :
            (nFIDs64 < INT_MIN) ? INT_MIN : (int) nFIDs64;
    if ((GIntBig) nFIDs != nFIDs64) {
        warning("ogrInfo: feature count overflow");
        INTEGER(vec1)[0]=NA_INTEGER;
        PROTECT(dFIDs=NEW_NUMERIC(1));
        pc++;
        NUMERIC_POINTER(dFIDs)[0] = (double) nFIDs64;
        setAttrib(vec1, install("dFIDs"), dFIDs);
    } else {
        // store number of FIDs
        INTEGER(vec1)[0]=nFIDs;
//.........这里部分代码省略.........

void ShpReader::VectorMetaDataInfo(OGRDataSource* OGRDataset, StudyControllerPtr studyController, VectorMapControllerPtr vectorMapController)
{	
	vFileMetaData* vHeader = vectorMapController->GetVectorMapModel()->GetVectorMetaData();
    OGRLayer *poLayer = OGRDataset->GetLayer( 0 );	
	char *originalWkt = NULL;	
	char *destinationWkt = NULL;

	OGREnvelope psEnvelope;
	poLayer->GetExtent( &psEnvelope );
	vHeader->originalExtents.x = psEnvelope.MinX;
	vHeader->originalExtents.dx = psEnvelope.MaxX;
	vHeader->originalExtents.y = psEnvelope.MinY; 
	vHeader->originalExtents.dy = psEnvelope.MaxY;
	Log::Inst().Write("Original Extents: ");
	Log::Inst().Write(_T("Lower, Left (x,y): ") +wxString(StringTools::ToStringW(vHeader->originalExtents.x, 2).c_str()) + _T(", ") 
		+ wxString(StringTools::ToStringW(vHeader->originalExtents.y, 2).c_str()));	
	Log::Inst().Write(_T("Upper, Right (x,y): ") +wxString(StringTools::ToStringW(vHeader->originalExtents.dx, 2).c_str()) + _T(", ") 
		+ wxString(StringTools::ToStringW(vHeader->originalExtents.dy, 2).c_str()));		

	vHeader->numFeatures= poLayer->GetFeatureCount();
	Log::Inst().Write("Number of Features: " + StringTools::ToString(vHeader->numFeatures));
	
	std::string type;
	OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn();
	if ( wkbFlatten( poFDefn->GetGeomType() ) == wkbPoint || wkbFlatten( poFDefn->GetGeomType() ) == wkbMultiPoint )
		 type="Point";
	else if ( wkbFlatten(poFDefn->GetGeomType() ) == wkbLineString || wkbFlatten(poFDefn->GetGeomType() ) == wkbMultiLineString )
		type="Polyline";
	else if ( wkbFlatten( poFDefn->GetGeomType() ) == wkbPolygon || wkbFlatten( poFDefn->GetGeomType() ) == wkbMultiPolygon )
		 type="Polygon";  	
	else
		type="Unknown";

	 vHeader->geometryType=type;
	 Log::Inst().Write("Geometry Type: "+ type);

	OGRSpatialReference* originalShpSR = poLayer->GetSpatialRef();	
	if (originalShpSR != NULL)
	{
		originalShpSR->exportToWkt(&originalWkt);	
		vHeader->originalProjection = std::string(originalWkt);
	}
	Log::Inst().Write("Original Projection: " + vHeader->originalProjection);
	Log::Inst().Write("");
	ProjectionToolPtr projTool = studyController->GetProjectionTool();
	if (projTool != NULL)
	{
		OGRSpatialReference *destinationCS= projTool->GetTargetCS ();
		destinationCS->exportToWkt(&destinationWkt);
		vHeader->destinationProjection = std::string(destinationWkt);
	}

	vHeader->currentExtents.x = vectorMapController->GetVectorMapModel()->GetVectorBoundary_MinX();
	vHeader->currentExtents.y = vectorMapController->GetVectorMapModel()->GetVectorBoundary_MinY();	

	vHeader->currentExtents.dx = vectorMapController->GetVectorMapModel()->GetVectorBoundary_MaxX();
	vHeader->currentExtents.dy = vectorMapController->GetVectorMapModel()->GetVectorBoundary_MaxY();
	//Log::Inst().Write("Upper, Right (x,y): " +
	//	StringTools::ToString(vHeader->currentExtents.dx) + ", " +
		//StringTools::ToString(vHeader->currentExtents.dy));
	
}

void DataCvt_Log::CalAsShp()
{
	double ptmptime_start = MPI_Wtime();
	OGRRegisterAll();
	//create a new shp to caculate
	CBaseOperate pbaseOperate;
	if(rankid==0)
	{
	     pbaseOperate.CreateCopyShp(m_DatasourceConStr.c_str(),m_PathOrTableName.c_str(),m_ResultPathOrTableName.c_str());
	}
	MPI_Barrier(MPI_COMM_WORLD);
	string pLayerName = m_ResultPathOrTableName;
	if(strcmp(m_DatasourceConStr.c_str(),"")==0)//input is shp file
	{
		m_DatasourceConStr = m_ResultPathOrTableName;
		string pShpPath = m_ResultPathOrTableName;
		string p_LayerName = GetFileNameOnly(pShpPath.c_str());
		int pindex = p_LayerName.find_first_of('.');
		pLayerName = p_LayerName.erase(pindex,p_LayerName.size()-1);
	}

	OGRDataSource* poDS = OGRSFDriverRegistrar::Open(m_DatasourceConStr.c_str(),TRUE); 
	OGRLayer* poLayer = poDS->GetLayerByName(pLayerName.c_str());
	OGRFeature * pFeature =poLayer->GetNextFeature();

	m_RowNum = poLayer->GetFeatureCount();
	//if -c argv is null, all field will be used
	m_ColNum= pFeature->GetFieldCount();
	if(c_CaculateCols.size()==0)
	{
		if(rankid==0)
			cout<<"!!!!!!!  No input columns ids, all fields in the input file will be used..."<<endl<<endl;
		for(int i=0;i<m_ColNum;i++)
		{
			c_CaculateCols.push_back(i);
		}
	}

	double ptmptime_end = MPI_Wtime();
	m_GdalIOInCal_Time=m_GdalIOInCal_Time+ptmptime_end-ptmptime_start;

	//assign number to each process
	int pnum_start,pnum_end;
	int pMyProcessNum=0; int pRemainder=m_RowNum%numproc;
	if(rankid<pRemainder)
	{
		pMyProcessNum = m_RowNum/numproc+1;
		pnum_start = rankid*pMyProcessNum;
		pnum_end = pnum_start+pMyProcessNum-1;
	}
	else
	{
		pMyProcessNum = m_RowNum/numproc;
		pnum_start = pRemainder*(pMyProcessNum+1)+(rankid-pRemainder)*pMyProcessNum;
		pnum_end = pnum_start+pMyProcessNum-1;
	}

	//postgis: fid begins from 1, not 0
	string pwhere="";
	if(strcmp(m_DatasourceConStr.substr(0,3).c_str(),"PG:")==0)//input is postgis
	{
		pwhere = "gid>="+toString(pnum_start+1)+" and gid<="+toString(pnum_end+1);
	}
	else//shpfile: fid begins from 0, not 1
	{
		pwhere = "fid>="+toString(pnum_start)+" and fid<="+toString(pnum_end);
	}

	poLayer->SetAttributeFilter(pwhere.c_str());
	pFeature = poLayer->GetNextFeature();
	while(pFeature!=NULL)
	{		
		for(int i=0;i<c_CaculateCols.size();i++)
		{			
			double tmp = pFeature->GetFieldAsDouble(c_CaculateCols[i]);
			if(tmp!=0)
				tmp = log(fabs(tmp));
			else
				tmp=0;
			double pstart = MPI_Wtime();
			pFeature->SetField(c_CaculateCols[i],tmp);
			poLayer->SetFeature(pFeature);
			double pend = MPI_Wtime();
			m_GdalIOInCal_Time = m_GdalIOInCal_Time+pend-pstart;
		}		
		pFeature = poLayer->GetNextFeature();
	}
	OGRDataSource::DestroyDataSource( poDS );
	double ptmptime_endall = MPI_Wtime();
	Mpi_Caculate_Time = ptmptime_endall-ptmptime_start-m_GdalIOInCal_Time;

	double tmp=Mpi_Caculate_Time;
	MPI_Reduce(&Mpi_Caculate_Time,&tmp,1,MPI_DOUBLE,MPI_MIN,0,MPI_COMM_WORLD);
	Mpi_Caculate_Time = tmp;
}