ogr_fromepsg.cpp

在linux环境下

/******************************************************************************
 * $Id: ogr_fromepsg.cpp,v 1.29 2003/06/23 14:48:42 warmerda Exp $
 *
 * Project:  OpenGIS Simple Features Reference Implementation
 * Purpose:  Generate an OGRSpatialReference object based on an EPSG
 *           PROJCS, or GEOGCS code.
 * Author:   Frank Warmerdam, [e-mail]
 *
 ******************************************************************************
 * Copyright (c) 2000, Frank Warmerdam
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ******************************************************************************
 *
 * $Log: ogr_fromepsg.cpp,v $
 * Revision 1.29  2003/06/23 14:48:42  warmerda
 * make OGREPSTDatumNameMassage() public
 *
 * Revision 1.28  2003/05/30 18:35:15  warmerda
 * fixed stateplane.csv error msg and clear() before SetLocalCS call
 *
 * Revision 1.27  2003/05/30 17:48:57  warmerda
 * Set GEOGCS|UNIT AUTHORITY information.
 *
 * Revision 1.26  2003/05/30 15:39:53  warmerda
 * Added override units capability for SetStatePlane()
 *
 * Revision 1.25  2003/05/20 19:03:24  warmerda
 * suggest GDAL_DATA in error message
 *
 * Revision 1.24  2003/05/20 18:43:31  warmerda
 * added error reporting
 *
 * Revision 1.23  2003/02/14 22:15:04  warmerda
 * expand tabs
 *
 * Revision 1.22  2003/02/14 16:29:02  warmerda
 * Added LAEA support
 *
 * Revision 1.21  2003/02/06 04:53:53  warmerda
 * fixed handling of empty fields to avoid emitting debug warning
 *
 * Revision 1.20  2003/02/05 17:14:11  warmerda
 * work around a bug in EPSG 6.2.2 with scale factor units for PCS 2935/2936
 *
 * Revision 1.19  2003/01/08 18:23:45  warmerda
 * ensure order is fixed up
 *
 * Revision 1.18  2002/12/16 17:06:32  warmerda
 * ensure that prime meridian offsets are applied to angular proj parms
 *
 * Revision 1.17  2002/12/14 22:59:14  warmerda
 * added Krovak in ESRI compatible way
 *
 * Revision 1.16  2002/12/13 06:35:45  warmerda
 * fixed rotatin sign in TOWGS84 from method 9607
 *
 * Revision 1.15  2002/12/02 00:07:50  warmerda
 * set angular units properly in GEOGCS
 *
 * Revision 1.14  2002/11/30 15:45:07  warmerda
 * be careful about how UOM codes are linear
 *
 * Revision 1.13  2002/11/29 20:49:10  warmerda
 * fixed GetTRFInfo filename corruption
 *
 * Revision 1.12  2002/11/28 15:59:58  warmerda
 * rewritten to use EPSG 6.2.2 tables
 *
 * Revision 1.11  2001/08/23 13:20:07  warmerda
 * fixed serious bug in translating LCC projections from EPSG to WKT
 *
 * Revision 1.10  2001/07/18 05:03:05  warmerda
 * added CPL_CVSID
 *
 * Revision 1.9  2001/03/15 03:19:23  warmerda
 * added polyconic (9818) support
 *
 * Revision 1.8  2001/03/14 21:33:46  warmerda
 * treaat 9717 (Lambert Near Conformal) like LCC1SP
 *
 * Revision 1.7  2001/01/19 21:10:46  warmerda
 * replaced tabs
 *
 * Revision 1.6  2000/12/07 03:46:04  danmo
 * Made papszDatumEquiv[] const
 *
 * Revision 1.5  2000/10/23 13:01:59  warmerda
 * Fixed OSRSetStatePlane
 *
 * Revision 1.4  2000/10/20 04:19:38  warmerda
 * added setstateplane
 *
 * Revision 1.3  2000/03/20 22:39:08  warmerda
 * Added C function.
 *
 * Revision 1.2  2000/03/20 14:59:05  warmerda
 * fixed semi major axis check
 *
 * Revision 1.1  2000/03/16 19:04:14  warmerda
 * New
 *
 */

#include "ogr_spatialref.h"
#include "cpl_csv.h"

CPL_CVSID("$Id: ogr_fromepsg.cpp,v 1.29 2003/06/23 14:48:42 warmerda Exp $");

#ifndef PI
#  define PI 3.14159265358979323846
#endif

static const char *papszDatumEquiv[] =
{
    "Militar_Geographische_Institut",
    "Militar_Geographische_Institute",
    "World_Geodetic_System_1984",
    "WGS_1984",
    "WGS_72_Transit_Broadcast_Ephemeris",
    "WGS_1972_Transit_Broadcast_Ephemeris",
    "World_Geodetic_System_1972",
    "WGS_1972",
    "European_Terrestrial_Reference_System_89",
    "European_Reference_System_1989",
    NULL
};

/************************************************************************/
/*                      OGREPSGDatumNameMassage()                       */
/*                                                                      */
/*      Massage an EPSG datum name into WMT format.  Also transform     */
/*      specific exception cases into WKT versions.                     */
/************************************************************************/

void OGREPSGDatumNameMassage( char ** ppszDatum )

{
    int         i, j;
    char        *pszDatum = *ppszDatum;

/* -------------------------------------------------------------------- */
/*      Translate non-alphanumeric values to underscores.               */
/* -------------------------------------------------------------------- */
    for( i = 0; pszDatum[i] != '\0'; i++ )
    {
        if( !(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z')
            && !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z')
            && !(pszDatum[i] >= '0' && pszDatum[i] <= '9') )
        {
            pszDatum[i] = '_';
        }
    }

/* -------------------------------------------------------------------- */
/*      Remove repeated and trailing underscores.                       */
/* -------------------------------------------------------------------- */
    for( i = 1, j = 0; pszDatum[i] != '\0'; i++ )
    {
        if( pszDatum[j] == '_' && pszDatum[i] == '_' )
            continue;

        pszDatum[++j] = pszDatum[i];
    }
    if( pszDatum[j] == '_' )
        pszDatum[j] = '\0';
    else
        pszDatum[j+1] = '\0';
    
/* -------------------------------------------------------------------- */
/*      Search for datum equivelences.  Specific massaged names get     */
/*      mapped to OpenGIS specified names.                              */
/* -------------------------------------------------------------------- */
    for( i = 0; papszDatumEquiv[i] != NULL; i += 2 )
    {
        if( EQUAL(*ppszDatum,papszDatumEquiv[i]) )
        {
            CPLFree( *ppszDatum );
            *ppszDatum = CPLStrdup( papszDatumEquiv[i+1] );
            break;
        }
    }
}

/************************************************************************/
/*                        EPSGAngleStringToDD()                         */
/*                                                                      */
/*      Convert an angle in the specified units to decimal degrees.     */
/************************************************************************/

static double
EPSGAngleStringToDD( const char * pszAngle, int nUOMAngle )

{
    double      dfAngle;
    
    if( nUOMAngle == 9110 )             /* DDD.MMSSsss */
    {
        char    *pszDecimal;
        
        dfAngle = ABS(atoi(pszAngle));
        pszDecimal = (char *) strchr(pszAngle,'.');
        if( pszDecimal != NULL && strlen(pszDecimal) > 1 )
        {
            char        szMinutes[3];
            char        szSeconds[64];

            szMinutes[0] = pszDecimal[1];
            if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' )
                szMinutes[1] = pszDecimal[2];
            else
                szMinutes[1] = '0';
            
            szMinutes[2] = '\0';
            dfAngle += atoi(szMinutes) / 60.0;

            if( strlen(pszDecimal) > 3 )
            {
                szSeconds[0] = pszDecimal[3];
                if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' )
                {
                    szSeconds[1] = pszDecimal[4];
                    szSeconds[2] = '.';
                    strcpy( szSeconds+3, pszDecimal + 5 );
                }
                else
                {
                    szSeconds[1] = '0';
                    szSeconds[2] = '\0';
                }
                dfAngle += atof(szSeconds) / 3600.0;
            }
        }

        if( pszAngle[0] == '-' )
            dfAngle *= -1;
    }
    else if( nUOMAngle == 9105 || nUOMAngle == 9106 )   /* grad */
    {
        dfAngle = 180 * (atof(pszAngle ) / 200);
    }
    else if( nUOMAngle == 9101 )                        /* radians */
    {
        dfAngle = 180 * (atof(pszAngle ) / PI);
    }
    else if( nUOMAngle == 9103 )                        /* arc-minute */
    {
        dfAngle = atof(pszAngle) / 60;
    }
    else if( nUOMAngle == 9104 )                        /* arc-second */
    {
        dfAngle = atof(pszAngle) / 3600;
    }
    else /* decimal degrees ... some cases missing but seeminly never used */
    {
        CPLAssert( nUOMAngle == 9102 || nUOMAngle == 0 );
        
        dfAngle = atof(pszAngle );
    }

    return( dfAngle );
}

/************************************************************************/
/*                        EPSGGetUOMAngleInfo()                         */
/************************************************************************/

int EPSGGetUOMAngleInfo( int nUOMAngleCode,
                         char **ppszUOMName,
                         double * pdfInDegrees )

{
    const char  *pszUOMName = NULL;
    double      dfInDegrees = 1.0;
    const char *pszFilename = CSVFilename( "unit_of_measure.csv" );
    char        szSearchKey[24];

    sprintf( szSearchKey, "%d", nUOMAngleCode );
    pszUOMName = CSVGetField( pszFilename,
                              "UOM_CODE", szSearchKey, CC_Integer,
                              "UNIT_OF_MEAS_NAME" );

/* -------------------------------------------------------------------- */
/*      If the file is found, read from there.  Note that FactorC is    */
/*      an empty field for any of the DMS style formats, and in this    */
/*      case we really want to return the default InDegrees value       */
/*      (1.0) from above.                                               */
/* -------------------------------------------------------------------- */
    if( pszUOMName != NULL )
    {
        double dfFactorB, dfFactorC;
        
        dfFactorB = 
            atof(CSVGetField( pszFilename,
                              "UOM_CODE", szSearchKey, CC_Integer,
                              "FACTOR_B" ));
        
        dfFactorC = 
            atof(CSVGetField( pszFilename,
                              "UOM_CODE", szSearchKey, CC_Integer,
                              "FACTOR_C" ));

        if( dfFactorC != 0.0 )
            dfInDegrees = (dfFactorB / dfFactorC) * (180.0 / PI);

        /* We do a special override of some of the DMS formats name */
        if( nUOMAngleCode == 9102 || nUOMAngleCode == 9107
            || nUOMAngleCode == 9108 || nUOMAngleCode == 9110 )
            pszUOMName = "degree";

        // For some reason, (FactorB) is not very precise in EPSG, use
        // a more exact form for grads.
        if( nUOMAngleCode == 9105 )
            dfInDegrees = 180.0 / 200.0;
    }

/* -------------------------------------------------------------------- */
/*      Otherwise handle a few well known units directly.               */
/* -------------------------------------------------------------------- */
    else
    {
        switch( nUOMAngleCode )
        {
          case 9101:
            pszUOMName = "radian";
            dfInDegrees = 180.0 / PI;
            break;
        
          case 9102:
          case 9107:
          case 9108:
          case 9110:
            pszUOMName = "degree";
            dfInDegrees = 1.0;
            break;

          case 9103:
            pszUOMName = "arc-minute";
            dfInDegrees = 1 / 60.0;
            break;

          case 9104:
            pszUOMName = "arc-second";
            dfInDegrees = 1 / 3600.0;
            break;
        
          case 9105:
            pszUOMName = "grad";
            dfInDegrees = 180.0 / 200.0;
            break;

          case 9106:
            pszUOMName = "gon";
            dfInDegrees = 180.0 / 200.0;
            break;
        
          case 9109:
            pszUOMName = "microradian";
            dfInDegrees = 180.0 / (3.14159265358979 * 1000000.0);
            break;

          default:
            return FALSE;
        }
    }

/* -------------------------------------------------------------------- */
/*      Return to caller.                                               */
/* -------------------------------------------------------------------- */
    if( ppszUOMName != NULL )
    {
        if( pszUOMName != NULL )
            *ppszUOMName = CPLStrdup( pszUOMName );
        else
            *ppszUOMName = NULL;
    }

    if( pdfInDegrees != NULL )
        *pdfInDegrees = dfInDegrees;