shape2ogr.cpp

来自「GIS系统支持库Geospatial Data Abstraction Libr」· C++ 代码 · 共 1,155 行 · 第 1/3 页

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

/* -------------------------------------------------------------------- */
/*      Polygon                                                         */
/*                                                                      */
/*      For now we assume the first ring is an outer ring, and          */
/*      everything else is an inner ring.  This must smarten up in      */
/*      the future.                                                     */
/* -------------------------------------------------------------------- */
    else if( psShape->nSHPType == SHPT_POLYGON
             || psShape->nSHPType == SHPT_POLYGONM
             || psShape->nSHPType == SHPT_POLYGONZ )
    {
        int     iRing;
        
	if ( psShape->nParts == 1 ) { /* Surely outer ring */
            OGRPolygon *poOGRPoly;
	    OGRLinearRing       *poRing;

	    poOGR = poOGRPoly = new OGRPolygon();
	    poRing = CreateLinearRing ( psShape, 0 );
	    poOGRPoly->addRingDirectly( poRing );
	} else {
	    int nOuter = 0; /* Number of outer rings */ 
	    int *direction; /* ring direction (1 CW,outer, -1 CCW, inner) */ 
	    int *outer;     /* list of outer rings */ 
	    int *outside;   /* outer ring index for inner rings, -1 for outer rings */ 

	    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;
	    }
	    
	    /* Find for each inner ring outer ring it is within */

	    /* Note: According to specification, rings can touch by one vertex, but not along segment. */
	    /* It may happen that inner ring touches the outer and PointInRing() retuirns 0 (outside) */
	    /* because of representation error. Because of that, cycle through all vertices of inner */
	    /* ring is used, until outer ring is found. Of course, this may fail as well, but probability */
	    /* is lower. We could use also points on inner segments or centroid. */ 
	    /* In theory, it may also happen, that vertex touching outer ring falls outside its outer */
	    /* ring and inside another (incorrect). */ 
	    /* The problem is how to find a point on inner ring which doesn't lie on outer ring. */
	    for( iRing = 0; iRing < psShape->nParts; iRing++ ) { /* cycle through inner rings */
		int  start, end;
		
		if ( direction[iRing] != -1 ) /* is not inner ring */
		    continue;
		
		RingStartEnd ( psShape, iRing, &start, &end );

		for ( int vert = start; vert < end; vert++ ) { /* cycle through inner ring's vertices */
		    double x, y;
		    
		    x = psShape->padfX[vert];
		    y = psShape->padfY[vert];

		    for( int oRing = 0; oRing < psShape->nParts; oRing++ ) { /* outer rings */
			int inside;
			    
			if ( direction[oRing] != 1 ) /* not outer */ 
			    continue;
			    
			inside = PointInRing ( psShape, oRing, x, y);

			if ( inside ) {
			    outside[iRing] = oRing;
			    break;
			}
		    }
		    if ( outside[iRing] >= 0 ) /* outer ring found */
			break;
	        }
	    }

            /* 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;
	        OGRLinearRing *poRing;

		/* 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 );
			poOGRPoly->addRingDirectly( poRing );
		    }
		}
	    } else { 
		OGRGeometryCollection *poOGRGeCo;

		poOGR = poOGRGeCo = new OGRMultiPolygon();

	        for ( iRing = 0; iRing < nOuter; iRing++ ) { /* outer rings */
		    int oRing; 

		    OGRPolygon    *poOGRPoly;
		    OGRLinearRing *poRing;

		    oRing = outer[iRing];

		    /* Outer */
		    poOGRPoly = new OGRPolygon();
		    poRing = CreateLinearRing ( psShape, oRing );
		    poOGRPoly->addRingDirectly( poRing );

		    /* Inner */
		    for( int iRing2 = 0; iRing2 < psShape->nParts; iRing2++ ) {
			if ( outside[iRing2] == oRing ) {
			    poRing = CreateLinearRing ( psShape, iRing2 );
			    poOGRPoly->addRingDirectly( poRing );
			}
		    }
                    poOGRGeCo->addGeometryDirectly (poOGRPoly);
		}
	    }
            CPLFree( direction );
            CPLFree( outer );
            CPLFree( outside );
	}
    }

/* -------------------------------------------------------------------- */
/*      Otherwise for now we just ignore the object.  Eventually we     */
/*      should implement multipatch.                                    */
/* -------------------------------------------------------------------- */
    else
    {
        if( psShape->nSHPType != SHPT_NULL )
            CPLDebug( "OGR", "Unsupported shape type in SHPReadOGRObject()" );
        /* nothing returned */
    }
    
/* -------------------------------------------------------------------- */
/*      Cleanup shape, and set feature id.                              */
/* -------------------------------------------------------------------- */
    SHPDestroyObject( psShape );

    return poOGR;
}

/************************************************************************/
/*                         SHPWriteOGRObject()                          */
/************************************************************************/

OGRErr SHPWriteOGRObject( SHPHandle hSHP, int iShape, OGRGeometry *poGeom )

{
/* ==================================================================== */
/*      Write "shape" with no geometry.                                 */
/* ==================================================================== */
    if( poGeom == NULL )
    {
        SHPObject       *psShape;

        psShape = SHPCreateSimpleObject( SHPT_NULL, 0, NULL, NULL, NULL );
        SHPWriteObject( hSHP, iShape, psShape );
        SHPDestroyObject( psShape );
    }

/* ==================================================================== */
/*      Write point geometry.                                           */
/* ==================================================================== */
    else if( hSHP->nShapeType == SHPT_POINT
             || hSHP->nShapeType == SHPT_POINTM
             || hSHP->nShapeType == SHPT_POINTZ )
    {
        SHPObject       *psShape;
        OGRPoint        *poPoint = (OGRPoint *) poGeom;
        double          dfX, dfY, dfZ = 0;

        if( poGeom->getGeometryType() != wkbPoint
            && poGeom->getGeometryType() != wkbPoint25D )        
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Attempt to write non-point (%s) geometry to"
                      " point shapefile.",
                      poGeom->getGeometryName() );

            return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
        }

        dfX = poPoint->getX();
        dfY = poPoint->getY();
        dfZ = poPoint->getZ();
        
        psShape = SHPCreateSimpleObject( hSHP->nShapeType, 1,
                                         &dfX, &dfY, &dfZ );
        SHPWriteObject( hSHP, iShape, psShape );
        SHPDestroyObject( psShape );
    }
/* ==================================================================== */
/*      MultiPoint.                                                     */
/* ==================================================================== */
    else if( hSHP->nShapeType == SHPT_MULTIPOINT
             || hSHP->nShapeType == SHPT_MULTIPOINTM
             || hSHP->nShapeType == SHPT_MULTIPOINTZ )
    {
        OGRMultiPoint   *poMP = (OGRMultiPoint *) poGeom;
        double          *padfX, *padfY, *padfZ;
        int             iPoint;
        SHPObject       *psShape;

        if( wkbFlatten(poGeom->getGeometryType()) != wkbMultiPoint )
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Attempt to write non-multipoint (%s) geometry to "
                      "multipoint shapefile.",
                      poGeom->getGeometryName() );

            return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
        }

        padfX = (double *) CPLMalloc(sizeof(double)*poMP->getNumGeometries());
        padfY = (double *) CPLMalloc(sizeof(double)*poMP->getNumGeometries());
        padfZ = (double *) CPLCalloc(sizeof(double),poMP->getNumGeometries());

        for( iPoint = 0; iPoint < poMP->getNumGeometries(); iPoint++ )
        {
            OGRPoint    *poPoint = (OGRPoint *) poMP->getGeometryRef(iPoint);
            
            padfX[iPoint] = poPoint->getX();
            padfY[iPoint] = poPoint->getY();
            padfZ[iPoint] = poPoint->getZ();
        }

        psShape = SHPCreateSimpleObject( hSHP->nShapeType,
                                         poMP->getNumGeometries(),
                                         padfX, padfY, padfZ );
        SHPWriteObject( hSHP, iShape, psShape );
        SHPDestroyObject( psShape );
        
        CPLFree( padfX );
        CPLFree( padfY );
        CPLFree( padfZ );
    }

/* ==================================================================== */
/*      Arcs from simple line strings.                                  */
/* ==================================================================== */
    else if( (hSHP->nShapeType == SHPT_ARC
              || hSHP->nShapeType == SHPT_ARCM
              || hSHP->nShapeType == SHPT_ARCZ)
             && wkbFlatten(poGeom->getGeometryType()) == wkbLineString )
    {
        OGRLineString   *poArc = (OGRLineString *) poGeom;
        double          *padfX, *padfY, *padfZ;
        int             iPoint;
        SHPObject       *psShape;

        if( poGeom->getGeometryType() != wkbLineString
            && poGeom->getGeometryType() != wkbLineString25D )
        {
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Attempt to write non-linestring (%s) geometry to "
                      "ARC type shapefile.",
                      poGeom->getGeometryName() );

            return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
        }

        padfX = (double *) CPLMalloc(sizeof(double)*poArc->getNumPoints());
        padfY = (double *) CPLMalloc(sizeof(double)*poArc->getNumPoints());
        padfZ = (double *) CPLCalloc(sizeof(double),poArc->getNumPoints());

        for( iPoint = 0; iPoint < poArc->getNumPoints(); iPoint++ )
        {
            padfX[iPoint] = poArc->getX( iPoint );
            padfY[iPoint] = poArc->getY( iPoint );
            padfZ[iPoint] = poArc->getZ( iPoint );
        }

        psShape = SHPCreateSimpleObject( hSHP->nShapeType,
                                         poArc->getNumPoints(),
                                         padfX, padfY, padfZ );
        SHPWriteObject( hSHP, iShape, psShape );
        SHPDestroyObject( psShape );
        
        CPLFree( padfX );
        CPLFree( padfY );
        CPLFree( padfZ );
    }
/* ==================================================================== */
/*      Arcs - Try to treat as MultiLineString.                         */
/* ==================================================================== */
    else if( hSHP->nShapeType == SHPT_ARC
             || hSHP->nShapeType == SHPT_ARCM
             || hSHP->nShapeType == SHPT_ARCZ )
    {
        OGRMultiLineString *poML;
        double          *padfX=NULL, *padfY=NULL, *padfZ=NULL;
        int             iGeom, iPoint, nPointCount = 0;
        SHPObject       *psShape;
        int             *panRingStart;

        poML = (OGRMultiLineString *) 
            OGRGeometryFactory::forceToMultiLineString( poGeom->clone() );

        if( wkbFlatten(poML->getGeometryType()) != wkbMultiLineString )
        {
            delete poML;
            CPLError( CE_Failure, CPLE_AppDefined,
                      "Attempt to write non-linestring (%s) geometry to "
                      "ARC type shapefile.",
                      poGeom->getGeometryName() );

            return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
        }

        panRingStart = (int *) 
            CPLMalloc(sizeof(int) * poML->getNumGeometries());

        for( iGeom = 0; iGeom < poML->getNumGeometries(); iGeom++ )
        {
            OGRLineString *poArc = (OGRLineString *)
                poML->getGeometryRef(iGeom);
            int nNewPoints = poArc->getNumPoints();

            panRingStart[iGeom] = nPointCount;

            padfX = (double *) 
                CPLRealloc( padfX, sizeof(double)*(nNewPoints+nPointCount) );
            padfY = (double *) 
                CPLRealloc( padfY, sizeof(double)*(nNewPoints+nPointCount) );
            padfZ = (double *) 
                CPLRealloc( padfZ, sizeof(double)*(nNewPoints+nPointCount) );

            for( iPoint = 0; iPoint < nNewPoints; iPoint++ )
            {
                padfX[nPointCount] = poArc->getX( iPoint );
                padfY[nPointCount] = poArc->getY( iPoint );
                padfZ[nPointCount] = poArc->getZ( iPoint );
                nPointCount++;
            }
        }
        
        psShape = SHPCreateObject( hSHP->nShapeType, iShape, 
                                   poML->getNumGeometries(), 
                                   panRingStart, NULL,