Project

General

Profile

1 2568 harris
/* utm2ll.c - convert Universal Transverse Mercator into latitude & longitude
2
**
3
** Partially based on code by Chuck Gantz <chuck.gantz@globalstar.com>.
4
**
5
** Copyright ? 2001 by Jef Poskanzer <jef@acme.com>.
6
** All rights reserved.
7
**
8
** Redistribution and use in source and binary forms, with or without
9
** modification, are permitted provided that the following conditions
10
** are met:
11
** 1. Redistributions of source code must retain the above copyright
12
**    notice, this list of conditions and the following disclaimer.
13
** 2. Redistributions in binary form must reproduce the above copyright
14
**    notice, this list of conditions and the following disclaimer in the
15
**    documentation and/or other materials provided with the distribution.
16
**
17
** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
18
** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19
** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20
** ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
21
** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22
** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
23
** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
24
** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
25
** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
26
** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
27
** SUCH DAMAGE.
28
*/
29
30
#include <stdio.h>
31
#include <stdlib.h>
32
#include <unistd.h>
33
#include <string.h>
34
#include <math.h>
35
36
#include "coords.h"
37
38
39
#define ABS(x) ((x)>=0?(x):(-x))
40
41
42
static char* argv0;
43
44
45
static void usage( void );
46
static void parse_error( char* argstr );
47
static void fprintll( FILE* fp, char* fmt, double lat, double lng );
48
49
50
int
51
main( int argc, char** argv )
52
    {
53
    char* format = "%ns %Ad  %ew %Od";
54
    int arglen, argn;
55
    char* argstr;
56
    double northing, easting;
57
    int zone;
58
    char letter[100];
59
    double x, y;
60
    double eccPrimeSquared;
61
    double e1;
62
    double N1, T1, C1, R1, D, M;
63
    double long_origin;
64
    double mu, phi1_rad;
65
    int northernHemisphere;	/* 1 for northern hemisphere, 0 for southern */
66
    double latitude, longitude;
67
68
    argv0 = argv[0];
69
    if ( argc >= 3 && strcmp( argv[1], "-f" ) == 0 )
70
	{
71
	format = argv[2];
72
	argv += 2;
73
	argc -= 2;
74
	}
75
    if ( argc < 2 )
76
	usage();
77
78
    /* Collect all args into a single string. */
79
    arglen = 0;
80
    for ( argn = 1; argn < argc; ++argn )
81
	arglen += strlen( argv[argn] + 1 );
82
    argstr = (char*) malloc( arglen );
83
    if ( argstr == (char*) 0 )
84
	{
85
	(void) fprintf( stderr, "%s: out of memory\n", argv0 );
86
	exit( 1 );
87
	}
88
    (void) strcpy( argstr, "" );
89
    for ( argn = 1; argn < argc; ++argn )
90
	{
91
	if ( argn > 1 )
92
	    (void) strcat( argstr, " " );
93
	(void) strcat( argstr, argv[argn] );
94
	}
95
96
    /* Trim leading and trailing blanks. */
97
    while ( argstr[0] == ' ' )
98
	(void) strcpy( argstr, &(argstr[1]) );
99
    while ( argstr[strlen(argstr)-1] == ' ' )
100
	argstr[strlen(argstr)-1] = '\0';
101
102
    /* Parse string into northing, easting, and zone. */
103
    if ( sscanf( argstr, "%lf %lf %d%[A-Z]", &northing, &easting, &zone, letter ) == 4 )
104
	;
105
    else if ( sscanf( argstr, "%lf %lf %d", &northing, &easting, &zone ) == 3 )
106
	(void) strcpy( letter, "S" );
107
    else
108
	parse_error( argstr );
109
110
    /* Now convert. */
111
    x = easting - 500000.0;	/* remove 500000 meter offset */
112
    y = northing;
113
    if ( ( *letter - 'N' ) >= 0 )
114
	northernHemisphere = 1;	/* northern hemisphere */
115
    else
116
	{
117
	northernHemisphere = 0;	/* southern hemisphere */
118
	y -= 10000000.0;	/* remove 1e7 meter offset */
119
	}
120
    long_origin = ( zone - 1 ) * 6 - 180 + 3;	/* +3 puts origin in middle of zone */
121
    eccPrimeSquared = EccentricitySquared / ( 1.0 - EccentricitySquared );
122
    e1 = ( 1.0 - sqrt( 1.0 - EccentricitySquared ) ) / ( 1.0 + sqrt( 1.0 - EccentricitySquared ) );
123
    M = y / K0;
124
    mu = M / ( EquatorialRadius * ( 1.0 - EccentricitySquared / 4 - 3 * EccentricitySquared * EccentricitySquared / 64 - 5 * EccentricitySquared * EccentricitySquared * EccentricitySquared / 256 ) );
125
    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 );
126
    N1 = EquatorialRadius / sqrt( 1.0 - EccentricitySquared * sin( phi1_rad ) * sin( phi1_rad ) );
127
    T1 = tan( phi1_rad ) * tan( phi1_rad );
128
    C1 = eccPrimeSquared * cos( phi1_rad ) * cos( phi1_rad );
129
    R1 = EquatorialRadius * ( 1.0 - EccentricitySquared ) / pow( 1.0 - EccentricitySquared * sin( phi1_rad ) * sin( phi1_rad ), 1.5 );
130
    D = x / ( N1 * K0 );
131
    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 );
132
    latitude = latitude * 180.0 / M_PI;
133
    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 );
134
    longitude = long_origin + longitude * 180.0 / M_PI;
135
136
    /* Show results. */
137
    fprintll( stdout, format, latitude, longitude );
138
    putchar( '\n' );
139
140
    /* All done. */
141
    exit( 0 );
142
    }
143
144
145
static void
146
usage( void )
147
    {
148
    (void) fprintf( stderr, "usage:  %s [-f format] UTM\n", argv0 );
149
    exit( 1 );
150
    }
151
152
153
static void
154
parse_error( char* argstr )
155
    {
156
    (void) fprintf(
157
	stderr, "%s: can't parse UTM coordinates '%s'\n", argv0, argstr );
158
    exit( 1 );
159
    }
160
161
162
static void
163
fprintll( FILE* fp, char* fmt, double lat, double lng )
164
    {
165
    char ns = lat >= 0.0 ? 'N' : 'S';
166
    char ew = lng >= 0.0 ? 'E' : 'W';
167
    double abs_lat = ABS( lat );
168
    double abs_lng = ABS( lng );
169
    int int_deg_lat = (int) lat;
170
    int int_deg_lng = (int) lng;
171
    int int_abs_deg_lat = ABS( int_deg_lat );
172
    int int_abs_deg_lng = ABS( int_deg_lng );
173
    double min_lat = ( abs_lat - int_abs_deg_lat ) * 60;
174
    double min_lng = ( abs_lng - int_abs_deg_lng ) * 60;
175
    int int_min_lat = (int) min_lat;
176
    int int_min_lng = (int) min_lng;
177
    double sec_lat = ( min_lat - int_min_lat ) * 60;
178
    double sec_lng = ( min_lng - int_min_lng ) * 60;
179
    int int_sec_lat = (int) sec_lat;
180
    int int_sec_lng = (int) sec_lng;
181
    char* cp;
182
    char c1, c2;
183
184
    for ( cp = fmt; *cp != '\0'; ++cp )
185
	{
186
	if ( *cp != '%' )
187
	    {
188
	    putc( *cp, fp );
189
	    continue;
190
	    }
191
	++cp;
192
	if ( *cp == '\0' )
193
	    {
194
	    putc( '%', fp );
195
	    break;
196
	    }
197
	if ( *cp == '%' )
198
	    {
199
	    putc( '%', fp );
200
	    continue;
201
	    }
202
	c1 = *cp;
203
	++cp;
204
	if ( *cp == '\0' )
205
	    {
206
	    putc( '%', fp );
207
	    putc( c1, fp );
208
	    break;
209
	    }
210
	c2 = *cp;
211
	if ( c1 == 'n' && c2 == 's' )
212
	    putc( ns, fp );
213
	else if ( c1 == 'e' && c2 == 'w' )
214
	    putc( ew, fp );
215
	else if ( c1 == 'a' && c2 == 'd' )
216
	    (void) fprintf( fp, "%0.8g", lat );
217
	else if ( c1 == 'o' && c2 == 'd' )
218
	    (void) fprintf( fp, "%0.8g", lng );
219
	else if ( c1 == 'A' && c2 == 'd' )
220
	    (void) fprintf( fp, "%0.8g", abs_lat );
221
	else if ( c1 == 'O' && c2 == 'd' )
222
	    (void) fprintf( fp, "%0.8g", abs_lng );
223
	else if ( c1 == 'a' && c2 == 'D' )
224
	    (void) fprintf( fp, "%d", int_deg_lat );
225
	else if ( c1 == 'o' && c2 == 'D' )
226
	    (void) fprintf( fp, "%d", int_deg_lng );
227
	else if ( c1 == 'A' && c2 == 'D' )
228
	    (void) fprintf( fp, "%d", int_abs_deg_lat );
229
	else if ( c1 == 'O' && c2 == 'D' )
230
	    (void) fprintf( fp, "%d", int_abs_deg_lng );
231
	else if ( c1 == 'a' && c2 == 'm' )
232
	    (void) fprintf( fp, "%0.6g", min_lat );
233
	else if ( c1 == 'o' && c2 == 'm' )
234
	    (void) fprintf( fp, "%0.6g", min_lng );
235
	else if ( c1 == 'a' && c2 == 'M' )
236
	    (void) fprintf( fp, "%d", int_min_lat );
237
	else if ( c1 == 'o' && c2 == 'M' )
238
	    (void) fprintf( fp, "%d", int_min_lng );
239
	else if ( c1 == 'a' && c2 == 's' )
240
	    (void) fprintf( fp, "%0.4g", sec_lat );
241
	else if ( c1 == 'o' && c2 == 's' )
242
	    (void) fprintf( fp, "%0.4g", sec_lng );
243
	else if ( c1 == 'a' && c2 == 'S' )
244
	    (void) fprintf( fp, "%d", int_sec_lat );
245
	else if ( c1 == 'o' && c2 == 'S' )
246
	    (void) fprintf( fp, "%d", int_sec_lng );
247
	else
248
	    {
249
	    putc( '%', fp );
250
	    putc( c1, fp );
251
	    putc( c2, fp );
252
	    }
253
	}
254
    }