mitab_spatialref.cpp

mitab,读取MapInfo的地图文件

        /*--------------------------------------------------------------
         * Polyconic
         *-------------------------------------------------------------*/
      case 27:
        m_poSpatialRef->SetPolyconic( sTABProj.adProjParams[1],
                                      sTABProj.adProjParams[0],
                                      sTABProj.adProjParams[2],
                                      sTABProj.adProjParams[3] );
        break;

        /*--------------------------------------------------------------
         * Cassini/Soldner
         *-------------------------------------------------------------*/
      case 30:
        m_poSpatialRef->SetCS( sTABProj.adProjParams[1],
                               sTABProj.adProjParams[0],
                               sTABProj.adProjParams[2],
                               sTABProj.adProjParams[3] );
        break;

      default:
        break;
    }

    /*-----------------------------------------------------------------
     * Collect units definition.
     *----------------------------------------------------------------*/
    if( sTABProj.nProjId != 1 && m_poSpatialRef->GetRoot() != NULL )
    {
        OGR_SRSNode     *poUnits = new OGR_SRSNode("UNIT");
        
        m_poSpatialRef->GetRoot()->AddChild(poUnits);

        poUnits->AddChild( new OGR_SRSNode( pszUnitsName ) );
        poUnits->AddChild( new OGR_SRSNode( pszUnitsConv ) );
    }

#ifdef MITAB_AFFINE_PARAMS  // See MITAB bug 1155
    /*-----------------------------------------------------------------
     * Collect affine definitions. (Added by Encom 2003)
     *----------------------------------------------------------------*/
    if( sTABProj.nAffineFlag==1 && m_poSpatialRef->GetRoot() != NULL )
    {
        m_poSpatialRef->nAffineFlag = 1;
        m_poSpatialRef->nAffineUnit = sTABProj.nAffineUnits;
        m_poSpatialRef->dAffineParamA = sTABProj.dAffineParamA;
        m_poSpatialRef->dAffineParamB = sTABProj.dAffineParamB;
        m_poSpatialRef->dAffineParamC = sTABProj.dAffineParamC;
        m_poSpatialRef->dAffineParamD = sTABProj.dAffineParamD;
        m_poSpatialRef->dAffineParamE = sTABProj.dAffineParamE;
        m_poSpatialRef->dAffineParamF = sTABProj.dAffineParamF;
    }
    else
    {
        m_poSpatialRef->nAffineFlag = 0;
    }
#endif // MITAB_AFFINE_PARAMS

    /*-----------------------------------------------------------------
     * Local (nonearth) coordinate systems have no Geographic relationship
     * so we just return from here. 
     *----------------------------------------------------------------*/
    if( sTABProj.nProjId == 0 )
        return m_poSpatialRef;

    /*-----------------------------------------------------------------
     * Set the datum.  We are only given the X, Y and Z shift for
     * the datum, so for now we just synthesize a name from this.
     * It would be better if we could lookup a name based on the shift.
     *
     * Since we have already encountered files in which adDatumParams[] values
     * were in the order of 1e-150 when they should have actually been zeros,
     * we will use an epsilon in our scan instead of looking for equality.
     *----------------------------------------------------------------*/
#define TAB_EQUAL(a, b) (((a)<(b) ? ((b)-(a)) : ((a)-(b))) < 1e-10)
    char        szDatumName[160];
    int         iDatumInfo;
    MapInfoDatumInfo *psDatumInfo = NULL;

    for( iDatumInfo = 0;
         asDatumInfoList[iDatumInfo].nMapInfoDatumID != -1;
         iDatumInfo++ )
    {
        psDatumInfo = asDatumInfoList + iDatumInfo;
        
        if( TAB_EQUAL(psDatumInfo->nEllipsoid, sTABProj.nEllipsoidId) &&
            ((sTABProj.nDatumId > 0 && 
              sTABProj.nDatumId == psDatumInfo->nMapInfoDatumID) ||
             (sTABProj.nDatumId <= 0
              && TAB_EQUAL(psDatumInfo->dfShiftX, sTABProj.dDatumShiftX)
              && TAB_EQUAL(psDatumInfo->dfShiftY, sTABProj.dDatumShiftY)
              && TAB_EQUAL(psDatumInfo->dfShiftZ, sTABProj.dDatumShiftZ)
              && TAB_EQUAL(psDatumInfo->dfDatumParm0,sTABProj.adDatumParams[0])
              && TAB_EQUAL(psDatumInfo->dfDatumParm1,sTABProj.adDatumParams[1])
              && TAB_EQUAL(psDatumInfo->dfDatumParm2,sTABProj.adDatumParams[2])
              && TAB_EQUAL(psDatumInfo->dfDatumParm3,sTABProj.adDatumParams[3])
              && TAB_EQUAL(psDatumInfo->dfDatumParm4,sTABProj.adDatumParams[4]))))
            break;

        psDatumInfo = NULL;
    }

    if( psDatumInfo == NULL )
    {
        if( sTABProj.adDatumParams[0] == 0.0
            && sTABProj.adDatumParams[1] == 0.0
            && sTABProj.adDatumParams[2] == 0.0
            && sTABProj.adDatumParams[3] == 0.0
            && sTABProj.adDatumParams[4] == 0.0 )
        {
            sprintf( szDatumName,
                     "MIF 999,%d,%.4g,%.4g,%.4g",
                     sTABProj.nEllipsoidId,
                     sTABProj.dDatumShiftX, 
                     sTABProj.dDatumShiftY, 
                     sTABProj.dDatumShiftZ );
        }
        else
        {
            sprintf( szDatumName,
                     "MIF 9999,%d,%.4g,%.4g,%.4g,%.15g,%.15g,%.15g,%.15g,%.15g",
                     sTABProj.nEllipsoidId,
                     sTABProj.dDatumShiftX, 
                     sTABProj.dDatumShiftY, 
                     sTABProj.dDatumShiftZ,
                     sTABProj.adDatumParams[0],
                     sTABProj.adDatumParams[1],
                     sTABProj.adDatumParams[2],
                     sTABProj.adDatumParams[3],
                     sTABProj.adDatumParams[4] );
        }
    }
    else if( strlen(psDatumInfo->pszOGCDatumName) > 0 )
    {
        strncpy( szDatumName, psDatumInfo->pszOGCDatumName,
                 sizeof(szDatumName) );
    }
    else
    {
        sprintf( szDatumName, "MIF %d", psDatumInfo->nMapInfoDatumID );
    }

    /*-----------------------------------------------------------------
     * Set the spheroid.
     *----------------------------------------------------------------*/
    double      dfSemiMajor=0.0, dfInvFlattening=0.0;
    const char *pszSpheroidName = NULL;

    for( int i = 0; asSpheroidInfoList[i].nMapInfoId != -1; i++ )
    {
        if( asSpheroidInfoList[i].nMapInfoId == sTABProj.nEllipsoidId )
        {
            dfSemiMajor = asSpheroidInfoList[i].dfA;
            dfInvFlattening = asSpheroidInfoList[i].dfInvFlattening;
            pszSpheroidName = asSpheroidInfoList[i].pszMapinfoName;
            break;
        }
    }

    // use WGS 84 if nothing is known.
    if( pszSpheroidName == NULL )
    {
        pszSpheroidName = "unknown";
        dfSemiMajor = 6378137.0;
        dfInvFlattening = 298.257223563;
    }

    /*-----------------------------------------------------------------
     * Set the prime meridian.
     *----------------------------------------------------------------*/
    double      dfPMOffset = 0.0;
    const char *pszPMName = "Greenwich";
    
    if( sTABProj.adDatumParams[4] != 0.0 )
    {
        dfPMOffset = sTABProj.adDatumParams[4];

        pszPMName = "non-Greenwich";
    }
                    
    /*-----------------------------------------------------------------
     * Create a GEOGCS definition.
     *----------------------------------------------------------------*/

    m_poSpatialRef->SetGeogCS( "unnamed",
                               szDatumName,
                               pszSpheroidName,
                               dfSemiMajor, dfInvFlattening,
                               pszPMName, dfPMOffset,
                               SRS_UA_DEGREE, atof(SRS_UA_DEGREE_CONV));

    if( psDatumInfo != NULL )
    {
        m_poSpatialRef->SetTOWGS84( psDatumInfo->dfShiftX, 
                                    psDatumInfo->dfShiftY,
                                    psDatumInfo->dfShiftZ,
                                    -psDatumInfo->dfDatumParm0, 
                                    -psDatumInfo->dfDatumParm1, 
                                    -psDatumInfo->dfDatumParm2, 
                                    psDatumInfo->dfDatumParm3 );
    }

    return m_poSpatialRef;
}

/**********************************************************************
 *                   TABFile::SetSpatialRef()
 *
 * Set the OGRSpatialReference for this dataset.
 * A reference to the OGRSpatialReference will be kept, and it will also
 * be converted into a TABProjInfo to be stored in the .MAP header.
 *
 * Returns 0 on success, and -1 on error.
 **********************************************************************/
int TABFile::SetSpatialRef(OGRSpatialReference *poSpatialRef)
{
    if (m_eAccessMode != TABWrite)
    {
        CPLError(CE_Failure, CPLE_NotSupported,
                 "SetSpatialRef() can be used only with Write access.");
        return -1;
    }

    if (m_poMAPFile == NULL )
    {
        CPLError(CE_Failure, CPLE_AssertionFailed,
                 "SetSpatialRef() failed: file has not been opened yet.");
        return -1;
    }

    if( poSpatialRef == NULL )
    {
        CPLError(CE_Failure, CPLE_AssertionFailed,
                 "SetSpatialRef() failed: Called with NULL poSpatialRef.");
        return -1;
    }

    /*-----------------------------------------------------------------
     * Keep a copy of the OGRSpatialReference...
     * Note: we have to take the reference count into account...
     *----------------------------------------------------------------*/
    if (m_poSpatialRef && m_poSpatialRef->Dereference() == 0)
        delete m_poSpatialRef;
    
    m_poSpatialRef = poSpatialRef->Clone();

    /*-----------------------------------------------------------------
     * Initialize TABProjInfo
     *----------------------------------------------------------------*/
    TABProjInfo     sTABProj;

    sTABProj.nProjId = 0;
    sTABProj.nEllipsoidId = 0; /* how will we set this? */
    sTABProj.nUnitsId = 7;
    sTABProj.adProjParams[0] = sTABProj.adProjParams[1] = 0.0;
    sTABProj.adProjParams[2] = sTABProj.adProjParams[3] = 0.0;
    sTABProj.adProjParams[4] = sTABProj.adProjParams[5] = 0.0;
    
    sTABProj.nDatumId = 0;
    sTABProj.dDatumShiftX = 0.0;
    sTABProj.dDatumShiftY = 0.0;
    sTABProj.dDatumShiftZ = 0.0;
    sTABProj.adDatumParams[0] = 0.0;
    sTABProj.adDatumParams[1] = 0.0;
    sTABProj.adDatumParams[2] = 0.0;
    sTABProj.adDatumParams[3] = 0.0;
    sTABProj.adDatumParams[4] = 0.0;

    // Encom 2003
    sTABProj.nAffineFlag   = 0;
    sTABProj.nAffineUnits  = 7;
    sTABProj.dAffineParamA = 0.0;
    sTABProj.dAffineParamB = 0.0;
    sTABProj.dAffineParamC = 0.0;
    sTABProj.dAffineParamD = 0.0;
    sTABProj.dAffineParamE = 0.0;
    sTABProj.dAffineParamF = 0.0;
    
    /*-----------------------------------------------------------------
     * Get the linear units and conversion.
     *----------------------------------------------------------------*/