shape2ogr.cpp

来自「支持各种栅格图像和矢量图像读取的库」· C++ 代码 · 共 1,356 行 · 第 1/4 页

/******************************************************************************
 * $Id: shape2ogr.cpp 10646 2007-01-18 02:38:10Z warmerdam $
 *
 * Project:  OpenGIS Simple Features Reference Implementation
 * Purpose:  Implements translation of Shapefile shapes into OGR
 *           representation.
 * Author:   Frank Warmerdam, [e-mail]
 *
 ******************************************************************************
 * Copyright (c) 1999,  Les Technologies SoftMap Inc.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/

#include "ogrshape.h"
#include "cpl_conv.h"

CPL_CVSID("$Id: shape2ogr.cpp 10646 2007-01-18 02:38:10Z warmerdam $");

static const double EPSILON = 1E-5;

/************************************************************************/
/*                            epsilonEqual()                            */
/************************************************************************/

static inline bool epsilonEqual(double a, double b, double eps) 
{
    return (::fabs(a - b) < eps);
}

/* A helper class which represents 2D bounding box */
class RingExtent
{
public:
    RingExtent();
    /* Calculates bounding box for given set of points */
    void calculate( int ringParts, double *ringX, double *ringY );
    /* Does this extent contains rhs*/
    bool contains( const RingExtent &rhs );
    /* Does this extent contains given point*/
    bool contains( double x, double y );
private:
    bool empty; /* extent is empty */
    double xMin, yMin; /* lower left vertex of the bounding box */
    double xMax, yMax; /* upper right vertex of the bounding box */
};

RingExtent::RingExtent():
    empty( true ),
    xMin( 0 ), //nan ?
    yMin( 0 ), 
    xMax( 0 ), 
    yMax( 0 )
{
}

void RingExtent::calculate( int ringParts, double *ringX, double *ringY )
{
    empty = ringParts <= 0;
    if ( !empty )
    {
        xMin = xMax = *ringX;
        yMin = yMax = *ringY;
        while( --ringParts > 0 )
        {
            if ( xMin > *ringX )
                xMin = *ringX;
            else if ( xMax < *ringX )
                xMax = *ringX;
            if ( yMin > *ringY )
               yMin = *ringY;
            else if ( yMax < *ringY )
                yMax = *ringY;
            ++ringX; 
            ++ringY; 
        }
    }
}

inline bool RingExtent::contains( double x, double y )
{
    if ( empty )
        return false;

    return xMin <= x && x <= xMax &&
        yMin <= y && y <= yMax;
}

inline bool RingExtent::contains( const RingExtent &extent )
{
    return contains( extent.xMin, extent.yMin ) &&
        contains( extent.xMax, extent.yMax );
}


/************************************************************************/
/*                        RingStartEnd                                  */
/*        set first and last vertex for given ring                      */
/************************************************************************/
void RingStartEnd ( SHPObject *psShape, int ring, int *start, int *end )
{
    if( psShape->panPartStart == NULL )
    {
	    *start = 0;
        *end = psShape->nVertices - 1;
    }
    else
    {
        if( ring == psShape->nParts - 1 )
            *end = psShape->nVertices - 1;
        else
            *end = psShape->panPartStart[ring+1] - 1;

        *start = psShape->panPartStart[ring];
    }
}
    
/************************************************************************/
/*                          RingDirection()                             */
/*                                                                      */
/*      return: 1 CW (clockwise)                                        */
/*             -1 CCW (counterclockwise)                                */
/*              0 error                                                 */
/************************************************************************/

int RingDirection ( SHPObject *Shape, int ring ) 
{
    int    i, start, end, v, next;
    double *x, *y, dx0, dy0, dx1, dy1, area2; 

    /* Note: first and last coordinate MUST be identical according to shapefile specification */
    if ( ring < 0 || ring >= Shape->nParts ) return ( 0 );

    x = Shape->padfX;
    y = Shape->padfY;

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

    if ( start == end ) return ( 0 ); // empty ring!?!?

    /* Find the lowest rightmost vertex */
    v = start;  
    for ( i = start + 1; i < end; i++ )
    {
        /* => v < end */
        if ( y[i] < y[v] || ( y[i] == y[v] && x[i] > x[v] ) )
        {
            v = i;
        }
    }
    
    /* Vertices may be duplicate, we have to go to nearest different in each direction */
    /* preceding */
    next = v - 1;
    while ( 1 )
    {
        if ( next < start ) 
        {
            next = end - 1; 
        }

        if( !epsilonEqual(x[next], x[v], EPSILON) 
            || !epsilonEqual(y[next], y[v], EPSILON) )
        {
            break;
        }

        if ( next == v ) /* So we cannot get into endless loop */
        {
            break;
        }

        next--;
    }
	    
    dx0 = x[next] - x[v];
    dy0 = y[next] - y[v];
    
    
    /* following */
    next = v + 1;
    while ( 1 )
    {
        if ( next >= end ) 
        {
            next = start; 
        }

        if ( !epsilonEqual(x[next], x[v], EPSILON) 
             || !epsilonEqual(y[next], y[v], EPSILON) )
        {
            break;
        }

        if ( next == v ) /* So we cannot get into endless loop */
        {
            break;
        }

        next++;
    }

    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] );