ogr_fromepsg.cpp

在linux环境下

#define PseudoStdParallelLat 8818
#define PseudoStdParallelScaleFactor 8819
#define FalseOriginLat       8821
#define FalseOriginLong      8822
#define StdParallel1Lat      8823
#define StdParallel2Lat      8824
#define FalseOriginEasting   8826
#define FalseOriginNorthing  8827
#define SphericalOriginLat   8828
#define SphericalOriginLong  8829
#define InitialLongitude     8830
#define ZoneWidth            8831

/************************************************************************/
/*                         EPSGGetProjTRFInfo()                         */
/*                                                                      */
/*      Transform a PROJECTION_TRF_CODE into a projection method,       */
/*      and a set of parameters.  The parameters identify will          */
/*      depend on the returned method, but they will all have been      */
/*      normalized into degrees and meters.                             */
/************************************************************************/

static int
EPSGGetProjTRFInfo( int nPCS, int * pnProjMethod,
                    int *panParmIds, double * padfProjParms )

{
    int         nProjMethod, i;
    double      adfProjParms[7];
    char        szTRFCode[16];
    char       *pszFilename = CPLStrdup(CSVFilename( "pcs.csv" ));

/* -------------------------------------------------------------------- */
/*      Get the proj method.  If this fails to return a meaningful      */
/*      number, then the whole function fails.                          */
/* -------------------------------------------------------------------- */
    sprintf( szTRFCode, "%d", nPCS );
    nProjMethod =
        atoi( CSVGetField( pszFilename,
                           "COORD_REF_SYS_CODE", szTRFCode, CC_Integer,
                           "COORD_OP_METHOD_CODE" ) );
    if( nProjMethod == 0 )
    {
        CPLFree( pszFilename );
        return FALSE;
    }

/* -------------------------------------------------------------------- */
/*      Get the parameters for this projection.                         */
/* -------------------------------------------------------------------- */

    for( i = 0; i < 7; i++ )
    {
        char    szParamUOMID[32], szParamValueID[32], szParamCodeID[32];
        char    *pszValue;
        int     nUOM;
        
        sprintf( szParamCodeID, "PARAMETER_CODE_%d", i+1 );
        sprintf( szParamUOMID, "PARAMETER_UOM_%d", i+1 );
        sprintf( szParamValueID, "PARAMETER_VALUE_%d", i+1 );

        if( panParmIds != NULL )
            panParmIds[i] = 
                atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szTRFCode,
                                  CC_Integer, szParamCodeID ));

        nUOM = atoi(CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szTRFCode,
                                 CC_Integer, szParamUOMID ));
        pszValue = CPLStrdup(
            CSVGetField( pszFilename, "COORD_REF_SYS_CODE", szTRFCode,
                         CC_Integer, szParamValueID ));

        // there is a bug in the EPSG 6.2.2 database for PCS 2935 and 2936
        // such that they have foot units for the scale factor.  Avoid this.
        if( (panParmIds[i] == NatOriginScaleFactor 
             || panParmIds[i] == InitialLineScaleFactor
             || panParmIds[i] == PseudoStdParallelScaleFactor) 
            && nUOM < 9200 )
            nUOM = 9201;

        if( nUOM >= 9100 && nUOM < 9200 )
            adfProjParms[i] = EPSGAngleStringToDD( pszValue, nUOM );
        else if( nUOM > 9000 && nUOM < 9100 )
        {
            double dfInMeters;

            if( !EPSGGetUOMLengthInfo( nUOM, NULL, &dfInMeters ) )
                dfInMeters = 1.0;
            adfProjParms[i] = atof(pszValue) * dfInMeters;
        }
        else if( EQUAL(pszValue,"") ) /* null field */
        {
            adfProjParms[i] = 0.0;
        }
        else /* really we should consider looking up other scaling factors */
        {
            if( nUOM != 9201 )
                CPLDebug( "OGR", "Non-unity scale factor units!" );
            adfProjParms[i] = atof(pszValue);
        }

        CPLFree( pszValue );
    }

/* -------------------------------------------------------------------- */
/*      Transfer requested data into passed variables.                  */
/* -------------------------------------------------------------------- */
    if( pnProjMethod != NULL )
        *pnProjMethod = nProjMethod;

    if( padfProjParms != NULL )
    {
        for( i = 0; i < 7; i++ )
            padfProjParms[i] = adfProjParms[i];
    }

    CPLFree( pszFilename );

    return TRUE;
}


/************************************************************************/
/*                           EPSGGetPCSInfo()                           */
/************************************************************************/

static int 
EPSGGetPCSInfo( int nPCSCode, char **ppszEPSGName,
                int *pnUOMLengthCode, int *pnUOMAngleCode,
                int *pnGeogCS, int *pnTRFCode )

{
    char        **papszRecord;
    char        szSearchKey[24];
    const char  *pszFilename = CSVFilename( "pcs.csv" );
    
/* -------------------------------------------------------------------- */
/*      Search the units database for this unit.  If we don't find      */
/*      it return failure.                                              */
/* -------------------------------------------------------------------- */
    sprintf( szSearchKey, "%d", nPCSCode );
    papszRecord = CSVScanFileByName( pszFilename, "COORD_REF_SYS_CODE",
                                     szSearchKey, CC_Integer );

    if( papszRecord == NULL )
        return FALSE;

/* -------------------------------------------------------------------- */
/*      Get the name, if requested.                                     */
/* -------------------------------------------------------------------- */
    if( ppszEPSGName != NULL )
    {
        *ppszEPSGName =
            CPLStrdup( CSLGetField( papszRecord,
                                    CSVGetFileFieldId(pszFilename,
                                                      "COORD_REF_SYS_NAME") ));
    }

/* -------------------------------------------------------------------- */
/*      Get the UOM Length code, if requested.                          */
/* -------------------------------------------------------------------- */
    if( pnUOMLengthCode != NULL )
    {
        const char      *pszValue;

        pszValue =
            CSLGetField( papszRecord,
                         CSVGetFileFieldId(pszFilename,"UOM_CODE"));
        if( atoi(pszValue) > 0 )
            *pnUOMLengthCode = atoi(pszValue);
        else
            *pnUOMLengthCode = 0;
    }

/* -------------------------------------------------------------------- */
/*      Get the UOM Angle code, if requested.                           */
/* -------------------------------------------------------------------- */
    if( pnUOMAngleCode != NULL )
    {
        const char      *pszValue;
        
        pszValue =
            CSLGetField( papszRecord,
                         CSVGetFileFieldId(pszFilename,"UOM_ANGLE_CODE") );
        
        if( atoi(pszValue) > 0 )
            *pnUOMAngleCode = atoi(pszValue);
        else
            *pnUOMAngleCode = 0;
    }

/* -------------------------------------------------------------------- */
/*      Get the GeogCS (Datum with PM) code, if requested.              */
/* -------------------------------------------------------------------- */
    if( pnGeogCS != NULL )
    {
        const char      *pszValue;

        pszValue =
            CSLGetField( papszRecord,
                         CSVGetFileFieldId(pszFilename,"SOURCE_GEOGCRS_CODE"));
        if( atoi(pszValue) > 0 )
            *pnGeogCS = atoi(pszValue);
        else
            *pnGeogCS = 0;
    }

/* -------------------------------------------------------------------- */
/*      Get the GeogCS (Datum with PM) code, if requested.              */
/* -------------------------------------------------------------------- */
    if( pnTRFCode != NULL )
    {
        const char      *pszValue;

        pszValue =
            CSLGetField( papszRecord,
                         CSVGetFileFieldId(pszFilename,"COORD_OP_CODE"));
                         
        
        if( atoi(pszValue) > 0 )
            *pnTRFCode = atoi(pszValue);
        else
            *pnTRFCode = 0;
    }

    return TRUE;
}

/************************************************************************/
/*                           SetEPSGGeogCS()                            */
/*                                                                      */
/*      FLAWS:                                                          */
/*       o Units are all hardcoded.                                     */
/*       o Axis hardcoded.                                              */
/************************************************************************/

static OGRErr SetEPSGGeogCS( OGRSpatialReference * poSRS, int nGeogCS )

{
    int  nDatumCode, nPMCode, nUOMAngle, nEllipsoidCode;
    char *pszGeogCSName = NULL, *pszDatumName = NULL, *pszEllipsoidName = NULL;
    char *pszPMName = NULL, *pszAngleName = NULL;
    double dfPMOffset, dfSemiMajor, dfInvFlattening, adfBursaTransform[7];
    double dfAngleInDegrees, dfAngleInRadians;

    if( !EPSGGetGCSInfo( nGeogCS, &pszGeogCSName,
                         &nDatumCode, &pszDatumName, 
                         &nPMCode, &nEllipsoidCode, &nUOMAngle ) )
        return OGRERR_UNSUPPORTED_SRS;

    if( !EPSGGetPMInfo( nPMCode, &pszPMName, &dfPMOffset ) )
        return OGRERR_UNSUPPORTED_SRS;

    OGREPSGDatumNameMassage( &pszDatumName );

    if( !EPSGGetEllipsoidInfo( nEllipsoidCode, &pszEllipsoidName, 
                               &dfSemiMajor, &dfInvFlattening ) )
        return OGRERR_UNSUPPORTED_SRS;

    if( !EPSGGetUOMAngleInfo( nUOMAngle, &pszAngleName, &dfAngleInDegrees ) )
    {
        pszAngleName = CPLStrdup("degree");
        dfAngleInDegrees = 1.0;
        nUOMAngle = -1;
    }

    if( dfAngleInDegrees == 1.0 )
        dfAngleInRadians = atof(SRS_UA_DEGREE_CONV);
    else
        dfAngleInRadians = atof(SRS_UA_DEGREE_CONV) * dfAngleInDegrees;

    poSRS->SetGeogCS( pszGeogCSName, pszDatumName, 
                      pszEllipsoidName, dfSemiMajor, dfInvFlattening,
                      pszPMName, dfPMOffset,
                      pszAngleName, dfAngleInRadians );

    if( EPSGGetWGS84Transform( nGeogCS, adfBursaTransform ) )
    {
        OGR_SRSNode     *poWGS84;
        char            szValue[48];

        poWGS84 = new OGR_SRSNode( "TOWGS84" );

        for( int iCoeff = 0; iCoeff < 7; iCoeff++ )
        {
            sprintf( szValue, "%g", adfBursaTransform[iCoeff] );
            poWGS84->AddChild( new OGR_SRSNode( szValue ) );
        }

        poSRS->GetAttrNode( "DATUM" )->AddChild( poWGS84 );
    }

    poSRS->SetAuthority( "GEOGCS", "EPSG", nGeogCS );
    poSRS->SetAuthority( "DATUM", "EPSG", nDatumCode );
    poSRS->SetAuthority( "SPHEROID", "EPSG", nEllipsoidCode );
    poSRS->SetAuthority( "PRIMEM", "EPSG", nPMCode );

    if( nUOMAngle > 0 )
        poSRS->SetAuthority( "GEOGCS|UNIT", "EPSG", nUOMAngle );

    CPLFree( pszAngleName );
    CPLFree( pszDatumName );
    CPLFree( pszEllipsoidName );
    CPLFree( pszGeogCSName );
    CPLFree( pszPMName );

    return OGRERR_NONE;
}

/************************************************************************/
/*                           OGR_FetchParm()                            */
/*                                                                      */
/*      Fetch a parameter from the parm list, based on it's EPSG        */
/*      parameter code.                                                 */
/************************************************************************/

static double OGR_FetchParm( double *padfProjParms, int *panParmIds, 
                             int nTargetId, double dfFromGreenwich )

{
    int i;
    double dfResult;

/* -------------------------------------------------------------------- */
/*      Set default in meters/degrees.                                  */
/* -------------------------------------------------------------------- */
    switch( nTargetId )
    {
      case NatOriginScaleFactor:
      case InitialLineScaleFactor:
      case PseudoStdParallelScaleFactor:
        dfResult = 1.0;
        break;
        
      case AngleRectifiedToSkewedGrid:
        dfResult = 90.0;
        break;

      default:
        dfResult = 0.0;
    }

/* -------------------------------------------------------------------- */
/*      Try to find actual value in parameter list.                     */
/* -------------------------------------------------------------------- */
    for( i = 0; i < 7; i++ )
    {
        if( panParmIds[i] == nTargetId )
        {
            dfResult = padfProjParms[i];
            break;
        }
    }

/* -------------------------------------------------------------------- */
/*      Correct longitude to be relative to greenwich                   */
/* -------------------------------------------------------------------- */
    switch( nTargetId )
    {
      case NatOriginLong:
      case ProjCenterLong:
      case FalseOriginLong:
      case SphericalOriginLong:
      case InitialLongitude:
        dfResult = dfResult + dfFromGreenwich;
        break;

      default:
        ;
    }

    return dfResult;
}

#define OGR_FP(x) OGR_FetchParm( adfProjParms, anParmIds, (x), \
                                 dfFromGreenwich )

/************************************************************************/
/*                           SetEPSGProjCS()                            */
/************************************************************************/

static OGRErr SetEPSGProjCS( OGRSpatialReference * poSRS, int nPCSCode )

{
    int         nGCSCode, nUOMAngleCode, nUOMLength, nTRFCode, nProjMethod;
    int         anParmIds[7];
    char        *pszPCSName = NULL, *pszUOMLengthName = NULL;
    double      adfProjParms[7], dfInMeters, dfFromGreenwich;
    OGRErr      nErr;
    OGR_SRSNode *poNode;

    if( !EPSGGetPCSInfo( nPCSCode, &pszPCSName, &nUOMLength, &nUOMAngleCode,
                         &nGCSCode, &nTRFCode ) )
        return OGRERR_UNSUPPORTED_SRS;