ogr_srs_esri.cpp

用于读取TAB、MIF、SHP文件的类

/* -------------------------------------------------------------------- */
/*      ArcGIS and related products now use a varient of Well Known     */
/*      Text.  Try to recognise this and ingest it.  WKT is usually     */
/*      all on one line, but we will accept multi-line formats and      */
/*      concatenate.                                                    */
/* -------------------------------------------------------------------- */
    if( EQUALN(papszPrj[0],"GEOGCS",6)
        || EQUALN(papszPrj[0],"PROJCS",6)
        || EQUALN(papszPrj[0],"LOCAL_CS",8) )
    {
        char    *pszWKT, *pszWKT2;
        OGRErr  eErr;
        int     i;

        pszWKT = CPLStrdup(papszPrj[0]);
        for( i = 1; papszPrj[i] != NULL; i++ )
        {
            pszWKT = (char *) 
                CPLRealloc(pszWKT,strlen(pszWKT)+strlen(papszPrj[i])+1);
            strcat( pszWKT, papszPrj[i] );
        }
        pszWKT2 = pszWKT;
        eErr = importFromWkt( &pszWKT2 );
        CPLFree( pszWKT );

        if( eErr == OGRERR_NONE )
            eErr = morphFromESRI();
        return eErr;
    }

/* -------------------------------------------------------------------- */
/*      Operate on the basis of the projection name.                    */
/* -------------------------------------------------------------------- */
    const char *pszProj = OSR_GDS( papszPrj, "Projection", NULL );

    if( pszProj == NULL )
    {
        CPLDebug( "OGR_ESRI", "Can't find Projection\n" );
        return OGRERR_CORRUPT_DATA;
    }

    else if( EQUAL(pszProj,"GEOGRAPHIC") )
    {
    }
    
    else if( EQUAL(pszProj,"utm") )
    {
        if( (int) OSR_GDV( papszPrj, "zone", 0.0 ) != 0 )
        {
            double      dfYShift = OSR_GDV( papszPrj, "Yshift", 0.0 );

            SetUTM( (int) OSR_GDV( papszPrj, "zone", 0.0 ),
                    dfYShift == 0.0 );
        }
        else
        {
            double      dfCentralMeridian, dfRefLat;
            int         nZone;

            dfCentralMeridian = OSR_GDV( papszPrj, "PARAM_1", 0.0 );
            dfRefLat = OSR_GDV( papszPrj, "PARAM_2", 0.0 );

            nZone = (int) ((dfCentralMeridian+183) / 6.0 + 0.0000001);
            SetUTM( nZone, dfRefLat >= 0.0 );
        }
    }

    else if( EQUAL(pszProj,"STATEPLANE") )
    {
        int nZone = (int) OSR_GDV( papszPrj, "zone", 0.0 );
        if( nZone != 0 )
            nZone = ESRIToUSGSZone( nZone );
        else
            nZone = (int) OSR_GDV( papszPrj, "fipszone", 0.0 );

        if( nZone != 0 )
        {
            if( EQUAL(OSR_GDS( papszPrj, "Datum", "NAD83"),"NAD27") )
                SetStatePlane( nZone, FALSE );
            else
                SetStatePlane( nZone, TRUE );
        }
    }

    else if( EQUAL(pszProj,"GREATBRITIAN_GRID") 
             || EQUAL(pszProj,"GREATBRITAIN_GRID") )
    {
        const char *pszWkt = 
            "PROJCS[\"OSGB 1936 / British National Grid\",GEOGCS[\"OSGB 1936\",DATUM[\"OSGB_1936\",SPHEROID[\"Airy 1830\",6377563.396,299.3249646]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",49],PARAMETER[\"central_meridian\",-2],PARAMETER[\"scale_factor\",0.999601272],PARAMETER[\"false_easting\",400000],PARAMETER[\"false_northing\",-100000],UNIT[\"metre\",1]]";

        importFromWkt( (char **) &pszWkt );
    }

    else if( EQUAL(pszProj,"ALBERS") )
    {
        SetACEA( OSR_GDV( papszPrj, "PARAM_1", 0.0 ), 
                 OSR_GDV( papszPrj, "PARAM_2", 0.0 ), 
                 OSR_GDV( papszPrj, "PARAM_4", 0.0 ), 
                 OSR_GDV( papszPrj, "PARAM_3", 0.0 ), 
                 OSR_GDV( papszPrj, "PARAM_5", 0.0 ), 
                 OSR_GDV( papszPrj, "PARAM_6", 0.0 ) );
    }

    else if( EQUAL(pszProj,"LAMBERT") )
    {
        SetLCC( OSR_GDV( papszPrj, "PARAM_1", 0.0 ),
                OSR_GDV( papszPrj, "PARAM_2", 0.0 ),
                OSR_GDV( papszPrj, "PARAM_4", 0.0 ),
                OSR_GDV( papszPrj, "PARAM_3", 0.0 ),
                OSR_GDV( papszPrj, "PARAM_5", 0.0 ),
                OSR_GDV( papszPrj, "PARAM_6", 0.0 ) );
    }

    else if( EQUAL(pszProj,"EQUIDISTANT_CONIC") )
    {
        int     nStdPCount = (int) OSR_GDV( papszPrj, "PARAM_1", 0.0 );

        if( nStdPCount == 1 )
        {
            SetEC( OSR_GDV( papszPrj, "PARAM_2", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_2", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_4", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_3", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_5", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_6", 0.0 ) );
        }
        else
        {
            SetEC( OSR_GDV( papszPrj, "PARAM_2", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_3", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_5", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_4", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_5", 0.0 ), 
                   OSR_GDV( papszPrj, "PARAM_7", 0.0 ) );
        }
    }

    else if( EQUAL(pszProj,"TRANSVERSE") )
    {
        SetTM( OSR_GDV( papszPrj, "PARAM_2", 0.0 ), 
               OSR_GDV( papszPrj, "PARAM_3", 0.0 ), 
               OSR_GDV( papszPrj, "PARAM_1", 0.0 ), 
               OSR_GDV( papszPrj, "PARAM_4", 0.0 ), 
               OSR_GDV( papszPrj, "PARAM_5", 0.0 ) );
    }

    else if( EQUAL(pszProj,"POLAR") )
    {
        SetPS( OSR_GDV( papszPrj, "PARAM_2", 0.0 ), 
               OSR_GDV( papszPrj, "PARAM_1", 0.0 ), 
               1.0,
               OSR_GDV( papszPrj, "PARAM_3", 0.0 ), 
               OSR_GDV( papszPrj, "PARAM_4", 0.0 ) );
    }

    else
    {
        CPLDebug( "OGR_ESRI", "Unsupported projection: %s", pszProj );
        SetLocalCS( pszProj );
    }

/* -------------------------------------------------------------------- */
/*      Try to translate the datum/spheroid.                            */
/* -------------------------------------------------------------------- */
    if( !IsLocal() && GetAttrNode( "GEOGCS" ) == NULL )
    {
        const char *pszDatum;

        pszDatum = OSR_GDS( papszPrj, "Datum", "");

        if( EQUAL(pszDatum,"NAD27") || EQUAL(pszDatum,"NAD83")
            || EQUAL(pszDatum,"WGS84") || EQUAL(pszDatum,"WGS72") )
        {
            SetWellKnownGeogCS( pszDatum );
        }
        else if( EQUAL( pszDatum, "EUR" ) )
        {
            SetWellKnownGeogCS( "EPSG:4230" );
        }
        else
        {
            const char *pszSpheroid;

            pszSpheroid = OSR_GDS( papszPrj, "Spheroid", "");
            
            if( EQUAL(pszSpheroid,"INT1909") 
                || EQUAL(pszSpheroid,"INTERNATIONAL1909") )
            {
                OGRSpatialReference oGCS;
                oGCS.importFromEPSG( 4022 );
                CopyGeogCSFrom( &oGCS );
            }
            else if( EQUAL(pszSpheroid,"AIRY") )
            {
                OGRSpatialReference oGCS;
                oGCS.importFromEPSG( 4001 );
                CopyGeogCSFrom( &oGCS );
            }
            else if( EQUAL(pszSpheroid,"CLARKE1866") )
            {
                OGRSpatialReference oGCS;
                oGCS.importFromEPSG( 4008 );
                CopyGeogCSFrom( &oGCS );
            }
            else if( EQUAL(pszSpheroid,"GRS80") )
            {
                OGRSpatialReference oGCS;
                oGCS.importFromEPSG( 4019 );
                CopyGeogCSFrom( &oGCS );
            }
            else if( EQUAL(pszSpheroid,"KRASOVSKY") 
                     || EQUAL(pszSpheroid,"KRASSOVSKY") )
            {
                OGRSpatialReference oGCS;
                oGCS.importFromEPSG( 4024 );
                CopyGeogCSFrom( &oGCS );
            }
            else
            {
                // If we don't know, default to WGS84 so there is something there.
                SetWellKnownGeogCS( "WGS84" );
            }
        }
    }

/* -------------------------------------------------------------------- */
/*      Linear units translation                                        */
/* -------------------------------------------------------------------- */
    if( IsLocal() || IsProjected() )
    {
        const char *pszValue;

        pszValue = OSR_GDS( papszPrj, "Units", NULL );
        if( pszValue == NULL )
            SetLinearUnitsAndUpdateParameters( SRS_UL_METER, 1.0 );
        else if( EQUAL(pszValue,"FEET") )
            SetLinearUnitsAndUpdateParameters( SRS_UL_FOOT, atof(SRS_UL_FOOT_CONV) );
        else
            SetLinearUnitsAndUpdateParameters( pszValue, 1.0 );
    }
    
    return OGRERR_NONE;
}

/************************************************************************/
/*                            morphToESRI()                             */
/************************************************************************/

/**
 * Convert in place from ESRI WKT format.
 *
 * The value notes of this coordinate system as modified in various manners
 * to adhere more closely to the WKT standard.  This mostly involves
 * translating a variety of ESRI names for projections, arguments and
 * datums to "standard" names, as defined by Adam Gawne-Cain's reference
 * translation of EPSG to WKT for the CT specification.
 *
 * This does the same as the C function OSRMorphToESRI().
 *
 * @return OGRERR_NONE unless something goes badly wrong.
 */

OGRErr OGRSpatialReference::morphToESRI()

{
    OGRErr      eErr;

/* -------------------------------------------------------------------- */
/*      Fixup ordering, missing linear units, etc.                      */
/* -------------------------------------------------------------------- */
    eErr = Fixup();
    if( eErr != OGRERR_NONE )
        return eErr;

/* -------------------------------------------------------------------- */
/*      Strip all CT parameters (AXIS, AUTHORITY, TOWGS84, etc).        */
/* -------------------------------------------------------------------- */
    eErr = StripCTParms();
    if( eErr != OGRERR_NONE )
        return eErr;

    if( GetRoot() == NULL )
        return OGRERR_NONE;

/* -------------------------------------------------------------------- */
/*      There is a special case for Hotine Oblique Mercator to split    */
/*      out the case with an angle to rectified grid.  Bug 423          */
/* -------------------------------------------------------------------- */
    const char *pszProjection = GetAttrValue("PROJECTION");
    
    if( pszProjection != NULL
        && EQUAL(pszProjection,SRS_PT_HOTINE_OBLIQUE_MERCATOR) 
        && fabs(GetProjParm(SRS_PP_AZIMUTH, 0.0 )-90) < 0.0001 
        && fabs(GetProjParm(SRS_PP_RECTIFIED_GRID_ANGLE, 0.0 )-90) < 0.0001 )
    {
        SetNode( "PROJCS|PROJECTION", 
                 "Hotine_Oblique_Mercator_Azimuth_Center" );

        /* ideally we should strip out of the rectified_grid_angle */
        pszProjection = GetAttrValue("PROJECTION");
    }

/* -------------------------------------------------------------------- */
/*      Polar_Stereographic maps to ESRI codes                          */
/*      Stereographic_South_Pole or Stereographic_North_Pole based      */
/*      on latitude.                                                    */
/* -------------------------------------------------------------------- */
    if( pszProjection != NULL
        && EQUAL(pszProjection,SRS_PT_POLAR_STEREOGRAPHIC) )
    {
        if( GetProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0 ) < 0.0 )
        {
            SetNode( "PROJCS|PROJECTION", 
                     "Stereographic_South_Pole" );
            pszProjection = GetAttrValue("PROJECTION");
        }
        else
        {
            SetNode( "PROJCS|PROJECTION", 
                     "Stereographic_North_Pole" );
            pszProjection = GetAttrValue("PROJECTION");
        }
    }

/* -------------------------------------------------------------------- */
/*      Translate PROJECTION keywords that are misnamed.                */
/* -------------------------------------------------------------------- */
    GetRoot()->applyRemapper( "PROJECTION", 
                              apszProjMapping+1, apszProjMapping, 2 );
    pszProjection = GetAttrValue("PROJECTION");

/* -------------------------------------------------------------------- */
/*      Translate DATUM keywords that are misnamed.                     */
/* -------------------------------------------------------------------- */
    InitDatumMappingTable();

    GetRoot()->applyRemapper( "DATUM", 
                              papszDatumMapping+2, papszDatumMapping+1, 3 );

/* -------------------------------------------------------------------- */
/*      Very specific handling for some well known geographic           */
/*      coordinate systems.                                             */
/* -------------------------------------------------------------------- */
    OGR_SRSNode *poGeogCS = GetAttrNode( "GEOGCS" );
    if( poGeogCS != NULL )
    {
        const char *pszGeogCSName = poGeogCS->GetChild(0)->GetValue();
        const char *pszAuthName = GetAuthorityName("GEOGCS");
        const char *pszUTMPrefix = NULL;