src/coords.c

   1 /*
   2 coords.c
   3 borrowed from:
   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.
  10
  11 */
  12 /* coords.h - include file for coords routines
  13 **
  14 ** Copyright © 2001 by Jef Poskanzer <jef@acme.com>.
  15 ** All rights reserved.
  16 **
  17 ** Redistribution and use in source and binary forms, with or without
  18 ** modification, are permitted provided that the following conditions
  19 ** are met:
  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.
  25 **
  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
  36 ** SUCH DAMAGE.
  37 */
  38 #ifdef HAVE_CONFIG_H
  39 #include "config.h"
  40 #endif
  41
  42 #ifdef HAVE_STDLIB_H
  43 #include <stdlib.h>
  44 #endif
  45 #ifdef HAVE_STRING_H
  46 #include <string.h>
  47 #endif
  48 #ifdef HAVE_MATH_H
  49 #include <math.h>
  50 #endif
  51
  52 #include "coords.h"
  53 #ifdef HAVE_VIKING
  54 #include "viking.h"
  55 #include "globals.h"
  56 #else
  57 #define DEG2RAD(x) ((x)*(M_PI/180))
  58 #define RAD2DEG(x) ((x)*(180/M_PI))
  59 #endif
  60 #include "degrees_converters.h"
  61 #include "misc/fpconv.h"
  62
  63 /**
  64  *
  65  */
  66 void a_coords_dtostr_buffer ( double d, char buffer[COORDS_STR_BUFFER_SIZE] )
  67 {
  68   int str_len = fpconv_dtoa(d, buffer);
  69   if ( str_len < COORDS_STR_BUFFER_SIZE )
  70     buffer[str_len] = '\0';
  71   else
  72     buffer[COORDS_STR_BUFFER_SIZE-1] = '\0';
  73 }
  74
  75 /**
  76  * Convert a double to a string WITHOUT LOCALE.
  77  *
  78  * Following GPX specifications, decimal values are xsd:decimal
  79  * So, they must use the period separator, not the localized one.
  80  *
  81  * The returned value must be freed by g_free.
  82  */
  83 char *a_coords_dtostr ( double d )
  84 {
  85   gchar *buffer = g_malloc(G_ASCII_DTOSTR_BUF_SIZE*sizeof(gchar));
  86   // Note that this doesn't necessary produce the shortest string
  87   //g_ascii_formatd (buffer, G_ASCII_DTOSTR_BUF_SIZE, "%.9lf", (gdouble)d);
  88   // Thus use third party code:
  89   a_coords_dtostr_buffer ( d, buffer );
  90   return buffer;
  91 }
  92
  93 #define PIOVER180 0.01745329252
  94
  95 #define K0 0.9996
  96
  97 /* WGS-84 */
  98 #define EquatorialRadius 6378137
  99 #define EccentricitySquared 0.00669438
 100
 101 static char coords_utm_letter( double latitude );
 102
 103 int a_coords_utm_equal( const struct UTM *utm1, const struct UTM *utm2 )
 104 {
 105   return ( utm1->easting == utm2->easting && utm1->northing == utm2->northing && utm1->zone == utm2->zone );
 106 }
 107
 108 double a_coords_utm_diff( const struct UTM *utm1, const struct UTM *utm2 )
 109 {
 110   static struct LatLon tmp1, tmp2;
 111   if ( utm1->zone == utm2->zone ) {
 112     return sqrt ( pow ( utm1->easting - utm2->easting, 2 ) + pow ( utm1->northing - utm2->northing, 2 ) );
 113   } else {
 114     a_coords_utm_to_latlon ( utm1, &tmp1 );
 115     a_coords_utm_to_latlon ( utm2, &tmp2 );
 116     return a_coords_latlon_diff ( &tmp1, &tmp2 );
 117   }
 118 }
 119
 120 double a_coords_latlon_diff ( const struct LatLon *ll1, const struct LatLon *ll2 )
 121 {
 122   static struct LatLon tmp1, tmp2;
 123   gdouble tmp3;
 124   tmp1.lat = ll1->lat * PIOVER180;
 125   tmp1.lon = ll1->lon * PIOVER180;
 126   tmp2.lat = ll2->lat * PIOVER180;
 127   tmp2.lon = ll2->lon * PIOVER180;
 128   tmp3 = EquatorialRadius * acos(sin(tmp1.lat)*sin(tmp2.lat)+cos(tmp1.lat)*cos(tmp2.lat)*cos(tmp1.lon-tmp2.lon));
 129   // For very small differences we can sometimes get NaN returned
 130   return isnan(tmp3)?0:tmp3;
 131 }
 132
 133 void a_coords_latlon_to_utm( const struct LatLon *latlon, struct UTM *utm )
 134     {
 135     double latitude;
 136     double longitude;
 137     double lat_rad, long_rad;
 138     double long_origin, long_origin_rad;
 139     double eccPrimeSquared;
 140     double N, T, C, A, M;
 141     int zone;
 142     double northing, easting;
 143
 144     longitude = latlon->lon;
 145     latitude = latlon->lat;
 146
 147     /* We want the longitude within -180..180. */
 148     if ( longitude < -180.0 )
 149         longitude += 360.0;
 150     if ( longitude > 180.0 )
 151         longitude -= 360.0;
 152
 153     /* Now convert. */
 154     lat_rad = DEG2RAD(latitude);
 155     long_rad = DEG2RAD(longitude);
 156     zone = (int) ( ( longitude + 180 ) / 6 ) + 1;
 157     if ( latitude >= 56.0 && latitude < 64.0 &&
 158          longitude >= 3.0 && longitude < 12.0 )
 159         zone = 32;
 160     /* Special zones for Svalbard. */
 161     if ( latitude >= 72.0 && latitude < 84.0 )
 162         {
 163         if      ( longitude >= 0.0  && longitude <  9.0 ) zone = 31;
 164         else if ( longitude >= 9.0  && longitude < 21.0 ) zone = 33;
 165         else if ( longitude >= 21.0 && longitude < 33.0 ) zone = 35;
 166         else if ( longitude >= 33.0 && longitude < 42.0 ) zone = 37;
 167         }
 168     long_origin = ( zone - 1 ) * 6 - 180 + 3;   /* +3 puts origin in middle of zone */
 169     long_origin_rad = DEG2RAD(long_origin);
 170     eccPrimeSquared = EccentricitySquared / ( 1.0 - EccentricitySquared );
 171     N = EquatorialRadius / sqrt( 1.0 - EccentricitySquared * sin( lat_rad ) * sin( lat_rad ) );
 172     T = tan( lat_rad ) * tan( lat_rad );
 173     C = eccPrimeSquared * cos( lat_rad ) * cos( lat_rad );
 174     A = cos( lat_rad ) * ( long_rad - long_origin_rad );
 175     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 ) );
 176     easting =
 177         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;
 178     northing =
 179         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 ) );
 180     if ( latitude < 0.0 )
 181         northing += 10000000.0;  /* 1e7 meter offset for southern hemisphere */
 182
 183     utm->northing = northing;
 184     utm->easting = easting;
 185     utm->zone = zone;
 186     utm->letter = coords_utm_letter( latitude );
 187
 188     /* All done. */
 189     }
 190
 191
 192 static char coords_utm_letter( double latitude )
 193     {
 194     /* This routine determines the correct UTM letter designator for the
 195     ** given latitude.  It returns 'Z' if the latitude is outside the UTM
 196     ** limits of 84N to 80S.
 197     */
 198     if ( latitude <= 84.0 && latitude >= 72.0 ) return 'X';
 199     else if ( latitude < 72.0 && latitude >= 64.0 ) return 'W';
 200     else if ( latitude < 64.0 && latitude >= 56.0 ) return 'V';
 201     else if ( latitude < 56.0 && latitude >= 48.0 ) return 'U';
 202     else if ( latitude < 48.0 && latitude >= 40.0 ) return 'T';
 203     else if ( latitude < 40.0 && latitude >= 32.0 ) return 'S';
 204     else if ( latitude < 32.0 && latitude >= 24.0 ) return 'R';
 205     else if ( latitude < 24.0 && latitude >= 16.0 ) return 'Q';
 206     else if ( latitude < 16.0 && latitude >= 8.0 ) return 'P';
 207     else if ( latitude <  8.0 && latitude >= 0.0 ) return 'N';
 208     else if ( latitude <  0.0 && latitude >= -8.0 ) return 'M';
 209     else if ( latitude < -8.0 && latitude >= -16.0 ) return 'L';
 210     else if ( latitude < -16.0 && latitude >= -24.0 ) return 'K';
 211     else if ( latitude < -24.0 && latitude >= -32.0 ) return 'J';
 212     else if ( latitude < -32.0 && latitude >= -40.0 ) return 'H';
 213     else if ( latitude < -40.0 && latitude >= -48.0 ) return 'G';
 214     else if ( latitude < -48.0 && latitude >= -56.0 ) return 'F';
 215     else if ( latitude < -56.0 && latitude >= -64.0 ) return 'E';
 216     else if ( latitude < -64.0 && latitude >= -72.0 ) return 'D';
 217     else if ( latitude < -72.0 && latitude >= -80.0 ) return 'C';
 218     else return 'Z';
 219     }
 220
 221
 222
 223 void a_coords_utm_to_latlon( const struct UTM *utm, struct LatLon *latlon )
 224     {
 225     double northing, easting;
 226     int zone;
 227     char letter[100];
 228     double x, y;
 229     double eccPrimeSquared;
 230     double e1;
 231     double N1, T1, C1, R1, D, M;
 232     double long_origin;
 233     double mu, phi1_rad;
 234     double latitude, longitude;
 235
 236     northing = utm->northing;
 237     easting = utm->easting;
 238     zone = utm->zone;
 239     letter[0] = utm->letter;
 240
 241     /* Now convert. */
 242     x = easting - 500000.0;     /* remove 500000 meter offset */
 243     y = northing;
 244     if ( ( *letter - 'N' ) < 0 ) {
 245       /* southern hemisphere */
 246       y -= 10000000.0;  /* remove 1e7 meter offset */
 247     }
 248
 249     long_origin = ( zone - 1 ) * 6 - 180 + 3;   /* +3 puts origin in middle of zone */
 250     eccPrimeSquared = EccentricitySquared / ( 1.0 - EccentricitySquared );
 251     e1 = ( 1.0 - sqrt( 1.0 - EccentricitySquared ) ) / ( 1.0 + sqrt( 1.0 - EccentricitySquared ) );
 252     M = y / K0;
 253     mu = M / ( EquatorialRadius * ( 1.0 - EccentricitySquared / 4 - 3 * EccentricitySquared * EccentricitySquared / 64 - 5 * EccentricitySquared * EccentricitySquared * EccentricitySquared / 256 ) );
 254     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 );
 255     N1 = EquatorialRadius / sqrt( 1.0 - EccentricitySquared * sin( phi1_rad ) * sin( phi1_rad ) );
 256     T1 = tan( phi1_rad ) * tan( phi1_rad );
 257     C1 = eccPrimeSquared * cos( phi1_rad ) * cos( phi1_rad );
 258     R1 = EquatorialRadius * ( 1.0 - EccentricitySquared ) / pow( 1.0 - EccentricitySquared * sin( phi1_rad ) * sin( phi1_rad ), 1.5 );
 259     D = x / ( N1 * K0 );
 260     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 );
 261     latitude = RAD2DEG(latitude);
 262     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 );
 263     longitude = long_origin + RAD2DEG(longitude);
 264
 265     /* Show results. */
 266
 267     latlon->lat = latitude;
 268     latlon->lon = longitude;
 269
 270     }
 271
 272 void a_coords_latlon_to_string ( const struct LatLon *latlon,
 273                                  gchar **lat,
 274                                  gchar **lon )
 275 {
 276   g_return_if_fail ( latlon != NULL );
 277 #ifdef HAVE_VIKING
 278   vik_degree_format_t format = a_vik_get_degree_format ();
 279
 280   switch (format) {
 281   case VIK_DEGREE_FORMAT_DDD:
 282     *lat = convert_lat_dec_to_ddd ( latlon->lat );
 283     *lon = convert_lon_dec_to_ddd ( latlon->lon );
 284     break;
 285   case VIK_DEGREE_FORMAT_DMM:
 286     *lat = convert_lat_dec_to_dmm ( latlon->lat );
 287     *lon = convert_lon_dec_to_dmm ( latlon->lon );
 288     break;
 289   case VIK_DEGREE_FORMAT_DMS:
 290     *lat = convert_lat_dec_to_dms ( latlon->lat );
 291     *lon = convert_lon_dec_to_dms ( latlon->lon );
 292     break;
 293   case VIK_DEGREE_FORMAT_RAW:
 294     *lat = g_strdup_printf ( "%.6f", latlon->lat );
 295     *lon = g_strdup_printf ( "%.6f", latlon->lon );
 296     break;
 297   default:
 298     g_critical("Houston, we've had a problem. format=%d", format);
 299   }
 300 #else
 301   *lat = convert_lat_dec_to_ddd ( latlon->lat );
 302   *lon = convert_lon_dec_to_ddd ( latlon->lon );
 303 #endif
 304 }