ogr_fromepsg.cpp

mitab,读取MapInfo的地图文件

 *
 * This method is the same as the C function OSRImportFromEPSG().
 *
 * @param nCode a GCS or PCS code from the horizontal coordinate system table.
 * 
 * @return OGRERR_NONE on success, or an error code on failure.
 */

OGRErr OGRSpatialReference::importFromEPSG( int nCode )

{
    OGRErr  eErr;

    bNormInfoSet = FALSE;

/* -------------------------------------------------------------------- */
/*      Clear any existing definition.                                  */
/* -------------------------------------------------------------------- */
    if( GetRoot() != NULL )
    {
        delete poRoot;
        poRoot = NULL;
    }

/* -------------------------------------------------------------------- */
/*      Verify that we can find the required filename(s).               */
/* -------------------------------------------------------------------- */
    if( CSVScanFileByName( CSVFilename( "gcs.csv" ),
                           "COORD_REF_SYS_CODE", 
                           "4269", CC_Integer ) == NULL )
    {
        CPLError( CE_Failure, CPLE_OpenFailed, 
                  "Unable to open EPSG support file %s.\n"
                  "Try setting the GDAL_DATA environment variable to point to the\n"
                  "directory containing EPSG csv files.", 
                  CSVFilename( "gcs.csv" ) );
        return OGRERR_FAILURE;
    }

/* -------------------------------------------------------------------- */
/*      Is this a GeogCS code?   this is inadequate as a criteria       */
/* -------------------------------------------------------------------- */
    if( EPSGGetGCSInfo( nCode, NULL, NULL, NULL, NULL, NULL, NULL ) )
        eErr = SetEPSGGeogCS( this, nCode );
    else
        eErr = SetEPSGProjCS( this, nCode );

/* -------------------------------------------------------------------- */
/*      If we get it as an unsupported code, try looking it up in       */
/*      the epsg.wkt coordinate system dictionary.                      */
/* -------------------------------------------------------------------- */
    if( eErr == OGRERR_UNSUPPORTED_SRS )
    {
        char szCode[32];
        sprintf( szCode, "%d", nCode );
        eErr = importFromDict( "epsg.wkt", szCode );
    }

/* -------------------------------------------------------------------- */
/*      If we get it as an unsupported code, try looking it up in       */
/*      the PROJ.4 support file(s).                                     */
/* -------------------------------------------------------------------- */
    static int bLoopingToProj4 = FALSE;

    if( eErr == OGRERR_UNSUPPORTED_SRS && !bLoopingToProj4 )
    {
        char szWrkDefn[100];

        sprintf( szWrkDefn, "+init=epsg:%d", nCode );
        bLoopingToProj4 = TRUE;
        CPLPushErrorHandler( CPLQuietErrorHandler );
        eErr = SetFromUserInput( szWrkDefn );
        CPLPopErrorHandler();
        bLoopingToProj4 = FALSE;

        if( eErr != OGRERR_NONE )
            eErr = OGRERR_UNSUPPORTED_SRS;
        else
        {
            if( IsProjected() )
                SetAuthority( "PROJCS", "EPSG", nCode );
            else if( IsGeographic() )
                SetAuthority( "GEOGCS", "EPSG", nCode );
        }
    }

    if( eErr == OGRERR_UNSUPPORTED_SRS )
    {
        CPLError( CE_Failure, CPLE_NotSupported,
                  "EPSG PCS/GCS code %d not found in EPSG support files.  Is this a valid\nEPSG coordinate system?", 
                  nCode );
    }

/* -------------------------------------------------------------------- */
/*      Make sure any peculiarities in the ordering are fixed up.       */
/* -------------------------------------------------------------------- */
    if( eErr == OGRERR_NONE )
        eErr = FixupOrdering();

    return eErr;
}

/************************************************************************/
/*                         OSRImportFromEPSG()                          */
/************************************************************************/

OGRErr CPL_STDCALL OSRImportFromEPSG( OGRSpatialReferenceH hSRS, int nCode )

{
    return ((OGRSpatialReference *) hSRS)->importFromEPSG( nCode );
}

/************************************************************************/
/*                           SetStatePlane()                            */
/************************************************************************/

/**
 * Set State Plane projection definition.
 *
 * This will attempt to generate a complete definition of a state plane
 * zone based on generating the entire SRS from the EPSG tables.  If the
 * EPSG tables are unavailable, it will produce a stubbed LOCAL_CS definition
 * and return OGRERR_FAILURE.
 *
 * This method is the same as the C function OSRSetStatePlaneWithUnits().
 *
 * @param nZone State plane zone number, in the USGS numbering scheme (as
 * dinstinct from the Arc/Info and Erdas numbering scheme. 
 *
 * @param bNAD83 TRUE if the NAD83 zone definition should be used or FALSE
 * if the NAD27 zone definition should be used.  
 *
 * @param pszOverrideUnitName Linear unit name to apply overriding the 
 * legal definition for this zone.
 *
 * @param dfOverrideUnit Linear unit conversion factor to apply overriding
 * the legal definition for this zone. 
 * 
 * @return OGRERR_NONE on success, or OGRERR_FAILURE on failure, mostly likely
 * due to the EPSG tables not being accessable. 
 */

OGRErr OGRSpatialReference::SetStatePlane( int nZone, int bNAD83,
                                           const char *pszOverrideUnitName,
                                           double dfOverrideUnit )

{
    int         nAdjustedId;
    int         nPCSCode;
    char        szID[32];

/* -------------------------------------------------------------------- */
/*      Get the index id from stateplane.csv.                           */
/* -------------------------------------------------------------------- */
    if( bNAD83 )
        nAdjustedId = nZone;
    else
        nAdjustedId = nZone + 10000;

/* -------------------------------------------------------------------- */
/*      Turn this into a PCS code.  We assume there will only be one    */
/*      PCS corresponding to each Proj_ code since the proj code        */
/*      already effectively indicates NAD27 or NAD83.                   */
/* -------------------------------------------------------------------- */
    sprintf( szID, "%d", nAdjustedId );
    nPCSCode =
        atoi( CSVGetField( CSVFilename( "stateplane.csv" ),
                           "ID", szID, CC_Integer,
                           "EPSG_PCS_CODE" ) );
    if( nPCSCode < 1 )
    {
        char    szName[128];
        static int bFailureReported = FALSE;

        if( !bFailureReported )
        {
            bFailureReported = TRUE;
            CPLError( CE_Warning, CPLE_OpenFailed, 
                      "Unable to find state plane zone in stateplane.csv,\n"
                      "likely because the GDAL data files cannot be found.  Using\n"
                      "incomplete definition of state plane zone.\n" );
        }

        Clear();
        if( bNAD83 )
        {
            sprintf( szName, "State Plane Zone %d / NAD83", nZone );
            SetLocalCS( szName );
            SetLinearUnits( SRS_UL_METER, 1.0 );
        }
        else
        {
            sprintf( szName, "State Plane Zone %d / NAD27", nZone );
            SetLocalCS( szName );
            SetLinearUnits( SRS_UL_US_FOOT, atof(SRS_UL_US_FOOT_CONV) );
        }

        return OGRERR_FAILURE;
    }

/* -------------------------------------------------------------------- */
/*      Define based on a full EPSG definition of the zone.             */
/* -------------------------------------------------------------------- */
    OGRErr eErr = importFromEPSG( nPCSCode );

    if( eErr != OGRERR_NONE )
        return eErr;

/* -------------------------------------------------------------------- */
/*      Apply units override if required.                               */
/*                                                                      */
/*      We will need to adjust the linear projection parameter to       */
/*      match the provided units, and clear the authority code.         */
/* -------------------------------------------------------------------- */
    if( dfOverrideUnit != 0.0 
        && fabs(dfOverrideUnit - GetLinearUnits()) > 0.0000000001 )
    {
        double dfFalseEasting = GetNormProjParm( SRS_PP_FALSE_EASTING );
        double dfFalseNorthing= GetNormProjParm( SRS_PP_FALSE_NORTHING);
        OGR_SRSNode *poPROJCS;

        SetLinearUnits( pszOverrideUnitName, dfOverrideUnit );
        
        SetNormProjParm( SRS_PP_FALSE_EASTING, dfFalseEasting );
        SetNormProjParm( SRS_PP_FALSE_NORTHING, dfFalseNorthing );

        poPROJCS = GetAttrNode( "PROJCS" );
        if( poPROJCS != NULL && poPROJCS->FindChild( "AUTHORITY" ) != -1 )
        {
            poPROJCS->DestroyChild( poPROJCS->FindChild( "AUTHORITY" ) );
        }
    }

    return OGRERR_NONE;
}

/************************************************************************/
/*                          OSRSetStatePlane()                          */
/************************************************************************/

OGRErr OSRSetStatePlane( OGRSpatialReferenceH hSRS, int nZone, int bNAD83 )

{
    return ((OGRSpatialReference *) hSRS)->SetStatePlane( nZone, bNAD83 );
}

/************************************************************************/
/*                     OSRSetStatePlaneWithUnits()                      */
/************************************************************************/

OGRErr OSRSetStatePlaneWithUnits( OGRSpatialReferenceH hSRS, 
                                  int nZone, int bNAD83,
                                  const char *pszOverrideUnitName,
                                  double dfOverrideUnit )

{
    return ((OGRSpatialReference *) hSRS)->SetStatePlane( nZone, bNAD83,
                                                          pszOverrideUnitName,
                                                          dfOverrideUnit );
}

/************************************************************************/
/*                          AutoIdentifyEPSG()                          */
/************************************************************************/

/**
 * Set EPSG authority info if possible.
 *
 * This method inspects a WKT definition, and adds EPSG authority nodes
 * where an aspect of the coordinate system can be easily and safely 
 * corresponded with an EPSG identifier.  In practice, this method will 
 * evolve over time.  In theory it can add authority nodes for any object
 * (ie. spheroid, datum, GEOGCS, units, and PROJCS) that could have an 
 * authority node.  Mostly this is useful to inserting appropriate 
 * PROJCS codes for common formulations (like UTM n WGS84). 
 *
 * If it success the OGRSpatialReference is updated in place, and the 
 * method return OGRERR_NONE.  If the method fails to identify the 
 * general coordinate system OGRERR_UNSUPPORTED_SRS is returned but no 
 * error message is posted via CPLError(). 
 *
 * This method is the same as the C function OSRAutoIdentifyEPSG().
 *
 * @return OGRERR_NONE or OGRERR_UNSUPPORTED_SRS.
 */

OGRErr OGRSpatialReference::AutoIdentifyEPSG()

{
/* -------------------------------------------------------------------- */
/*      Is this a UTM coordinate system with a common GEOGCS?           */
/* -------------------------------------------------------------------- */
    int nZone, bNorth;
    if( (nZone = GetUTMZone( &bNorth )) != 0 
        && GetAuthorityCode( "PROJCS") == NULL )
    {
        const char *pszAuthName, *pszAuthCode;

        pszAuthName = GetAuthorityName( "PROJCS|GEOGCS" );
        pszAuthCode = GetAuthorityCode( "PROJCS|GEOGCS" );

        if( pszAuthName == NULL ||  pszAuthCode == NULL )
        {
            /* don't exactly recognise datum */
        }
        else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4326 )
        { // WGS84
            if( bNorth ) 
                SetAuthority( "PROJCS", "EPSG", 32600 + nZone );
            else
                SetAuthority( "PROJCS", "EPSG", 32700 + nZone );
        }
        else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4267 
                 && nZone >= 3 && nZone <= 22 && bNorth )
            SetAuthority( "PROJCS", "EPSG", 26700 + nZone ); // NAD27
        else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4269
                 && nZone >= 3 && nZone <= 23 && bNorth )
            SetAuthority( "PROJCS", "EPSG", 26900 + nZone ); // NAD83
        else if( EQUAL(pszAuthName,"EPSG") && atoi(pszAuthCode) == 4322 )
        { // WGS72
            if( bNorth ) 
                SetAuthority( "PROJCS", "EPSG", 32200 + nZone );
            else
                SetAuthority( "PROJCS", "EPSG", 32300 + nZone );
        }
    }

/* -------------------------------------------------------------------- */
/*      Return.                                                         */
/* ---