git.sesse.net Git - tcxmerge/blob - LatLong-UTMconversion.cpp

   1 //LatLong- UTM conversion.cpp
   2 //Lat Long - UTM, UTM - Lat Long conversions
   3
   4 #include <math.h>
   5 #include <stdio.h>
   6 #include <stdlib.h>
   7 #include "constants.h"
   8 #include "LatLong-UTMconversion.h"
   9
  10
  11 /*Reference ellipsoids derived from Peter H. Dana's website-
  12 http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.html
  13 Department of Geography, University of Texas at Austin
  14 Internet: pdana@mail.utexas.edu
  15 3/22/95
  16
  17 Source
  18 Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System
  19 1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency
  20 */
  21
  22
  23
  24 void LLtoUTM(int ReferenceEllipsoid, const double Lat, const double Long,
  25                          double &UTMNorthing, double &UTMEasting, char* UTMZone)
  26 {
  27 //converts lat/long to UTM coords.  Equations from USGS Bulletin 1532
  28 //East Longitudes are positive, West longitudes are negative.
  29 //North latitudes are positive, South latitudes are negative
  30 //Lat and Long are in decimal degrees
  31         //Written by Chuck Gantz- chuck.gantz@globalstar.com
  32
  33         double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius;
  34         double eccSquared = ellipsoid[ReferenceEllipsoid].eccentricitySquared;
  35         double k0 = 0.9996;
  36
  37         double LongOrigin;
  38         double eccPrimeSquared;
  39         double N, T, C, A, M;
  40
  41 //Make sure the longitude is between -180.00 .. 179.9
  42         double LongTemp = (Long+180)-int((Long+180)/360)*360-180; // -180.00 .. 179.9;
  43
  44         double LatRad = Lat*deg2rad;
  45         double LongRad = LongTemp*deg2rad;
  46         double LongOriginRad;
  47         int    ZoneNumber;
  48
  49         ZoneNumber = int((LongTemp + 180)/6) + 1;
  50
  51         if( Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0 )
  52                 ZoneNumber = 32;
  53
  54   // Special zones for Svalbard
  55         if( Lat >= 72.0 && Lat < 84.0 )
  56         {
  57           if(      LongTemp >= 0.0  && LongTemp <  9.0 ) ZoneNumber = 31;
  58           else if( LongTemp >= 9.0  && LongTemp < 21.0 ) ZoneNumber = 33;
  59           else if( LongTemp >= 21.0 && LongTemp < 33.0 ) ZoneNumber = 35;
  60           else if( LongTemp >= 33.0 && LongTemp < 42.0 ) ZoneNumber = 37;
  61          }
  62         LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
  63         LongOriginRad = LongOrigin * deg2rad;
  64
  65         //compute the UTM Zone from the latitude and longitude
  66         sprintf(UTMZone, "%d%c", ZoneNumber, UTMLetterDesignator(Lat));
  67
  68         eccPrimeSquared = (eccSquared)/(1-eccSquared);
  69
  70         N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad));
  71         T = tan(LatRad)*tan(LatRad);
  72         C = eccPrimeSquared*cos(LatRad)*cos(LatRad);
  73         A = cos(LatRad)*(LongRad-LongOriginRad);
  74
  75         M = a*((1       - eccSquared/4          - 3*eccSquared*eccSquared/64    - 5*eccSquared*eccSquared*eccSquared/256)*LatRad
  76                                 - (3*eccSquared/8       + 3*eccSquared*eccSquared/32    + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad)
  77                                                                         + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad)
  78                                                                         - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad));
  79
  80         UTMEasting = (double)(k0*N*(A+(1-T+C)*A*A*A/6
  81                                         + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120)
  82                                         + 500000.0);
  83
  84         UTMNorthing = (double)(k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24
  85                                  + (61-58*T+T*T+600*C-330*eccPrimeSquared)*A*A*A*A*A*A/720)));
  86         if(Lat < 0)
  87                 UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
  88 }
  89
  90 char UTMLetterDesignator(double Lat)
  91 {
  92 //This routine determines the correct UTM letter designator for the given latitude
  93 //returns 'Z' if latitude is outside the UTM limits of 84N to 80S
  94         //Written by Chuck Gantz- chuck.gantz@globalstar.com
  95         char LetterDesignator;
  96
  97         if((84 >= Lat) && (Lat >= 72)) LetterDesignator = 'X';
  98         else if((72 > Lat) && (Lat >= 64)) LetterDesignator = 'W';
  99         else if((64 > Lat) && (Lat >= 56)) LetterDesignator = 'V';
 100         else if((56 > Lat) && (Lat >= 48)) LetterDesignator = 'U';
 101         else if((48 > Lat) && (Lat >= 40)) LetterDesignator = 'T';
 102         else if((40 > Lat) && (Lat >= 32)) LetterDesignator = 'S';
 103         else if((32 > Lat) && (Lat >= 24)) LetterDesignator = 'R';
 104         else if((24 > Lat) && (Lat >= 16)) LetterDesignator = 'Q';
 105         else if((16 > Lat) && (Lat >= 8)) LetterDesignator = 'P';
 106         else if(( 8 > Lat) && (Lat >= 0)) LetterDesignator = 'N';
 107         else if(( 0 > Lat) && (Lat >= -8)) LetterDesignator = 'M';
 108         else if((-8> Lat) && (Lat >= -16)) LetterDesignator = 'L';
 109         else if((-16 > Lat) && (Lat >= -24)) LetterDesignator = 'K';
 110         else if((-24 > Lat) && (Lat >= -32)) LetterDesignator = 'J';
 111         else if((-32 > Lat) && (Lat >= -40)) LetterDesignator = 'H';
 112         else if((-40 > Lat) && (Lat >= -48)) LetterDesignator = 'G';
 113         else if((-48 > Lat) && (Lat >= -56)) LetterDesignator = 'F';
 114         else if((-56 > Lat) && (Lat >= -64)) LetterDesignator = 'E';
 115         else if((-64 > Lat) && (Lat >= -72)) LetterDesignator = 'D';
 116         else if((-72 > Lat) && (Lat >= -80)) LetterDesignator = 'C';
 117         else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits
 118
 119         return LetterDesignator;
 120 }
 121
 122
 123 void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone,
 124                           double& Lat,  double& Long )
 125 {
 126 //converts UTM coords to lat/long.  Equations from USGS Bulletin 1532
 127 //East Longitudes are positive, West longitudes are negative.
 128 //North latitudes are positive, South latitudes are negative
 129 //Lat and Long are in decimal degrees.
 130         //Written by Chuck Gantz- chuck.gantz@globalstar.com
 131
 132         double k0 = 0.9996;
 133         double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius;
 134         double eccSquared = ellipsoid[ReferenceEllipsoid].eccentricitySquared;
 135         double eccPrimeSquared;
 136         double e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared));
 137         double N1, T1, C1, R1, D, M;
 138         double LongOrigin;
 139         double mu, phi1, phi1Rad;
 140         double x, y;
 141         int ZoneNumber;
 142         char* ZoneLetter;
 143         int NorthernHemisphere; //1 for northern hemispher, 0 for southern
 144
 145         x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
 146         y = UTMNorthing;
 147
 148         ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10);
 149         if((*ZoneLetter - 'N') >= 0)
 150                 NorthernHemisphere = 1;//point is in northern hemisphere
 151         else
 152         {
 153                 NorthernHemisphere = 0;//point is in southern hemisphere
 154                 y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
 155         }
 156
 157         LongOrigin = (ZoneNumber - 1)*6 - 180 + 3;  //+3 puts origin in middle of zone
 158
 159         eccPrimeSquared = (eccSquared)/(1-eccSquared);
 160
 161         M = y / k0;
 162         mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256));
 163
 164         phi1Rad = mu    + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu)
 165                                 + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu)
 166                                 +(151*e1*e1*e1/96)*sin(6*mu);
 167         phi1 = phi1Rad*rad2deg;
 168
 169         N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad));
 170         T1 = tan(phi1Rad)*tan(phi1Rad);
 171         C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad);
 172         R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5);
 173         D = x/(N1*k0);
 174
 175         Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24
 176                                         +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720);
 177         Lat = Lat * rad2deg;
 178
 179         Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1)
 180                                         *D*D*D*D*D/120)/cos(phi1Rad);
 181         Long = LongOrigin + Long * rad2deg;
 182
 183 }