shape2ogr.cpp

用于读取TAB、MIF、SHP文件的类

    dx1 = x[next] - x[v];
    dy1 = y[next] - y[v];

    area2 = dx1 * dy0 - dx0 * dy1;
    
    if ( area2 > 0 )      /* CCW */
	return -1;
    else if ( area2 < 0 )  /* CW */
	return 1;
    
    return 0; /* error */
}

/************************************************************************/
/*                          PointInRing()                               */
/*                                                                      */
/*      return: 1 point is inside the ring                              */
/*              0 point is outside the ring                             */
/*                                                                      */
/*              for point exactly on the boundary it returns 0 or 1     */
/************************************************************************/

int PointInRing ( SHPObject *Shape, int ring, double x, double y ) 
{
    int    i, start, end, c = 0;
    double *xp, *yp; 

    if ( ring < 0 || ring >= Shape->nParts ) return ( 0 );

    xp = Shape->padfX;
    yp = Shape->padfY;

    RingStartEnd ( Shape, ring, &start, &end );

    /* This code was originaly written by Randolph Franklin, here it is slightly modified. */
    for (i = start; i < end; i++) {
        if ( ( ( ( yp[i] <= y ) && ( y < yp[i+1] ) ) || 
               ( ( yp[i+1] <= y ) && ( y < yp[i] ) ) ) &&
             ( x < (xp[i+1] - xp[i]) * (y - yp[i]) / (yp[i+1] - yp[i]) + xp[i] )
	     ) 
	{
            c = !c;
	}
    }
    return c;
}

/************************************************************************/
/*                          RingInRing()                                */
/*                                                                      */
/* Checks point by point using PointInRing if oRing contains iRing      */
/*                                                                      */
/*      return: 1 oRing contains iRing                                  */
/*              0 iRing is not contained by oRing                       */
/*                                                                      */
/*              for point exactly on the boundary it returns 0 or 1     */
/************************************************************************/
int RingInRing( SHPObject *Shape, int oRing, int iRing )
{
    int iRingStart, iRingEnd;
    RingStartEnd ( Shape, iRing, &iRingStart, &iRingEnd );
    while( iRingStart < iRingEnd )
    {
        if ( !PointInRing( Shape, oRing, Shape->padfX[iRingStart], Shape->padfY[iRingStart] ) )
            return 0;
        ++iRingStart;
    }
    return 1;
}
    
/************************************************************************/
/*                        CreateLinearRing                              */
/*                                                                      */
/************************************************************************/
OGRLinearRing * CreateLinearRing ( SHPObject *psShape, int ring )
{
    OGRLinearRing *poRing;
    int nRingStart, nRingEnd, nRingPoints;

    poRing = new OGRLinearRing();

    RingStartEnd ( psShape, ring, &nRingStart, &nRingEnd );

    nRingPoints = nRingEnd - nRingStart + 1;

    poRing->setPoints( nRingPoints, psShape->padfX + nRingStart, 
	       psShape->padfY + nRingStart, psShape->padfZ + nRingStart );

    return ( poRing );
}
    
/************************************************************************/
/*                          SHPReadOGRObject()                          */
/*                                                                      */
/*      Read an item in a shapefile, and translate to OGR geometry      */
/*      representation.                                                 */
/************************************************************************/

OGRGeometry *SHPReadOGRObject( SHPHandle hSHP, int iShape )
{
    // CPLDebug( "Shape", "SHPReadOGRObject( iShape=%d )\n", iShape );

    SHPObject   *psShape;
    OGRGeometry *poOGR = NULL;

    psShape = SHPReadObject( hSHP, iShape );

    if( psShape == NULL )
    {
        return NULL;
    }

/* -------------------------------------------------------------------- */
/*      Point.                                                          */
/* -------------------------------------------------------------------- */
    else if( psShape->nSHPType == SHPT_POINT
             || psShape->nSHPType == SHPT_POINTM
             || psShape->nSHPType == SHPT_POINTZ )
    {
        poOGR = new OGRPoint( psShape->padfX[0], psShape->padfY[0],
                              psShape->padfZ[0] );

        if( psShape->nSHPType == SHPT_POINT )
        {
            poOGR->setCoordinateDimension( 2 );
        }
    }

/* -------------------------------------------------------------------- */
/*      Multipoint.                                                     */
/* -------------------------------------------------------------------- */
    else if( psShape->nSHPType == SHPT_MULTIPOINT
             || psShape->nSHPType == SHPT_MULTIPOINTM
             || psShape->nSHPType == SHPT_MULTIPOINTZ )
    {
        OGRMultiPoint *poOGRMPoint = new OGRMultiPoint();
        int             i;

        for( i = 0; i < psShape->nVertices; i++ )
        {
            OGRPoint    *poPoint;

            poPoint = new OGRPoint( psShape->padfX[i], psShape->padfY[i],
                                    psShape->padfZ[i] );
            
            poOGRMPoint->addGeometry( poPoint );

            delete poPoint;
        }
        
        poOGR = poOGRMPoint;

        if( psShape->nSHPType == SHPT_MULTIPOINT )
            poOGR->setCoordinateDimension( 2 );
    }

/* -------------------------------------------------------------------- */
/*      Arc (LineString)                                                */
/*                                                                      */
/*      I am ignoring parts though they can apply to arcs as well.      */
/* -------------------------------------------------------------------- */
    else if( psShape->nSHPType == SHPT_ARC
             || psShape->nSHPType == SHPT_ARCM
             || psShape->nSHPType == SHPT_ARCZ )
    {
        if( psShape->nParts == 1 )
        {
            OGRLineString *poOGRLine = new OGRLineString();

            poOGRLine->setPoints( psShape->nVertices,
                                  psShape->padfX, psShape->padfY, psShape->padfZ );
        
            poOGR = poOGRLine;
        }
        else
        {
            int iRing;
            OGRMultiLineString *poOGRMulti;
        
            poOGR = poOGRMulti = new OGRMultiLineString();
            
            for( iRing = 0; iRing < psShape->nParts; iRing++ )
            {
                OGRLineString   *poLine;
                int     nRingPoints;
                int     nRingStart;

                poLine = new OGRLineString();

                if( psShape->panPartStart == NULL )
                {
                    nRingPoints = psShape->nVertices;
                    nRingStart = 0;
                }
                else
                {

                    if( iRing == psShape->nParts - 1 )
                        nRingPoints =
                            psShape->nVertices - psShape->panPartStart[iRing];
                    else
                        nRingPoints = psShape->panPartStart[iRing+1]
                            - psShape->panPartStart[iRing];
                    nRingStart = psShape->panPartStart[iRing];
                }
            
                poLine->setPoints( nRingPoints, 
                                   psShape->padfX + nRingStart,
                                   psShape->padfY + nRingStart,
                                   psShape->padfZ + nRingStart );

                poOGRMulti->addGeometryDirectly( poLine );
            }
        }

        if( psShape->nSHPType == SHPT_ARC )
            poOGR->setCoordinateDimension( 2 );
    }

/* -------------------------------------------------------------------- */
/*      Polygon                                                         */
/*                                                                      */
/* As for now Z coordinate is not handled correctly                     */
/* -------------------------------------------------------------------- */
    else if( psShape->nSHPType == SHPT_POLYGON
             || psShape->nSHPType == SHPT_POLYGONM
             || psShape->nSHPType == SHPT_POLYGONZ )
    {
        int iRing;
        
        //CPLDebug( "Shape", "Shape type: polygon with nParts=%d \n", psShape->nParts );

        if ( psShape->nParts == 1 )
        {
            /* Surely outer ring */
            OGRPolygon *poOGRPoly = NULL;
            OGRLinearRing *poRing = NULL;

            poOGR = poOGRPoly = new OGRPolygon();
            poRing = CreateLinearRing ( psShape, 0 );
            poOGRPoly->addRingDirectly( poRing );
        }
        else
        {
            /* Multipart polygon. */

            int nOuter = 0;         /* Number of outer rings */ 
            int *direction = NULL;  /* ring direction (1 CW,outer, -1 CCW, inner) */ 
            int *outer = NULL;      /* list of outer rings */ 
                
            /* Outer ring index for inner rings, -1 for outer rings */ 
            int *outside = NULL;   

            direction = (int *) CPLMalloc( psShape->nParts * sizeof(int) );
            outer = (int *) CPLMalloc( psShape->nParts * sizeof(int) );
            outside = (int *) CPLMalloc( psShape->nParts * sizeof(int) );
            
            /* First distinguish outer from inner rings */
            for( iRing = 0; iRing < psShape->nParts; iRing++ )
            {
                direction[iRing] = RingDirection ( psShape, iRing);
                if ( direction[iRing] == 1 )
                {
                    outer[nOuter] = iRing;
                    nOuter++;
                }
                outside[iRing] = -1;
            }

            if ( nOuter == 1 )
            { 
                /* One outer ring?
                 * we do not need to check anything as we assume that shapefile is valid
                 * ie. if there is one outer ring then all inner rings are inside.
                 */
                for( iRing = 0; iRing < psShape->nParts; iRing++ )
                {
                    if ( direction[iRing] == -1 )
                    {
                        outside[iRing] = outer[ 0 ];
                    }
                }
            }
            else
            {
                /* Calculate ring extents */
                RingExtent *extents = new RingExtent[ psShape->nParts ];

                for( iRing = 0; iRing < psShape->nParts; iRing++ )
                {
                    int start = 0;
                    int end = 0;
                    RingStartEnd ( psShape, iRing, &start, &end );
                    extents[ iRing ].calculate( end-start, psShape->padfX + start, psShape->padfY + start );
                }

                /* Outer ring index */
                int oRing;

                /* This loops tries to assign inner to outer rings.
                 * The fundamental assumption is that if the extent of an inner ring is contained
                 * in the extent of exactly one outer ring then this inner ring is inside the outer ring.
                 * Otherwise that shapefile is broken.
                 *
                 * XXX - mloskot - This loop does incorrect classification, see Bug 1217
                 */
                for( iRing = 0; iRing < psShape->nParts; iRing++ )
                {
                    if ( direction[ iRing ] != 1 )
                    {
                        /* Inner ring */

                        /* Find the outer ring for iRing */
                        for( oRing = 0; oRing < nOuter; oRing++ )
                        {
                            if ( extents[ outer[ oRing ] ].contains( extents[ iRing ] ) )
                            {
                                /* The outer ring contains the inner ring, we have to execute
                                 * some additional tests.
                                 */ 
                                if ( outside[ iRing ] == -1 )
                                {
                                    /* this is the first outer ring, assume it containes this inner iRing */
                                    outside[ iRing ] = outer[ oRing ];
                                }
                                else
                                { 
                                    /* This is another outer ring, which contains iRing, have to distinguish
                                     * between this and previous outer rings using RingInPoint.
                                     * Here is a place for some optimization as sometimes we may check
                                     * the same ring many times.
                                     */
                                    if ( !RingInRing( psShape, outside[ iRing ], iRing ) )
                                    {
                                        outside[ iRing ] = outer[ oRing ];
                                    }
                                }                            
                            }
                        }
                    }
                }

                /* Destroy ring extents object. */
                delete[] extents;
            }

            /* The inner rings are preliminary sorted, but:
             * - Some of them are not assigned to any outer ring. Probably broken shapefile
             * - Some inner rings are assigned to outer rings using extents only.
             * Additional geometry tests are ommited here due to perfomance penalties.
             * Promote any unassigned inner rings to be outside rings.
             */
            for ( iRing = 0; iRing < psShape->nParts; iRing++ )
            {
                /* Outer rings */
                if( direction[iRing] != 1 && outside[iRing] == -1 )
                {
                    direction[iRing] = 1; /* this isn't exactly true! */
                    outer[nOuter++] = iRing;
                }
            }

            if ( nOuter == 1 )
            {
                /* One outer ring and more inner rings */
                OGRPolygon    *poOGRPoly = NULL;
                OGRLinearRing *poRing = NULL;

                /* Outer */
                poOGR = poOGRPoly = new OGRPolygon();
                poRing = CreateLinearRing ( psShape, outer[0] );
                poOGRPoly->addRingDirectly( poRing );

                /* Inner */
                for( iRing = 0; iRing < psShape->nParts; iRing++ )
                {
                    if ( direction[iRing] == -1 )
                    {
                        poRing = CreateLinearRing ( psShape, iRing );