geo_normalize.c

GIS系统支持库Geospatial Data Abstraction Library代码.GDAL is a translator library for raster geospatial dat

/******************************************************************************
 * $Id: geo_normalize.c,v 1.42 2005/02/17 01:21:38 fwarmerdam Exp $
 *
 * Project:  libgeotiff
 * Purpose:  Code to normalize PCS and other composite codes in a GeoTIFF file.
 * Author:   Frank Warmerdam, [e-mail]
 *
 ******************************************************************************
 * Copyright (c) 1999, 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: geo_normalize.c,v $
 * Revision 1.42  2005/02/17 01:21:38  fwarmerdam
 * fixed handling of ProjFalseOrigin{Easting,Northing}GeoKey
 *
 * Revision 1.41  2004/12/01 22:06:42  fwarmerdam
 * bug 698: GTIFGetGCSInfo should not fail on missing pm if pm info not req.
 *
 * Revision 1.40  2004/07/09 17:27:37  warmerda
 * Added 9122 as an alias for simple degrees.
 *
 * Revision 1.39  2004/06/07 12:57:13  warmerda
 * fallback to using gdal_datum.csv if datum.csv not found
 *
 * Revision 1.38  2004/03/19 12:20:40  dron
 * Initialize projection parameters in GTIFFetchProjParms() before using.
 *
 * Revision 1.37  2003/07/08 17:31:30  warmerda
 * cleanup various warnings
 *
 * Revision 1.36  2003/01/28 18:31:58  warmerda
 * Default dfInDegrees in GTIFAngleToDD().
 *
 * Revision 1.35  2003/01/15 04:39:16  warmerda
 * Added GTIFDeaccessCSV
 *
 * Revision 1.34  2003/01/15 03:37:40  warmerda
 * added GTIFFreeMemory()
 *
 * Revision 1.33  2002/12/05 19:21:01  warmerda
 * fixed dfInDegrees to actually be in degrees, not radians!
 *
 * Revision 1.32  2002/11/30 16:01:11  warmerda
 * fixed some problems in GTIFGetUOMAngleInfo
 *
 * Revision 1.31  2002/11/30 15:44:35  warmerda
 * fixed GetCTParms EPSG code mappings
 *
 * Revision 1.30  2002/11/28 22:27:42  warmerda
 * preliminary upgrade to EPSG 6.2.2 tables
 *
 * Revision 1.29  2002/06/19 03:51:15  warmerda
 * migrated cpl_csv.h into cpl_serv.h
 *
 * Revision 1.28  2002/01/03 21:28:25  warmerda
 * call CSVDeaccess(NULL) at end of GTIFPrintDefn()
 *
 * Revision 1.27  2001/04/17 13:41:10  warmerda
 * fix memory leaks in GTIFPrintDefn()
 *
 * Revision 1.26  2001/04/17 13:23:07  warmerda
 * added support for reading custom ellipsoid definitions
 *
 * Revision 1.25  2001/03/05 04:55:26  warmerda
 * CVSDeaccess at end of GTIFGetDefn to avoid file leak
 *
 * Revision 1.24  2001/03/05 03:26:29  warmerda
 * fixed memory leaks in GTIFPrintDefn()
 *
 * Revision 1.23  2001/02/23 13:49:48  warmerda
 * Fixed GTIFPrintDefn() to use fprintf( fp ), instead of printf().
 *
 * Revision 1.22  2000/10/13 14:30:57  warmerda
 * fixed LCC parm order when parameters read directly from geotiff file
 *
 * Revision 1.21  2000/09/15 19:30:14  warmerda
 * report units of linear proj parms
 *
 * Revision 1.20  2000/09/15 18:21:07  warmerda
 * Fixed order of parameters for LCC 2SP.  When parameters
 * were read from EPSG CSV files the standard parallels and origin
 * were mixed up.  This affects alot of state plane zones!
 *
 * Revision 1.19  2000/06/09 14:05:43  warmerda
 * added default knowledge of NAD27/NAD83/WGS72/WGS84
 *
 * Revision 1.18  1999/12/10 21:28:12  warmerda
 * fixed Stereographic to look for ProjCenterLat/Long
 *
 * Revision 1.17  1999/12/10 20:06:58  warmerda
 * fixed up scale geokey used for a couple of projections
 *
 * Revision 1.16  1999/12/10 19:50:21  warmerda
 * Added EquidistantConic support, fixed return of StdParallel2GeoKey for
 * LCC2, and Albers.
 *
 * Revision 1.15  1999/12/10 19:39:26  warmerda
 * Fixed bug setting the false northing for files with
 * ProjCenterNorthingGeoKey set in GTIFGetDefn().
 *
 * Revision 1.14  1999/09/17 14:58:37  warmerda
 * Added ProjRectifiedGridAngleGeoKey(3096) and support for it's
 * use with Oblique Mercator in geo_normalize.c.
 *
 * Revision 1.13  1999/09/17 00:55:26  warmerda
 * added GTIFGetUOMAngleInfo(), and UOMAngle in GTIFDefn
 *
 * Revision 1.12  1999/09/15 18:51:31  warmerda
 * Map 9808 to TM South Oriented, not TM Modified Alaska.
 *
 * Revision 1.11  1999/09/15 16:44:06  warmerda
 * Change meter to metre to match EPSG database in GTIFGetUOMLengthInfo()
 * shortcut.
 *
 * Revision 1.10  1999/09/15 16:35:15  warmerda
 * Fixed the fractions of second handling properly in GTIFAngleStringToDD().
 *
 * Revision 1.9  1999/09/15 14:24:17  warmerda
 * Fixed serious bug in geo_normalize.c with translation of
 * DD.MMSSsss values.  Return value was seriously off if any
 * fraction of a second was included in the string.
 *
 * Revision 1.8  1999/07/13 03:12:52  warmerda
 * Make scale a parameter of CT_Stereographic.
 *
 * Revision 1.7  1999/05/04 03:13:22  warmerda
 * fixed a serious bug in parsing DMSmmss.sss values, and a bug in forming DMS strings
 *
 * Revision 1.6  1999/05/03 17:50:31  warmerda
 * avoid warnings on IRIX
 *
 * Revision 1.5  1999/04/28 20:04:51  warmerda
 * Added doxygen style documentation.
 * Use GTIFPCSToMapSys() and related functions to partially normalize
 * projections when we don't have the CSV files.
 *
 * Revision 1.4  1999/03/18 21:34:59  geotiff
 * added GTIFDecToDMS
 *
 * Revision 1.3  1999/03/17 19:53:15  geotiff
 * sys includes moved to cpl_serv.h
 *
 * Revision 1.2  1999/03/10 18:24:06  geotiff
 * corrected to use int'
 *
 * Revision 1.1  1999/03/09 15:57:04  geotiff
 * New
 *
 * Revision 1.4  1999/03/03 02:29:38  warmerda
 * Define PI if not already defined.
 *
 * Revision 1.3  1999/03/02 21:10:57  warmerda
 * added lots of projections
 *
 * Revision 1.2  1999/02/24 16:24:15  warmerda
 * Continuing to evolve
 *
 * Revision 1.1  1999/02/22 18:51:08  warmerda
 * New
 *
 */
 
#include "cpl_serv.h"
#include "geo_tiffp.h"
#include "geovalues.h"
#include "geo_normalize.h"

#ifndef KvUserDefined
#  define KvUserDefined 32767
#endif

#ifndef PI
#  define PI 3.14159265358979323846
#endif

/* EPSG Codes for projection parameters.  Unfortunately, these bear no
   relationship to the GeoTIFF codes even though the names are so similar. */

#define EPSGNatOriginLat         8801
#define EPSGNatOriginLong        8802
#define EPSGNatOriginScaleFactor 8805
#define EPSGFalseEasting         8806
#define EPSGFalseNorthing        8807
#define EPSGProjCenterLat        8811
#define EPSGProjCenterLong       8812
#define EPSGAzimuth              8813
#define EPSGAngleRectifiedToSkewedGrid 8814
#define EPSGInitialLineScaleFactor 8815
#define EPSGProjCenterEasting    8816
#define EPSGProjCenterNorthing   8817
#define EPSGPseudoStdParallelLat 8818
#define EPSGPseudoStdParallelScaleFactor 8819
#define EPSGFalseOriginLat       8821
#define EPSGFalseOriginLong      8822
#define EPSGStdParallel1Lat      8823
#define EPSGStdParallel2Lat      8824
#define EPSGFalseOriginEasting   8826
#define EPSGFalseOriginNorthing  8827
#define EPSGSphericalOriginLat   8828
#define EPSGSphericalOriginLong  8829
#define EPSGInitialLongitude     8830
#define EPSGZoneWidth            8831

/************************************************************************/
/*                           GTIFGetPCSInfo()                           */
/************************************************************************/

int GTIFGetPCSInfo( int nPCSCode, char **ppszEPSGName, 
                    short *pnProjOp, short *pnUOMLengthCode, 
                    short *pnGeogCS )

{
    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 = (short) atoi(pszValue);
        else
            *pnUOMLengthCode = KvUserDefined;
    }

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

        pszValue =
            CSLGetField( papszRecord,
                         CSVGetFileFieldId(pszFilename,"COORD_OP_CODE"));
        if( atoi(pszValue) > 0 )
            *pnProjOp = (short) atoi(pszValue);
        else
            *pnUOMLengthCode = KvUserDefined;
    }

/* -------------------------------------------------------------------- */
/*      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 = (short) atoi(pszValue);
        else
            *pnGeogCS = KvUserDefined;
    }

    return TRUE;
}

/************************************************************************/
/*                           GTIFAngleToDD()                            */
/*                                                                      */
/*      Convert a numeric angle to decimal degress.                     */
/************************************************************************/

double GTIFAngleToDD( double dfAngle, int nUOMAngle )

{
    if( nUOMAngle == 9110 )		/* DDD.MMSSsss */
    {
        char	szAngleString[32];

        sprintf( szAngleString, "%12.7f", dfAngle );
        dfAngle = GTIFAngleStringToDD( szAngleString, nUOMAngle );
    }
    else
    {
        double		dfInDegrees = 1.0;
        
        GTIFGetUOMAngleInfo( nUOMAngle, NULL, &dfInDegrees );
        dfAngle = dfAngle * dfInDegrees;
    }

    return( dfAngle );
}

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

double GTIFAngleStringToDD( const char * pszAngle, int nUOMAngle )

{
    double	dfAngle;