C++ OGRLinearRing类代码示例

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

    EarlySetConfigOptions(nArgc, papszArgv);

    OGRRegisterAll();

/* -------------------------------------------------------------------- */
/*      Processing command line arguments.                              */
/* -------------------------------------------------------------------- */
    nArgc = OGRGeneralCmdLineProcessor(nArgc, &papszArgv, 0);

    if( nArgc < 1 )
        exit(-nArgc);

    char *pszWHERE = nullptr;
    const char *pszDataSource = nullptr;
    char **papszLayers = nullptr;
    OGRGeometry *poSpatialFilter = nullptr;
    int nRepeatCount = 1;
    bool bAllLayers = false;
    char *pszSQLStatement = nullptr;
    const char *pszDialect = nullptr;
    int nRet = 0;
    const char* pszGeomField = nullptr;
    char **papszOpenOptions = nullptr;
    char **papszExtraMDDomains = nullptr;
    bool bListMDD = false;
    bool bShowMetadata = true;
    bool bFeatureCount = true;
    bool bExtent = true;
    bool bDatasetGetNextFeature = false;
    bool bReadOnly = false;
    bool bUpdate = false;
    const char* pszWKTFormat = "WKT2";

    for( int iArg = 1; iArg < nArgc; iArg++ )
    {
        if( EQUAL(papszArgv[iArg], "--utility_version") )
        {
            printf("%s was compiled against GDAL %s and "
                   "is running against GDAL %s\n",
                   papszArgv[0], GDAL_RELEASE_NAME,
                   GDALVersionInfo("RELEASE_NAME"));
            CSLDestroy(papszArgv);
            return 0;
        }
        else if( EQUAL(papszArgv[iArg], "--help") )
        {
            Usage();
        }
        else if( EQUAL(papszArgv[iArg], "-ro") )
        {
            bReadOnly = true;
        }
        else if( EQUAL(papszArgv[iArg], "-update") )
        {
            bUpdate = true;
        }
        else if( EQUAL(papszArgv[iArg], "-q") ||
                 EQUAL(papszArgv[iArg], "-quiet"))
        {
            bVerbose = false;
        }
        else if( EQUAL(papszArgv[iArg], "-qq") )
        {
            /* Undocumented: mainly only useful for AFL testing */
            bVerbose = false;
            bSuperQuiet = true;
        }
        else if( EQUAL(papszArgv[iArg], "-fid") )
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(1);
            nFetchFID = CPLAtoGIntBig(papszArgv[++iArg]);
        }
        else if( EQUAL(papszArgv[iArg], "-spat") )
        {
            CHECK_HAS_ENOUGH_ADDITIONAL_ARGS(4);

            OGRLinearRing oRing;
            oRing.addPoint(CPLAtof(papszArgv[iArg+1]),
                           CPLAtof(papszArgv[iArg+2]));
            oRing.addPoint(CPLAtof(papszArgv[iArg+1]),
                           CPLAtof(papszArgv[iArg+4]));
            oRing.addPoint(CPLAtof(papszArgv[iArg+3]),
                           CPLAtof(papszArgv[iArg+4]));
            oRing.addPoint(CPLAtof(papszArgv[iArg+3]),
                           CPLAtof(papszArgv[iArg+2]));
            oRing.addPoint(CPLAtof(papszArgv[iArg+1]),
                           CPLAtof(papszArgv[iArg+2]));

            poSpatialFilter = new OGRPolygon();
            static_cast<OGRPolygon *>(poSpatialFilter)->addRing(&oRing);
            iArg += 4;
        }
        else if( EQUAL(papszArgv[iArg], "-geomfield") )
        {
//.........这里部分代码省略.........

int OGRPolygonLabelPoint(OGRPolygon *poPoly, OGRPoint *poLabelPoint)
{
    double        slope;
    OGRRawPoint   point1, point2;
    int           i, j, k, nfound;
    double        x, y, *xintersect, temp;
    double        hi_y, lo_y;
    int           wrong_order, n;
    double        len, max_len=0;
    double        skip;
    OGREnvelope   oEnv;

    if (poPoly == NULL)
        return OGRERR_FAILURE;

    poPoly->getEnvelope(&oEnv);

    poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);
    poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);

    //if(get_centroid(p, lp, &miny, &maxy) == -1) return(-1);

    if(OGRIntersectPointPolygon(poLabelPoint, poPoly) == TRUE) /* cool, done */
        return OGRERR_NONE;

    /* do it the hard way - scanline */

    skip = (oEnv.MaxY - oEnv.MinY)/NUM_SCANLINES;

    n=0;
    for(j=0; j<OGR_NUM_RINGS(poPoly); j++) 
    {
        /* count total number of points */
        n += OGR_GET_RING(poPoly, j)->getNumPoints();
    }

    xintersect = (double *)calloc(n, sizeof(double));
    if (xintersect == NULL)
        return OGRERR_FAILURE;

    for(k=1; k<=NUM_SCANLINES; k++) 
    { 
        /* sample the shape in the y direction */
    
        y = oEnv.MaxY - k*skip; 

        /* need to find a y that won't intersect any vertices exactly */  
        hi_y = y - 1; /* first initializing lo_y, hi_y to be any 2 pnts on either side of lp->y */
        lo_y = y + 1;
        for(j=0; j<OGR_NUM_RINGS(poPoly); j++) 
        {
            OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);

            if((lo_y < y) && (hi_y >= y)) 
                break; /* already initialized */
            for(i=0; i < poRing->getNumPoints(); i++) 
            {   
                if((lo_y < y) && (hi_y >= y)) 
                    break; /* already initialized */
                if(poRing->getY(i) < y)
                    lo_y = poRing->getY(i);
                if(poRing->getY(i) >= y)
                    hi_y = poRing->getY(i);
            }
        }

        n=0;
        for(j=0; j<OGR_NUM_RINGS(poPoly); j++) 
        {
            OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);

            for(i=0; i < poRing->getNumPoints(); i++) 
            {
                if((poRing->getY(i) < y) && 
                   ((y - poRing->getY(i)) < (y - lo_y)))
                    lo_y = poRing->getY(i);
                if((poRing->getY(i) >= y) && 
                   ((poRing->getY(i) - y) < (hi_y - y)))
                    hi_y = poRing->getY(i);
            }      
        }

        if(lo_y == hi_y) 
            return OGRERR_FAILURE;
        else  
            y = (hi_y + lo_y)/2.0;    
    
        nfound = 0;
        for(j=0; j<OGR_NUM_RINGS(poPoly); j++)   /* for each line */
        {
            OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);
            point1.x = poRing->getX(poRing->getNumPoints()-1);
            point1.y = poRing->getY(poRing->getNumPoints()-1);

            for(i=0; i < poRing->getNumPoints(); i++) 
            {
                point2.x = poRing->getX(i);
                point2.y = poRing->getY(i);
        
                if(EDGE_CHECK(point1.y, y, point2.y) == CLIP_MIDDLE) 
//.........这里部分代码省略.........

OGRGeometry *OGRIngresLayer::TranslateGeometry( const char *pszGeom )

{
    OGRGeometry *poGeom = NULL;

/* -------------------------------------------------------------------- */
/*      Parse the tuple list into an array of x/y vertices.  The        */
/*      input may look like "(2,3)" or "((2,3),(4,5),...)".  Extra      */
/*      spaces may occur between tokens.                                */
/* -------------------------------------------------------------------- */
    double *padfXY = NULL;
    int    nVertMax = 0, nVertCount = 0;
    int    nDepth = 0;
    const char *pszNext = pszGeom;

    while( *pszNext != '\0' )
    {
        while( *pszNext == ' ' )
            pszNext++;

        if( *pszNext == '(' )
        {
            pszNext++;
            nDepth++;
            continue;
        }

        if( *pszNext == ')' )
        {
            pszNext++;
            CPLAssert( nDepth == 1 );
            nDepth--;
            break;
        }

        if( *pszNext == ',' )
        {
            pszNext++;
            CPLAssert( nDepth == 1 );
            continue;
        }

        if( nVertCount == nVertMax )
        {
            nVertMax = nVertMax * 2 + 1;
            padfXY = (double *) 
                CPLRealloc(padfXY, sizeof(double) * nVertMax * 2 );
        }

        if( !ParseXY( &pszNext, padfXY + nVertCount*2 ) )
        {
            CPLDebug( "INGRES", "Error parsing geometry: %s", 
                      pszGeom );
            CPLFree( padfXY );
            return NULL;
        }
        
        CPLAssert( *pszNext == ')' );
        nVertCount++;
        pszNext++;
        nDepth--;

        while( *pszNext == ' ' )
            pszNext++;
    }

    CPLAssert( nDepth == 0 );

/* -------------------------------------------------------------------- */
/*      Handle Box/IBox.                                                */
/* -------------------------------------------------------------------- */
    if( EQUAL(osIngresGeomType,"BOX")
        || EQUAL(osIngresGeomType,"IBOX") )
    {
        CPLAssert( nVertCount == 2 );

        OGRLinearRing *poRing = new OGRLinearRing();
        poRing->addPoint( padfXY[0], padfXY[1] );
        poRing->addPoint( padfXY[2], padfXY[1] );
        poRing->addPoint( padfXY[2], padfXY[3] );
        poRing->addPoint( padfXY[0], padfXY[3] );
        poRing->addPoint( padfXY[0], padfXY[1] );

        OGRPolygon *poPolygon = new OGRPolygon();
        poPolygon->addRingDirectly( poRing );

        poGeom = poPolygon;
    }

/* -------------------------------------------------------------------- */
/*      Handle Point/IPoint                                             */
/* -------------------------------------------------------------------- */
    else if( EQUAL(osIngresGeomType,"POINT")
             || EQUAL(osIngresGeomType,"IPOINT") )
    {
        CPLAssert( nVertCount == 1 );

        poGeom = new OGRPoint( padfXY[0], padfXY[1] );
    }

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

OGRGeometryH OGRBuildPolygonFromEdges( OGRGeometryH hLines,
                                       int bBestEffort,
                                       int bAutoClose,
                                       double dfTolerance,
                                       OGRErr * peErr )

{
    int         bSuccess = TRUE;
    OGRGeometryCollection *poLines = (OGRGeometryCollection *) hLines;
    OGRPolygon  *poPolygon = new OGRPolygon();

    (void) bBestEffort;

/* -------------------------------------------------------------------- */
/*      Setup array of line markers indicating if they have been        */
/*      added to a ring yet.                                            */
/* -------------------------------------------------------------------- */
    int         nEdges = poLines->getNumGeometries();
    int         *panEdgeConsumed, nRemainingEdges = nEdges;

    panEdgeConsumed = (int *) CPLCalloc(sizeof(int),nEdges);

/* ==================================================================== */
/*      Loop generating rings.                                          */
/* ==================================================================== */
    while( nRemainingEdges > 0 )
    {
        int             iEdge;
        OGRLineString   *poLine;

/* -------------------------------------------------------------------- */
/*      Find the first unconsumed edge.                                 */
/* -------------------------------------------------------------------- */
        for( iEdge = 0; panEdgeConsumed[iEdge]; iEdge++ ) {}

        poLine = (OGRLineString *) poLines->getGeometryRef(iEdge);

/* -------------------------------------------------------------------- */
/*      Start a new ring, copying in the current line directly          */
/* -------------------------------------------------------------------- */
        OGRLinearRing   *poRing = new OGRLinearRing();

        AddEdgeToRing( poRing, poLine, FALSE );

        panEdgeConsumed[iEdge] = TRUE;
        nRemainingEdges--;

/* ==================================================================== */
/*      Loop adding edges to this ring until we make a whole pass       */
/*      within finding anything to add.                                 */
/* ==================================================================== */
        int             bWorkDone = TRUE;
        double          dfBestDist = dfTolerance;

        while( !CheckPoints(poRing,0,poRing,poRing->getNumPoints()-1,NULL)
               && nRemainingEdges > 0
               && bWorkDone )
        {
            int         iBestEdge = -1, bReverse = FALSE;

            bWorkDone = FALSE;
            dfBestDist = dfTolerance;

            // We consider linking the end to the beginning.  If this is
            // closer than any other option we will just close the loop.

            //CheckPoints(poRing,0,poRing,poRing->getNumPoints()-1,&dfBestDist);

            // Find unused edge with end point closest to our loose end.
            for( iEdge = 0; iEdge < nEdges; iEdge++ )
            {
                if( panEdgeConsumed[iEdge] )
                    continue;

                poLine = (OGRLineString *) poLines->getGeometryRef(iEdge);

                if( CheckPoints(poLine,0,poRing,poRing->getNumPoints()-1,
                                &dfBestDist) )
                {
                    iBestEdge = iEdge;
                    bReverse = FALSE;
                }
                if( CheckPoints(poLine,poLine->getNumPoints()-1,
                                poRing,poRing->getNumPoints()-1,
                                &dfBestDist) )
                {
                    iBestEdge = iEdge;
                    bReverse = TRUE;
                }
            }

            // We found one within tolerance - add it.
            if( iBestEdge != -1 )
            {
                poLine = (OGRLineString *)
                    poLines->getGeometryRef(iBestEdge);

                AddEdgeToRing( poRing, poLine, bReverse );

                panEdgeConsumed[iBestEdge] = TRUE;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
					/* XXX - Should those pointers be checked against NULL? */ 
					assert(NULL != poField);
					assert(NULL != poFieldCur);

					if( poField->GetType() != poFieldCur->GetType() 
						|| poField->GetWidth() != poFieldCur->GetWidth() 
						|| poField->GetPrecision() != poFieldCur->GetPrecision() 
						|| !EQUAL( poField->GetNameRef(), poFieldCur->GetNameRef() ) )
					{
						fprintf( stderr, "Schema of attributes of layer %s of %s does not match ... skipping it.\n",
                                 poLayer->GetLayerDefn()->GetName(), papszArgv[nFirstSourceDataset]);
                        if (bFirstWarningForNonMatchingAttributes)
                        {
                            fprintf( stderr, "Note : you can override this behaviour with -accept_different_schemas option\n"
                                             "but this may result in a tileindex incompatible with MapServer\n");
                            bFirstWarningForNonMatchingAttributes = FALSE;
                        }
                        bSkip = TRUE; 
                        break;
					}
				}
                
                if (bSkip)
                    continue;
			}


/* -------------------------------------------------------------------- */
/*      Get layer extents, and create a corresponding polygon           */
/*      geometry.                                                       */
/* -------------------------------------------------------------------- */
            OGREnvelope sExtents;
            OGRPolygon oRegion;
            OGRLinearRing oRing;

            if( poLayer->GetExtent( &sExtents, TRUE ) != OGRERR_NONE )
            {
                fprintf( stderr, "GetExtent() failed on layer %s of %s, skipping.\n", 
                        poLayer->GetLayerDefn()->GetName(), 
                        papszArgv[nFirstSourceDataset] );
                continue;
            }
            
            oRing.addPoint( sExtents.MinX, sExtents.MinY );
            oRing.addPoint( sExtents.MinX, sExtents.MaxY );
            oRing.addPoint( sExtents.MaxX, sExtents.MaxY );
            oRing.addPoint( sExtents.MaxX, sExtents.MinY );
            oRing.addPoint( sExtents.MinX, sExtents.MinY );

            oRegion.addRing( &oRing );

/* -------------------------------------------------------------------- */
/*      Add layer to tileindex.                                         */
/* -------------------------------------------------------------------- */
            char        szLocation[5000];
            OGRFeature  oTileFeat( poDstLayer->GetLayerDefn() );

            sprintf( szLocation, "%s,%d", 
                     fileNameToWrite, iLayer );
            oTileFeat.SetGeometry( &oRegion );
            oTileFeat.SetField( iTileIndexField, szLocation );

            if( poDstLayer->CreateFeature( &oTileFeat ) != OGRERR_NONE )
            {
                fprintf( stderr, "Failed to create feature on tile index ... terminating." );
                GDALClose( (GDALDatasetH) poDstDS );

//.........这里部分代码省略.........
        else if( poFeatureDefn->GetGeomType() == wkbLineString25D
                 || (wkbFlatten(poFeatureDefn->GetGeomType()) == wkbUnknown
                     && aoVertices.size() > 1) )
        {
            // We should likely be applying ringstart to break things into
            // a multilinestring in some cases.
            if( aoVertices.size() > 1 )
            {
                OGRLineString *poLS = new OGRLineString();

                poLS->setNumPoints( static_cast<int>(aoVertices.size()) );

                for( unsigned int i = 0; i < aoVertices.size(); i++ )
                    poLS->setPoint( i,
                                    aoVertices[i].x,
                                    aoVertices[i].y,
                                    aoVertices[i].z );
                if (poSRS)
                    poLS->assignSpatialReference(poSRS);

                poFeature->SetGeometryDirectly( poLS );
            }
            else
            {
                // report issue?
            }
        }

/* -------------------------------------------------------------------- */
/*      Polygon - Currently we have no way to recognise if we are       */
/*      dealing with a multipolygon when we have more than one          */
/*      ring.  Also, PCIDSK allows the rings to be in arbitrary         */
/*      order, not necessarily outside first which we are not yet       */
/*      ready to address in the following code.                         */
/* -------------------------------------------------------------------- */
        else if( poFeatureDefn->GetGeomType() == wkbPolygon25D )
        {
            std::vector<PCIDSK::int32> anRingStart;
            OGRPolygon *poPoly = new OGRPolygon();

            if( iRingStartField != -1 )
                anRingStart = aoFields[iRingStartField].GetValueCountedInt();

            for( unsigned int iRing = 0; iRing < anRingStart.size()+1; iRing++ )
            {
                int iStartVertex;
                if( iRing == 0 )
                    iStartVertex = 0;
                else
                    iStartVertex = anRingStart[iRing-1];

                int iEndVertex;
                if( iRing == anRingStart.size() )
                    iEndVertex = static_cast<int>(aoVertices.size()) - 1;
                else
                    iEndVertex = anRingStart[iRing] - 1;

                int iVertex;
                OGRLinearRing *poRing = new OGRLinearRing();
                poRing->setNumPoints( iEndVertex - iStartVertex + 1 );
                for( iVertex = iStartVertex; iVertex <= iEndVertex; iVertex++ )
                {
                    poRing->setPoint( iVertex - iStartVertex,
                                      aoVertices[iVertex].x,
                                      aoVertices[iVertex].y,
                                      aoVertices[iVertex].z );
                }

                poPoly->addRingDirectly( poRing );
            }

            if (poSRS)
                poPoly->assignSpatialReference(poSRS);

            poFeature->SetGeometryDirectly( poPoly );
        }
    }

/* -------------------------------------------------------------------- */
/*      Trap exceptions and report as CPL errors.                       */
/* -------------------------------------------------------------------- */
    catch( PCIDSK::PCIDSKException ex )
    {
        delete poFeature;
        CPLError( CE_Failure, CPLE_AppDefined,
                  "%s", ex.what() );
        return NULL;
    }
    catch(...)
    {
        delete poFeature;
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Non-PCIDSK exception trapped." );
        return NULL;
    }

    m_nFeaturesRead++;

    return poFeature;
}

static int TestSpatialFilter( OGRLayer *poLayer )

{
    int bRet = TRUE;
    OGRFeature  *poFeature, *poTargetFeature;
    OGRPolygon  oInclusiveFilter, oExclusiveFilter;
    OGRLinearRing oRing;
    OGREnvelope sEnvelope;
    int         nInclusiveCount;

/* -------------------------------------------------------------------- */
/*      Read the target feature.                                        */
/* -------------------------------------------------------------------- */
    poLayer->ResetReading();
    poTargetFeature = poLayer->GetNextFeature();

    if( poTargetFeature == NULL )
    {
        printf( "INFO: Skipping Spatial Filter test for %s.\n"
                "      No features in layer.\n",
                poLayer->GetName() );
        return bRet;
    }

    if( poTargetFeature->GetGeometryRef() == NULL )
    {
        printf( "INFO: Skipping Spatial Filter test for %s,\n"
                "      target feature has no geometry.\n",
                poTargetFeature->GetDefnRef()->GetName() );
        OGRFeature::DestroyFeature(poTargetFeature);
        return bRet;
    }

    poTargetFeature->GetGeometryRef()->getEnvelope( &sEnvelope );

/* -------------------------------------------------------------------- */
/*      Construct inclusive filter.                                     */
/* -------------------------------------------------------------------- */
    
    oRing.setPoint( 0, sEnvelope.MinX - 20.0, sEnvelope.MinY - 20.0 );
    oRing.setPoint( 1, sEnvelope.MinX - 20.0, sEnvelope.MaxY + 10.0 );
    oRing.setPoint( 2, sEnvelope.MaxX + 10.0, sEnvelope.MaxY + 10.0 );
    oRing.setPoint( 3, sEnvelope.MaxX + 10.0, sEnvelope.MinY - 20.0 );
    oRing.setPoint( 4, sEnvelope.MinX - 20.0, sEnvelope.MinY - 20.0 );
    
    oInclusiveFilter.addRing( &oRing );

    poLayer->SetSpatialFilter( &oInclusiveFilter );

/* -------------------------------------------------------------------- */
/*      Verify that we can find the target feature.                     */
/* -------------------------------------------------------------------- */
    poLayer->ResetReading();

    while( (poFeature = poLayer->GetNextFeature()) != NULL )
    {
        if( poFeature->Equal(poTargetFeature) )
        {
            OGRFeature::DestroyFeature(poFeature);
            break;
        }
        else
            OGRFeature::DestroyFeature(poFeature);
    }

    if( poFeature == NULL )
    {
        bRet = FALSE;
        printf( "ERROR: Spatial filter eliminated a feature unexpectedly!\n");
    }
    else if( bVerbose )
    {
        printf( "INFO: Spatial filter inclusion seems to work.\n" );
    }

    nInclusiveCount = poLayer->GetFeatureCount();

/* -------------------------------------------------------------------- */
/*      Construct exclusive filter.                                     */
/* -------------------------------------------------------------------- */
    oRing.setPoint( 0, sEnvelope.MinX - 20.0, sEnvelope.MinY - 20.0 );
    oRing.setPoint( 1, sEnvelope.MinX - 10.0, sEnvelope.MinY - 20.0 );
    oRing.setPoint( 2, sEnvelope.MinX - 10.0, sEnvelope.MinY - 10.0 );
    oRing.setPoint( 3, sEnvelope.MinX - 20.0, sEnvelope.MinY - 10.0 );
    oRing.setPoint( 4, sEnvelope.MinX - 20.0, sEnvelope.MinY - 20.0 );
    
    oExclusiveFilter.addRing( &oRing );

    poLayer->SetSpatialFilter( &oExclusiveFilter );

/* -------------------------------------------------------------------- */
/*      Verify that we can find the target feature.                     */
/* -------------------------------------------------------------------- */
    poLayer->ResetReading();

    while( (poFeature = poLayer->GetNextFeature()) != NULL )
    {
        if( poFeature->Equal(poTargetFeature) )
        {
            OGRFeature::DestroyFeature(poFeature);
//.........这里部分代码省略.........

OGRFeature *OGROGDILayer::GetNextRawFeature()
{
/* -------------------------------------------------------------------- */
/*      Retrieve object from OGDI server and create new feature         */
/* -------------------------------------------------------------------- */
    ecs_Result *psResult = cln_GetNextObject(m_nClientID);
    if (! ECSSUCCESS(psResult))
    {
        if( ECSERROR( psResult ) &&
            (psResult->message == nullptr ||
             strstr(psResult->message, "End of selection") == nullptr) )
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Access to next object of layer '%s' failed: %s\n",
                      m_pszOGDILayerName,
                      psResult->message ? psResult->message : "(no error string)" );
        }
        // We probably reached EOF... keep track of shape count.
        m_nTotalShapeCount = m_iNextShapeId - m_nFilteredOutShapes;
        return nullptr;
    }

    OGRFeature *poFeature = new OGRFeature(m_poFeatureDefn);

    poFeature->SetFID( m_iNextShapeId++ );
    m_nFeaturesRead++;

/* -------------------------------------------------------------------- */
/*      Process geometry                                                */
/* -------------------------------------------------------------------- */
    if (m_eFamily == Point)
    {
        ecs_Point *psPoint = &(ECSGEOM(psResult).point);
        OGRPoint *poOGRPoint = new OGRPoint(psPoint->c.x, psPoint->c.y);

        poOGRPoint->assignSpatialReference(m_poSpatialRef);
        poFeature->SetGeometryDirectly(poOGRPoint);
    }
    else if (m_eFamily == Line)
    {
        ecs_Line        *psLine = &(ECSGEOM(psResult).line);
        OGRLineString   *poOGRLine = new OGRLineString();

        poOGRLine->setNumPoints( psLine->c.c_len );

        for( int i = 0; i < (int) psLine->c.c_len; i++ )
        {
            poOGRLine->setPoint(i, psLine->c.c_val[i].x, psLine->c.c_val[i].y);
        }

        poOGRLine->assignSpatialReference(m_poSpatialRef);
        poFeature->SetGeometryDirectly(poOGRLine);
    }
    else if (m_eFamily == Area)
    {
        ecs_Area        *psArea = &(ECSGEOM(psResult).area);
        OGRPolygon      *poOGRPolygon = new OGRPolygon();

        for( int iRing = 0; iRing < (int) psArea->ring.ring_len; iRing++ )
        {
            ecs_FeatureRing     *psRing = &(psArea->ring.ring_val[iRing]);
            OGRLinearRing       *poOGRRing = new OGRLinearRing();

            poOGRRing->setNumPoints( psRing->c.c_len );

            for( int i = 0; i < (int) psRing->c.c_len; i++ )
            {
                poOGRRing->setPoint(i, psRing->c.c_val[i].x,
                                    psRing->c.c_val[i].y);
            }
            poOGRPolygon->addRingDirectly(poOGRRing);
        }

        // __TODO__
        // When OGR supports polygon centroids then we should carry them here

        poOGRPolygon->assignSpatialReference(m_poSpatialRef);
        poFeature->SetGeometryDirectly(poOGRPolygon);
    }
    else if (m_eFamily == Text)
    {
        // __TODO__
        // For now text is treated as a point and string is lost
        //
        ecs_Text       *psText = &(ECSGEOM(psResult).text);
        OGRPoint        *poOGRPoint = new OGRPoint(psText->c.x, psText->c.y);

        poOGRPoint->assignSpatialReference(m_poSpatialRef);
        poFeature->SetGeometryDirectly(poOGRPoint);
    }
    else
    {
        CPLAssert(false);
    }

/* -------------------------------------------------------------------- */
/*      Set attributes                                                  */
/* -------------------------------------------------------------------- */
    char *pszAttrList = ECSOBJECTATTR(psResult);

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

OGRErr 
OGROCIWritableLayer::TranslateElementGroup( OGRGeometry *poGeometry )

{
    switch( wkbFlatten(poGeometry->getGeometryType()) )
    {
      case wkbPoint:
      {
          OGRPoint *poPoint = (OGRPoint *) poGeometry;

          PushElemInfo( nOrdinalCount+1, 1, 1 );

          PushOrdinal( poPoint->getX() );
          PushOrdinal( poPoint->getY() );
          if( nDimension == 3 )
              PushOrdinal( poPoint->getZ() );

          return OGRERR_NONE;
      }

      case wkbLineString:
      {
          OGRLineString *poLine = (OGRLineString *) poGeometry;
          int  iVert;
          
          PushElemInfo( nOrdinalCount+1, 2, 1 );

          for( iVert = 0; iVert < poLine->getNumPoints(); iVert++ )
          {
              PushOrdinal( poLine->getX(iVert) );
              PushOrdinal( poLine->getY(iVert) );
              if( nDimension == 3 )
                  PushOrdinal( poLine->getZ(iVert) );
          }
          return OGRERR_NONE;
      }

      case wkbPolygon:
      {
          OGRPolygon *poPoly = (OGRPolygon *) poGeometry;
          int iRing;

          for( iRing = -1; iRing < poPoly->getNumInteriorRings(); iRing++ )
          {
              OGRLinearRing *poRing;
              int            iVert;

              if( iRing == -1 )
                  poRing = poPoly->getExteriorRing();
              else
                  poRing = poPoly->getInteriorRing(iRing);

              if( iRing == -1 )
                  PushElemInfo( nOrdinalCount+1, 1003, 1 );
              else
                  PushElemInfo( nOrdinalCount+1, 2003, 1 );

              if( (iRing == -1 && poRing->isClockwise())
                  || (iRing != -1 && !poRing->isClockwise()) )
              {
                  for( iVert = poRing->getNumPoints()-1; iVert >= 0; iVert-- )
                  {
                      PushOrdinal( poRing->getX(iVert) );
                      PushOrdinal( poRing->getY(iVert) );
                      if( nDimension == 3 )
                          PushOrdinal( poRing->getZ(iVert) );
                  }
              }
              else
              {
                  for( iVert = 0; iVert < poRing->getNumPoints(); iVert++ )
                  {
                      PushOrdinal( poRing->getX(iVert) );
                      PushOrdinal( poRing->getY(iVert) );
                      if( nDimension == 3 )
                          PushOrdinal( poRing->getZ(iVert) );
                  }
              }
          }

          return OGRERR_NONE;
      }

      default:
      {
          return OGRERR_FAILURE;
      }
    }
}

int
NBHeightMapper::loadShapeFile(const std::string& file) {
#ifdef HAVE_GDAL
#if GDAL_VERSION_MAJOR < 2
    OGRRegisterAll();
    OGRDataSource* ds = OGRSFDriverRegistrar::Open(file.c_str(), FALSE);
#else
    GDALAllRegister();
    GDALDataset* ds = (GDALDataset*)GDALOpenEx(file.c_str(), GDAL_OF_VECTOR | GA_ReadOnly, NULL, NULL, NULL);
#endif
    if (ds == NULL) {
        throw ProcessError("Could not open shape file '" + file + "'.");
    }

    // begin file parsing
    OGRLayer* layer = ds->GetLayer(0);
    layer->ResetReading();

    // triangle coordinates are stored in WGS84 and later matched with network coordinates in WGS84
    // build coordinate transformation
    OGRSpatialReference* sr_src = layer->GetSpatialRef();
    OGRSpatialReference sr_dest;
    sr_dest.SetWellKnownGeogCS("WGS84");
    OGRCoordinateTransformation* toWGS84 = OGRCreateCoordinateTransformation(sr_src, &sr_dest);
    if (toWGS84 == 0) {
        WRITE_WARNING("Could not create geocoordinates converter; check whether proj.4 is installed.");
    }

    int numFeatures = 0;
    OGRFeature* feature;
    layer->ResetReading();
    while ((feature = layer->GetNextFeature()) != NULL) {
        OGRGeometry* geom = feature->GetGeometryRef();
        assert(geom != 0);

        // @todo gracefull handling of shapefiles with unexpected contents or any error handling for that matter
        assert(std::string(geom->getGeometryName()) == std::string("POLYGON"));
        // try transform to wgs84
        geom->transform(toWGS84);
        OGRLinearRing* cgeom = ((OGRPolygon*) geom)->getExteriorRing();
        // assume TIN with with 4 points and point0 == point3
        assert(cgeom->getNumPoints() == 4);
        PositionVector corners;
        for (int j = 0; j < 3; j++) {
            Position pos((double) cgeom->getX(j), (double) cgeom->getY(j), (double) cgeom->getZ(j));
            corners.push_back(pos);
            myBoundary.add(pos);
        }
        addTriangle(corners);
        numFeatures++;

        /*
        OGRwkbGeometryType gtype = geom->getGeometryType();
        switch (gtype) {
            case wkbPolygon: {
                break;
            }
            case wkbPoint: {
                WRITE_WARNING("got wkbPoint");
                break;
            }
            case wkbLineString: {
                WRITE_WARNING("got wkbLineString");
                break;
            }
            case wkbMultiPoint: {
                WRITE_WARNING("got wkbMultiPoint");
                break;
            }
            case wkbMultiLineString: {
                WRITE_WARNING("got wkbMultiLineString");
                break;
            }
            case wkbMultiPolygon: {
                WRITE_WARNING("got wkbMultiPolygon");
                break;
            }
            default:
                WRITE_WARNING("Unsupported shape type occurred");
            break;
        }
        */
        OGRFeature::DestroyFeature(feature);
    }
#if GDAL_VERSION_MAJOR < 2
    OGRDataSource::DestroyDataSource(ds);
#else
    GDALClose(ds);
#endif
    OCTDestroyCoordinateTransformation(toWGS84);
    OGRCleanupAll();
    return numFeatures;
#else
    UNUSED_PARAMETER(file);
    WRITE_ERROR("Cannot load shape file since SUMO was compiled without GDAL support.");
    return 0;
#endif
}

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

            if( STARTS_WITH_CI(osIngresGeomType, "I") )
                osPoint.Printf( "(%d,%d)",
                                (int) floor(poLS->getX(i)),
                                (int) floor(poLS->getY(i)) );
            else
                osPoint.Printf( "(%.15g,%.15g)",
                                poLS->getX(i), poLS->getY(i) );

            if( osPoint == osLastPoint )
            {
                CPLDebug( "INGRES",
                          "Dropping duplicate point in linestring(2).");
                continue;
            }
            osLastPoint = osPoint;

            if( osRetGeomText.size() > 1 )
                osRetGeomText += "," + osPoint;
            else
                osRetGeomText += osPoint;
        }
        osRetGeomText += ")";

        return OGRERR_NONE;
    }

/* -------------------------------------------------------------------- */
/*      Polygon                                                         */
/* -------------------------------------------------------------------- */
    if( wkbFlatten(poGeom->getGeometryType()) == wkbPolygon )
    {
        OGRPolygon *poPoly = (OGRPolygon *) poGeom;
        OGRLinearRing *poLS = poPoly->getExteriorRing();
        int i, nPoints;

        if( poLS == NULL )
            return OGRERR_FAILURE;

        if( poPoly->getNumInteriorRings() > 0 )
        {
            CPLError( CE_Warning, CPLE_AppDefined,
                      "%d inner rings discarded from polygon being converted\n"
                      "to old ingres spatial data type '%s'.",
                      poPoly->getNumInteriorRings(),
                      osIngresGeomType.c_str() );
        }

        if( EQUAL(osIngresGeomType,"POLYGON")
            && poLS->getNumPoints() > 124 )
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Attempt to place %d vertex linestring in %s field.",
                      poLS->getNumPoints(),
                      osIngresGeomType.c_str() );
            return OGRERR_FAILURE;
        }
        else if( EQUAL(osIngresGeomType,"IPOLYGON")
                 && poLS->getNumPoints() > 248 )
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Attempt to place %d vertex linestring in %s field.",
                      poLS->getNumPoints(),
                      osIngresGeomType.c_str() );
            return OGRERR_FAILURE;
        }

//.........这里部分代码省略.........
                {
                    OGRLineString *poLine = new OGRLineString;
                    int nVerticesInThisPart;

                    if( i == nParts-1 )
                        nVerticesInThisPart = nPoints - panPartStart[i];
                    else
                        nVerticesInThisPart = 
                            panPartStart[i+1] - panPartStart[i];

                    poLine->setPoints( nVerticesInThisPart, 
                                       padfX + panPartStart[i], 
                                       padfY + panPartStart[i], 
                                       padfZ + panPartStart[i] );

                    poMulti->addGeometryDirectly( poLine );
                }
            }
        } /* ARC */

/* -------------------------------------------------------------------- */
/*      Polygon                                                         */
/* -------------------------------------------------------------------- */
        else if(    nSHPType == SHPT_POLYGON
                 || nSHPType == SHPT_POLYGONZ
                 || nSHPType == SHPT_POLYGONM
                 || nSHPType == SHPT_POLYGONZM )
        {
            OGRPolygon *poMulti = new OGRPolygon;
            *ppoGeom = poMulti;

            for( i = 0; i < nParts; i++ )
            {
                OGRLinearRing *poRing = new OGRLinearRing;
                int nVerticesInThisPart;

                if( i == nParts-1 )
                    nVerticesInThisPart = nPoints - panPartStart[i];
                else
                    nVerticesInThisPart = 
                        panPartStart[i+1] - panPartStart[i];

                poRing->setPoints( nVerticesInThisPart, 
                                   padfX + panPartStart[i], 
                                   padfY + panPartStart[i], 
                                   padfZ + panPartStart[i] );

                poMulti->addRingDirectly( poRing );
            }
        } /* polygon */

/* -------------------------------------------------------------------- */
/*      Multipatch                                                      */
/* -------------------------------------------------------------------- */
        else if( bIsMultiPatch )
        {
            /* return to this later */
        } 

        CPLFree( panPartStart );
        CPLFree( padfX );
        CPLFree( padfY );
        CPLFree( padfZ );

        if( !bHasZ )
            (*ppoGeom)->setCoordinateDimension( 2 );

OGRFeature *OGRBNALayer::BuildFeatureFromBNARecord (BNARecord* record, long fid)
{
    OGRFeature  *poFeature = new OGRFeature( poFeatureDefn );
    for( int i = 0; i < nIDs; i++ )
    {
        poFeature->SetField( i, record->ids[i] ? record->ids[i] : "");
    }
    poFeature->SetFID( fid );
    if (bnaFeatureType == BNA_POINT)
    {
        poFeature->SetGeometryDirectly(
            new OGRPoint( record->tabCoords[0][0], record->tabCoords[0][1] ) );
    }
    else if (bnaFeatureType == BNA_POLYLINE)
    {
        OGRLineString* lineString = new OGRLineString ();
        lineString->setCoordinateDimension(2);
        lineString->setNumPoints(record->nCoords);
        for( int i = 0; i < record->nCoords; i++)
        {
            lineString->setPoint( i, record->tabCoords[i][0],
                                  record->tabCoords[i][1] );
        }
        poFeature->SetGeometryDirectly(lineString);
    }
    else if (bnaFeatureType == BNA_POLYGON)
    {
        double firstX = record->tabCoords[0][0];
        double firstY = record->tabCoords[0][1];
        int isFirstPolygon = 1;
        double secondaryFirstX = 0, secondaryFirstY = 0;

        OGRLinearRing* ring = new OGRLinearRing ();
        ring->setCoordinateDimension(2);
        ring->addPoint(record->tabCoords[0][0], record->tabCoords[0][1] );

        /* record->nCoords is really a safe upper bound */
        int nbPolygons = 0;
        OGRPolygon** tabPolygons =
            static_cast<OGRPolygon**>(
                CPLMalloc( record->nCoords * sizeof(OGRPolygon*) ) );

        int i = 1;
        for( ; i < record->nCoords; i++ )
        {
            ring->addPoint(record->tabCoords[i][0], record->tabCoords[i][1] );
            if (isFirstPolygon == 1 &&
                record->tabCoords[i][0] == firstX &&
                record->tabCoords[i][1] == firstY)
            {
                OGRPolygon* polygon = new OGRPolygon ();
                polygon->addRingDirectly(ring);
                tabPolygons[nbPolygons] = polygon;
                nbPolygons++;

                if (i == record->nCoords - 1)
                {
                    break;
                }

                isFirstPolygon = 0;

                i++;
                secondaryFirstX = record->tabCoords[i][0];
                secondaryFirstY = record->tabCoords[i][1];
                ring = new OGRLinearRing ();
                ring->setCoordinateDimension(2);
                ring->addPoint( record->tabCoords[i][0],
                                record->tabCoords[i][1] );
            }
            else if (isFirstPolygon == 0 &&
                    record->tabCoords[i][0] == secondaryFirstX &&
                    record->tabCoords[i][1] == secondaryFirstY)
            {

                OGRPolygon* polygon = new OGRPolygon ();
                polygon->addRingDirectly(ring);
                tabPolygons[nbPolygons] = polygon;
                nbPolygons++;

                if (i < record->nCoords - 1)
                {
                    // After the closing of a subpolygon, the first
                    // coordinates of the first polygon should be
                    // recalled... in theory
                    if (record->tabCoords[i+1][0] == firstX
                        && record->tabCoords[i+1][1] == firstY)
                    {
                        if (i + 1 == record->nCoords - 1)
                            break;
                        i ++;
                    }
                    else
                    {
#if 0
                        CPLError(CE_Warning, CPLE_AppDefined,
                                 "Geometry of polygon of fid %d starting at "
                                 "line %d is not strictly conformant. "
                                 "Trying to go on...\n",
                                 fid,
//.........这里部分代码省略.........

OGRErr OGRBNALayer::ICreateFeature( OGRFeature *poFeature )

{
    OGRGeometry     *poGeom = poFeature->GetGeometryRef();
    char eol[3];
    const char* partialEol = (poDS->GetMultiLine()) ? eol : poDS->GetCoordinateSeparator();

    if (poGeom == NULL || poGeom->IsEmpty() )
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "OGR BNA driver cannot write features with empty geometries.");
        return OGRERR_FAILURE;
    }

    if (poDS->GetUseCRLF())
    {
        eol[0] = 13;
        eol[1] = 10;
        eol[2] = 0;
    }
    else
    {
        eol[0] = 10;
        eol[1] = 0;
    }

    if ( ! bWriter )
    {
        return OGRERR_FAILURE;
    }

    if( poFeature->GetFID() == OGRNullFID )
        poFeature->SetFID( nFeatures++ );

    VSILFILE* fp = poDS->GetOutputFP();
    int nbPairPerLine = poDS->GetNbPairPerLine();

    switch( poGeom->getGeometryType() )
    {
        case wkbPoint:
        case wkbPoint25D:
        {
            OGRPoint* point = (OGRPoint*)poGeom;
            WriteFeatureAttributes(fp, poFeature);
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "1"));
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", partialEol));
            WriteCoord(fp, point->getX(), point->getY());
            CPL_IGNORE_RET_VAL(VSIFPrintfL( fp, "%s", eol));
            break;
        }

        case wkbPolygon:
        case wkbPolygon25D:
        {
            OGRPolygon* polygon = (OGRPolygon*)poGeom;
            OGRLinearRing* ring = polygon->getExteriorRing();
            if (ring == NULL)
            {
                return OGRERR_FAILURE;
            }

            double firstX = ring->getX(0);
            double firstY = ring->getY(0);
            int nBNAPoints = ring->getNumPoints();
            int is_ellipse = FALSE;

            /* This code tries to detect an ellipse in a polygon geometry */
            /* This will only work presumably on ellipses already read from a BNA file */
            /* Mostly a BNA to BNA feature... */
            if (poDS->GetEllipsesAsEllipses() &&
                polygon->getNumInteriorRings() == 0 &&
                nBNAPoints == 361)
            {
                double oppositeX = ring->getX(180);
                double oppositeY = ring->getY(180);
                double quarterX = ring->getX(90);
                double quarterY = ring->getY(90);
                double antiquarterX = ring->getX(270);
                double antiquarterY = ring->getY(270);
                double center1X = 0.5*(firstX + oppositeX);
                double center1Y = 0.5*(firstY + oppositeY);
                double center2X = 0.5*(quarterX + antiquarterX);
                double center2Y = 0.5*(quarterY + antiquarterY);
                if (fabs(center1X - center2X) < 1e-5 && fabs(center1Y - center2Y) < 1e-5 &&
                    fabs(oppositeY - firstY) < 1e-5 &&
                    fabs(quarterX - antiquarterX) < 1e-5)
                {
                    double major_radius = fabs(firstX - center1X);
                    double minor_radius = fabs(quarterY - center1Y);
                    is_ellipse = TRUE;
                    for(int i=0;i<360;i++)
                    {
                        if (!(fabs(center1X + major_radius * cos(i * (M_PI / 180)) - ring->getX(i)) < 1e-5 &&
                              fabs(center1Y + minor_radius * sin(i * (M_PI / 180)) - ring->getY(i)) < 1e-5))
                        {
                            is_ellipse = FALSE;
                            break;
                        }
                    }
                    if ( is_ellipse == TRUE )
//.........这里部分代码省略.........

OGRGeometry *OGRGRASSLayer::GetFeatureGeometry ( long nFeatureId, int *cat )
{
    CPLDebug ( "GRASS", "OGRGRASSLayer::GetFeatureGeometry nFeatureId = %d", nFeatureId );

    int cidx = paFeatureIndex[(int)nFeatureId];

    int type, id;
    Vect_cidx_get_cat_by_index ( poMap, iLayerIndex, cidx, cat, &type, &id );

    //CPLDebug ( "GRASS", "cat = %d type = %d id = %d", *cat, type, id );

    OGRGeometry *poOGR = NULL;
    int bIs3D = Vect_is_3d(poMap);

    switch ( type ) {
	case GV_POINT:
        {
	    Vect_read_line ( poMap, poPoints, poCats, id);
            if (bIs3D)
                poOGR = new OGRPoint( poPoints->x[0], poPoints->y[0], poPoints->z[0] );
            else
                poOGR = new OGRPoint( poPoints->x[0], poPoints->y[0] );
        }
        break;
	    
	case GV_LINE:
	case GV_BOUNDARY:
        {
	    Vect_read_line ( poMap, poPoints, poCats, id);
	    OGRLineString *poOGRLine = new OGRLineString();
            if (bIs3D)
                poOGRLine->setPoints( poPoints->n_points, 
                                      poPoints->x, poPoints->y, poPoints->z );
            else
                poOGRLine->setPoints( poPoints->n_points, 
                                      poPoints->x, poPoints->y );

            poOGR = poOGRLine;
        }
        break;

	case GV_AREA:
        {
	    Vect_get_area_points ( poMap, id, poPoints );
	    
	    OGRPolygon 		*poOGRPoly;
	    poOGRPoly = new OGRPolygon();

	    OGRLinearRing       *poRing;
	    poRing = new OGRLinearRing();
            if (bIs3D)
                poRing->setPoints( poPoints->n_points,
                                poPoints->x, poPoints->y, poPoints->z );
            else
                poRing->setPoints( poPoints->n_points,
                                poPoints->x, poPoints->y ); 

	    poOGRPoly->addRingDirectly( poRing );

	    // Islands
	    int nisles = Vect_get_area_num_isles ( poMap, id );
	    for ( int i = 0; i < nisles; i++ ) {
		int isle =  Vect_get_area_isle ( poMap, id, i );
		Vect_get_isle_points ( poMap, isle, poPoints );

		poRing = new OGRLinearRing();
                if (bIs3D)
                    poRing->setPoints( poPoints->n_points,
                                    poPoints->x, poPoints->y, poPoints->z );
                else
                    poRing->setPoints( poPoints->n_points,
                                    poPoints->x, poPoints->y );

		poOGRPoly->addRingDirectly( poRing );
	    }
	    
	    poOGR = poOGRPoly;
        }   
        break;

	default: // Should not happen
        {
	    CPLError( CE_Failure, CPLE_AppDefined, "Unknown GRASS feature type.");
	    return NULL;
        }
    }
	    
    return poOGR;
}