geo_normalize.c

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

    
    if( nUOMAngle == 9110 )		/* DDD.MMSSsss */
    {
        char	*pszDecimal;
        
        dfAngle = ABS(atoi(pszAngle));
        pszDecimal = 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 == KvUserDefined
                   || nUOMAngle == 0 );
        
        dfAngle = atof(pszAngle );
    }

    return( dfAngle );
}

/************************************************************************/
/*                           GTIFGetGCSInfo()                           */
/*                                                                      */
/*      Fetch the datum, and prime meridian related to a particular     */
/*      GCS.                                                            */
/************************************************************************/

int GTIFGetGCSInfo( int nGCSCode, char ** ppszName,
                    short * pnDatum, short * pnPM, short *pnUOMAngle )

{
    char	szSearchKey[24];
    int		nDatum, nPM, nUOMAngle;

/* -------------------------------------------------------------------- */
/*      Search the database for the corresponding datum code.           */
/* -------------------------------------------------------------------- */
    sprintf( szSearchKey, "%d", nGCSCode );

    nDatum = atoi(CSVGetField( CSVFilename("gcs.csv" ),
                               "COORD_REF_SYS_CODE", szSearchKey, CC_Integer,
                               "DATUM_CODE" ) );

/* -------------------------------------------------------------------- */
/*      Handle some "well known" GCS codes directly if the table        */
/*      wasn't found.                                                   */
/* -------------------------------------------------------------------- */
    if( nDatum < 1 )
    {
        const char * pszName = NULL;
        nPM = PM_Greenwich;
        nUOMAngle = Angular_DMS_Hemisphere; 
        if( nGCSCode == GCS_NAD27 )
        {
            nDatum = Datum_North_American_Datum_1927;
            pszName = "NAD27";
        }
        else if( nGCSCode == GCS_NAD83 )
        {
            nDatum = Datum_North_American_Datum_1983;
            pszName = "NAD83";
        }
        else if( nGCSCode == GCS_WGS_84 )
        {
            nDatum = Datum_WGS84;
            pszName = "WGS 84";
        }
        else if( nGCSCode == GCS_WGS_72 )
        {
            nDatum = Datum_WGS72;
            pszName = "WGS 82";
        }
        else
            return FALSE;

        if( ppszName != NULL )
            *ppszName = CPLStrdup( pszName );
        if( pnDatum != NULL )
            *pnDatum = (short) nDatum;
        if( pnPM != NULL )
            *pnPM = (short) nPM;
        if( pnUOMAngle != NULL )
            *pnUOMAngle = (short) nUOMAngle;

        return TRUE;
    }

    if( pnDatum != NULL )
        *pnDatum = (short) nDatum;
    
/* -------------------------------------------------------------------- */
/*      Get the PM.                                                     */
/* -------------------------------------------------------------------- */
    if( pnPM != NULL )
    {
        nPM = atoi(CSVGetField( CSVFilename("gcs.csv" ),
                            "COORD_REF_SYS_CODE", szSearchKey, CC_Integer,
                            "PRIME_MERIDIAN_CODE" ) );

        if( nPM < 1 )
            return FALSE;

        *pnPM = (short) nPM;
    }

/* -------------------------------------------------------------------- */
/*      Get the angular units.                                          */
/* -------------------------------------------------------------------- */
    nUOMAngle = atoi(CSVGetField( CSVFilename("gcs.csv" ),
                                  "COORD_REF_SYS_CODE",szSearchKey, CC_Integer,
                                  "UOM_CODE" ) );

    if( nUOMAngle < 1 )
        return FALSE;

    if( pnUOMAngle != NULL )
        *pnUOMAngle = (short) nUOMAngle;

/* -------------------------------------------------------------------- */
/*      Get the name, if requested.                                     */
/* -------------------------------------------------------------------- */
    if( ppszName != NULL )
        *ppszName =
            CPLStrdup(CSVGetField( CSVFilename("gcs.csv" ),
                                   "COORD_REF_SYS_CODE",szSearchKey,CC_Integer,
                                   "COORD_REF_SYS_NAME" ));
    
    return( TRUE );
}

/************************************************************************/
/*                        GTIFGetEllipsoidInfo()                        */
/*                                                                      */
/*      Fetch info about an ellipsoid.  Axes are always returned in     */
/*      meters.  SemiMajor computed based on inverse flattening         */
/*      where that is provided.                                         */
/************************************************************************/

int GTIFGetEllipsoidInfo( int nEllipseCode, char ** ppszName,
                          double * pdfSemiMajor, double * pdfSemiMinor )

{
    char	szSearchKey[24];
    double	dfSemiMajor, dfToMeters = 1.0;
    int		nUOMLength;
    
/* -------------------------------------------------------------------- */
/*      Get the semi major axis.                                        */
/* -------------------------------------------------------------------- */
    sprintf( szSearchKey, "%d", nEllipseCode );

    dfSemiMajor =
        atof(CSVGetField( CSVFilename("ellipsoid.csv" ),
                          "ELLIPSOID_CODE", szSearchKey, CC_Integer,
                          "SEMI_MAJOR_AXIS" ) );

/* -------------------------------------------------------------------- */
/*      Try some well known ellipsoids.                                 */
/* -------------------------------------------------------------------- */
    if( dfSemiMajor == 0.0 )
    {
        double     dfInvFlattening, dfSemiMinor;
        const char *pszName = NULL;
        
        if( nEllipseCode == Ellipse_Clarke_1866 )
        {
            pszName = "Clarke 1866";
            dfSemiMajor = 6378206.4;
            dfSemiMinor = 6356583.8;
            dfInvFlattening = 0.0;
        }
        else if( nEllipseCode == Ellipse_GRS_1980 )
        {
            pszName = "GRS 1980";
            dfSemiMajor = 6378137.0;
            dfSemiMinor = 0.0;
            dfInvFlattening = 298.257222101;
        }
        else if( nEllipseCode == Ellipse_WGS_84 )
        {
            pszName = "WGS 84";
            dfSemiMajor = 6378137.0;
            dfSemiMinor = 0.0;
            dfInvFlattening = 298.257223563;
        }
        else if( nEllipseCode == 7043 )
        {
            pszName = "WGS 72";
            dfSemiMajor = 6378135.0;
            dfSemiMinor = 0.0;
            dfInvFlattening = 298.26;
        }
        else
            return FALSE;

        if( dfSemiMinor == 0.0 )
            dfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening);

        if( pdfSemiMinor != NULL )
            *pdfSemiMinor = dfSemiMinor;
        if( pdfSemiMajor != NULL )
            *pdfSemiMajor = dfSemiMajor;
        if( ppszName != NULL )
            *ppszName = CPLStrdup( pszName );

        return TRUE;
    }

/* -------------------------------------------------------------------- */
/*	Get the translation factor into meters.				*/
/* -------------------------------------------------------------------- */
    nUOMLength = atoi(CSVGetField( CSVFilename("ellipsoid.csv" ),
                                   "ELLIPSOID_CODE", szSearchKey, CC_Integer,
                                   "UOM_CODE" ));
    GTIFGetUOMLengthInfo( nUOMLength, NULL, &dfToMeters );

    dfSemiMajor *= dfToMeters;
    
    if( pdfSemiMajor != NULL )
        *pdfSemiMajor = dfSemiMajor;
    
/* -------------------------------------------------------------------- */
/*      Get the semi-minor if requested.  If the Semi-minor axis        */
/*      isn't available, compute it based on the inverse flattening.    */
/* -------------------------------------------------------------------- */
    if( pdfSemiMinor != NULL )
    {
        *pdfSemiMinor =
            atof(CSVGetField( CSVFilename("ellipsoid.csv" ),
                              "ELLIPSOID_CODE", szSearchKey, CC_Integer,
                              "SEMI_MINOR_AXIS" )) * dfToMeters;

        if( *pdfSemiMinor == 0.0 )
        {
            double	dfInvFlattening;
            
            dfInvFlattening = 
                atof(CSVGetField( CSVFilename("ellipsoid.csv" ),
                                  "ELLIPSOID_CODE", szSearchKey, CC_Integer,
                                  "INV_FLATTENING" ));
            *pdfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening);
        }
    }

/* -------------------------------------------------------------------- */
/*      Get the name, if requested.                                     */
/* -------------------------------------------------------------------- */
    if( ppszName != NULL )
        *ppszName =
            CPLStrdup(CSVGetField( CSVFilename("ellipsoid.csv" ),
                                   "ELLIPSOID_CODE", szSearchKey, CC_Integer,
                                   "ELLIPSOID_NAME" ));
    
    return( TRUE );
}

/************************************************************************/
/*                           GTIFGetPMInfo()                            */
/*                                                                      */
/*      Get the offset between a given prime meridian and Greenwich     */
/*      in degrees.                                                     */
/************************************************************************/

int GTIFGetPMInfo( int nPMCode, char ** ppszName, double *pdfOffset )

{
    char	szSearchKey[24];
    int		nUOMAngle;
    const char *pszFilename = CSVFilename("prime_meridian.csv");

/* -------------------------------------------------------------------- */
/*      Use a special short cut for Greenwich, since it is so common.   */
/* -------------------------------------------------------------------- */
    if( nPMCode == PM_Greenwich )
    {
        if( pdfOffset != NULL )
            *pdfOffset = 0.0;
        if( ppszName != NULL )
            *ppszName = CPLStrdup( "Greenwich" );
        return TRUE;
    }

/* -------------------------------------------------------------------- */
/*      Search the database for the corresponding datum code.           */
/* -------------------------------------------------------------------- */
    sprintf( szSearchKey, "%d", nPMCode );

    nUOMAngle =
        atoi(CSVGetField( pszFilename, 
                          "PRIME_MERIDIAN_CODE", szSearchKey, CC_Integer,
                          "UOM_CODE" ) );
    if( nUOMAngle < 1 )
        return FALSE;

/* -------------------------------------------------------------------- */