Metacat

/* utm2ll.c - convert Universal Transverse Mercator into latitude & longitude

**

** Partially based on code by Chuck Gantz <chuck.gantz@globalstar.com>.

**

** Copyright ? 2001 by Jef Poskanzer <jef@acme.com>.

** All rights reserved.

**

** Redistribution and use in source and binary forms, with or without

** modification, are permitted provided that the following conditions

** are met:

** 1. Redistributions of source code must retain the above copyright

**    notice, this list of conditions and the following disclaimer.

** 2. Redistributions in binary form must reproduce the above copyright

**    notice, this list of conditions and the following disclaimer in the

**    documentation and/or other materials provided with the distribution.

**

** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND

** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE

** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE

** ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE

** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL

** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS

** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY

** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF

** SUCH DAMAGE.

*/

#include <stdio.h>

#include <stdlib.h>

#include <unistd.h>

#include <string.h>

#include <math.h>

#include "coords.h"

#define ABS(x) ((x)>=0?(x):(-x))

static char* argv0;

static void usage( void );

static void parse_error( char* argstr );

static void fprintll( FILE* fp, char* fmt, double lat, double lng );

int

main( int argc, char** argv )

    {

    char* format = "%ns %Ad  %ew %Od";

    int arglen, argn;

    char* argstr;

    double northing, easting;

    int zone;

    char letter[100];

    double x, y;

    double eccPrimeSquared;

    double e1;

    double N1, T1, C1, R1, D, M;

    double long_origin;

    double mu, phi1_rad;

    int northernHemisphere;	/* 1 for northern hemisphere, 0 for southern */

    double latitude, longitude;

    argv0 = argv[0];

    if ( argc >= 3 && strcmp( argv[1], "-f" ) == 0 )

	{

	format = argv[2];

	argv += 2;

	argc -= 2;

	}

    if ( argc < 2 )

	usage();

    /* Collect all args into a single string. */

    arglen = 0;

    for ( argn = 1; argn < argc; ++argn )

	arglen += strlen( argv[argn] + 1 );

    argstr = (char*) malloc( arglen );

    if ( argstr == (char*) 0 )

	{

	(void) fprintf( stderr, "%s: out of memory\n", argv0 );

	exit( 1 );

	}

    (void) strcpy( argstr, "" );

    for ( argn = 1; argn < argc; ++argn )

	{

	if ( argn > 1 )

	    (void) strcat( argstr, " " );

	(void) strcat( argstr, argv[argn] );

	}

    /* Trim leading and trailing blanks. */

    while ( argstr[0] == ' ' )

	(void) strcpy( argstr, &(argstr[1]) );

    while ( argstr[strlen(argstr)-1] == ' ' )

	argstr[strlen(argstr)-1] = '\0';

    /* Parse string into northing, easting, and zone. */

    if ( sscanf( argstr, "%lf %lf %d%[A-Z]", &northing, &easting, &zone, letter ) == 4 )

	;

    else if ( sscanf( argstr, "%lf %lf %d", &northing, &easting, &zone ) == 3 )

	(void) strcpy( letter, "S" );

    else

	parse_error( argstr );

    /* Now convert. */

    x = easting - 500000.0;	/* remove 500000 meter offset */

    y = northing;

    if ( ( *letter - 'N' ) >= 0 )

	northernHemisphere = 1;	/* northern hemisphere */

    else

	{

	northernHemisphere = 0;	/* southern hemisphere */

	y -= 10000000.0;	/* remove 1e7 meter offset */

	}

    long_origin = ( zone - 1 ) * 6 - 180 + 3;	/* +3 puts origin in middle of zone */

    eccPrimeSquared = EccentricitySquared / ( 1.0 - EccentricitySquared );

    e1 = ( 1.0 - sqrt( 1.0 - EccentricitySquared ) ) / ( 1.0 + sqrt( 1.0 - EccentricitySquared ) );

    M = y / K0;

    mu = M / ( EquatorialRadius * ( 1.0 - EccentricitySquared / 4 - 3 * EccentricitySquared * EccentricitySquared / 64 - 5 * EccentricitySquared * EccentricitySquared * EccentricitySquared / 256 ) );

    phi1_rad = mu + ( 3 * e1 / 2 - 27 * e1 * e1 * e1 / 32 )* sin( 2 * mu ) + ( 21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32 ) * sin( 4 * mu ) + ( 151 * e1 * e1 * e1 / 96 ) * sin( 6 *mu );

    N1 = EquatorialRadius / sqrt( 1.0 - EccentricitySquared * sin( phi1_rad ) * sin( phi1_rad ) );

    T1 = tan( phi1_rad ) * tan( phi1_rad );

    C1 = eccPrimeSquared * cos( phi1_rad ) * cos( phi1_rad );

    R1 = EquatorialRadius * ( 1.0 - EccentricitySquared ) / pow( 1.0 - EccentricitySquared * sin( phi1_rad ) * sin( phi1_rad ), 1.5 );

    D = x / ( N1 * K0 );

    latitude = phi1_rad - ( N1 * tan( phi1_rad ) / R1 ) * ( D * D / 2 -( 5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * eccPrimeSquared ) * D * D * D * D / 24 + ( 61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * eccPrimeSquared - 3 * C1 * C1 ) * D * D * D * D * D * D / 720 );

    latitude = latitude * 180.0 / M_PI;

    longitude = ( D - ( 1 + 2 * T1 + C1 ) * D * D * D / 6 + ( 5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * eccPrimeSquared + 24 * T1 * T1 ) * D * D * D * D * D / 120 ) / cos( phi1_rad );

    longitude = long_origin + longitude * 180.0 / M_PI;

    /* Show results. */

    fprintll( stdout, format, latitude, longitude );

    putchar( '\n' );

    /* All done. */

    exit( 0 );

    }

static void

usage( void )

    {

    (void) fprintf( stderr, "usage:  %s [-f format] UTM\n", argv0 );

    exit( 1 );

    }

static void

parse_error( char* argstr )

    {

    (void) fprintf(

	stderr, "%s: can't parse UTM coordinates '%s'\n", argv0, argstr );

    exit( 1 );

    }

static void

fprintll( FILE* fp, char* fmt, double lat, double lng )

    {

    char ns = lat >= 0.0 ? 'N' : 'S';

    char ew = lng >= 0.0 ? 'E' : 'W';

    double abs_lat = ABS( lat );

    double abs_lng = ABS( lng );

    int int_deg_lat = (int) lat;

    int int_deg_lng = (int) lng;

    int int_abs_deg_lat = ABS( int_deg_lat );

    int int_abs_deg_lng = ABS( int_deg_lng );

    double min_lat = ( abs_lat - int_abs_deg_lat ) * 60;

    double min_lng = ( abs_lng - int_abs_deg_lng ) * 60;

    int int_min_lat = (int) min_lat;

    int int_min_lng = (int) min_lng;

    double sec_lat = ( min_lat - int_min_lat ) * 60;

    double sec_lng = ( min_lng - int_min_lng ) * 60;

    int int_sec_lat = (int) sec_lat;

    int int_sec_lng = (int) sec_lng;

    char* cp;

    char c1, c2;

    for ( cp = fmt; *cp != '\0'; ++cp )

	{

	if ( *cp != '%' )

	    {

	    putc( *cp, fp );

	    continue;

	    }

	++cp;

	if ( *cp == '\0' )

	    {

	    putc( '%', fp );

	    break;

	    }

	if ( *cp == '%' )

	    {

	    putc( '%', fp );

	    continue;

	    }

	c1 = *cp;

	++cp;

	if ( *cp == '\0' )

	    {

	    putc( '%', fp );

	    putc( c1, fp );

	    break;

	    }

	c2 = *cp;

	if ( c1 == 'n' && c2 == 's' )

	    putc( ns, fp );

	else if ( c1 == 'e' && c2 == 'w' )

	    putc( ew, fp );

	else if ( c1 == 'a' && c2 == 'd' )

	    (void) fprintf( fp, "%0.8g", lat );

	else if ( c1 == 'o' && c2 == 'd' )

	    (void) fprintf( fp, "%0.8g", lng );

	else if ( c1 == 'A' && c2 == 'd' )

	    (void) fprintf( fp, "%0.8g", abs_lat );

	else if ( c1 == 'O' && c2 == 'd' )

	    (void) fprintf( fp, "%0.8g", abs_lng );

	else if ( c1 == 'a' && c2 == 'D' )

	    (void) fprintf( fp, "%d", int_deg_lat );

	else if ( c1 == 'o' && c2 == 'D' )

	    (void) fprintf( fp, "%d", int_deg_lng );

	else if ( c1 == 'A' && c2 == 'D' )

	    (void) fprintf( fp, "%d", int_abs_deg_lat );

	else if ( c1 == 'O' && c2 == 'D' )

	    (void) fprintf( fp, "%d", int_abs_deg_lng );

	else if ( c1 == 'a' && c2 == 'm' )

	    (void) fprintf( fp, "%0.6g", min_lat );

	else if ( c1 == 'o' && c2 == 'm' )

	    (void) fprintf( fp, "%0.6g", min_lng );

	else if ( c1 == 'a' && c2 == 'M' )

	    (void) fprintf( fp, "%d", int_min_lat );

	else if ( c1 == 'o' && c2 == 'M' )

	    (void) fprintf( fp, "%d", int_min_lng );

	else if ( c1 == 'a' && c2 == 's' )

	    (void) fprintf( fp, "%0.4g", sec_lat );

	else if ( c1 == 'o' && c2 == 's' )

	    (void) fprintf( fp, "%0.4g", sec_lng );

	else if ( c1 == 'a' && c2 == 'S' )

	    (void) fprintf( fp, "%d", int_sec_lat );

	else if ( c1 == 'o' && c2 == 'S' )

	    (void) fprintf( fp, "%d", int_sec_lng );

	else

	    {

	    putc( '%', fp );

	    putc( c1, fp );

	    putc( c2, fp );

	    }

	}

    }