4 http://acme.com/software/coords/
5 I (Evan Battaglia <viking@greentorch.org>) have only made some small changes such as
6 renaming functions and defining LatLon and UTM structs.
7 2004-02-10 -- I also added a function of my own -- a_coords_utm_diff() -- that I felt belonged in coords.c
8 2004-02-21 -- I also added a_coords_utm_equal().
9 2005-11-23 -- Added a_coords_dtostr() for lack of a better place.
12 /* coords.h - include file for coords routines
14 ** Copyright © 2001 by Jef Poskanzer <jef@acme.com>.
15 ** All rights reserved.
17 ** Redistribution and use in source and binary forms, with or without
18 ** modification, are permitted provided that the following conditions
20 ** 1. Redistributions of source code must retain the above copyright
21 ** notice, this list of conditions and the following disclaimer.
22 ** 2. Redistributions in binary form must reproduce the above copyright
23 ** notice, this list of conditions and the following disclaimer in the
24 ** documentation and/or other materials provided with the distribution.
26 ** THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
27 ** ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 ** IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
29 ** ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
30 ** FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
31 ** DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
32 ** OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
33 ** HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
34 ** LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
35 ** OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
46 #define M_PI 3.14159265358979
50 * Convert a double to a string WITHOUT LOCALE.
52 * Following GPX specifications, decimal values are xsd:decimal
53 * So, they must use the period separator, not the localized one.
55 * The returned value must be freed by g_free.
57 char *a_coords_dtostr ( double d
)
59 /* In order to ignore locale, we do all the stuff manually */
60 double integer
, decimal
;
63 /* 6 decimals are sufficient (~0,1m) */
64 /* Cf. http://www.tbs-sct.gc.ca/rpm-gbi/guides/Latlong_f.asp */
65 decimal
= d
- integer
;
66 decimal
= decimal
* 1000000;
67 decimal
= trunc ( decimal
);
68 decimal
= fabs ( decimal
);
71 return g_strdup_printf ( "%g.%06g", integer
, decimal
);
74 #define PIOVER180 0.01745329252
79 #define EquatorialRadius 6378137
80 #define EccentricitySquared 0.00669438
82 static char coords_utm_letter( double latitude
);
84 int a_coords_utm_equal( const struct UTM
*utm1
, const struct UTM
*utm2
)
86 return ( utm1
->easting
== utm2
->easting
&& utm1
->northing
== utm2
->northing
&& utm1
->zone
== utm2
->zone
);
89 double a_coords_utm_diff( const struct UTM
*utm1
, const struct UTM
*utm2
)
91 static struct LatLon tmp1
, tmp2
;
92 if ( utm1
->zone
== utm2
->zone
) {
93 return sqrt ( pow ( utm1
->easting
- utm2
->easting
, 2 ) + pow ( utm1
->northing
- utm2
->northing
, 2 ) );
95 a_coords_utm_to_latlon ( utm1
, &tmp1
);
96 a_coords_utm_to_latlon ( utm2
, &tmp2
);
97 return a_coords_latlon_diff ( &tmp1
, &tmp2
);
101 double a_coords_latlon_diff ( const struct LatLon
*ll1
, const struct LatLon
*ll2
)
103 static struct LatLon tmp1
, tmp2
;
104 tmp1
.lat
= ll1
->lat
* PIOVER180
;
105 tmp1
.lon
= ll1
->lon
* PIOVER180
;
106 tmp2
.lat
= ll2
->lat
* PIOVER180
;
107 tmp2
.lon
= ll2
->lon
* PIOVER180
;
108 return EquatorialRadius
* acos(sin(tmp1
.lat
)*sin(tmp2
.lat
)+cos(tmp1
.lat
)*cos(tmp2
.lat
)*cos(tmp1
.lon
-tmp2
.lon
));
111 void a_coords_latlon_to_utm( const struct LatLon
*latlon
, struct UTM
*utm
)
115 double lat_rad
, long_rad
;
116 double long_origin
, long_origin_rad
;
117 double eccPrimeSquared
;
118 double N
, T
, C
, A
, M
;
120 double northing
, easting
;
122 longitude
= latlon
->lon
;
123 latitude
= latlon
->lat
;
125 /* We want the longitude within -180..180. */
126 if ( longitude
< -180.0 )
128 if ( longitude
> 180.0 )
132 lat_rad
= latitude
* M_PI
/ 180.0;
133 long_rad
= longitude
* M_PI
/ 180.0;
134 zone
= (int) ( ( longitude
+ 180 ) / 6 ) + 1;
135 if ( latitude
>= 56.0 && latitude
< 64.0 &&
136 longitude
>= 3.0 && longitude
< 12.0 )
138 /* Special zones for Svalbard. */
139 if ( latitude
>= 72.0 && latitude
< 84.0 )
141 if ( longitude
>= 0.0 && longitude
< 9.0 ) zone
= 31;
142 else if ( longitude
>= 9.0 && longitude
< 21.0 ) zone
= 33;
143 else if ( longitude
>= 21.0 && longitude
< 33.0 ) zone
= 35;
144 else if ( longitude
>= 33.0 && longitude
< 42.0 ) zone
= 37;
146 long_origin
= ( zone
- 1 ) * 6 - 180 + 3; /* +3 puts origin in middle of zone */
147 long_origin_rad
= long_origin
* M_PI
/ 180.0;
148 eccPrimeSquared
= EccentricitySquared
/ ( 1.0 - EccentricitySquared
);
149 N
= EquatorialRadius
/ sqrt( 1.0 - EccentricitySquared
* sin( lat_rad
) * sin( lat_rad
) );
150 T
= tan( lat_rad
) * tan( lat_rad
);
151 C
= eccPrimeSquared
* cos( lat_rad
) * cos( lat_rad
);
152 A
= cos( lat_rad
) * ( long_rad
- long_origin_rad
);
153 M
= EquatorialRadius
* ( ( 1.0 - EccentricitySquared
/ 4 - 3 * EccentricitySquared
* EccentricitySquared
/ 64 - 5 * EccentricitySquared
* EccentricitySquared
* EccentricitySquared
/ 256 ) * lat_rad
- ( 3 * EccentricitySquared
/ 8 + 3 * EccentricitySquared
* EccentricitySquared
/ 32 + 45 * EccentricitySquared
* EccentricitySquared
* EccentricitySquared
/ 1024 ) * sin( 2 * lat_rad
) + ( 15 * EccentricitySquared
* EccentricitySquared
/ 256 + 45 * EccentricitySquared
* EccentricitySquared
* EccentricitySquared
/ 1024 ) * sin( 4 * lat_rad
) - ( 35 * EccentricitySquared
* EccentricitySquared
* EccentricitySquared
/ 3072 ) * sin( 6 * lat_rad
) );
155 K0
* N
* ( A
+ ( 1 - T
+ C
) * A
* A
* A
/ 6 + ( 5 - 18 * T
+ T
* T
+ 72 * C
- 58 * eccPrimeSquared
) * A
* A
* A
* A
* A
/ 120 ) + 500000.0;
157 K0
* ( M
+ N
* tan( lat_rad
) * ( A
* A
/ 2 + ( 5 - T
+ 9 * C
+ 4 * C
* C
) * A
* A
* A
* A
/ 24 + ( 61 - 58 * T
+ T
* T
+ 600 * C
- 330 * eccPrimeSquared
) * A
* A
* A
* A
* A
* A
/ 720 ) );
158 if ( latitude
< 0.0 )
159 northing
+= 10000000.0; /* 1e7 meter offset for southern hemisphere */
161 utm
->northing
= northing
;
162 utm
->easting
= easting
;
164 utm
->letter
= coords_utm_letter( latitude
);
170 static char coords_utm_letter( double latitude
)
172 /* This routine determines the correct UTM letter designator for the
173 ** given latitude. It returns 'Z' if the latitude is outside the UTM
174 ** limits of 84N to 80S.
176 if ( latitude
<= 84.0 && latitude
>= 72.0 ) return 'X';
177 else if ( latitude
< 72.0 && latitude
>= 64.0 ) return 'W';
178 else if ( latitude
< 64.0 && latitude
>= 56.0 ) return 'V';
179 else if ( latitude
< 56.0 && latitude
>= 48.0 ) return 'U';
180 else if ( latitude
< 48.0 && latitude
>= 40.0 ) return 'T';
181 else if ( latitude
< 40.0 && latitude
>= 32.0 ) return 'S';
182 else if ( latitude
< 32.0 && latitude
>= 24.0 ) return 'R';
183 else if ( latitude
< 24.0 && latitude
>= 16.0 ) return 'Q';
184 else if ( latitude
< 16.0 && latitude
>= 8.0 ) return 'P';
185 else if ( latitude
< 8.0 && latitude
>= 0.0 ) return 'N';
186 else if ( latitude
< 0.0 && latitude
>= -8.0 ) return 'M';
187 else if ( latitude
< -8.0 && latitude
>= -16.0 ) return 'L';
188 else if ( latitude
< -16.0 && latitude
>= -24.0 ) return 'K';
189 else if ( latitude
< -24.0 && latitude
>= -32.0 ) return 'J';
190 else if ( latitude
< -32.0 && latitude
>= -40.0 ) return 'H';
191 else if ( latitude
< -40.0 && latitude
>= -48.0 ) return 'G';
192 else if ( latitude
< -48.0 && latitude
>= -56.0 ) return 'F';
193 else if ( latitude
< -56.0 && latitude
>= -64.0 ) return 'E';
194 else if ( latitude
< -64.0 && latitude
>= -72.0 ) return 'D';
195 else if ( latitude
< -72.0 && latitude
>= -80.0 ) return 'C';
201 void a_coords_utm_to_latlon( const struct UTM
*utm
, struct LatLon
*latlon
)
203 double northing
, easting
;
207 double eccPrimeSquared
;
209 double N1
, T1
, C1
, R1
, D
, M
;
212 int northernHemisphere
; /* 1 for northern hemisphere, 0 for southern */
213 double latitude
, longitude
;
215 northing
= utm
->northing
;
216 easting
= utm
->easting
;
218 letter
[0] = utm
->letter
;
221 x
= easting
- 500000.0; /* remove 500000 meter offset */
223 if ( ( *letter
- 'N' ) >= 0 )
224 northernHemisphere
= 1; /* northern hemisphere */
227 northernHemisphere
= 0; /* southern hemisphere */
228 y
-= 10000000.0; /* remove 1e7 meter offset */
230 long_origin
= ( zone
- 1 ) * 6 - 180 + 3; /* +3 puts origin in middle of zone */
231 eccPrimeSquared
= EccentricitySquared
/ ( 1.0 - EccentricitySquared
);
232 e1
= ( 1.0 - sqrt( 1.0 - EccentricitySquared
) ) / ( 1.0 + sqrt( 1.0 - EccentricitySquared
) );
234 mu
= M
/ ( EquatorialRadius
* ( 1.0 - EccentricitySquared
/ 4 - 3 * EccentricitySquared
* EccentricitySquared
/ 64 - 5 * EccentricitySquared
* EccentricitySquared
* EccentricitySquared
/ 256 ) );
235 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
);
236 N1
= EquatorialRadius
/ sqrt( 1.0 - EccentricitySquared
* sin( phi1_rad
) * sin( phi1_rad
) );
237 T1
= tan( phi1_rad
) * tan( phi1_rad
);
238 C1
= eccPrimeSquared
* cos( phi1_rad
) * cos( phi1_rad
);
239 R1
= EquatorialRadius
* ( 1.0 - EccentricitySquared
) / pow( 1.0 - EccentricitySquared
* sin( phi1_rad
) * sin( phi1_rad
), 1.5 );
241 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 );
242 latitude
= latitude
* 180.0 / M_PI
;
243 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
);
244 longitude
= long_origin
+ longitude
* 180.0 / M_PI
;
248 latlon
->lat
= latitude
;
249 latlon
->lon
= longitude
;