ogr_srs_proj4.cpp

mitab,读取MapInfo的地图文件

    }

    else if( EQUAL(pszProj,"krovak") )
    {
        SetKrovak( OSR_GDV( papszNV, "lat_0", 0.0 ), 
                   OSR_GDV( papszNV, "lon_0", 0.0 ), 
                   OSR_GDV( papszNV, "alpha", 0.0 ), 
                   0.0, // pseudo_standard_parallel_1
                   OSR_GDV( papszNV, "k", 1.0 ), 
                   OSR_GDV( papszNV, "x_0", 0.0 ), 
                   OSR_GDV( papszNV, "y_0", 0.0 ) );
    }

    else if( EQUAL(pszProj,"tpeqd") )
    {
        SetTPED( OSR_GDV( papszNV, "lat_1", 0.0 ), 
                 OSR_GDV( papszNV, "lon_1", 0.0 ), 
                 OSR_GDV( papszNV, "lat_2", 0.0 ), 
                 OSR_GDV( papszNV, "lon_2", 0.0 ), 
                 OSR_GDV( papszNV, "x_0", 0.0 ), 
                 OSR_GDV( papszNV, "y_0", 0.0 ) );
    }

    else
    {
        CPLDebug( "OGR_PROJ4", "Unsupported projection: %s", pszProj );
        CSLDestroy( papszNV );
        return OGRERR_CORRUPT_DATA;
    }

/* -------------------------------------------------------------------- */
/*      Try to translate the datum.                                     */
/* -------------------------------------------------------------------- */
    const char *pszValue;
    int  bFullyDefined = FALSE;

    pszValue = CSLFetchNameValue(papszNV, "datum");
    if( pszValue == NULL )
    {
        /* do nothing */
    }
    else if( (EQUAL(pszValue,"NAD27") || EQUAL(pszValue,"NAD83")
              || EQUAL(pszValue,"WGS84") || EQUAL(pszValue,"WGS72"))
             && dfFromGreenwich == 0.0 )
    {
        SetWellKnownGeogCS( pszValue );
        bFullyDefined = TRUE;
    }
    else if( EQUAL(pszValue,"potsdam") )
    {
        OGRSpatialReference oGCS;
        oGCS.importFromEPSG( 4314 );
        CopyGeogCSFrom( &oGCS );
        bFullyDefined = TRUE;
    }
    else if( EQUAL(pszValue,"nzgd49") )
    {
        OGRSpatialReference oGCS;
        oGCS.importFromEPSG( 4272 );
        CopyGeogCSFrom( &oGCS );
        bFullyDefined = TRUE;
    }
    else
    {
        /* we don't recognise the datum, and ignore it */
    }

/* -------------------------------------------------------------------- */
/*      Set the ellipsoid information.                                   */
/* -------------------------------------------------------------------- */
    double dfSemiMajor, dfInvFlattening, dfSemiMinor;

    pszValue = CSLFetchNameValue(papszNV, "ellps");
    if( pszValue != NULL && !bFullyDefined )
    {
        for( i = 0; ogr_pj_ellps[i] != NULL; i += 4 )
        {
            if( !EQUAL(ogr_pj_ellps[i],pszValue) )
                continue;

            CPLAssert( EQUALN(ogr_pj_ellps[i+1],"a=",2) );
            
            dfSemiMajor = atof(ogr_pj_ellps[i+1]+2);
            if( EQUALN(ogr_pj_ellps[i+2],"rf=",3) )
                dfInvFlattening = atof(ogr_pj_ellps[i+2]+3);
            else
            {
                CPLAssert( EQUALN(ogr_pj_ellps[i+2],"b=",2) );
                dfSemiMinor = atof(ogr_pj_ellps[i+2]+2);
                
                if( ABS(dfSemiMajor/dfSemiMinor) - 1.0 < 0.0000000000001 )
                    dfInvFlattening = 0.0;
                else
                    dfInvFlattening = -1.0 / (dfSemiMinor/dfSemiMajor - 1.0);
            }
            
            SetGeogCS( ogr_pj_ellps[i+3], "unknown", ogr_pj_ellps[i], 
                       dfSemiMajor, dfInvFlattening,
                       pszPM, dfFromGreenwich );

            bFullyDefined = TRUE;
            break;
        }
    }

    if( !bFullyDefined )
    {
        dfSemiMajor = OSR_GDV( papszNV, "a", 0.0 );
        if( dfSemiMajor == 0.0 )
        {
            dfSemiMajor = OSR_GDV( papszNV, "R", 0.0 );
            if( dfSemiMajor != 0.0 )
            {
                dfSemiMinor = -1.0;
                dfInvFlattening = 0.0;
            }
            else
            {
                CPLDebug( "OGR_PROJ4", "Can't find ellipse definition, default to WGS84:\n%s", 
                          pszProj4 );
                
                dfSemiMajor = SRS_WGS84_SEMIMAJOR;
                dfSemiMinor = -1.0;
                dfInvFlattening = SRS_WGS84_INVFLATTENING;
            }
        }
        else
        {
            dfSemiMinor = OSR_GDV( papszNV, "b", -1.0 );
            dfInvFlattening = OSR_GDV( papszNV, "rf", -1.0 );
        }
        
        if( dfSemiMinor == -1.0 && dfInvFlattening == -1.0 )
        {
            CPLDebug( "OGR_PROJ4", "Can't find ellipse definition in:\n%s", 
                      pszProj4 );
            CSLDestroy( papszNV );
            return OGRERR_UNSUPPORTED_SRS;
        }

        if( dfInvFlattening == -1.0 )
        {
            if( ABS(dfSemiMajor/dfSemiMinor) - 1.0 < 0.0000000000001 )
                dfInvFlattening = 0.0;
            else
                dfInvFlattening = -1.0 / (dfSemiMinor/dfSemiMajor - 1.0);
        }
        
        SetGeogCS( "unnamed ellipse", "unknown", "unnamed",
                   dfSemiMajor, dfInvFlattening,
                   pszPM, dfFromGreenwich );
        
        bFullyDefined = TRUE;
    }

/* -------------------------------------------------------------------- */
/*      Handle TOWGS84 conversion.                                      */
/* -------------------------------------------------------------------- */
    pszValue = CSLFetchNameValue(papszNV, "towgs84");
    if(pszValue!=NULL)
    {
        char **papszToWGS84 = CSLTokenizeStringComplex( pszValue, ",", 
                                                        FALSE, TRUE );

        if( CSLCount(papszToWGS84) >= 7 )
            SetTOWGS84( atof(papszToWGS84[0]), 
                        atof(papszToWGS84[1]), 
                        atof(papszToWGS84[2]), 
                        atof(papszToWGS84[3]), 
                        atof(papszToWGS84[4]), 
                        atof(papszToWGS84[5]), 
                        atof(papszToWGS84[6]) );
        else if( CSLCount(papszToWGS84) >= 3 )
            SetTOWGS84( atof(papszToWGS84[0]), 
                        atof(papszToWGS84[1]), 
                        atof(papszToWGS84[2]) );
        else
            CPLError( CE_Warning, CPLE_AppDefined, 
                      "Seemingly corrupt +towgs84 option (%s), ignoring.", 
                      pszValue );
                        
        CSLDestroy(papszToWGS84);
    }

/* -------------------------------------------------------------------- */
/*      Linear units translation                                        */
/* -------------------------------------------------------------------- */
    if( IsProjected() || IsLocal() )
    {
        pszValue = CSLFetchNameValue(papszNV, "to_meter");

        if( pszValue != NULL && atof(pszValue) > 0.0 )
        {
            SetLinearUnits( "unknown", atof(pszValue) );
        }
        else if( (pszValue = CSLFetchNameValue(papszNV, "units")) != NULL )
        {
            if( EQUAL(pszValue,"meter" ) || EQUAL(pszValue,"m") )
                SetLinearUnits( SRS_UL_METER, 1.0 );
            else if( EQUAL(pszValue,"us-ft" ) )
                SetLinearUnits( SRS_UL_US_FOOT, atof(SRS_UL_US_FOOT_CONV) );
            else if( EQUAL(pszValue,"ft" ) )
                SetLinearUnits( SRS_UL_FOOT, atof(SRS_UL_FOOT_CONV) );
            else if( EQUAL(pszValue,"yd" ) )
                SetLinearUnits( pszValue, 0.9144 );
            else if( EQUAL(pszValue,"us-yd" ) )
                SetLinearUnits( pszValue, 0.914401828803658 );
            else // This case is untranslatable.  Should add all proj.4 unts
                SetLinearUnits( pszValue, 1.0 );
        }
    }

/* -------------------------------------------------------------------- */
/*      Adjust linear parameters into PROJCS units if the linear        */
/*      units are not meters.                                           */
/* -------------------------------------------------------------------- */
    if( GetLinearUnits() != 1.0 && IsProjected() )
    {
        OGR_SRSNode *poPROJCS = GetAttrNode( "PROJCS" );
        int  i;

        for( i = 0; i < poPROJCS->GetChildCount(); i++ )
        {
            OGR_SRSNode *poParm = poPROJCS->GetChild(i);
            if( !EQUAL(poParm->GetValue(),"PARAMETER") 
                || poParm->GetChildCount() != 2 )
                continue;

            const char *pszParmName = poParm->GetChild(0)->GetValue();

            if( IsLinearParameter(pszParmName) )
                SetNormProjParm(pszParmName,GetProjParm(pszParmName));
        }        
    }


/* -------------------------------------------------------------------- */
/*      do we want to insert a PROJ.4 EXTENSION item?                   */
/* -------------------------------------------------------------------- */
    if( strstr(pszProj4,"wktext") != NULL )
        SetExtension( GetRoot()->GetValue(), "PROJ4", pszProj4 );
        
    CSLDestroy( papszNV );
    
    return OGRERR_NONE;
}


/************************************************************************/
/*                          OSRExportToProj4()                          */
/************************************************************************/

OGRErr CPL_STDCALL OSRExportToProj4( OGRSpatialReferenceH hSRS, 
                                     char ** ppszReturn )

{
    *ppszReturn = NULL;

    return ((OGRSpatialReference *) hSRS)->exportToProj4( ppszReturn );
}

/************************************************************************/
/*                           exportToProj4()                            */
/************************************************************************/

/**
 * Export coordinate system in PROJ.4 format.
 *
 * Converts the loaded coordinate reference system into PROJ.4 format
 * to the extent possible.  The string returned in ppszProj4 should be
 * deallocated by the caller with CPLFree() when no longer needed.
 *
 * LOCAL_CS coordinate systems are not translatable.  An empty string
 * will be returned along with OGRERR_NONE.  
 *
 * This method is the equivelent of the C function OSRExportToProj4().
 *
 * @param ppszProj4 pointer to which dynamically allocated PROJ.4 definition 
 * will be assigned. 
 *
 * @return OGRERR_NONE on success or an error code on failure. 
 */

OGRErr OGRSpatialReference::exportToProj4( char ** ppszProj4 ) const

{
    char        szProj4[512];
    const char *pszProjection = GetAttrValue("PROJECTION");

    szProj4[0] = '\0';

    if( GetRoot() == NULL )
    {
        *ppszProj4 = CPLStrdup("");
        CPLError( CE_Failure, CPLE_NotSupported,
                  "No translation an empty SRS to PROJ.4 format is known.", 
                  pszProjection );
        return OGRERR_UNSUPPORTED_SRS;
    }

/* -------------------------------------------------------------------- */
/*      Do we have a PROJ.4 override definition?                        */
/* -------------------------------------------------------------------- */
    const char *pszPredefProj4 = GetExtension( GetRoot()->GetValue(), 
                                               "PROJ4", NULL );
    if( pszPredefProj4 != NULL )
    {
        *ppszProj4 = CPLStrdup( pszPredefProj4 );
        return OGRERR_NONE;
    }

/* -------------------------------------------------------------------- */
/*      Get the prime meridian info.                                    */
/* -------------------------------------------------------------------- */
    const OGR_SRSNode *poPRIMEM = GetAttrNode( "PRIMEM" );
    double dfFromGreenwich = 0.0;

    if( poPRIMEM != NULL && poPRIMEM->GetChildCount() >= 2 
        && atof(poPRIMEM->GetChild(1)->GetValue()) != 0.0 )
    {
        dfFromGreenwich = atof(poPRIMEM->GetChild(1)->GetValue());
    }

/* ==================================================================== */
/*      Handle the projection definition.                               */
/* ==================================================================== */

    if( pszProjection == NULL && IsGeographic() )
    {
        sprintf( szProj4+strlen(szProj4), "+proj=longlat " );