s57reader.cpp

支持各种栅格图像和矢量图像读取的库

        if( dfEndAngle - dfStartAngle > 360.0 )
        {
            double dfTempAngle;

            dfTempAngle = dfStartAngle;
            dfStartAngle = dfEndAngle;
            dfEndAngle = dfTempAngle;

            while( dfEndAngle < dfStartAngle )
                dfStartAngle -= 360.0;
        }
    }

    dfRadius = sqrt( (dfCenterX - dfStartX) * (dfCenterX - dfStartX)
                     + (dfCenterY - dfStartY) * (dfCenterY - dfStartY) );
    
    return S57StrokeArcToOGRGeometry_Angles( dfCenterX, dfCenterY, 
                                             dfRadius, 
                                             dfStartAngle, dfEndAngle,
                                             nVertexCount );
}

/************************************************************************/
/*                             FetchLine()                              */
/************************************************************************/

int S57Reader::FetchLine( DDFRecord *poSRecord, 
                          int iStartVertex, int iDirection,
                          OGRLineString *poLine )

{
    int             nVCount;
    DDFField        *poSG2D = poSRecord->FindField( "SG2D" );
    DDFField        *poAR2D = poSRecord->FindField( "AR2D" );
    DDFSubfieldDefn *poXCOO=NULL, *poYCOO=NULL;
    int bStandardFormat = TRUE;

    if( poSG2D == NULL && poAR2D != NULL )
        poSG2D = poAR2D;

/* -------------------------------------------------------------------- */
/*      Get some basic definitions.                                     */
/* -------------------------------------------------------------------- */
    if( poSG2D != NULL )
    {
        poXCOO = poSG2D->GetFieldDefn()->FindSubfieldDefn("XCOO");
        poYCOO = poSG2D->GetFieldDefn()->FindSubfieldDefn("YCOO");

        if( poXCOO == NULL || poYCOO == NULL )
        {
            CPLDebug( "S57", "XCOO or YCOO are NULL" );
            return FALSE;
        }

        nVCount = poSG2D->GetRepeatCount();
    }
    else
    {
        return TRUE;
    }

/* -------------------------------------------------------------------- */
/*      It is legitimate to have zero vertices for line segments        */
/*      that just have the start and end node (bug 840).                */
/* -------------------------------------------------------------------- */
    if( nVCount == 0 )
        return TRUE;
 
/* -------------------------------------------------------------------- */
/*      Make sure out line is long enough to hold all the vertices      */
/*      we will apply.                                                  */
/* -------------------------------------------------------------------- */
    int nVBase;

    if( iDirection < 0 )
        nVBase = iStartVertex + nVCount;
    else
        nVBase = iStartVertex;

    if( poLine->getNumPoints() < iStartVertex + nVCount )
        poLine->setNumPoints( iStartVertex + nVCount );
        
/* -------------------------------------------------------------------- */
/*      Are the SG2D and XCOO/YCOO definitions in the form we expect?   */
/* -------------------------------------------------------------------- */
    if( poSG2D->GetFieldDefn()->GetSubfieldCount() != 2 )
        bStandardFormat = FALSE;

    if( !EQUAL(poXCOO->GetFormat(),"b24") 
        || !EQUAL(poYCOO->GetFormat(),"b24") )
        bStandardFormat = FALSE;

/* -------------------------------------------------------------------- */
/*      Collect the vertices:                                           */
/*                                                                      */
/*      This approach assumes that the data is LSB organized int32      */
/*      binary data as per the specification.  We avoid lots of         */
/*      extra calls to low level DDF methods as they are quite          */
/*      expensive.                                                      */
/* -------------------------------------------------------------------- */
    if( bStandardFormat )
    {
        const char  *pachData;
        int         nBytesRemaining;

        pachData = poSG2D->GetSubfieldData(poYCOO,&nBytesRemaining,0);
        
        for( int i = 0; i < nVCount; i++ )
        {
            double      dfX, dfY;
            GInt32      nXCOO, nYCOO;

            memcpy( &nYCOO, pachData, 4 );
            pachData += 4;
            memcpy( &nXCOO, pachData, 4 );
            pachData += 4;

#ifdef CPL_MSB
            CPL_SWAP32PTR( &nXCOO );
            CPL_SWAP32PTR( &nYCOO );
#endif
            dfX = nXCOO / (double) nCOMF;
            dfY = nYCOO / (double) nCOMF;

            poLine->setPoint( nVBase, dfX, dfY );

            nVBase += iDirection;
        }
    }

/* -------------------------------------------------------------------- */
/*      Collect the vertices:                                           */
/*                                                                      */
/*      The generic case where we use low level but expensive DDF       */
/*      methods to get the data.  This should work even if some         */
/*      things are changed about the SG2D fields such as making them    */
/*      floating point or a different byte order.                       */
/* -------------------------------------------------------------------- */
    else
    {
        for( int i = 0; i < nVCount; i++ )
        {
            double      dfX, dfY;
            const char  *pachData;
            int         nBytesRemaining;

            pachData = poSG2D->GetSubfieldData(poXCOO,&nBytesRemaining,i);
                
            dfX = poXCOO->ExtractIntData(pachData,nBytesRemaining,NULL)
                / (double) nCOMF;

            pachData = poSG2D->GetSubfieldData(poYCOO,&nBytesRemaining,i);

            dfY = poXCOO->ExtractIntData(pachData,nBytesRemaining,NULL)
                / (double) nCOMF;
                
            poLine->setPoint( nVBase, dfX, dfY );

            nVBase += iDirection;
        }
    }

/* -------------------------------------------------------------------- */
/*      If this is actually an arc, turn the start, end and center      */
/*      of rotation into a "stroked" arc linestring.                    */
/* -------------------------------------------------------------------- */
    if( poAR2D != NULL && poLine->getNumPoints() >= 3 )
    {
        OGRLineString *poArc;
        int i, iLast = poLine->getNumPoints() - 1;
        
        poArc = S57StrokeArcToOGRGeometry_Points( 
            poLine->getX(iLast-0), poLine->getY(iLast-0), 
            poLine->getX(iLast-1), poLine->getY(iLast-1),
            poLine->getX(iLast-2), poLine->getY(iLast-2),
            30 );

        if( poArc != NULL )
        {
            for( i = 0; i < poArc->getNumPoints(); i++ )
                poLine->setPoint( iLast-2+i, poArc->getX(i), poArc->getY(i) );

            delete poArc;
        }
    }

    return TRUE;
}

/************************************************************************/
/*                       AssemblePointGeometry()                        */
/************************************************************************/

void S57Reader::AssemblePointGeometry( DDFRecord * poFRecord,
                                       OGRFeature * poFeature )

{
    DDFField    *poFSPT;
    int         nRCNM, nRCID;

/* -------------------------------------------------------------------- */
/*      Feature the spatial record containing the point.                */
/* -------------------------------------------------------------------- */
    poFSPT = poFRecord->FindField( "FSPT" );
    if( poFSPT == NULL )
        return;

    if( poFSPT->GetRepeatCount() != 1 )
    {
#ifdef DEBUG
        fprintf( stderr, 
                 "Point features with other than one spatial linkage.\n" );
        poFRecord->Dump( stderr );
#endif
        CPLDebug( "S57", 
           "Point feature encountered with other than one spatial linkage." );
    }
        
    nRCID = ParseName( poFSPT, 0, &nRCNM );

    double      dfX = 0.0, dfY = 0.0, dfZ = 0.0;

    if( !FetchPoint( nRCNM, nRCID, &dfX, &dfY, &dfZ ) )
    {
        CPLAssert( FALSE );
        return;
    }

    if( dfZ == 0.0 )
        poFeature->SetGeometryDirectly( new OGRPoint( dfX, dfY ) );
    else
        poFeature->SetGeometryDirectly( new OGRPoint( dfX, dfY, dfZ ) );
}

/************************************************************************/
/*                      AssembleSoundingGeometry()                      */
/************************************************************************/

void S57Reader::AssembleSoundingGeometry( DDFRecord * poFRecord,
                                          OGRFeature * poFeature )

{
    DDFField    *poFSPT;
    int         nRCNM, nRCID;
    DDFRecord   *poSRecord;
    

/* -------------------------------------------------------------------- */
/*      Feature the spatial record containing the point.                */
/* -------------------------------------------------------------------- */
    poFSPT = poFRecord->FindField( "FSPT" );
    if( poFSPT == NULL )
        return;

    CPLAssert( poFSPT->GetRepeatCount() == 1 );
        
    nRCID = ParseName( poFSPT, 0, &nRCNM );

    if( nRCNM == RCNM_VI )
        poSRecord = oVI_Index.FindRecord( nRCID );
    else
        poSRecord = oVC_Index.FindRecord( nRCID );

    if( poSRecord == NULL )
        return;

/* -------------------------------------------------------------------- */
/*      Extract vertices.                                               */
/* -------------------------------------------------------------------- */
    OGRMultiPoint       *poMP = new OGRMultiPoint();
    DDFField            *poField;
    int                 nPointCount, i, nBytesLeft;
    DDFSubfieldDefn    *poXCOO, *poYCOO, *poVE3D;
    const char         *pachData;

    poField = poSRecord->FindField( "SG2D" );
    if( poField == NULL )
        poField = poSRecord->FindField( "SG3D" );
    if( poField == NULL )
        return;

    poXCOO = poField->GetFieldDefn()->FindSubfieldDefn( "XCOO" );
    poYCOO = poField->GetFieldDefn()->FindSubfieldDefn( "YCOO" );
    poVE3D = poField->GetFieldDefn()->FindSubfieldDefn( "VE3D" );

    nPointCount = poField->GetRepeatCount();

    pachData = poField->GetData();
    nBytesLeft = poField->GetDataSize();

    for( i = 0; i < nPointCount; i++ )
    {
        double          dfX, dfY, dfZ = 0.0;
        int             nBytesConsumed;

        dfY = poYCOO->ExtractIntData( pachData, nBytesLeft,
                                      &nBytesConsumed ) / (double) nCOMF;
        nBytesLeft -= nBytesConsumed;
        pachData += nBytesConsumed;

        dfX = poXCOO->ExtractIntData( pachData, nBytesLeft,
                                      &nBytesConsumed ) / (double) nCOMF;
        nBytesLeft -= nBytesConsumed;
        pachData += nBytesConsumed;
        
        if( poVE3D != NULL )
        {
            dfZ = poYCOO->ExtractIntData( pachData, nBytesLeft,
                                          &nBytesConsumed ) / (double) nSOMF;
            nBytesLeft -= nBytesConsumed;
            pachData += nBytesConsumed;
        }

        poMP->addGeometryDirectly( new OGRPoint( dfX, dfY, dfZ ) );
    }

    poFeature->SetGeometryDirectly( poMP );
}

/************************************************************************/
/*                        AssembleLineGeometry()                        */
/************************************************************************/

void S57Reader::AssembleLineGeometry( DDFRecord * poFRecord,
                                      OGRFeature * poFeature )

{
    DDFField    *poFSPT;
    int         nEdgeCount;
    OGRLineString *poLine = new OGRLineString();

/* -------------------------------------------------------------------- */
/*      Find the FSPT field.                                            */
/* -------------------------------------------------------------------- */
    poFSPT = poFRecord->FindField( "FSPT" );
    if( poFSPT == NULL )
        return;

    nEdgeCount = poFSPT->GetRepeatCount();

/* ==================================================================== */
/*      Loop collecting edges.                                          */
/* ==================================================================== */
    for( int iEdge = 0; iEdge < nEdgeCount; iEdge++ )
    {
        DDFRecord       *poSRecord;
        int             nRCID;

/* -------------------------------------------------------------------- */
/*      Find the spatial record for this edge.                          */
/* -------------------------------------------------------------------- */
        nRCID = ParseName( poFSPT, iEdge );

        poSRecord = oVE_Index.FindRecord( nRCID );
        if( poSRecord == NULL )
        {
            CPLError( CE_Warning, CPLE_AppDefined,
                      "Couldn't find spatial record %d.\n"
                      "Feature OBJL=%s, RCID=%d may have corrupt or"
                      "missing geometry.",
                      nRCID,
                      poFeature->GetDefnRef()->GetName(),