C++ OGRDataSource::GetLayerByName方法代码示例

bool V2vProj::Compute(const data::VectorBarral * barrel)
{
    OGRDataSource * poSourceDs = VectorOpen(barrel->GetSrcDataSource().c_str(),
                                            GA_ReadOnly);
    ON_SCOPE_EXIT([&]() {OGRDataSource::DestroyDataSource(poSourceDs); });
    OGRDataSource * poOutputDs = VectorOpen(barrel->GetDstDataSource().c_str(),
                                            GA_Update);
    ON_SCOPE_EXIT([&]() {OGRDataSource::DestroyDataSource(poOutputDs); });
    OGRLayer * poSrcLayer = poSourceDs->GetLayerByName(
                                barrel->GetSrcLayer().c_str());
    OGRLayer * poDstLayer = poOutputDs->GetLayerByName(
                                barrel->GetDstLayer().c_str());
    OGRSpatialReference * poSourceSRS = poSrcLayer->GetSpatialRef();
    OGRCoordinateTransformation * poCT = poCT = OGRCreateCoordinateTransformation(
            poSourceSRS, m_ogrSr);
    OGRFeatureDefn * poDstFeatureDefn = poDstLayer->GetLayerDefn();
    auto features = barrel->GetFeatures();
    std::for_each(begin(features), end(features)
    , [&](int fid) {
        poSrcLayer->GetFeature(fid);
        OGRFeature * poDstFeature = OGRFeature::CreateFeature(poDstFeatureDefn);
        ON_SCOPE_EXIT([&]() {OGRFeature::DestroyFeature(poDstFeature); });
        poDstFeature->SetFrom(poSrcLayer->GetFeature(fid));
        OGRGeometry * poDstGeometry = poDstFeature->GetGeometryRef();
        OGRGeometry * poReprojectedGeom = OGRGeometryFactory::transformWithOptions(
                                              poDstGeometry, poCT, NULL);
        poDstFeature->SetGeometryDirectly(poReprojectedGeom);
        poDstLayer->CreateFeature(poDstFeature);
    });
    return true;
}

void OGRDataSourceWithTransaction::RemapLayers()
{
    std::set<OGRLayerWithTransaction*>::iterator oIter = m_oSetLayers.begin();
    for(; oIter != m_oSetLayers.end(); ++oIter )
    {
        OGRLayerWithTransaction* poWrappedLayer = *oIter;
        if( m_poBaseDataSource == NULL )
            poWrappedLayer->m_poDecoratedLayer = NULL;
        else
        {
            poWrappedLayer->m_poDecoratedLayer =
                m_poBaseDataSource->GetLayerByName(poWrappedLayer->GetName());
        }
    }
    m_oMapLayers.clear();
}

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

JNIEXPORT jlong JNICALL Java_es_gva_cit_jogr_OGRDataSource_getLayerByNameNat
  (JNIEnv *env, jobject obj, jlong cPtr, jstring name){
  	  	
  	OGRDataSource 	*ds = (OGRDataSource *) 0 ;
  	const char 		*layername;
  	OGRLayer		*layer;
  	long			ptro_layer=-1;
  	
  	ds = *(OGRDataSource **)&cPtr;
  	layername = env->GetStringUTFChars( name, 0 );
  	if(ds!=NULL){
	  	layer=ds->GetLayerByName(layername);
  		if(layer!=NULL)ptro_layer = (long)&(*layer);
  	}
  	env->ReleaseStringUTFChars( name, layername );
  	return (jlong)ptro_layer;
  }

    MyOGRHandler(const std::string& filename) :
        m_data_source(NULL),
        m_layer_point(NULL),
        m_filter(true),
        m_tohstore() {
        OGRRegisterAll();

        OGRSFDriver* driver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName("PostgreSQL");
        if (driver == NULL) {
            std::cerr << "PostgreSQL OGR driver not available.\n";
            exit(1);
        }

        // using COPY is much faster than INSERT
        CPLSetConfigOption("PG_USE_COPY", "YES");
        const char* options[] = { NULL };
        m_data_source = driver->CreateDataSource(filename.c_str(), const_cast<char**>(options));
        if (m_data_source == NULL) {
            std::cerr << "Database open failed.\n";
            exit(1);
        }

        // OGR can't create a table with hstore column, so we do it ourselves here
        OGRLayer* dummy = m_data_source->ExecuteSQL("CREATE TABLE nodes (id VARCHAR, tags hstore);", NULL, NULL);
        if (dummy) {
            m_data_source->ReleaseResultSet(dummy);
        }
        dummy = m_data_source->ExecuteSQL("SELECT AddGeometryColumn('nodes', 'geom', 4326, 'POINT', 2);", NULL, NULL);
        if (dummy) {
            m_data_source->ReleaseResultSet(dummy);
        }

        m_layer_point = m_data_source->GetLayerByName("nodes");
        if (!m_layer_point) {
            std::cerr << "Something went wrong setting up the 'nodes' layer.\n";
            exit(1);
        }

        // using transactions makes this much faster than without
        m_layer_point->StartTransaction();

        m_filter.add(false, "created_by");
        m_filter.add(false, "odbl");
    }

OGRLayer    *OGRDataSourceWithTransaction::GetLayerByName(const char *pszName)
{
    if( !m_poBaseDataSource ) return NULL;
    return WrapLayer(m_poBaseDataSource->GetLayerByName(pszName));
}

OGRLayer * OGRDataSource::ExecuteSQL( const char *pszStatement,
                                      OGRGeometry *poSpatialFilter,
                                      const char *pszDialect )

{
    const char *pszError;
    swq_select *psSelectInfo = NULL;

    (void) pszDialect;

/* -------------------------------------------------------------------- */
/*      Handle CREATE INDEX statements specially.                       */
/* -------------------------------------------------------------------- */
    if( EQUALN(pszStatement,"CREATE INDEX",12) )
    {
        ProcessSQLCreateIndex( pszStatement );
        return NULL;
    }
    
/* -------------------------------------------------------------------- */
/*      Handle DROP INDEX statements specially.                         */
/* -------------------------------------------------------------------- */
    if( EQUALN(pszStatement,"DROP INDEX",10) )
    {
        ProcessSQLDropIndex( pszStatement );
        return NULL;
    }
    
/* -------------------------------------------------------------------- */
/*      Preparse the SQL statement.                                     */
/* -------------------------------------------------------------------- */
    pszError = swq_select_preparse( pszStatement, &psSelectInfo );
    if( pszError != NULL )
    {
        CPLError( CE_Failure, CPLE_AppDefined, 
                  "SQL: %s", pszError );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Validate that all the source tables are recognised, count       */
/*      fields.                                                         */
/* -------------------------------------------------------------------- */
    int  nFieldCount = 0, iTable;

    for( iTable = 0; iTable < psSelectInfo->table_count; iTable++ )
    {
        swq_table_def *psTableDef = psSelectInfo->table_defs + iTable;
        OGRLayer *poSrcLayer;
        OGRDataSource *poTableDS = this;

        if( psTableDef->data_source != NULL )
        {
            poTableDS = (OGRDataSource *) 
                OGROpenShared( psTableDef->data_source, FALSE, NULL );
            if( poTableDS == NULL )
            {
                if( strlen(CPLGetLastErrorMsg()) == 0 )
                    CPLError( CE_Failure, CPLE_AppDefined, 
                              "Unable to open secondary datasource\n"
                              "`%s' required by JOIN.",
                              psTableDef->data_source );

                swq_select_free( psSelectInfo );
                return NULL;
            }

            // This drops explicit reference, but leave it open for use by
            // code in ogr_gensql.cpp
            poTableDS->Dereference();
        }

        poSrcLayer = poTableDS->GetLayerByName( psTableDef->table_name );

        if( poSrcLayer == NULL )
        {
            CPLError( CE_Failure, CPLE_AppDefined, 
                      "SELECT from table %s failed, no such table/featureclass.",
                      psTableDef->table_name );
            swq_select_free( psSelectInfo );
            return NULL;
        }

        nFieldCount += poSrcLayer->GetLayerDefn()->GetFieldCount();
    }
    
/* -------------------------------------------------------------------- */
/*      Build the field list for all indicated tables.                  */
/* -------------------------------------------------------------------- */
    swq_field_list sFieldList;
    int            nFIDIndex = 0;

    memset( &sFieldList, 0, sizeof(sFieldList) );
    sFieldList.table_count = psSelectInfo->table_count;
    sFieldList.table_defs = psSelectInfo->table_defs;

    sFieldList.count = 0;
    sFieldList.names = (char **) CPLMalloc( sizeof(char *) * (nFieldCount+1) );
    sFieldList.types = (swq_field_type *)  
        CPLMalloc( sizeof(swq_field_type) * (nFieldCount+1) );
//.........这里部分代码省略.........

//.........这里部分代码省略.........
/*      Analysze the statement to determine which tables will be used.  */
/* -------------------------------------------------------------------- */
    std::set<LayerDesc> oSetLayers;
    std::set<CPLString> oSetSpatialIndex;
    CPLString osModifiedSQL;
    OGR2SQLITEGetPotentialLayerNames(pszStatement, oSetLayers,
                                     oSetSpatialIndex, osModifiedSQL);
    std::set<LayerDesc>::iterator oIter = oSetLayers.begin();

    if( strcmp(pszStatement, osModifiedSQL.c_str()) != 0 )
        CPLDebug("OGR", "Modified SQL: %s", osModifiedSQL.c_str());
    pszStatement = osModifiedSQL.c_str(); /* do not use it anymore */

    int bFoundOGRStyle = ( osModifiedSQL.ifind("OGR_STYLE") != std::string::npos );

/* -------------------------------------------------------------------- */
/*      For each of those tables, create a Virtual Table.               */
/* -------------------------------------------------------------------- */
    for(; oIter != oSetLayers.end(); ++oIter)
    {
        const LayerDesc& oLayerDesc = *oIter;
        /*CPLDebug("OGR", "Layer desc : %s, %s, %s, %s",
                 oLayerDesc.osOriginalStr.c_str(),
                 oLayerDesc.osSubstitutedName.c_str(),
                 oLayerDesc.osDSName.c_str(),
                 oLayerDesc.osLayerName.c_str());*/

        CPLString osSQL;
        OGRLayer* poLayer = NULL;
        CPLString osTableName;
        int nExtraDS;
        if( oLayerDesc.osDSName.size() == 0 )
        {
            poLayer = poDS->GetLayerByName(oLayerDesc.osLayerName);
            /* Might be a false positive (unlikely) */
            if( poLayer == NULL )
                continue;

            osTableName = oLayerDesc.osLayerName;

            nExtraDS = -1;

            osSQL.Printf("CREATE VIRTUAL TABLE \"%s\" USING VirtualOGR(%d,'%s',%d)",
                         OGRSQLiteEscapeName(osTableName).c_str(),
                         nExtraDS,
                         OGRSQLiteEscape(osTableName).c_str(),
                         bFoundOGRStyle);
        }
        else
        {
            OGRDataSource* poOtherDS = (OGRDataSource* )
                OGROpen(oLayerDesc.osDSName, FALSE, NULL);
            if( poOtherDS == NULL )
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Cannot open datasource '%s'",
                         oLayerDesc.osDSName.c_str() );
                delete poSQLiteDS;
                VSIUnlink(pszTmpDBName);
                CPLFree(pszTmpDBName);
                return NULL;
            }
            
            poLayer = poOtherDS->GetLayerByName(oLayerDesc.osLayerName);
            if( poLayer == NULL )
            {

void lasclip(std::string &outfile, std::string &shapefile,
		std::string &layername, std::vector<std::string> &files,
		std::set<int> &classes, bool quiet) {

	if (outfile.empty())
		g_argerr("An output file is required.");
	if (shapefile.empty())
		g_argerr("A shape file is required.");
	if (files.size() == 0)
		g_argerr("At least one input file is required.");
	if (classes.size() == 0)
		g_warn("No classes specified, matching all classes.");

	/* Attempt to open and load geometries from the shape file. */
	OGRRegisterAll();
	OGRLayer *layer;
	OGRFeature *feat;
	OGRGeometry *og;
	OGRwkbGeometryType type;
	gg::GeometryCollection *geomColl;
	gg::Geometry *geom;

	OGRDataSource *ds = OGRSFDriverRegistrar::Open(shapefile.c_str(), FALSE);
	if (ds == nullptr)
		g_runerr("Couldn't open shapefile.");
	if (layername.empty()) {
		layer = ds->GetLayer(0);
	} else {
		layer = ds->GetLayerByName(layername.c_str());
	}
	if (layer == nullptr)
		g_runerr("Couldn't get layer.");

	type = layer->GetGeomType();
	if (type != wkbPolygon && type != wkbMultiPolygon)
		g_runerr("Geometry must be polygon or multipolygon.");

	const GEOSContextHandle_t gctx = OGRGeometry::createGEOSContext();
	const gg::GeometryFactory *gf = gg::GeometryFactory::getDefaultInstance();
	const gg::CoordinateSequenceFactory *cf =
			gf->getCoordinateSequenceFactory();
	std::vector<gg::Geometry *> geoms;

	while ((feat = layer->GetNextFeature()) != NULL) {
		og = feat->GetGeometryRef();
		geom = (gg::Geometry *) og->exportToGEOS(gctx);
		geoms.push_back(geom);
	}

	GDALClose(ds);

	if (geoms.size() == 0)
		g_runerr("No geometries were found.");

	/* The geometry collection is used for checking whether a las file intersects
	 the region of interest. */
	geomColl = gf->createGeometryCollection(geoms);
	const gg::Envelope *env = geomColl->getEnvelopeInternal();
	Bounds cbounds(env->getMinX(), env->getMinY(), env->getMaxX(),
			env->getMaxY());

	/* Loop over files and figure out which ones are relevant. */
	liblas::ReaderFactory rf;
	liblas::Header *dsth = nullptr;
	std::vector<unsigned int> indices;

	for (unsigned int i = 0; i < files.size(); ++i) {

		std::ifstream in(files[i].c_str(), std::ios::in | std::ios::binary);
		liblas::Reader r = rf.CreateWithStream(in);
		liblas::Header h = r.GetHeader();

		if (i == 0)
			dsth = new liblas::Header(h);

		std::vector<gg::Coordinate> coords;
		coords.push_back(gg::Coordinate(h.GetMinX(), h.GetMinY()));
		coords.push_back(gg::Coordinate(h.GetMaxX(), h.GetMinY()));
		coords.push_back(gg::Coordinate(h.GetMaxX(), h.GetMaxY()));
		coords.push_back(gg::Coordinate(h.GetMinX(), h.GetMaxY()));
		coords.push_back(gg::Coordinate(h.GetMinX(), h.GetMinY()));

		gg::CoordinateSequence *cs = cf->create(&coords);
		gg::LinearRing *lr = gf->createLinearRing(cs);
		gg::Polygon *bounds = gf->createPolygon(lr, NULL);

		if (bounds->intersects(geomColl))
			indices.push_back(i);

		in.close();
	}

	if (indices.size() == 0)
		g_runerr("No files matched the given bounds.");

	std::ofstream out(outfile, std::ios::out | std::ios::binary);
	liblas::WriterFactory wf;
	liblas::Writer w(out, *dsth);
	liblas::Header::RecordsByReturnArray recs;
	int count = 0;
//.........这里部分代码省略.........

OGRLayer * OGRDataSource::ExecuteSQL( const char *pszStatement,
                                      OGRGeometry *poSpatialFilter,
                                      const char *pszDialect )

{
    const char *pszError;
    swq_select *psSelectInfo = NULL;

    (void) pszDialect;

    swq_field_list sFieldList;
    int            nFIDIndex = 0;
    OGRGenSQLResultsLayer *poResults = NULL;

    memset( &sFieldList, 0, sizeof(sFieldList) );

/* -------------------------------------------------------------------- */
/*      Handle CREATE INDEX statements specially.                       */
/* -------------------------------------------------------------------- */
    if( EQUALN(pszStatement,"CREATE INDEX",12) )
    {
        ProcessSQLCreateIndex( pszStatement );
        return NULL;
    }
    
/* -------------------------------------------------------------------- */
/*      Handle DROP INDEX statements specially.                         */
/* -------------------------------------------------------------------- */
    if( EQUALN(pszStatement,"DROP INDEX",10) )
    {
        ProcessSQLDropIndex( pszStatement );
        return NULL;
    }
    
/* -------------------------------------------------------------------- */
/*      Preparse the SQL statement.                                     */
/* -------------------------------------------------------------------- */
    pszError = swq_select_preparse( pszStatement, &psSelectInfo );
    if( pszError != NULL )
    {
        CPLError( CE_Failure, CPLE_AppDefined, 
                  "SQL: %s", pszError );
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Validate that all the source tables are recognised, count       */
/*      fields.                                                         */
/* -------------------------------------------------------------------- */
    int  nFieldCount = 0, iTable, iField;
    int  iEDS;
    int  nExtraDSCount = 0;
    OGRDataSource** papoExtraDS = NULL;
    OGRSFDriverRegistrar *poReg=OGRSFDriverRegistrar::GetRegistrar();

    for( iTable = 0; iTable < psSelectInfo->table_count; iTable++ )
    {
        swq_table_def *psTableDef = psSelectInfo->table_defs + iTable;
        OGRLayer *poSrcLayer;
        OGRDataSource *poTableDS = this;

        if( psTableDef->data_source != NULL )
        {
            poTableDS = (OGRDataSource *) 
                OGROpenShared( psTableDef->data_source, FALSE, NULL );
            if( poTableDS == NULL )
            {
                if( strlen(CPLGetLastErrorMsg()) == 0 )
                    CPLError( CE_Failure, CPLE_AppDefined, 
                              "Unable to open secondary datasource\n"
                              "`%s' required by JOIN.",
                              psTableDef->data_source );

                swq_select_free( psSelectInfo );
                goto end;
            }

            /* Keep in an array to release at the end of this function */
            papoExtraDS = (OGRDataSource** )CPLRealloc(papoExtraDS,
                               sizeof(OGRDataSource*) * (nExtraDSCount + 1));
            papoExtraDS[nExtraDSCount++] = poTableDS;
        }

        poSrcLayer = poTableDS->GetLayerByName( psTableDef->table_name );

        if( poSrcLayer == NULL )
        {
            CPLError( CE_Failure, CPLE_AppDefined, 
                      "SELECT from table %s failed, no such table/featureclass.",
                      psTableDef->table_name );
            swq_select_free( psSelectInfo );
            goto end;
        }

        nFieldCount += poSrcLayer->GetLayerDefn()->GetFieldCount();
    }
    
/* -------------------------------------------------------------------- */
/*      Build the field list for all indicated tables.                  */
/* -------------------------------------------------------------------- */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
/* -------------------------------------------------------------------- */
    std::set<LayerDesc> oSetLayers;
    std::set<CPLString> oSetSpatialIndex;
    CPLString osModifiedSQL;
    OGR2SQLITEGetPotentialLayerNames(pszStatement, oSetLayers,
                                     oSetSpatialIndex, osModifiedSQL);
    std::set<LayerDesc>::iterator oIter = oSetLayers.begin();

    if( strcmp(pszStatement, osModifiedSQL.c_str()) != 0 )
        CPLDebug("OGR", "Modified SQL: %s", osModifiedSQL.c_str());
    pszStatement = osModifiedSQL.c_str(); /* do not use it anymore */

    int bFoundOGRStyle = ( osModifiedSQL.ifind("OGR_STYLE") != std::string::npos );

/* -------------------------------------------------------------------- */
/*      For each of those tables, create a Virtual Table.               */
/* -------------------------------------------------------------------- */
    OGRLayer* poSingleSrcLayer = NULL;
    for(; oIter != oSetLayers.end(); ++oIter)
    {
        const LayerDesc& oLayerDesc = *oIter;
        /*CPLDebug("OGR", "Layer desc : %s, %s, %s, %s",
                 oLayerDesc.osOriginalStr.c_str(),
                 oLayerDesc.osSubstitutedName.c_str(),
                 oLayerDesc.osDSName.c_str(),
                 oLayerDesc.osLayerName.c_str());*/

        CPLString osSQL;
        OGRLayer* poLayer = NULL;
        CPLString osTableName;
        int nExtraDS;
        if( oLayerDesc.osDSName.size() == 0 )
        {
            poLayer = poDS->GetLayerByName(oLayerDesc.osLayerName);
            /* Might be a false positive (unlikely) */
            if( poLayer == NULL )
                continue;

            osTableName = oLayerDesc.osLayerName;

            nExtraDS = -1;
        }
        else
        {
            OGRDataSource* poOtherDS = (OGRDataSource* )
                OGROpen(oLayerDesc.osDSName, FALSE, NULL);
            if( poOtherDS == NULL )
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Cannot open datasource '%s'",
                         oLayerDesc.osDSName.c_str() );
                delete poSQLiteDS;
                VSIUnlink(pszTmpDBName);
                CPLFree(pszTmpDBName);
                return NULL;
            }

            poLayer = poOtherDS->GetLayerByName(oLayerDesc.osLayerName);
            if( poLayer == NULL )
            {
                CPLError(CE_Failure, CPLE_AppDefined,
                         "Cannot find layer '%s' in '%s'",
                         oLayerDesc.osLayerName.c_str(),
                         oLayerDesc.osDSName.c_str() );
                delete poOtherDS;
                delete poSQLiteDS;

void DrawShape::featuredetect(int XPoint,int YPoint)
{
	if(!sfile.empty() && XPoint>=(x_pos_shift-increase_width/2)*6.40 && XPoint<=(x_pos_shift-increase_width/2)*6.40+(100+increase_width)*6.4 && YPoint>=(y_pos_shift-increase_height/2)*4.8 && YPoint<=(y_pos_shift-increase_height/2)*4.8+(100+increase_height)*4.8 )//check if point is inside the view area 
	{


		OGRDataSource       *poDS;
		string dfile = sfile  + ".shp";
		string shp = "g_4326/" + sfile + ".shp";
		poDS = OGRSFDriverRegistrar::Open(shp.c_str(), FALSE );  //comment till here
		if( poDS == NULL )
		{
			WApplication::instance()->doJavaScript("alert('Open failed')");
			cout << "Open failed";
		}
		OGRLayer  *poLayer;
		poLayer = poDS->GetLayerByName( sfile.c_str() ); // comment here 
		OGRFeature *poFeature;
		OGREnvelope *   	 psExtent = new OGREnvelope();
		poLayer->GetExtent(psExtent);

		double xMin = psExtent->MinX;
		double yMin = psExtent->MinY;
		double xMax = psExtent->MaxX;
		double yMax = psExtent->MaxY;

		stringstream strm;
		strm << xMin;
		string exp;
		double scaleFactor;
		string bound_string = strm.str();
		int size = bound_string.size();
		size_t found=bound_string.find("+");
		if (found!=string::npos) {
			exp =  bound_string.substr(found+1,size);     
		}

		if(exp.empty()) {

			stringstream strExtent;
			strExtent << yMin;
			bound_string = strExtent.str();
			size = bound_string.size();
			found = bound_string.find("+");
			if(found!=string::npos) {
				exp = bound_string.substr(found+1,size); 
			}
		}
		if(exp.empty()) {
			stringstream strExtent;
			strExtent << xMax;
			bound_string = strExtent.str();
			size = bound_string.size();
			found = bound_string.find("+");
			if(found!=string::npos) {
				exp = bound_string.substr(found+1,size); 
			}
		}
		if(exp.empty()) {
			stringstream strExtent;
			strExtent << yMax;
			bound_string = strExtent.str();
			size = bound_string.size();
			found = bound_string.find("+");
			if(found!=string::npos) {
				exp = bound_string.substr(found+1,size); 
			}
		}
		if(!exp.empty()) {
			int exponent  = boost::lexical_cast<int>(exp);

			exponent-=3;
			scaleFactor = pow (10,exponent);
		}
		else
		{
			scaleFactor = 1;
		}

		xMin/=scaleFactor;
		xMax/=scaleFactor;
		yMin/=scaleFactor;
		yMax/=scaleFactor;

		double gWidth = xMax - xMin;
		double gHeight = yMax - yMin;
		double widthFactor = 1;
		double pwidth = abs(gWidth-gHeight);
		double s = gWidth - gHeight;
		if(s<0.16)
			gWidth = gHeight + 0.16;

		double scaledX=(XPoint*(gWidth* widthFactor*(100-increase_width)/100)/640 + (xMin+(-x_pos_shift+increase_width/2)/100*gWidth* widthFactor))*scaleFactor;
		double scaledY=(YPoint*(gHeight*widthFactor*(-1+increase_height/100))/480 + (yMax+(y_pos_shift-increase_height/2)/100*gHeight*widthFactor))*scaleFactor;
		OGRPoint *poPoint =new OGRPoint();
		poPoint->setX(scaledX);
		poPoint->setY(scaledY);


		poLayer->ResetReading();
//.........这里部分代码省略.........

int main(int argc, char **argv)
{
  // Get data from ogr
  OGRRegisterAll();
  std::cout << "Opening: " << argv[1] << std::endl;

  OGRDataSource *shp = OGRSFDriverRegistrar::Open(argv[1], FALSE);
  IsValid(shp, "Error opening file.");

  std::cout << "Shape contains " << shp->GetLayerCount() << " layers." << std::endl;
  OGRLayer *layer = shp->GetLayerByName(argv[2]);
  IsValid(layer, "Couldn't grab layer");

  OGRSpatialReference *srcSRS = NULL;
  srcSRS = layer->GetSpatialRef();

  // Set up writing
  const char *kDriverName = "ESRI Shapefile";
  OGRSFDriver *shpDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(kDriverName);
  IsValid(shpDriver, "Couldn't grab the shapefile driver.");

  IsValid(argv[3], "Please provide a output shp.");
  std::cout << "Writing to: " << argv[3] << std::endl;
  OGRDataSource *shpOut = shpDriver->CreateDataSource(argv[3], NULL);
  IsValid(shpOut, "Couldn't open output file");

  OGRLayer *outLayer = shpOut->CreateLayer(layer->GetName(), srcSRS, wkbMultiLineString, NULL);
  IsValid(outLayer, "Couldn't create an output layer");

  // copy over the fields from the source file
  OGRFeatureDefn *source = layer->GetLayerDefn();
  for(int i=0; i < source->GetFieldCount(); i++){
    OGRFieldDefn *field = source->GetFieldDefn(i);
    if(outLayer->CreateField(field) != OGRERR_NONE) {
      std::cout << "Couldn't make layer" << std::endl; exit(1);
    };
  }

  // Loop through features and grab the hull and put it into CGAL then
  // skeletonize the points
  OGRFeature *feature;
  int count = 0;
  while((feature = layer->GetNextFeature()) != NULL)
  {
    OGRMultiPolygon *geometry = dynamic_cast<OGRMultiPolygon *>(OGRGeometryFactory::forceToMultiPolygon(feature->GetGeometryRef()));
    IsValid(geometry, "No geometry.");

    OGRFeature *outFeature = OGRFeature::CreateFeature(outLayer->GetLayerDefn());
    IsValid(outFeature, "Couldn't make a feature.");

    for(int i=0; i < source->GetFieldCount(); i++){
      OGRField *field = feature->GetRawFieldRef(i);
      outFeature->SetField(i, field);
    }

    OGRGeometry* line = NULL;
    for(int i=0; i < geometry->getNumGeometries(); i++){
      OGRGeometry* segment = BuildMultiLine(geometry->getGeometryRef(i));
      if(segment != NULL){
        if(line == NULL) { line = new OGRLineString; }
        OGRGeometry* tmp = line->Union(segment);
        if(tmp != NULL){
          delete line;
          line = tmp;
        }
        delete segment;
      }
    }
    outFeature->SetGeometry(line);
    if(outLayer->CreateFeature(outFeature) != OGRERR_NONE){
      std::cout << "Couldn't create feature." << std::endl;
      exit(1);
    }

    // clean up
    OGRFeature::DestroyFeature(outFeature);
    std::cout << std::endl << ++count << std::endl;
  }

  // cleanup
  OGRDataSource::DestroyDataSource(shp);
  OGRDataSource::DestroyDataSource(shpOut);
  return 0;
}

void readQuadIndex(int series, QString file, QString layerName, Projection *pj)
{
  OGRDataSource       *ds = OGRSFDriverRegistrar::Open(file.toLatin1().data(), false);
  if (!ds) {
    fprintf(stderr, "Could not open quad index '%s'.\n", file.toLatin1().data());
    exit(-1);
  }

  OGRLayer  *layer;
  layer = ds->GetLayerByName(layerName.toLatin1().data());
  if (!layer) {
    fprintf(stderr, "Could not read layer '%s'.\n", layerName.toLatin1().data());
    exit(-1);
  }

  OGRSpatialReference *srs = layer->GetSpatialRef();
  if (!srs) {
    fprintf(stderr, "Missing spatial reference for layer '%s'.\n", 
            layerName.toLatin1().data());
    exit(-1);
  }
  char *proj = NULL;
  srs->exportToProj4(&proj);
  if (!proj) {
    fprintf(stderr, "Error computing PROJ4 spatial reference for layer '%s'.\n", 
            layerName.toLatin1().data());
    exit(-1);
  }
  Projection pjIndex(proj);
  CPLFree(proj);


  layer->ResetReading();
  OGRFeatureDefn *def = layer->GetLayerDefn();

  int idFieldNr = def->GetFieldIndex("ID");
  int nameFieldNr = def->GetFieldIndex("NAME");
  if (idFieldNr < 0 || nameFieldNr < 0) {
    fprintf(stderr, "Missing index layer fields.\n");
    exit(-1);
  }


  OGRFeature *f;
  while ((f = layer->GetNextFeature()) != NULL) {
    QString id(f->GetFieldAsString(idFieldNr));
    QString name(f->GetFieldAsString(nameFieldNr));

    //    printf("Quad id: %s; name: %s\n", id.toLatin1().data(), name.toLatin1().data());

    OGRGeometry *g;
    g = f->GetGeometryRef();
    if (g != NULL && wkbFlatten(g->getGeometryType()) == wkbPolygon) {
      OGRPolygon *p = (OGRPolygon *)g;
      OGRLinearRing *r = p->getExteriorRing();
      if (!r) {
        fprintf(stderr, "Quad has no exterior polygon ring %s\n", 
                id.toLatin1().data());
        continue;
      }

      int size = r->getNumPoints();
      QPolygonF boundary;
      for (int i = 0; i < size; i++) {
        OGRPoint p;
        r->getPoint(i, &p);
        boundary << QPointF(p.getX(), p.getY());
      }

      QPolygonF projBoundary = pj->transformFrom(&pjIndex, boundary);
      
      quads[id] = Quad(series, id, name, projBoundary);
    } else {
      fprintf(stderr, "Missing or invalid geometry for quad %s\n", 
              id.toLatin1().data());
    }
  }

  OGRDataSource::DestroyDataSource(ds);
}

static OGRGeometryCollection* LoadGeometry( const char* pszDS,
                                            const char* pszSQL,
                                            const char* pszLyr,
                                            const char* pszWhere )
{
    OGRDataSource       *poDS;
    OGRLayer            *poLyr;
    OGRFeature          *poFeat;
    OGRGeometryCollection *poGeom = NULL;
        
    poDS = OGRSFDriverRegistrar::Open( pszDS, FALSE );
    if ( poDS == NULL )
        return NULL;

    if ( pszSQL != NULL )
        poLyr = poDS->ExecuteSQL( pszSQL, NULL, NULL ); 
    else if ( pszLyr != NULL )
        poLyr = poDS->GetLayerByName( pszLyr );
    else
        poLyr = poDS->GetLayer(0);
        
    if ( poLyr == NULL )
    {
        fprintf( stderr,
            "FAILURE: Failed to identify source layer from datasource.\n" );
        OGRDataSource::DestroyDataSource( poDS );
        return NULL;
    }
    
    if ( pszWhere )
        poLyr->SetAttributeFilter( pszWhere );
        
    while ( (poFeat = poLyr->GetNextFeature()) != NULL )
    {
        OGRGeometry* poSrcGeom = poFeat->GetGeometryRef();
        if ( poSrcGeom )
        {
            OGRwkbGeometryType eType =
                wkbFlatten( poSrcGeom->getGeometryType() );
            
            if ( poGeom == NULL )
                poGeom = new OGRMultiPolygon();

            if ( eType == wkbPolygon )
                poGeom->addGeometry( poSrcGeom );
            else if ( eType == wkbMultiPolygon )
            {
                int iGeom;
                int nGeomCount =
                    ((OGRMultiPolygon *)poSrcGeom)->getNumGeometries();

                for ( iGeom = 0; iGeom < nGeomCount; iGeom++ )
                {
                    poGeom->addGeometry(
                        ((OGRMultiPolygon *)poSrcGeom)->getGeometryRef(iGeom) );
                }
            }
            else
            {
                fprintf( stderr, "FAILURE: Geometry not of polygon type.\n" );
                OGRGeometryFactory::destroyGeometry( poGeom );
                OGRFeature::DestroyFeature( poFeat );
                if ( pszSQL != NULL )
                    poDS->ReleaseResultSet( poLyr );
                OGRDataSource::DestroyDataSource( poDS );
                return NULL;
            }
        }
    
        OGRFeature::DestroyFeature( poFeat );
    }
    
    if( pszSQL != NULL )
        poDS->ReleaseResultSet( poLyr );
    OGRDataSource::DestroyDataSource( poDS );
    
    return poGeom;
}