shape2ogr.cpp

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

/******************************************************************************
 * $Id: shape2ogr.cpp,v 1.48 2006/09/21 20:58:47 fwarmerdam Exp $
 *
 * 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.
 ******************************************************************************
 *
 * $Log: shape2ogr.cpp,v $
 * Revision 1.48  2006/09/21 20:58:47  fwarmerdam
 * Comment out noisy debug message.
 *
 * Revision 1.47  2006/09/02 19:56:20  mloskot
 * Re-enabled reading NULL geometries from shapefile
 * (was broken in shape2ogr.cpp, rev 1.46)
 *
 * Revision 1.46  2006/08/28 14:00:02  mloskot
 * Added stronger test of Shapefile reading failures, e.g. truncated files.
 * The problem was discovered by Tim Sutton and reported here https://svn.qgis.org/trac/ticket/200
 *
 * Revision 1.45  2006/07/20 12:50:52  mloskot
 * Refactored shape2ogr.cpp file. Added comments related to Bug 1217 (not fixed yet, see the bug report).
 *
 * Revision 1.44  2006/06/30 09:20:20  mloskot
 * Fixed null dates issue reported by Aaron Koning. Fixed zero-fill of year value in OGRFeature::GetFieldAsString. Added assertions and NULL pointer tests.
 *
 * Revision 1.43  2006/06/21 20:43:02  fwarmerdam
 * support datasets without dbf: bug 1211
 *
 * Revision 1.42  2006/04/02 18:46:54  fwarmerdam
 * updated date support
 *
 * Revision 1.41  2006/02/15 04:26:55  fwarmerdam
 * added date support
 *
 * Revision 1.40  2006/01/05 02:15:39  fwarmerdam
 * implement DeleteFeature support
 *
 * Revision 1.39  2005/10/12 14:50:28  fwarmerdam
 * recover gracefully from POLYGON EMPTY on write
 *
 * Revision 1.38  2005/09/30 23:53:35  mbrudka
 * Fixed bug 940: http://bugzilla.remotesensing.org/show_bug.cgi?id=940
 *
 * Revision 1.37  2005/09/21 00:53:20  fwarmerdam
 * fixup OGRFeatureDefn and OGRSpatialReference refcount handling
 *
 * Revision 1.36  2005/09/20 21:27:25  mbrudka
 * Some optimizations in reading multi-ring polygons.
 *
 * Revision 1.35  2005/08/04 15:31:09  fwarmerdam
 * Added additional patch for computing ring direction in cases of
 * vertices that are very near or coincident.
 *
 * Revision 1.34  2005/08/04 15:10:46  fwarmerdam
 * Committed fix to avoid infinite loop on degenerate shapefiles.  Patch
 * provided by Baumann Konstantin.  See Var2_STRASSE.shp.
 *
 * Revision 1.33  2005/07/20 01:44:25  fwarmerdam
 * try and preserve geometry dimension properly
 *
 * Revision 1.32  2004/07/06 19:35:51  warmerda
 * Default to using FTString for unrecognised field types like logical.
 *
 * Revision 1.31  2004/01/16 22:40:40  warmerda
 * any inner rings not assigned to a polygo are promoted to be outer rings
 *
 * Revision 1.30  2003/12/11 19:19:40  warmerda
 * Fixed leak of geometry objects when building multipolygons.
 *
 * Revision 1.29  2003/11/09 18:54:53  warmerda
 * added auto-FID field creation if new file has no fields on first feature wrt
 *
 * Revision 1.28  2003/10/18 19:01:10  warmerda
 * added Radims patch to recognise multipolygons on read properly - bug 213
 *
 * Revision 1.27  2003/06/10 14:44:11  warmerda
 * Added support for writing shapes with NULL geometries.
 *
 * Revision 1.26  2003/05/27 21:39:53  warmerda
 * added support for writing MULTILINESTRINGs as ARCs
 *
 * Revision 1.25  2003/05/21 04:03:54  warmerda
 * expand tabs
 *
 * Revision 1.24  2003/02/25 14:35:11  warmerda
 * Incorporated more correct support for reading multi-part arcs as
 * MultiLineStrings.
 *
 * Revision 1.23  2002/10/04 04:31:40  warmerda
 * Fixes for writing 3D multipoints.
 *
 * Revision 1.22  2002/10/04 04:30:46  warmerda
 * Fixed some bugs in 3D support for writing polygons.
 *
 * Revision 1.21  2002/08/15 15:10:01  warmerda
 * Fixed problem when GetFeature() called for a non-existant id.
 * http://bugzilla.remotesensing.org/show_bug.cgi?id=175
 *
 * Revision 1.20  2002/08/14 23:27:25  warmerda
 * Fixed reading of Z coordinate for SHPT_ARCZ type.
 *
 * Revision 1.19  2002/04/17 21:49:02  warmerda
 * Fixed a bug writing arcs with Z coordinates.
 *
 * Revision 1.18  2002/04/17 20:03:27  warmerda
 * Try to preserve Z values on read.
 *
 * Revision 1.17  2002/04/17 20:01:39  warmerda
 * Ensure Z coordinate preserved on read.
 *
 * Revision 1.16  2002/04/16 16:21:28  warmerda
 * Ensure polygons written with correct winding.
 *
 * Revision 1.15  2002/03/27 21:04:38  warmerda
 * Added support for reading, and creating lone .dbf files for wkbNone geometry
 * layers.  Added support for creating a single .shp file instead of a directory
 * if a path ending in .shp is passed to the data source create method.
 *
 * Revision 1.14  2001/12/06 18:14:22  warmerda
 * handle case where there is no DBF file
 *
 * Revision 1.13  2001/09/27 14:53:31  warmerda
 * added proper 3D support for all element types
 *
 * Revision 1.12  2001/09/04 15:39:07  warmerda
 * tightened up NULL attribute handling
 *
 * Revision 1.11  2001/07/18 04:55:16  warmerda
 * added CPL_CSVID
 *
 * Revision 1.10  2001/02/06 14:30:15  warmerda
 * don't crash on attempt to write geometryless features
 *
 * Revision 1.9  2000/11/29 14:09:36  warmerda
 * Extended shapefile OGR driver to support MULTIPOLYGON objects.
 *
 * Revision 1.8  2000/01/26 22:05:45  warmerda
 * fixed bug with 2d/3d arcs
 *
 * Revision 1.7  1999/12/23 14:51:07  warmerda
 * Improved support in face of 3D geometries.
 *
 * Revision 1.6  1999/11/04 21:18:01  warmerda
 * fix bug with shapeid for new shapes
 *
 * Revision 1.5  1999/07/27 01:52:04  warmerda
 * added support for writing polygon/multipoint/arc files
 *
 * Revision 1.4  1999/07/27 00:52:40  warmerda
 * added write methods
 *
 * Revision 1.3  1999/07/26 13:59:25  warmerda
 * added feature writing api
 *
 * Revision 1.2  1999/07/12 15:39:18  warmerda
 * added setting of shape geometry
 *
 * Revision 1.1  1999/07/05 18:58:07  warmerda
 * New
 *
 * Revision 1.3  1999/06/11 19:22:13  warmerda
 * added OGRFeature and OGRFeatureDefn related translation
 *
 * Revision 1.2  1999/05/31 20:41:49  warmerda
 * added multipoint support
 *
 * Revision 1.1  1999/05/31 19:26:03  warmerda
 * New
 */

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

CPL_CVSID("$Id: shape2ogr.cpp,v 1.48 2006/09/21 20:58:47 fwarmerdam Exp $");

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++;
    }