/****************************************************************************** * $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, warmerda@home.com * ****************************************************************************** * 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] );