Jpp 19.3.0
the software that should make you happy
Loading...
Searching...
No Matches
LatLong-UTMconversion.hh File Reference
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

Go to the source code of this file.

Classes

class  Ellipsoid
 

Functions

char UTMLetterDesignator (double Lat)
 
void LLtoUTM (int ReferenceEllipsoid, const double Lat, const double Long, double &UTMNorthing, double &UTMEasting, char *UTMZone)
 
void UTMtoLL (int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char *UTMZone, double &Lat, double &Long)
 

Variables

static Ellipsoid ellipsoid []
 

Function Documentation

◆ UTMLetterDesignator()

char UTMLetterDesignator ( double Lat)
inline

Definition at line 74 of file LatLong-UTMconversion.hh.

74 {
75 //This routine determines the correct UTM letter designator for the given latitude
76 //returns 'Z' if latitude is outside the UTM limits of 84N to 80S
77 //Written by Chuck Gantz- chuck.gantz@globalstar.com
78 char LetterDesignator;
79
80 if ((84 >= Lat) && (Lat >= 72)) LetterDesignator = 'X';
81 else if ((72 > Lat) && (Lat >= 64)) LetterDesignator = 'W';
82 else if ((64 > Lat) && (Lat >= 56)) LetterDesignator = 'V';
83 else if ((56 > Lat) && (Lat >= 48)) LetterDesignator = 'U';
84 else if ((48 > Lat) && (Lat >= 40)) LetterDesignator = 'T';
85 else if ((40 > Lat) && (Lat >= 32)) LetterDesignator = 'S';
86 else if ((32 > Lat) && (Lat >= 24)) LetterDesignator = 'R';
87 else if ((24 > Lat) && (Lat >= 16)) LetterDesignator = 'Q';
88 else if ((16 > Lat) && (Lat >= 8)) LetterDesignator = 'P';
89 else if ((8 > Lat) && (Lat >= 0)) LetterDesignator = 'N';
90 else if ((0 > Lat) && (Lat >= -8)) LetterDesignator = 'M';
91 else if ((-8 > Lat) && (Lat >= -16)) LetterDesignator = 'L';
92 else if ((-16 > Lat) && (Lat >= -24)) LetterDesignator = 'K';
93 else if ((-24 > Lat) && (Lat >= -32)) LetterDesignator = 'J';
94 else if ((-32 > Lat) && (Lat >= -40)) LetterDesignator = 'H';
95 else if ((-40 > Lat) && (Lat >= -48)) LetterDesignator = 'G';
96 else if ((-48 > Lat) && (Lat >= -56)) LetterDesignator = 'F';
97 else if ((-56 > Lat) && (Lat >= -64)) LetterDesignator = 'E';
98 else if ((-64 > Lat) && (Lat >= -72)) LetterDesignator = 'D';
99 else if ((-72 > Lat) && (Lat >= -80)) LetterDesignator = 'C';
100 else LetterDesignator = 'Z'; //This is here as an error flag to show that the Latitude is outside the UTM limits
101
102 return LetterDesignator;
103}

◆ LLtoUTM()

void LLtoUTM ( int ReferenceEllipsoid,
const double Lat,
const double Long,
double & UTMNorthing,
double & UTMEasting,
char * UTMZone )
inline

Definition at line 107 of file LatLong-UTMconversion.hh.

108 {
109 //converts lat/long to UTM coords. Equations from USGS Bulletin 1532
110 //East Longitudes are positive, West longitudes are negative.
111 //North latitudes are positive, South latitudes are negative
112 //Lat and Long are in decimal degrees
113 //Written by Chuck Gantz- chuck.gantz@globalstar.com
114
115 const double PI = 3.14159265;
116 //const double FOURTHPI = PI / 4;
117 const double deg2rad = PI / 180;
118 //const double rad2deg = 180.0 / PI;
119
120 double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius;
121 double eccSquared = ellipsoid[ReferenceEllipsoid].eccentricitySquared;
122 double k0 = 0.9996;
123
124 double LongOrigin;
125 double eccPrimeSquared;
126 double N, T, C, A, M;
127
128 //Make sure the longitude is between -180.00 .. 179.9
129 double LongTemp = (Long + 180) - int((Long + 180) / 360)*360 - 180; // -180.00 .. 179.9;
130
131 double LatRad = Lat*deg2rad;
132 double LongRad = LongTemp*deg2rad;
133 double LongOriginRad;
134 int ZoneNumber;
135
136 ZoneNumber = int((LongTemp + 180) / 6) + 1;
137
138 if (Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0)
139 ZoneNumber = 32;
140
141 // Special zones for Svalbard
142 if (Lat >= 72.0 && Lat < 84.0) {
143 if (LongTemp >= 0.0 && LongTemp < 9.0) ZoneNumber = 31;
144 else if (LongTemp >= 9.0 && LongTemp < 21.0) ZoneNumber = 33;
145 else if (LongTemp >= 21.0 && LongTemp < 33.0) ZoneNumber = 35;
146 else if (LongTemp >= 33.0 && LongTemp < 42.0) ZoneNumber = 37;
147 }
148 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
149 LongOriginRad = LongOrigin * deg2rad;
150
151 //compute the UTM Zone from the latitude and longitude
152 sprintf(UTMZone, "%d%c", ZoneNumber, UTMLetterDesignator(Lat));
153
154 eccPrimeSquared = (eccSquared) / (1 - eccSquared);
155
156 N = a / sqrt(1 - eccSquared * sin(LatRad) * sin(LatRad));
157 T = tan(LatRad) * tan(LatRad);
158 C = eccPrimeSquared * cos(LatRad) * cos(LatRad);
159 A = cos(LatRad)*(LongRad - LongOriginRad);
160
161 M = a * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatRad
162 - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(2 * LatRad)
163 + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(4 * LatRad)
164 - (35 * eccSquared * eccSquared * eccSquared / 3072) * sin(6 * LatRad));
165
166 UTMEasting = (double) (k0 * N * (A + (1 - T + C) * A * A * A / 6
167 + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120)
168 + 500000.0);
169
170 UTMNorthing = (double) (k0 * (M + N * tan(LatRad)*(A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
171 + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A * A * A * A * A * A / 720)));
172 if (Lat < 0)
173 UTMNorthing += 10000000.0; //10000000 meter offset for southern hemisphere
174}
static Ellipsoid ellipsoid[]
char UTMLetterDesignator(double Lat)
const double a

◆ UTMtoLL()

void UTMtoLL ( int ReferenceEllipsoid,
const double UTMNorthing,
const double UTMEasting,
const char * UTMZone,
double & Lat,
double & Long )
inline

Definition at line 178 of file LatLong-UTMconversion.hh.

179 {
180 //converts UTM coords to lat/long. Equations from USGS Bulletin 1532
181 //East Longitudes are positive, West longitudes are negative.
182 //North latitudes are positive, South latitudes are negative
183 //Lat and Long are in decimal degrees.
184 //Written by Chuck Gantz- chuck.gantz@globalstar.com
185
186 const double PI = 3.14159265;
187 //const double FOURTHPI = PI / 4;
188 //const double deg2rad = PI / 180;
189 const double rad2deg = 180.0 / PI;
190
191 double k0 = 0.9996;
192 double a = ellipsoid[ReferenceEllipsoid].EquatorialRadius;
193 double eccSquared = ellipsoid[ReferenceEllipsoid].eccentricitySquared;
194 double eccPrimeSquared;
195 double e1 = (1 - sqrt(1 - eccSquared)) / (1 + sqrt(1 - eccSquared));
196 double N1, T1, C1, R1, D, M;
197 double LongOrigin;
198 double mu, phi1Rad;
199 //double phi1;
200 double x, y;
201 int ZoneNumber;
202 char* ZoneLetter;
203 // int NorthernHemisphere; //1 for northern hemispher, 0 for southern
204
205 x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
206 y = UTMNorthing;
207
208 ZoneNumber = strtoul(UTMZone, &ZoneLetter, 10);
209 /*
210 if((*ZoneLetter - 'N') >= 0)
211 NorthernHemisphere = 1;//point is in northern hemisphere
212 else
213 {
214 NorthernHemisphere = 0;//point is in southern hemisphere
215 y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
216 }
217 */
218 if ((*ZoneLetter - 'N') < 0) {
219 y -= 10000000.0; //remove 10,000,000 meter offset used for southern hemisphere
220 }
221
222 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; //+3 puts origin in middle of zone
223
224 eccPrimeSquared = (eccSquared) / (1 - eccSquared);
225
226 M = y / k0;
227 mu = M / (a * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256));
228
229 phi1Rad = mu + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * sin(2 * mu)
230 + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * sin(4 * mu)
231 +(151 * e1 * e1 * e1 / 96) * sin(6 * mu);
232 // phi1 = phi1Rad*rad2deg;
233
234 N1 = a / sqrt(1 - eccSquared * sin(phi1Rad) * sin(phi1Rad));
235 T1 = tan(phi1Rad) * tan(phi1Rad);
236 C1 = eccPrimeSquared * cos(phi1Rad) * cos(phi1Rad);
237 R1 = a * (1 - eccSquared) / pow(1 - eccSquared * sin(phi1Rad) * sin(phi1Rad), 1.5);
238 D = x / (N1 * k0);
239
240 Lat = phi1Rad - (N1 * tan(phi1Rad) / R1)*(D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * eccPrimeSquared) * D * D * D * D / 24
241 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * eccPrimeSquared - 3 * C1 * C1) * D * D * D * D * D * D / 720);
242 Lat = Lat * rad2deg;
243
244 Long = (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * eccPrimeSquared + 24 * T1 * T1)
245 * D * D * D * D * D / 120) / cos(phi1Rad);
246 Long = LongOrigin + Long * rad2deg;
247
248}
#define R1(x)
T pow(const T &x, const double y)
Power .
Definition JMath.hh:97

Variable Documentation

◆ ellipsoid

Ellipsoid ellipsoid[]
static
Initial value:
={
Ellipsoid(-1, "Placeholder", 0, 0),
Ellipsoid(1, "Airy", 6377563, 0.00667054),
Ellipsoid(2, "Australian National", 6378160, 0.006694542),
Ellipsoid(3, "Bessel 1841", 6377397, 0.006674372),
Ellipsoid(4, "Bessel 1841 (Nambia) ", 6377484, 0.006674372),
Ellipsoid(5, "Clarke 1866", 6378206, 0.006768658),
Ellipsoid(6, "Clarke 1880", 6378249, 0.006803511),
Ellipsoid(7, "Everest", 6377276, 0.006637847),
Ellipsoid(8, "Fischer 1960 (Mercury) ", 6378166, 0.006693422),
Ellipsoid(9, "Fischer 1968", 6378150, 0.006693422),
Ellipsoid(10, "GRS 1967", 6378160, 0.006694605),
Ellipsoid(11, "GRS 1980", 6378137, 0.00669438),
Ellipsoid(12, "Helmert 1906", 6378200, 0.006693422),
Ellipsoid(13, "Hough", 6378270, 0.00672267),
Ellipsoid(14, "International", 6378388, 0.00672267),
Ellipsoid(15, "Krassovsky", 6378245, 0.006693422),
Ellipsoid(16, "Modified Airy", 6377340, 0.00667054),
Ellipsoid(17, "Modified Everest", 6377304, 0.006637847),
Ellipsoid(18, "Modified Fischer 1960", 6378155, 0.006693422),
Ellipsoid(19, "South American 1969", 6378160, 0.006694542),
Ellipsoid(20, "WGS 60", 6378165, 0.006693422),
Ellipsoid(21, "WGS 66", 6378145, 0.006694542),
Ellipsoid(22, "WGS-72", 6378135, 0.006694318),
Ellipsoid(23, "WGS-84", 6378137, 0.00669438)
}

Definition at line 45 of file LatLong-UTMconversion.hh.

45 {// id, Ellipsoid name, Equatorial Radius, square of eccentricity
46 Ellipsoid(-1, "Placeholder", 0, 0), //placeholder only, To allow array indices to match id numbers
47 Ellipsoid(1, "Airy", 6377563, 0.00667054),
48 Ellipsoid(2, "Australian National", 6378160, 0.006694542),
49 Ellipsoid(3, "Bessel 1841", 6377397, 0.006674372),
50 Ellipsoid(4, "Bessel 1841 (Nambia) ", 6377484, 0.006674372),
51 Ellipsoid(5, "Clarke 1866", 6378206, 0.006768658),
52 Ellipsoid(6, "Clarke 1880", 6378249, 0.006803511),
53 Ellipsoid(7, "Everest", 6377276, 0.006637847),
54 Ellipsoid(8, "Fischer 1960 (Mercury) ", 6378166, 0.006693422),
55 Ellipsoid(9, "Fischer 1968", 6378150, 0.006693422),
56 Ellipsoid(10, "GRS 1967", 6378160, 0.006694605),
57 Ellipsoid(11, "GRS 1980", 6378137, 0.00669438),
58 Ellipsoid(12, "Helmert 1906", 6378200, 0.006693422),
59 Ellipsoid(13, "Hough", 6378270, 0.00672267),
60 Ellipsoid(14, "International", 6378388, 0.00672267),
61 Ellipsoid(15, "Krassovsky", 6378245, 0.006693422),
62 Ellipsoid(16, "Modified Airy", 6377340, 0.00667054),
63 Ellipsoid(17, "Modified Everest", 6377304, 0.006637847),
64 Ellipsoid(18, "Modified Fischer 1960", 6378155, 0.006693422),
65 Ellipsoid(19, "South American 1969", 6378160, 0.006694542),
66 Ellipsoid(20, "WGS 60", 6378165, 0.006693422),
67 Ellipsoid(21, "WGS 66", 6378145, 0.006694542),
68 Ellipsoid(22, "WGS-72", 6378135, 0.006694318),
69 Ellipsoid(23, "WGS-84", 6378137, 0.00669438)
70};