Rev 1212 Rev 1269
Line 22... Line 22...
22 22
23 this version retains the core algorithms as static functions 23 this version retains the core algorithms as static functions
24   24  
25 */ 25 */
26   26  
27 #include "stdafx.h" 27 //#include "stdafx.h"
28 #include "Common.h" 28 //#include "Common.h"
29 #include "osgb.h" 29 //#include "osgb.h"
30 #include "utm.h" 30 #include "utm.h"
31 #include <math.h> 31 #include <math.h>
32   32  
33 // Some constants used by these functions. 33 // Some constants used by these functions.
34 static const double fe = 500000.0; 34 static const double fe = 500000.0;
Line 115... Line 115...
115 // and download MADTRAN.ZIP 115 // and download MADTRAN.ZIP
116 // - If the UTM zone is out of range, the y-zone character is set to the 116 // - If the UTM zone is out of range, the y-zone character is set to the
117 // asterisk character ('*'). 117 // asterisk character ('*').
118 //======================================================================= 118 //=======================================================================
119   119  
120 void LatLonToUtm (double a, double f, int& utmXZone, char& utmYZone, 120 void LatLonToUtm (double a, double f, int* utmXZone, char* utmYZone,
121 double& easting, double& northing, double lat, double lon) 121 double* easting, double* northing, double lat, double lon)
122 { 122 {
123 double recf; 123 double recf;
124 double b; 124 double b;
125 double eSquared; 125 double eSquared;
126 double e2Squared; 126 double e2Squared;
Line 140... Line 140...
140 double tmd; 140 double tmd;
141 double t1, t2, t3, t6, t7; 141 double t1, t2, t3, t6, t7;
142 double nfn; 142 double nfn;
143   143  
144 if (lon <= 0.0) { 144 if (lon <= 0.0) {
145 utmXZone = 30 + (int)(lon / 6.0); 145 *utmXZone = 30 + (int)(lon / 6.0);
146 } else { 146 } else {
147 utmXZone = 31 + (int)(lon / 6.0); 147 *utmXZone = 31 + (int)(lon / 6.0);
148 } 148 }
149 if (lat < 84.0 && lat >= 72.0) { 149 if (lat < 84.0 && lat >= 72.0) {
150 // Special case: zone X is 12 degrees from north to south, not 8. 150 // Special case: zone X is 12 degrees from north to south, not 8.
151 utmYZone = cArray[19]; 151 *utmYZone = cArray[19];
152 } else { 152 } else {
153 utmYZone = cArray[(int)((lat + 80.0) / 8.0)]; 153 *utmYZone = cArray[(int)((lat + 80.0) / 8.0)];
154 } 154 }
155 if (lat >= 84.0 || lat < -80.0) { 155 if (lat >= 84.0 || lat < -80.0) {
156 // Invalid coordinate; the vertical zone is set to the invalid 156 // Invalid coordinate; the vertical zone is set to the invalid
157 // character. 157 // character.
158 utmYZone = '*'; 158 *utmYZone = '*';
159 } 159 }
160   160  
-   161 double latRad = lat * (M_PI/180);
-   162 double lonRad = lon * (M_PI/180);
161 double latRad = lat * deg2rad; 163 // double latRad = lat * deg2rad;
162 double lonRad = lon * deg2rad; 164 // double lonRad = lon * deg2rad;
163 recf = 1.0 / f; 165 recf = 1.0 / f;
164 b = a * (recf - 1.0) / recf; 166 b = a * (recf - 1.0) / recf;
165 eSquared = CalculateESquared (a, b); 167 eSquared = CalculateESquared (a, b);
166 e2Squared = CalculateE2Squared (a, b); 168 e2Squared = CalculateE2Squared (a, b);
167 tn = (a - b) / (a + b); 169 tn = (a - b) / (a + b);
Line 173... Line 175...
173 cp = 15.0 * a * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn) 175 cp = 15.0 * a * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn)
174 - (tn * tn * tn * tn * tn)) / 4.0) / 16.0; 176 - (tn * tn * tn * tn * tn)) / 4.0) / 16.0;
175 dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0 177 dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0
176 * (tn * tn * tn * tn * tn) / 16.0) / 48.0; 178 * (tn * tn * tn * tn * tn) / 16.0) / 48.0;
177 ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0; 179 ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0;
178 olam = (utmXZone * 6 - 183) * deg2rad; 180 olam = (*utmXZone * 6 - 183) * (M_PI/180);
179 dlam = lonRad - olam; 181 dlam = lonRad - olam;
180 s = sin (latRad); 182 s = sin (latRad);
181 c = cos (latRad); 183 c = cos (latRad);
182 t = s / c; 184 t = s / c;
183 eta = e2Squared * (c * c); 185 eta = e2Squared * (c * c);
Line 186... Line 188...
186 t1 = tmd * ok; 188 t1 = tmd * ok;
187 t2 = sn * s * c * ok / 2.0; 189 t2 = sn * s * c * ok / 2.0;
188 t3 = sn * s * (c * c * c) * ok * (5.0 - (t * t) + 9.0 * eta + 4.0 190 t3 = sn * s * (c * c * c) * ok * (5.0 - (t * t) + 9.0 * eta + 4.0
189 * (eta * eta)) / 24.0; 191 * (eta * eta)) / 24.0;
190 if (latRad < 0.0) nfn = 10000000.0; else nfn = 0; 192 if (latRad < 0.0) nfn = 10000000.0; else nfn = 0;
191 northing = nfn + t1 + (dlam * dlam) * t2 + (dlam * dlam * dlam 193 *northing = nfn + t1 + (dlam * dlam) * t2 + (dlam * dlam * dlam
192 * dlam) * t3 + (dlam * dlam * dlam * dlam * dlam * dlam) + 0.5; 194 * dlam) * t3 + (dlam * dlam * dlam * dlam * dlam * dlam) + 0.5;
193 t6 = sn * c * ok; 195 t6 = sn * c * ok;
194 t7 = sn * (c * c * c) * (1.0 - (t * t) + eta) / 6.0; 196 t7 = sn * (c * c * c) * (1.0 - (t * t) + eta) / 6.0;
195 easting = fe + dlam * t6 + (dlam * dlam * dlam) * t7 + 0.5; 197 *easting = fe + dlam * t6 + (dlam * dlam * dlam) * t7 + 0.5;
196 if (northing >= 9999999.0) northing = 9999999.0; 198 if (*northing >= 9999999.0) *northing = 9999999.0;
197 } 199 }
198   200  
199 //======================================================================= 201 //=======================================================================
200 // Purpose: 202 // Purpose:
201 // This function converts the specified lat/lon coordinate to a UTM 203 // This function converts the specified lat/lon coordinate to a UTM
202 // coordinate in the WGS84 datum. (See the comment block for the 204 // coordinate in the WGS84 datum. (See the comment block for the
203 // LatLonToUtm() member function.) 205 // LatLonToUtm() member function.)
204 //======================================================================= 206 //=======================================================================
205   207  
206 void LatLonToUtmWGS84 (int& utmXZone, char& utmYZone, 208 void LatLonToUtmWGS84 (int* utmXZone, char* utmYZone,
207 double& easting, double& northing, double lat, double lon) 209 double* easting, double* northing, double lat, double lon)
208 { 210 {
209 LatLonToUtm (6378137.0, 1 / 298.257223563, utmXZone, utmYZone, 211 LatLonToUtm (6378137.0, 1 / 298.257223563, utmXZone, utmYZone,
210 easting, northing, lat, lon); 212 easting, northing, lat, lon);
211 } 213 }
212   -  
213   -  
214   -  
215   -  
216 //======================================================================= -  
217 // Purpose: -  
218 // This function converts the specified UTM coordinate to a lat/lon -  
219 // coordinate. -  
220 // Pre: -  
221 // - utmXZone must be between 1 and 60, inclusive. -  
222 // - utmYZone must be one of: CDEFGHJKLMNPQRSTUVWX -  
223 // Parameters: -  
224 // double a: -  
225 // Ellipsoid semi-major axis, in meters. (For WGS84 datum, use 6378137.0) -  
226 // double f: -  
227 // Ellipsoid flattening. (For WGS84 datum, use 1 / 298.257223563) -  
228 // int utmXZone: -  
229 // The horizontal zone number of the UTM coordinate. -  
230 // char utmYZone: -  
231 // The vertical zone letter of the UTM coordinate. -  
232 // double easting, double northing: -  
233 // The UTM coordinate to convert. -  
234 // double& lat: -  
235 // Upon exit, lat contains the latitude. -  
236 // double& lon: -  
237 // Upon exit, lon contains the longitude. -  
238 // Notes: -  
239 // The code in this function is a C conversion of some of the source code -  
240 // from the Mapping Datum Transformation Software (MADTRAN) program, written -  
241 // in PowerBasic. To get the source code for MADTRAN, go to: -  
242 // -  
243 // http://164.214.2.59/publications/guides/MADTRAN/index.html -  
244 // -  
245 // and download MADTRAN.ZIP -  
246 //======================================================================= -  
247   -  
248 void UtmToLatLon (double a, double f, int utmXZone, char utmYZone, -  
249 double easting, double northing, double& lat, double& lon) -  
250 { -  
251 double recf; -  
252 double b; -  
253 double eSquared; -  
254 double e2Squared; -  
255 double tn; -  
256 double ap; -  
257 double bp; -  
258 double cp; -  
259 double dp; -  
260 double ep; -  
261 double nfn; -  
262 double tmd; -  
263 double sr; -  
264 double sn; -  
265 double ftphi; -  
266 double s; -  
267 double c; -  
268 double t; -  
269 double eta; -  
270 double de; -  
271 double dlam; -  
272 double olam; -  
273   -  
274 recf = 1.0 / f; -  
275 b = a * (recf - 1) / recf; -  
276 eSquared = CalculateESquared (a, b); -  
277 e2Squared = CalculateE2Squared (a, b); -  
278 tn = (a - b) / (a + b); -  
279 ap = a * (1.0 - tn + 5.0 * ((tn * tn) - (tn * tn * tn)) / 4.0 + 81.0 * -  
280 ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 64.0); -  
281 bp = 3.0 * a * (tn - (tn * tn) + 7.0 * ((tn * tn * tn) -  
282 - (tn * tn * tn * tn)) / 8.0 + 55.0 * (tn * tn * tn * tn * tn) / 64.0) -  
283 / 2.0; -  
284 cp = 15.0 * a * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn) -  
285 - (tn * tn * tn * tn * tn)) / 4.0) / 16.0; -  
286 dp = 35.0 * a * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0 -  
287 * (tn * tn * tn * tn * tn) / 16.0) / 48.0; -  
288 ep = 315.0 * a * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0; -  
289 if ((utmYZone <= 'M' && utmYZone >= 'C') -  
290 || (utmYZone <= 'm' && utmYZone >= 'c')) { -  
291 nfn = 10000000.0; -  
292 } else { -  
293 nfn = 0; -  
294 } -  
295 tmd = (northing - nfn) / ok; -  
296 sr = sphsr (a, eSquared, 0.0); -  
297 ftphi = tmd / sr; -  
298 double t10, t11, t14, t15; -  
299 for (int i = 0; i < 5; i++) { -  
300 t10 = sphtmd (ap, bp, cp, dp, ep, ftphi); -  
301 sr = sphsr (a, eSquared, ftphi); -  
302 ftphi = ftphi + (tmd - t10) / sr; -  
303 } -  
304 sr = sphsr (a, eSquared, ftphi); -  
305 sn = sphsn (a, eSquared, ftphi); -  
306 s = sin (ftphi); -  
307 c = cos (ftphi); -  
308 t = s / c; -  
309 eta = e2Squared * (c * c); -  
310 de = easting - fe; -  
311 t10 = t / (2.0 * sr * sn * (ok * ok)); -  
312 t11 = t * (5.0 + 3.0 * (t * t) + eta - 4.0 * (eta * eta) - 9.0 * (t * t) -  
313 * eta) / (24.0 * sr * (sn * sn * sn) * (ok * ok * ok * ok)); -  
314 lat = ftphi - (de * de) * t10 + (de * de * de * de) * t11; -  
315 t14 = 1.0 / (sn * c * ok); -  
316 t15 = (1.0 + 2.0 * (t * t) + eta) / (6 * (sn * sn * sn) * c -  
317 * (ok * ok * ok)); -  
318 dlam = de * t14 - (de * de * de) * t15; -  
319 olam = (utmXZone * 6 - 183.0) * deg2rad; -  
320 lon = olam + dlam; -  
321 lon *= rad2deg; -  
322 lat *= rad2deg; -  
323 } -  
324   -  
325 //======================================================================= -  
326 // Purpose: -  
327 // This function converts the specified UTM coordinate to a lat/lon -  
328 // coordinate in the WGS84 datum. (See the comment block for the -  
329 // UtmToLatLon() member function. -  
330 //======================================================================= -  
331   -  
332 void UtmToLatLonWGS84 (int utmXZone, char utmYZone, double easting, -  
333 double northing, double& lat, double& lon) -  
334 { -  
335 UtmToLatLon (6378137.0, 1 / 298.257223563, utmXZone, utmYZone, -  
336 easting, northing, lat, lon); -  
337 } -  
338   -  
339 //======================================================================= -  
340 /** -  
341 @func Build a position string -  
342 @parm target. must be 30 characters or longer. -  
343 */ -  
344 //======================================================================= -  
345   -  
346 void CUtmPoint::GetString(TCHAR *position) const -  
347 { -  
348 _stprintf(position, -  
349 _T("%02d%c %06d %07d"), -  
350 m_xzone, m_yzone, -  
351 (int)m_easting, -  
352 (int)m_northing); -  
353 } -  
354   -  
355 //======================================================================= -  
356 /** -  
357 @func get the position of a UTM point -  
358 @parm point out -  
359 */ -  
360 //======================================================================= -  
361   -  
362 void CUtmPoint::ToPosition(CPosition &pos) const -  
363 { -  
364 double lat,lon; -  
365 UtmToLatLonWGS84(m_xzone,m_yzone,m_easting,m_northing, -  
366 lat,lon); -  
367 pos.Clear(); -  
368 pos.SetLatitude(lat); -  
369 pos.SetLongitude(lon); -  
370 } -  
371   -  
372 //======================================================================= -  
373 /** -  
374 @func turn a position into a UTM point -  
375 @parm position -  
376 @rdesc true if it was in range -  
377 */ -  
378 //======================================================================= -  
379   -  
380 bool CUtmPoint::FromPosition(const CPosition &pos) -  
381 { -  
382 Clear(); -  
383 if(!IsPositionInUtmSpace(pos)) -  
384 return false; -  
385 LatLonToUtmWGS84(m_xzone,m_yzone,m_easting,m_northing, -  
386 pos.GetLatitude(), -  
387 pos.GetLongitude()); -  
388 return true; -  
389 } -