mitab_spatialref.cpp

GIS系统支持库Geospatial Data Abstraction Library代码.GDAL is a translator library for raster geospatial dat

         *-------------------------------------------------------------*/
      case 16:
        m_poSpatialRef->SetSinusoidal( sTABProj.adProjParams[0],
                                       0.0, 0.0 );
        break;

        /*--------------------------------------------------------------
         * Gall Stereographic
         *-------------------------------------------------------------*/
      case 17:
        m_poSpatialRef->SetGS( sTABProj.adProjParams[0], 0.0, 0.0 );
        break;
        
        /*--------------------------------------------------------------
         * New Zealand Map Grid
         *-------------------------------------------------------------*/
      case 18:
        m_poSpatialRef->SetNZMG( sTABProj.adProjParams[1],
                                 sTABProj.adProjParams[0],
                                 sTABProj.adProjParams[2] * dfConv,
                                 sTABProj.adProjParams[3] * dfConv );
        break;

        /*--------------------------------------------------------------
         * Lambert Conic Conformal (Belgium)
         *-------------------------------------------------------------*/
      case 19:
        m_poSpatialRef->SetLCCB( sTABProj.adProjParams[2],
                                 sTABProj.adProjParams[3],
                                 sTABProj.adProjParams[1],
                                 sTABProj.adProjParams[0],
                                 sTABProj.adProjParams[4] * dfConv,
                                 sTABProj.adProjParams[5] * dfConv );
        break;

        /*--------------------------------------------------------------
         * Stereographic
         *-------------------------------------------------------------*/
      case 20:
      case 31: /* this is called Double Stereographic, whats the diff? */
        m_poSpatialRef->SetStereographic( sTABProj.adProjParams[1],
                                          sTABProj.adProjParams[0],
                                          sTABProj.adProjParams[2],
                                          sTABProj.adProjParams[3] * dfConv,
                                          sTABProj.adProjParams[4] * dfConv );
        break;

        /*--------------------------------------------------------------
         * Swiss Oblique Mercator / Cylindrical
         *-------------------------------------------------------------*/
      case 25:
        m_poSpatialRef->SetSOC( sTABProj.adProjParams[1],
                                sTABProj.adProjParams[0],
                                sTABProj.adProjParams[2] * dfConv,
                                sTABProj.adProjParams[3] * dfConv );
        break;

        /*--------------------------------------------------------------
         * Regional Mercator (regular mercator with a latitude).
         *-------------------------------------------------------------*/
      case 26:
        m_poSpatialRef->SetMercator( sTABProj.adProjParams[1],
                                     sTABProj.adProjParams[0],
                                     1.0, 0.0, 0.0 );
        break;

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

        /*--------------------------------------------------------------
         * Cassini/Soldner
         *-------------------------------------------------------------*/
      case 30:
        m_poSpatialRef->SetCS( sTABProj.adProjParams[1],
                               sTABProj.adProjParams[0],
                               sTABProj.adProjParams[2] * dfConv,
                               sTABProj.adProjParams[3] * dfConv );
        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 ) );
    }

    /*-----------------------------------------------------------------
     * 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)
            && 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.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;
    
    /*-----------------------------------------------------------------
     * Get the linear units and conversion.
     *----------------------------------------------------------------*/
    char        *pszLinearUnits;
    double      dfLinearConv;

    dfLinearConv = poSpatialRef->GetLinearUnits( &pszLinearUnits );
    if( dfLinearConv == 0.0 )
        dfLinearConv = 1.0;

    /*-----------------------------------------------------------------
     * Transform the projection and projection parameters.
     *----------------------------------------------------------------*/
    const char *pszProjection = poSpatialRef->GetAttrValue("PROJECTION");
    double      *parms = sTABProj.adProjParams;

    if( pszProjection == NULL && poSpatialRef->GetAttrNode("LOCAL_CS") != NULL)
    {
        /* nonearth */
        sTABProj.nProjId = 0;
    }

    else if( pszProjection == NULL )
    {
        sTABProj.nProjId = 1;
    }

    else if( EQUAL(pszProjection,SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
    {
        sTABProj.nProjId = 9;
        parms[0] = poSpatialRef->GetProjParm(SRS_PP_LONGITUDE_OF_CENTER,0.0);
        parms[1] = poSpatialRef->GetProjParm(SRS_PP_LATITUDE_OF_CENTER,0.0);
        parms[2] = poSpatialRef->GetProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
        parms[3] = poSpatialRef->GetProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
        parms[4] = poSpatialRef->GetProjParm(SRS_PP_FALSE_EASTING,0.0) / dfLinearConv;
        parms[5] = poSpatialRef->GetProjParm(SRS_PP_FALSE_NORTHING,0.0) / dfLinearConv;
    }

    else if( EQUAL(pszProjection,SRS_PT_AZIMUTHAL_EQUIDISTANT) )
    {
        sTABProj.nProjId = 5;
        parms[0] = poSpatialRef->GetProjParm(SRS_PP_LONGITUDE_OF_CENTER,0.0);
        parms[1] = poSpatialRef->GetProjParm(SRS_PP_LATITUDE_OF_CENTER,0.0);
        parms[2] = 90.0;