ntffilereader.cpp

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

            }
            else
            {
                psAttDesc->poCodeList = poCodeList;
            }
        }

/* -------------------------------------------------------------------- */
/*      Handle database header record.                                  */
/* -------------------------------------------------------------------- */
        else if( poRecord->GetType() == NRT_DHR )
        {
            int         iChar;
            pszProduct = CPLStrdup(poRecord->GetField(3,22));
            for( iChar = strlen(pszProduct)-1;
                 iChar > 0 && pszProduct[iChar] == ' ';
                 pszProduct[iChar--] = '\0' ) {}

            pszPVName = CPLStrdup(poRecord->GetField(76+3,76+22));
            for( iChar = strlen(pszPVName)-1;
                 iChar > 0 && pszPVName[iChar] == ' ';
                 pszPVName[iChar--] = '\0' ) {}

        }

        delete poRecord;
    }

/* -------------------------------------------------------------------- */
/*      Did we fall off the end without finding what we were looking    */
/*      for?                                                            */
/* -------------------------------------------------------------------- */
    if( poRecord->GetType() == NRT_VTR )
    {
        delete poRecord;
        CPLError( CE_Failure, CPLE_AppDefined,
                  "Cound not find section header record in %s.\n", 
                  pszFilename );
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Classify the product type.                                      */
/* -------------------------------------------------------------------- */
    if( EQUALN(pszProduct,"LAND-LINE",9) && atof(pszPVName+5) < 1.3 )
        nProduct = NPC_LANDLINE;
    else if( EQUALN(pszProduct,"LAND-LINE",9) )
        nProduct = NPC_LANDLINE99;
    else if( EQUAL(pszProduct,"OS_LANDRANGER_CONT") ) // Panorama
        nProduct = NPC_LANDRANGER_CONT;
    else if( EQUAL(pszProduct,"L-F_PROFILE_CON") ) // Panorama
        nProduct = NPC_LANDFORM_PROFILE_CONT;
    else if( EQUALN(pszProduct,"Strategi",8) )
        nProduct = NPC_STRATEGI;
    else if( EQUALN(pszProduct,"Meridian_02",11) )
        nProduct = NPC_MERIDIAN2;
    else if( EQUALN(pszProduct,"Meridian_01",11) )
        nProduct = NPC_MERIDIAN;
    else if( EQUAL(pszProduct,NTF_BOUNDARYLINE) 
             && EQUALN(pszPVName,"A10N_FC",7) )
        nProduct = NPC_BOUNDARYLINE;
    else if( EQUAL(pszProduct,NTF_BOUNDARYLINE) 
             && EQUALN(pszPVName,"A20N_FC",7) )
        nProduct = NPC_BL2000;
    else if( EQUALN(pszProduct,"BaseData.GB",11) )
        nProduct = NPC_BASEDATA;
    else if( EQUALN(pszProduct,"OSCAR_ASSET",11) )
        nProduct = NPC_OSCAR_ASSET;
    else if( EQUALN(pszProduct,"OSCAR_TRAFF",11) )
        nProduct = NPC_OSCAR_TRAFFIC;
    else if( EQUALN(pszProduct,"OSCAR_ROUTE",11) )
        nProduct = NPC_OSCAR_ROUTE;
    else if( EQUALN(pszProduct,"OSCAR_NETWO",11) )
        nProduct = NPC_OSCAR_NETWORK;
    else if( EQUALN(pszProduct,"ADDRESS_POI",11) )
        nProduct = NPC_ADDRESS_POINT;
    else if( EQUALN(pszProduct,"CODE_POINT",10) )
    {
        if( GetAttDesc( "RH" ) == NULL )
            nProduct = NPC_CODE_POINT;
        else
            nProduct = NPC_CODE_POINT_PLUS;
    }
    else if( EQUALN(pszProduct,"OS_LANDRANGER_DTM",17) )
        nProduct = NPC_LANDRANGER_DTM;
    else if( EQUALN(pszProduct,"L-F_PROFILE_DTM",15) )
        nProduct = NPC_LANDFORM_PROFILE_DTM;
    else if( EQUALN(pszProduct,"NEXTMap Britian DTM",19) )
        nProduct = NPC_LANDFORM_PROFILE_DTM; // Treat as landform

    if( poDS->GetOption("FORCE_GENERIC") != NULL
        && !EQUAL(poDS->GetOption("FORCE_GENERIC"),"OFF") )
        nProduct = NPC_UNKNOWN;

    // No point in caching lines if there are no polygons.
    if( nProduct != NPC_BOUNDARYLINE && nProduct != NPC_BL2000 )
        bCacheLines = FALSE;
    
/* -------------------------------------------------------------------- */
/*      Handle the section header record.                               */
/* -------------------------------------------------------------------- */
    nSavedFeatureId = nBaseFeatureId;
    nStartPos = VSIFTell(fp);
    
    pszTileName = CPLStrdup(poRecord->GetField(3,12));        // SECT_REF
    while( pszTileName[strlen(pszTileName)-1] == ' ' )
        pszTileName[strlen(pszTileName)-1] = '\0';

    nCoordWidth = atoi(poRecord->GetField(15,19));            // XYLEN
    if( nCoordWidth == 0 )
        nCoordWidth = 10;
    
    nZWidth = atoi(poRecord->GetField(31,35));                // ZLEN
    if( nZWidth == 0 )
        nZWidth = 10;
    
    dfXYMult = atoi(poRecord->GetField(21,30)) / 1000.0;      // XY_MULT
    dfXOrigin = atoi(poRecord->GetField(47,56));
    dfYOrigin = atoi(poRecord->GetField(57,66));
    dfTileXSize = atoi(poRecord->GetField(23+74,32+74));
    dfTileYSize = atoi(poRecord->GetField(33+74,42+74));
    dfZMult = atoi(poRecord->GetField(37,46)) / 1000.0;

/* -------------------------------------------------------------------- */
/*      Setup scale and transformation factor for text height.          */
/* -------------------------------------------------------------------- */
    if( poRecord->GetLength() >= 187 )
        dfScale = atoi(poRecord->GetField(148+31,148+39));
    else if( nProduct == NPC_STRATEGI )
        dfScale = 250000;
    else if( nProduct == NPC_MERIDIAN || nProduct == NPC_MERIDIAN2 )
        dfScale = 100000;
    else if( nProduct == NPC_LANDFORM_PROFILE_CONT )
        dfScale = 10000;
    else if( nProduct == NPC_LANDRANGER_CONT )
        dfScale = 50000;
    else if( nProduct == NPC_OSCAR_ASSET
             || nProduct == NPC_OSCAR_TRAFFIC
             || nProduct == NPC_OSCAR_NETWORK
             || nProduct == NPC_OSCAR_ROUTE )
        dfScale = 10000;
    else if( nProduct == NPC_BASEDATA )
        dfScale = 625000;
    else if( nProduct == NPC_BOUNDARYLINE )
        dfScale = 10000;
    else
        dfScale = 10000;

    if( dfScale != 0.0 )
        dfPaperToGround = dfScale / 1000.0;
    else
        dfPaperToGround = 0.0;

    delete poRecord;

/* -------------------------------------------------------------------- */
/*      Ensure we have appropriate layers defined.                      */
/* -------------------------------------------------------------------- */
    CPLErrorReset();

    if( !IsRasterProduct() )
        EstablishLayers();
    else
        EstablishRasterAccess();
    
    return CPLGetLastErrorType() != CE_Failure;
}

/************************************************************************/
/*                            DumpReadable()                            */
/************************************************************************/

void NTFFileReader::DumpReadable( FILE *fpLog )

{
    fprintf( fpLog, "Tile Name = %s\n", pszTileName );
    fprintf( fpLog, "Product = %s\n", pszProduct );
    fprintf( fpLog, "NTFLevel = %d\n", nNTFLevel );
    fprintf( fpLog, "XYLEN = %d\n", nCoordWidth );
    fprintf( fpLog, "XY_MULT = %g\n", dfXYMult );
    fprintf( fpLog, "X_ORIG = %g\n", dfXOrigin );
    fprintf( fpLog, "Y_ORIG = %g\n", dfYOrigin ); 
    fprintf( fpLog, "XMAX = %g\n", dfTileXSize );
    fprintf( fpLog, "YMAX = %g\n", dfTileYSize );
}

/************************************************************************/
/*                          ProcessGeometry()                           */
/*                                                                      */
/*      Drop duplicate vertices from line strings ... they mess up      */
/*      FME's polygon handling sometimes.                               */
/************************************************************************/

OGRGeometry *NTFFileReader::ProcessGeometry( NTFRecord * poRecord,
                                             int * pnGeomId )

{
    int            nGType, nNumCoord;
    OGRGeometry    *poGeometry = NULL;

    if( poRecord->GetType() == NRT_GEOMETRY3D )
        return ProcessGeometry3D( poRecord, pnGeomId );

    else if( poRecord->GetType() != NRT_GEOMETRY )
        return NULL;

    nGType = atoi(poRecord->GetField(9,9));            // GTYPE
    nNumCoord = atoi(poRecord->GetField(10,13));       // NUM_COORD
    if( pnGeomId != NULL )
        *pnGeomId = atoi(poRecord->GetField(3,8));     // GEOM_ID

/* -------------------------------------------------------------------- */
/*      Point                                                           */
/* -------------------------------------------------------------------- */
    if( nGType == 1 )
    {
        double      dfX, dfY;
        
        dfX = atoi(poRecord->GetField(14,14+GetXYLen()-1)) * GetXYMult() 
            + GetXOrigin();
        dfY = atoi(poRecord->GetField(14+GetXYLen(),14+GetXYLen()*2-1))
            * GetXYMult() + GetYOrigin();
      
        poGeometry = new OGRPoint( dfX, dfY );
    }
    
/* -------------------------------------------------------------------- */
/*      Line (or arc)                                                   */
/* -------------------------------------------------------------------- */
    else if( nGType == 2 || nGType == 3 || nGType == 4 )
    {
        OGRLineString      *poLine = new OGRLineString;
        double             dfX, dfY, dfXLast=0.0, dfYLast=0.0;
        int                iCoord, nOutCount = 0;

        poGeometry = poLine;
        poLine->setNumPoints( nNumCoord );
        for( iCoord = 0; iCoord < nNumCoord; iCoord++ )
        {
            int            iStart = 14 + iCoord * (GetXYLen()*2+1);

            dfX = atoi(poRecord->GetField(iStart+0,
                                          iStart+GetXYLen()-1)) 
                * GetXYMult() + GetXOrigin();
            dfY = atoi(poRecord->GetField(iStart+GetXYLen(),
                                          iStart+GetXYLen()*2-1)) 
                * GetXYMult() + GetYOrigin();

            if( iCoord == 0 )
            {
                dfXLast = dfX;
                dfYLast = dfY;
                poLine->setPoint( nOutCount++, dfX, dfY );
            }
            else if( dfXLast != dfX || dfYLast != dfY )
            {
                dfXLast = dfX;
                dfYLast = dfY;
                poLine->setPoint( nOutCount++, dfX, dfY );
            }
        }
        poLine->setNumPoints( nOutCount );

        CacheAddByGeomId( atoi(poRecord->GetField(3,8)), poLine );
    }

/* -------------------------------------------------------------------- */
/*      Arc defined by three points on the arc.                         */
/* -------------------------------------------------------------------- */
    else if( nGType == 5 && nNumCoord == 3 )
    {
        double  adfX[3], adfY[3];
        int     iCoord;

        for( iCoord = 0; iCoord < nNumCoord; iCoord++ )
        {
            int            iStart = 14 + iCoord * (GetXYLen()*2+1);

            adfX[iCoord] = atoi(poRecord->GetField(iStart+0,
                                                  iStart+GetXYLen()-1)) 
                * GetXYMult() + GetXOrigin();
            adfY[iCoord] = atoi(poRecord->GetField(iStart+GetXYLen(),
                                                  iStart+GetXYLen()*2-1)) 
                * GetXYMult() + GetYOrigin();
        }

        poGeometry = NTFStrokeArcToOGRGeometry_Points( adfX[0], adfY[0], 
                                                       adfX[1], adfY[1],
                                                       adfX[2], adfY[2], 72 );
    }

/* -------------------------------------------------------------------- */
/*      Circle                                                          */
/* -------------------------------------------------------------------- */
    else if( nGType == 7 )
    {
        double  dfCenterX, dfCenterY, dfArcX, dfArcY, dfRadius;
        int     iCenterStart = 14;
        int     iArcStart = 14 + 2 * GetXYLen() + 1;

        dfCenterX = atoi(poRecord->GetField(iCenterStart,
                                            iCenterStart+GetXYLen()-1))
            * GetXYMult() + GetXOrigin();
        dfCenterY = atoi(poRecord->GetField(iCenterStart+GetXYLen(),
                                            iCenterStart+GetXYLen()*2-1))
            * GetXYMult() + GetYOrigin();
        
        dfArcX = atoi(poRecord->GetField(iArcStart,
                                         iArcStart+GetXYLen()-1))
            * GetXYMult() + GetXOrigin();
        dfArcY = atoi(poRecord->GetField(iArcStart+GetXYLen(),
                                         iArcStart+GetXYLen()*2-1))
            * GetXYMult() + GetYOrigin();

        dfRadius = sqrt( (dfCenterX - dfArcX) * (dfCenterX - dfArcX)
                         + (dfCenterY - dfArcY) * (dfCenterY - dfArcY) );

        poGeometry = NTFStrokeArcToOGRGeometry_Angles( dfCenterX, dfCenterY,
                                                       dfRadius, 
                                                       0.0, 360.0, 
                                                       72 );
    }

    else
    {
        fprintf( stderr, "GType = %d\n", nGType );
        CPLAssert( FALSE );
    }

    if( poGeometry != NULL )
        poGeometry->assignSpatialReference( poDS->GetSpatialRef() );

    return poGeometry;
}

/************************************************************************/
/*                         ProcessGeometry3D()                          */
/************************************************************************/

OGRGeometry *NTFFileReader::ProcessGeometry3D( NTFRecord * poRecord,
                                               int * pnGeomId )

{
    int            nGType, nNumCoord;
    OGRGeometry    *poGeometry = NULL;

    if( poRecord->GetType() != NRT_GEOMETRY3D )
        return NULL;

    nGType = atoi(poRecord->GetField(9,9));            // GTYPE
    nNumCoord = atoi(poRecord->GetField(10,13));       // NUM_COORD
    if( pnGeomId != NULL )
        *pnGeomId = atoi(poRecord->GetField(3,8));     // GEOM_ID

    if( nGType == 1 )
    {
        double      dfX, dfY, dfZ;
        
        dfX = atoi(poRecord->GetField(14,14+GetXYLen()-1)) * GetXYMult() 
            + GetXOrigin();
        dfY = atoi(poRecord->GetField(14+GetXYLen(),14+GetXYLen()*2-1))
            * GetXYMult() + GetYOrigin();
        dfZ = atoi(poRecord->GetField(14+1+2*GetXYLen(),
                                      14+1+2*GetXYLen()+nZWidth-1)) * dfZMult;