vpgl_datum_conversion.cxx
Go to the documentation of this file.
1 #include <cmath>
2 #include "vpgl_datum_conversion.h"
3 //******************************************************************************
4 //:
5 // \file
6 //
7 // This file is a subset of routines that will be used
8 // to create and list the datum coordinate shift formulas. Included are
9 // the following examples of datum shifts, demonstrating the functionality
10 // of these routines.
11 //
12 // These datum shifts have the following interesting property:
13 //
14 // TO GET:
15 // WGS84 from NAD27 -> compute d_lat, d_lon, d_H using NAD27 lat/lon,
16 // THEN
17 // lat(WGS84) = lat(NAD27)+d_lat
18 // lon(WGS84) = lon(NAD27)+d_lon
19 // elev(WGS84) = elev(NAD27)+d_H
20 //
21 // TO GET:
22 // NAD27 from WGS84 -> compute d_lat, d_lon, d_H using WGS84 lat/lon,
23 // THEN
24 // lat(NAD27) = lat(WGS84)-d_lat
25 // lon(NAD27) = lon(WGS84)-d_lon <- note "-" signs
26 // elev(NAD27) = elev(WGS84)-d_H
27 //
28 #ifdef _MSC_VER
29 # include <vcl_msvc_warnings.h>
30 #endif
31 #include <vnl/vnl_math.h>
32 
33 #define degree_to_rad (vnl_math::pi_over_180) // Degree to rad conv.
34 #define dcos(x) std::cos((x)*vnl_math::pi_over_180)
35 #define dsin(x) std::sin((x)*vnl_math::pi_over_180)
36 #define EPSILON 1.0e-12
37 
38 
39 double ipow(double x, int i)
40 {
41  double y = 1.0;
42  if (i<0) { i *= -1; x = y/x; }
43  while (i) { if (i%2) y*=x; i/=2; x*=x; }
44  return y;
45 }
46 
47 
48 //*************************************************************************
49 // FUNCTION: wgs72_to_wgs84_deltas()
50 // Computes deltas that are either added or
51 // subtracted instead of added to the input lat, lon, elev.
52 //*************************************************************************
53 
54 void wgs72_to_wgs84_deltas
55  (double phi, //!< input lat (degrees)
56  double *delta_phi, //!< lat shift (degrees)
57  double *delta_lamda, //!< lon shift (degrees)
58  double *delta_hgt) //!< elev shift (meters)
59 {
60  // Information found from "Supplement to Department of Defense
61  // world geodetic systems 1984 Technical Report". (pg 25-8)
62 
63  double delta_f, a, delta_a, delta_r; // constants
64 
65  //****************************** Instructions: ******************************
66  //
67  // To obtain WGS 84 Coordinates:
68  // Add the delta_phi, delta_lamda, and delta_height changes calculated
69  // using WGS 72 Coordinates to the WGS 72 Coordinates (phi, lamda, height)
70  // given.
71  //
72  // NOTE:Latitude is positive North, and Longitude is positive East (0 to 360)
73  //***************************************************************************
74 
75  // Constant definition
76  delta_f = .3121057e-7;
77  a = 6378135; // meters
78  delta_a = 2.000; // meters
79  delta_r = 1.400; // meters
80 
81  // First compute delta_phi and delta_lamda in arc seconds
82 
83  *delta_phi = (4.5 * dcos(phi))/(a * dsin(1/3600.0)) +
84  (delta_f * dsin(2.0 * phi)) / (dsin(1/3600.0));
85 
86  *delta_lamda = .554;
87 
88  *delta_hgt = 4.5 * dsin(phi) + a * delta_f * (dsin(phi) * dsin(phi))
89  - delta_a + delta_r;
90 
91  // Convert to decimal degrees
92  *delta_phi /= 3600.0;
93  *delta_lamda /= 3600.0;
94 }
95 
96 //*************************************************************************
97 // FUNCTION: wgs72_to_wgs84()
98 //*************************************************************************
99 
100 void wgs72_to_wgs84
101  (
102  double phi, //!< input lat, lon, elev coord (degrees)
103  double lamda,
104  double height,
105  double *wgs84_phi, //!< lat shift (degrees)
106  double *wgs84_lamda, //!< lon shift (degrees)
107  double *wgs84_hgt) //!< elev shift (meters)
108 {
109  double delta_phi, delta_lamda, delta_hgt;
110 
111  wgs72_to_wgs84_deltas(phi, &delta_phi, &delta_lamda, &delta_hgt);
112 
113  *wgs84_phi = phi + delta_phi;
114  *wgs84_lamda = lamda + delta_lamda;
115  *wgs84_hgt = height + delta_hgt;
116 }
117 
118 //*************************************************************************
119 // FUNCTION: wgs84_to_wgs72()
120 //*************************************************************************
121 
122 void wgs84_to_wgs72
123  (
124  double phi, //!< input lat, lon, elev coord (degrees)
125  double lamda,
126  double height,
127  double *wgs72_phi, //!< lat shift (degrees)
128  double *wgs72_lamda, //!< lon shift (degrees)
129  double *wgs72_hgt) //!< elev shift (meters)
130 {
131  double delta_phi, delta_lamda, delta_hgt;
132 
133  wgs72_to_wgs84_deltas(phi, &delta_phi, &delta_lamda, &delta_hgt);
134 
135  *wgs72_phi = phi - delta_phi;
136  *wgs72_lamda = lamda - delta_lamda;
137  *wgs72_hgt = height - delta_hgt;
138 }
139 
140 //*************************************************************************
141 // FUNCTION: nad27m_to_wgs84_deltas()
142 // Area: Mexico/Central America
143 //*************************************************************************
144 
145 void nad27m_to_wgs84_deltas
146  (double phi, //!< input lat (degrees)
147  double lamda, //!< input lon (degrees)
148  double /* height */, //!< input elv (meters)
149  double *delta_phi, //!< lat shift (arc sec)
150  double *delta_lamda, //!< lon shift (arc sec)
151  double *delta_hgt) //!< elev shift (meters)
152 {
153  // Information found from "Supplement to Department of Defense
154  // world geodetic systems 1984 Technical Report. (pg 20-36)
155 
156 
157  double U,V, K;
158 
159  K = .05235988;
160 
161  // Convert the longitude from -180 -> +180 to 0 -> 360
162  if (lamda < 0)
163  lamda += 360.0;
164 
165  // deta_phi and delta_lamda are in arc seconds
166 
167  U = K * (phi - 20);
168  V = K * (lamda -290);
169 
170  *delta_phi = .67775 + 10.02259*U - 18.72391*V + 12.66341*ipow(U,2) - 32.47686*ipow(V,2) - 52.70768*ipow(U,3) - 13.86527*U*ipow(V,2)
171  - 16.82734*ipow(V,3) - 59.54646*V*ipow(U,3) + 120.64887*ipow(U,5) + 7.23362*ipow(U,2)*ipow(V,3) + 95.89131*ipow(U,5)*V -2.89651*U*ipow(V,5)
172  - 1.23778*ipow(V,6) + 70.19213*ipow(U,8) + .84596*ipow(U,2)*ipow(V,6) + .14848*ipow(V,8) - 3.83786*ipow(U,3)*ipow(V,6) + .10405*U*ipow(V,9)
173  + 30.09795*ipow(U,8)*ipow(V,3) + 12.93585*ipow(U,6)*ipow(V,6) - .58747*ipow(U,3)*ipow(V,9) + .7485*ipow(U,4)*ipow(V,9) +
174  6.32462*ipow(U,8)*ipow(V,7) - .12736*ipow(U,6)*ipow(V,9) - .59993*ipow(U,9)*ipow(V,9);
175 
176  *delta_lamda = -.33643 - 3.11777*V -5.16881*ipow(U,2) - 3.17318*ipow(V,2) - 34.59331*ipow(U,3) + .97215*U*ipow(V,2) - .5818*ipow(V,3) +
177  223.78114*ipow(U,5) + 358.07266*ipow(U,6) + 270.68064*ipow(U,7) - 87.99549*ipow(U,6)*V - .09789*U*ipow(V,6) + 636.41982*ipow(U,7)*V -
178  729.91514*ipow(U,6)*ipow(V,2) - .60122*ipow(U,3)*ipow(V,5) - 375.38863*ipow(U,6)*ipow(V,3) - 6.97538*ipow(U,4)*ipow(V,5) -
179  .00138*ipow(V,9) - 363.45977*ipow(U,9)*V + 50.19955*ipow(U,8)*ipow(V,2) - 12.01575*ipow(U,5)*ipow(V,6) - .4314*ipow(U,6)*ipow(V,9) +
180  .43907*ipow(U,9)*ipow(V,9);
181 
182  *delta_hgt = -20.013 + 12.815*U - 8.084*V +74.122*ipow(U,2)*V +39.523*ipow(U,2)*ipow(V,2) - 23.158*ipow(U,4)*V + 194.444*ipow(U,6) +
183  36.797*ipow(U,5)*V + .181*U*ipow(V,5) - 292.155*ipow(U,8) + 3.749*ipow(U,5)*ipow(V,5) +2.537*ipow(U,8)*ipow(V,6);
184 
185  // Convert to decimal degrees
186  *delta_phi /= 3600.0;
187  *delta_lamda /= 3600.0;
188 }
189 
190 
191 //*************************************************************************
192 // FUNCTION: nad27m_to_wgs84()
193 //*************************************************************************
194 
195 void nad27m_to_wgs84
196  (double phi, //!< input lat, lon, elev coord (degrees)
197  double lamda,
198  double height,
199  double *wgs84_phi, //!< lat shift (degrees)
200  double *wgs84_lamda, //!< lon shift (degrees)
201  double *wgs84_hgt) //!< elev shift (meters)
202 {
203  double delta_phi, delta_lamda, delta_hgt;
204 
205  nad27m_to_wgs84_deltas(phi, lamda, height,
206  &delta_phi, &delta_lamda, &delta_hgt);
207 
208  *wgs84_phi = phi + delta_phi;
209  *wgs84_lamda = lamda + delta_lamda;
210  *wgs84_hgt = height + delta_hgt;
211 }
212 
213 
214 //*************************************************************************
215 // FUNCTION: wgs84_to_nad27m()
216 //*************************************************************************
217 
218 void wgs84_to_nad27m
219  (double phi, //!< input lat, lon, elev coord (degrees)
220  double lamda,
221  double height,
222  double *nad27m_phi, //!< lat shift (degrees)
223  double *nad27m_lamda, //!< lon shift (degrees)
224  double *nad27m_hgt) //!< elev shift (meters)
225 {
226  double delta_phi, delta_lamda, delta_hgt;
227 
228  nad27m_to_wgs84_deltas(phi, lamda, height,
229  &delta_phi, &delta_lamda, &delta_hgt);
230 
231  *nad27m_phi = phi - delta_phi;
232  *nad27m_lamda = lamda - delta_lamda;
233  *nad27m_hgt = height - delta_hgt;
234 }
235 
236 
237 //*************************************************************************
238 // FUNCTION: nad27n_to_wgs84_deltas()
239 //*************************************************************************
240 
241 void nad27n_to_wgs84_deltas
242  (double phi, //!< input lat (degrees)
243  double lamda, //!< input lon (degrees)
244  double /* height */, //!< input elv (meters)
245  double *delta_phi, //!< lat shift (arc sec)
246  double *delta_lamda, //!< lon shift (arc sec)
247  double *delta_hgt) //!< elev shift (meters)
248 {
249  // Information found from "Supplement to Department of Defense
250  // world geodetic systems 1984 Technical Report. (pg 19-2)
251 
252  double U,V, K;
253 
254  K = .05235988;
255 
256  // Convert the longitude from -180 -> +180 to 0 -> 360
257  if (lamda < 0)
258  lamda += 360.0;
259 
260  U = K * (phi - 37.0);
261  V = K * (lamda -265.0);
262 
263  // deta_phi and delta_lamda are in arc seconds
264 
265  *delta_phi = .16984 - .76173*U + .09585*V +1.09919*ipow(U,2) - 4.57801*ipow(U,3) - 1.13239*ipow(U,2)*V + .49831*ipow(V,3)
266  - .98399*ipow(U,3)*V + .12415*U*ipow(V,3) + .1145*ipow(V,4) +27.05396*ipow(U,5) + 2.03449*ipow(U,4)*V +.73357*ipow(U,2)*ipow(V,3)
267  - .37548*ipow(V,5) - .14197*ipow(V,6) -59.96555*ipow(U,7) + .07439*ipow(V,7) -4.76082*ipow(U,8) + .03385*ipow(V,8) +
268  49.0432*ipow(U,9) - 1.30575*ipow(U,6)*ipow(V,3) - .07653*ipow(U,3)*ipow(V,9) + .08646*ipow(U,4)*ipow(V,9);
269 
270  *delta_lamda = -.88437 + 2.05061*V + .26361*ipow(U,2) - .76804*U*V + .13374*ipow(V,2) - 1.31974*ipow(U,3) - .52162*ipow(U,2)*V -
271  1.05853*U*ipow(V,2) -.49211*ipow(U,2)*ipow(V,2) + 2.17204*U*ipow(V,3) -.06004*ipow(V,4) +.30139*ipow(U,4)*V + 1.88585*U*ipow(V,4)
272  - .81162*U*ipow(V,5) -.05183*ipow(V,6) -.96723*U*ipow(V,6) -.12948*ipow(U,3)*ipow(V,5) + 3.41827*ipow(U,9) -.44507*ipow(U,8)*V
273  +.18882*U*ipow(V,8) -.01444*ipow(V,9) +.04794*U*ipow(V,9) -.59013*ipow(U,9)*ipow(V,3);
274 
275  *delta_hgt = -36.526 +3.9*U -4.723*V -21.553*ipow(U,2) +7.294*U*V +8.886*ipow(V,2) -8.44*ipow(U,2)*V -2.93*U*ipow(V,2) +
276  56.937*ipow(U,4) -58.756*ipow(U,3)*V -4.061*ipow(V,4) +4.447*ipow(U,4)*V +4.903*ipow(U,2)*ipow(V,3) -55.873*ipow(U,6)
277  +212.005*ipow(U,5)*V + 3.081*ipow(V,6) -254.511*ipow(U,7)*V -.756*ipow(V,8) +30.654*ipow(U,8)*V -.122*U*ipow(V,9);
278 
279  // Convert to decimal degrees
280  *delta_phi /= 3600.0;
281  *delta_lamda /= 3600.0;
282 }
283 
284 
285 //*************************************************************************
286 // FUNCTION: nad27n_to_wgs84()
287 //*************************************************************************
288 
289 void nad27n_to_wgs84
290  (double phi, //!< input lat, lon, elev coord (degrees)
291  double lamda,
292  double height,
293  double *wgs84_phi, //!< lat shift (degrees)
294  double *wgs84_lamda, //!< lon shift (degrees)
295  double *wgs84_hgt) //!< elev shift (meters)
296 {
297  double delta_phi, delta_lamda, delta_hgt;
298 
299  nad27n_to_wgs84_deltas(phi, lamda, height,
300  &delta_phi, &delta_lamda, &delta_hgt);
301 
302  *wgs84_phi = phi + delta_phi;
303  *wgs84_lamda = lamda + delta_lamda;
304  *wgs84_hgt = height + delta_hgt;
305 }
306 
307 //*************************************************************************
308 // FUNCTION: wgs84_to_nad27n()
309 //*************************************************************************
310 
311 void wgs84_to_nad27n
312  (double phi, //!< input lat, lon, elev coord (degrees)
313  double lamda,
314  double height,
315  double *nad27n_phi, //!< lat shift (degrees)
316  double *nad27n_lamda, //!< lon shift (degrees)
317  double *nad27n_hgt) //!< elev shift (meters)
318 {
319  double delta_phi, delta_lamda, delta_hgt;
320 
321  nad27n_to_wgs84_deltas(phi, lamda, height,
322  &delta_phi, &delta_lamda, &delta_hgt);
323 
324  *nad27n_phi = phi - delta_phi;
325  *nad27n_lamda = lamda - delta_lamda;
326  *nad27n_hgt = height - delta_hgt;
327 }
328 
329 //*************************************************************************
330 // ROUTINE: nad27_wgs84_alternate
331 //
332 // REVISION: October 8, 1990 - Original
333 // REVISION: May 15, 1992 - Rajiv Gupta
334 //
335 // This program is a test program designed to convert a
336 // latitude/longitude coordinate in NAD27 to the WGS84
337 // coordinate system.
338 //*************************************************************************
339 
340 
341 //*************************************************************************
342 // FUNCTION: nad27n_to_wgs84_alternate()
343 // AUTHOR: Mark Young
344 //*************************************************************************
345 
346 void nad27n_to_wgs84_alternate
347  (double nad27_lat, double nad27_lon, double nad27_el,
348  double *wgs84_lat, double *wgs84_lon, double *wgs84_el)
349 {
350  double u,v,k;
351 
352  double a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14;
353  double a15,a16,a17,a18,a19,a20,a21,a22,a23;
354 
355  double b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14;
356  double b15,b16,b17,b18,b19,b20,b21,b22,b23;
357 
358  double c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14;
359  double c15,c16,c17,c18,c19,c20;
360 
361  double delta_lat;
362  double delta_lon;
363  double delta_H;
364 #ifdef DEBUG_DATUM
365  double N_height;
366 #endif
367 
368  double delta_lat_p;
369  double delta_lon_p;
370 
371 #ifdef DEBUG_DATUM
372  double prin_lat_rad;
373  double ft_lat;
374  double ft_lon;
375  double c_error;
376 #endif
377 
378 #if 0
379  int lon_int;
380  double lon_int_fl;
381  double lon_frac;
382 #endif
383 
384 #ifdef DEBUG_DATUM
385  double new_lat;
386  double new_lon;
387 #endif
388 
389  double prin_lat; // Est. lat of principal point
390  double prin_lon; // Est. lon of principal point
391 
392 #ifdef DEBUG_DATUM
393  // Display the program title.
394 
395  std::cout << "\n NAD27 to WGS84 Conversion routine\n\n"
396  << " Enter latitude and longitude : ";
397  std::cin >> nad27_lat >> nad27_lon;
398 
399 #endif
400 
401  prin_lat = nad27_lat;
402  prin_lon = nad27_lon;
403 
404 #ifdef DEBUG_DATUM
405  prin_lat_rad = prin_lat * degree_to_rad;
406 #endif
407 
408 //lon_int = (int) prin_lon;
409 //lon_int_fl = (double)lon_int;
410 //lon_frac = std::fabs(prin_lon - lon_int);
411 
412  if (prin_lon < 0)
413  prin_lon = 360.0 + prin_lon;
414 // prin_lon = 360.0 + lon_int_fl + lon_frac;
415 
416  k = .05235988;
417  u = k*(prin_lat - 37);
418  v = k*(prin_lon - 265);
419 
420  // Initialize latitude constants
421 
422  a1 = 0.16984;
423  a2 = -0.76173;
424  a3 = 0.09585;
425  a4 = 1.09919;
426  a5 = -4.57801;
427  a6 = -1.13239;
428  a7 = 0.49831;
429  a8 = -.98399;
430  a9 = 0.12415;
431  a10 = 0.11450;
432  a11 = 27.05396;
433  a12 = 2.03449;
434  a13 = 0.73357;
435  a14 = -0.37548;
436  a15 = -0.14197;
437  a16 = -59.96555;
438  a17 = 0.07439;
439  a18 = -4.76082;
440  a19 = 0.03385;
441  a20 = 49.04320;
442  a21 = -1.30575;
443  a22 = -0.07653;
444  a23 = 0.08646;
445 
446  delta_lat = a1 + a2*u + a3*v + a4*u*u + a5*u*u*u + a6*u*u*v
447  + a7*v*v*v + a8*u*u*u*v + a9*u*v*v*v + a10*v*v*v*v
448  + a11*u*u*u*u*u + a12*u*u*u*u*v + a13*u*u*v*v*v
449  + a14*v*v*v*v*v + a15*v*v*v*v*v*v + a16*u*u*u*u*u*u*u
450  + a17*v*v*v*v*v*v*v + a18*u*u*u*u*u*u*u*u
451  + a19*v*v*v*v*v*v*v*v + a20*u*u*u*u*u*u*u*u*u
452  + a21*u*u*u*u*u*u*v*v*v + a22*u*u*u*v*v*v*v*v*v*v*v*v
453  + a23*u*u*u*u*v*v*v*v*v*v*v*v*v;
454 
455  // Initialize longitude constants
456 
457  b1 = -0.88437;
458  b2 = 2.05061;
459  b3 = 0.26361;
460  b4 = -0.76804;
461  b5 = 0.13374;
462  b6 = -1.31974;
463  b7 = -0.52162;
464  b8 = -1.05853;
465  b9 = -0.49211;
466  b10 = 2.17204;
467  b11 = -0.06004;
468  b12 = 0.30139;
469  b13 = 1.88585;
470  b14 = -0.81162;
471  b15 = -0.05183;
472  b16 = -0.96723;
473  b17 = -.12948;
474  b18 = 3.41827;
475  b19 = -0.44507;
476  b20 = 0.18882;
477  b21 = -0.01444;
478  b22 = 0.04794;
479  b23 = -0.59013;
480 
481  delta_lon = b1 + b2*v + b3*u*u + b4*u*v + b5*v*v + b6*u*u*u
482  + b7*u*u*v + b8*u*v*v + b9*u*u*v*v + b10*u*v*v*v
483  + b11*v*v*v*v + b12*u*u*u*u*v + b13*u*v*v*v*v
484  + b14*u*v*v*v*v*v + b15*v*v*v*v*v*v
485  + b16*u*v*v*v*v*v*v + b17*u*u*u*v*v*v*v*v
486  + b18*u*u*u*u*u*u*u*u*u + b19*u*u*u*u*u*u*u*u*v
487  + b20*u*v*v*v*v*v*v*v*v + b21*v*v*v*v*v*v*v*v*v
488  + b22*u*v*v*v*v*v*v*v*v*v + b23*u*u*u*u*u*u*u*u*u*v*v*v;
489 
490 
491  c1 = -36.526;
492  c2 = 3.900;
493  c3 = -4.723;
494  c4 = -21.553;
495  c5 = 7.294;
496  c6 = 8.886;
497  c7 = -8.440;
498  c8 = -2.930;
499  c9 = 56.937;
500  c10 = -58.756;
501  c11 = -4.061;
502  c12 = 4.447;
503  c13 = 4.903;
504  c14 = -55.873;
505  c15 = 212.005;
506  c16 = 3.081;
507  c17 = -254.511;
508  c18 = -0.756;
509  c19 = 30.654;
510  c20 = -0.122;
511 
512  delta_H = c1 + c2*u + c3*v + c4*u*u + c5*u*v + c6*v*v + c7*u*u*v
513  + c8*u*v*v + c9*u*u*u*u + c10*u*u*u*v + c11*v*v*v*v
514  + c12*u*u*u*u*v + c13*u*u*v*v*v + c14*u*u*u*u*u*u
515  + c15*u*u*u*u*u*v + c16*v*v*v*v*v*v + c17*u*u*u*u*u*u*u*v
516  + c18*v*v*v*v*v*v*v*v + c19*u*u*u*u*u*u*u*u*v
517  + c20*v*v*v*v*v*v*v*v*v;
518 
519 
520 #ifdef DEBUG_DATUM
521  double d1 = 5.068;
522  double d2 = -11.570;
523  double d3 = -8.574;
524  double d4 = 27.839;
525  double d5 = -51.911;
526  double d6 = 29.496;
527  double d7 = 28.343;
528  double d8 = 24.481;
529  double d9 = 11.424;
530  double d10 = 132.550;
531  double d11 = -110.232;
532  double d12 = 41.018;
533  double d13 = -64.953;
534  double d14 = -128.293;
535  double d15 = 51.241;
536  double d16 = -4.326;
537  double d17 = 104.097;
538  double d18 = -128.031;
539  double d19 = 110.694;
540  double d20 = 36.330;
541  double d21 = 243.149;
542  double d22 = -15.790;
543  double d23 = -38.043;
544  double d24 = -40.277;
545  double d25 = 2.746;
546  double d26 = -7.321;
547  double d27 = -394.404;
548  double d28 = -927.540;
549  double d29 = 63.390;
550  double d30 = 10.626;
551  double d31 = -0.520;
552  double d32 = -117.207;
553  double d33 = 16.352;
554 
555  N_height = d1 + d2*u + d3*v + d4*u*u + d5*u*v + d6*v*v + d7*u*u*u
556  + d8*u*v*v + d9*v*v*v + d10*u*u*u*v + d11*u*u*v*v
557  + d12*u*v*v*v + d13*v*v*v*v + d14*u*u*u*v*v + d15*u*v*v*v*v
558  + d16*v*v*v*v*v + d17*u*u*u*u*v*v + d18*u*u*u*v*v*v
559  + d19*u*u*v*v*v*v + d20*v*v*v*v*v*v + d21*u*u*u*u*u*u*v
560  + d22*u*u*v*v*v*v*v + d23*u*v*v*v*v*v*v + d24*u*u*v*v*v*v*v*v
561  + d25*u*v*v*v*v*v*v*v + d26*v*v*v*v*v*v*v*v
562  + d27*u*u*u*u*u*u*u*u*u + d28*u*u*u*u*u*u*u*u*v
563  + d29*u*u*u*u*v*v*v*v*v + d30*u*v*v*v*v*v*v*v*v
564  + d31*u*v*v*v*v*v*v*v*v*v + d32*u*u*u*u*u*u*u*u*v*v*v*v
565  + d33*u*u*u*u*u*v*v*v*v*v*v*v*v;
566 #endif
567 
568  delta_lat_p = delta_lat; // NOTE: delta_lat is in arc seconds
569  delta_lon_p = delta_lon; // NOTE: delta_lon is in arc seconds
570 
571 #ifdef DEBUG_DATUM
572  ft_lat = 101.0 * delta_lat_p;
573  ft_lon = 101.0 * std::cos(prin_lat_rad) * delta_lon_p;
574 
575  c_error = (double)std::sqrt(ft_lat*ft_lat + ft_lon*ft_lon);
576 #endif
577 
578  *wgs84_lat = nad27_lat + (delta_lat_p/3600.0);
579  *wgs84_lon = nad27_lon + (delta_lon_p/3600.0);
580  *wgs84_el = nad27_el + delta_H;
581 
582 #ifdef DEBUG_DATUM
583  new_lat = *wgs84_lat;
584  new_lon = *wgs84_lon;
585  std::cout << "\n d_lat = " << delta_lat
586  << "\n d_lon = " << delta_lon
587  << "\n d_H = " << delta_H
588  << "\n N = " << N_height
589  << "\n circular error = " << c_error
590  << "\n New Lat = " << new_lat << ", New Lon =\n" << new_lon
591  << '\n';
592 #endif
593 }
594 
595 //*************************************************************************
596 // FUNCTION: wgs84_to_nad27n_alternate()
597 // AUTHOR: Mark Young
598 // This is a copy of the above routine in which the deltas are
599 // subtracted instead of added to the input lat, lon, elev.
600 //*************************************************************************
601 
602 void wgs84_to_nad27n_alternate
603  (double wgs84_lat, double wgs84_lon, double wgs84_el,
604  double *nad27n_lat, double *nad27n_lon, double *nad27n_el)
605 {
606  double u,v,k;
607 
608  double a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14;
609  double a15,a16,a17,a18,a19,a20,a21,a22,a23;
610 
611  double b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14;
612  double b15,b16,b17,b18,b19,b20,b21,b22,b23;
613 
614  double c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14;
615  double c15,c16,c17,c18,c19,c20;
616 
617  double delta_lat;
618  double delta_lon;
619  double delta_H;
620 #ifdef DEBUG_DATUM
621  double N_height;
622 #endif
623 
624  double delta_lat_p;
625  double delta_lon_p;
626 
627 #ifdef DEBUG_DATUM
628  double prin_lat_rad;
629  double ft_lat;
630  double ft_lon;
631  double c_error;
632 #endif
633 
634 #if 0
635  int lon_int;
636  double lon_int_fl;
637  double lon_frac;
638 #endif
639 
640 #ifdef DEBUG_DATUM
641  double new_lat;
642  double new_lon;
643 #endif
644 
645  double prin_lat; // Est. lat of principal point
646  double prin_lon; // Est. lon of principal point
647 
648 #ifdef DEBUG_DATUM
649  // Display the program title.
650 
651  std::cout << "\n wgs84 to NAD27N Conversion routine\n\n";
652  " Enter latitude, longitude ";
653  std::cin >> wgs84_lat >> wgs84_lon;
654 
655 #endif
656 
657  prin_lat = wgs84_lat;
658  prin_lon = wgs84_lon;
659 
660 #ifdef DEBUG_DATUM
661  prin_lat_rad = prin_lat * degree_to_rad;
662 #endif
663 
664 //lon_int = (int) prin_lon;
665 //lon_int_fl = (double)lon_int;
666 //lon_frac = std::fabs(prin_lon - lon_int);
667 
668  if (prin_lon < 0)
669  prin_lon = 360.0 + prin_lon;
670 // prin_lon = 360.0 + lon_int_fl + lon_frac;
671 
672  k = .05235988;
673  u = k*(prin_lat - 37);
674  v = k*(prin_lon - 265);
675 
676  // Initialize latitude constants
677 
678  a1 = 0.16984;
679  a2 = -0.76173;
680  a3 = 0.09585;
681  a4 = 1.09919;
682  a5 = -4.57801;
683  a6 = -1.13239;
684  a7 = 0.49831;
685  a8 = -.98399;
686  a9 = 0.12415;
687  a10 = 0.11450;
688  a11 = 27.05396;
689  a12 = 2.03449;
690  a13 = 0.73357;
691  a14 = -0.37548;
692  a15 = -0.14197;
693  a16 = -59.96555;
694  a17 = 0.07439;
695  a18 = -4.76082;
696  a19 = 0.03385;
697  a20 = 49.04320;
698  a21 = -1.30575;
699  a22 = -0.07653;
700  a23 = 0.08646;
701 
702  delta_lat = a1 + a2*u + a3*v + a4*u*u + a5*u*u*u + a6*u*u*v
703  + a7*v*v*v + a8*u*u*u*v + a9*u*v*v*v + a10*v*v*v*v
704  + a11*u*u*u*u*u + a12*u*u*u*u*v + a13*u*u*v*v*v
705  + a14*v*v*v*v*v + a15*v*v*v*v*v*v + a16*u*u*u*u*u*u*u
706  + a17*v*v*v*v*v*v*v + a18*u*u*u*u*u*u*u*u
707  + a19*v*v*v*v*v*v*v*v + a20*u*u*u*u*u*u*u*u*u
708  + a21*u*u*u*u*u*u*v*v*v + a22*u*u*u*v*v*v*v*v*v*v*v*v
709  + a23*u*u*u*u*v*v*v*v*v*v*v*v*v;
710 
711  // Initialize longitude constants
712 
713  b1 = -0.88437;
714  b2 = 2.05061;
715  b3 = 0.26361;
716  b4 = -0.76804;
717  b5 = 0.13374;
718  b6 = -1.31974;
719  b7 = -0.52162;
720  b8 = -1.05853;
721  b9 = -0.49211;
722  b10 = 2.17204;
723  b11 = -0.06004;
724  b12 = 0.30139;
725  b13 = 1.88585;
726  b14 = -0.81162;
727  b15 = -0.05183;
728  b16 = -0.96723;
729  b17 = -.12948;
730  b18 = 3.41827;
731  b19 = -0.44507;
732  b20 = 0.18882;
733  b21 = -0.01444;
734  b22 = 0.04794;
735  b23 = -0.59013;
736 
737  delta_lon = b1 + b2*v + b3*u*u + b4*u*v + b5*v*v + b6*u*u*u
738  + b7*u*u*v + b8*u*v*v + b9*u*u*v*v + b10*u*v*v*v
739  + b11*v*v*v*v + b12*u*u*u*u*v + b13*u*v*v*v*v
740  + b14*u*v*v*v*v*v + b15*v*v*v*v*v*v
741  + b16*u*v*v*v*v*v*v + b17*u*u*u*v*v*v*v*v
742  + b18*u*u*u*u*u*u*u*u*u + b19*u*u*u*u*u*u*u*u*v
743  + b20*u*v*v*v*v*v*v*v*v + b21*v*v*v*v*v*v*v*v*v
744  + b22*u*v*v*v*v*v*v*v*v*v + b23*u*u*u*u*u*u*u*u*u*v*v*v;
745 
746 
747  c1 = -36.526;
748  c2 = 3.900;
749  c3 = -4.723;
750  c4 = -21.553;
751  c5 = 7.294;
752  c6 = 8.886;
753  c7 = -8.440;
754  c8 = -2.930;
755  c9 = 56.937;
756  c10 = -58.756;
757  c11 = -4.061;
758  c12 = 4.447;
759  c13 = 4.903;
760  c14 = -55.873;
761  c15 = 212.005;
762  c16 = 3.081;
763  c17 = -254.511;
764  c18 = -0.756;
765  c19 = 30.654;
766  c20 = -0.122;
767 
768  delta_H = c1 + c2*u + c3*v + c4*u*u + c5*u*v + c6*v*v + c7*u*u*v
769  + c8*u*v*v + c9*u*u*u*u + c10*u*u*u*v + c11*v*v*v*v
770  + c12*u*u*u*u*v + c13*u*u*v*v*v + c14*u*u*u*u*u*u
771  + c15*u*u*u*u*u*v + c16*v*v*v*v*v*v + c17*u*u*u*u*u*u*u*v
772  + c18*v*v*v*v*v*v*v*v + c19*u*u*u*u*u*u*u*u*v
773  + c20*v*v*v*v*v*v*v*v*v;
774 
775 #ifdef DEBUG_DATUM
776  d1 = 5.068;
777  d2 = -11.570;
778  d3 = -8.574;
779  d4 = 27.839;
780  d5 = -51.911;
781  d6 = 29.496;
782  d7 = 28.343;
783  d8 = 24.481;
784  d9 = 11.424;
785  d10 = 132.550;
786  d11 = -110.232;
787  d12 = 41.018;
788  d13 = -64.953;
789  d14 = -128.293;
790  d15 = 51.241;
791  d16 = -4.326;
792  d17 = 104.097;
793  d18 = -128.031;
794  d19 = 110.694;
795  d20 = 36.330;
796  d21 = 243.149;
797  d22 = -15.790;
798  d23 = -38.043;
799  d24 = -40.277;
800  d25 = 2.746;
801  d26 = -7.321;
802  d27 = -394.404;
803  d28 = -927.540;
804  d29 = 63.390;
805  d30 = 10.626;
806  d31 = -0.520;
807  d32 = -117.207;
808  d33 = 16.352;
809 
810  N_height = d1 + d2*u + d3*v + d4*u*u + d5*u*v + d6*v*v + d7*u*u*u
811  + d8*u*v*v + d9*v*v*v + d10*u*u*u*v + d11*u*u*v*v
812  + d12*u*v*v*v + d13*v*v*v*v + d14*u*u*u*v*v + d15*u*v*v*v*v
813  + d16*v*v*v*v*v + d17*u*u*u*u*v*v + d18*u*u*u*v*v*v
814  + d19*u*u*v*v*v*v + d20*v*v*v*v*v*v + d21*u*u*u*u*u*u*v
815  + d22*u*u*v*v*v*v*v + d23*u*v*v*v*v*v*v + d24*u*u*v*v*v*v*v*v
816  + d25*u*v*v*v*v*v*v*v + d26*v*v*v*v*v*v*v*v
817  + d27*u*u*u*u*u*u*u*u*u + d28*u*u*u*u*u*u*u*u*v
818  + d29*u*u*u*u*v*v*v*v*v + d30*u*v*v*v*v*v*v*v*v
819  + d31*u*v*v*v*v*v*v*v*v*v + d32*u*u*u*u*u*u*u*u*v*v*v*v
820  + d33*u*u*u*u*u*v*v*v*v*v*v*v*v;
821 #endif
822 
823  delta_lat_p = delta_lat; // NOTE: delta_lat is in arc seconds
824  delta_lon_p = delta_lon; // NOTE: delta_lon is in arc seconds
825 
826 #ifdef DEBUG_DATUM
827  ft_lat = 101.0 * delta_lat_p;
828  ft_lon = 101.0 * std::cos(prin_lat_rad) * delta_lon_p;
829 
830  c_error = (double)std::sqrt(ft_lat*ft_lat + ft_lon*ft_lon);
831 #endif
832 
833  // SUBTRACT the deltas
834  *nad27n_lat = wgs84_lat - (delta_lat_p/3600.0);
835  *nad27n_lon = wgs84_lon - (delta_lon_p/3600.0);
836  *nad27n_el = wgs84_el - delta_H;
837 
838 #ifdef DEBUG_DATUM
839  new_lat = *nad27n_lat;
840  new_lon = *nad27n_lon;
841  std::cout << "\n d_lat = " << delta_lat
842  << "\n d_lon = " << delta_lon
843  << "\n d_H = " << delta_H
844  << "\n N = " << N_height
845  << "\n circular error = " << c_error
846  << "\n New Lat = " << new_lat << ", New Lon =\n" << new_lon
847  << '\n';
848 #endif
849 }
850 
851 
852 //**********************************************************************
853 // GEO_DETIC2CENTRIC(geodetic_lat, A, B):
854 //**********************************************************************
855 double geo_detic2centric
856  (double geodetic_lat, //!< geodetic latitude of input point
857  double A,
858  double B) //!< Major and minor axes of earth
859 {
860  return(std::atan((B/A)*(B/A) * std::tan(geodetic_lat)));
861 }
862 
863 //**********************************************************************
864 // GEO_CENTRIC2DETIC(geocentric_lat, A, B):
865 //**********************************************************************
866 double geo_centric2detic
867  (double geocentric_lat, //!< geocentric latitude of input point
868  double A,
869  double B) //!< Major and minor axes of earth
870 {
871  return(std::atan((A/B)*(A/B) * std::tan(geocentric_lat)));
872 }
873 
874 
875 //**********************************************************************
876 // GRS_to_latlong (x, y, z, lat, lon, el, A, B):
877 // Computes, from the GRS coordinates of a point, its the lat, long,
878 // and elevation in a GEODETIC coordinate system.
879 // It used an iterative method based on an alternate way
880 // of computing x, y, z from lat, lon, el based on the equations:
881 //
882 // $N = A/\sqrt{1 - e^2 \sin(\mbox{lat})^2}$
883 // $x = (N + \mbox{el}) \cos(\mbox{lat})\cos(\mbox{lon})$
884 // $y = (N + \mbox{el}) \cos(\mbox{lat})\sin(\mbox{lon})$
885 // $z = (N (1-e^2) + \mbox{el}) \sin(\mbox{lat})$
886 //
887 // Nov. 28, 2000 - JLM
888 // The input x, y, z, values are assumed to be in meters
889 // The input elevation and A, B values are assumed to be in meters
890 // The output lat and lon values are in radians
891 //
892 //**********************************************************************
893 void GRS_to_latlong
894  (double x, double y, double z, //!< Input GRS coords
895  double *geodetic_lat,
896  double *lon,
897  double *el, //!< output coordinates of point
898  double A,
899  double B) //!< Major and minor axes of earth
900 {
901  double new_lat; //!< used in iteration
902  double N; //!< Local distance to earth center
903  double xy_dist; //!< dist in x-y plane
904  double ee; //!< eccentricity square
905 
906  xy_dist = std::sqrt(x*x + y*y);
907  ee = 1 - (B/A)*(B/A);
908 
909  // Compute geocentric = geodetic longitude from the
910  // dot-product of (x, y, 0) and (1, 0, 0) and assign
911  // proper sign to it based on y.
912 
913  *lon = std::acos(x/xy_dist);
914  if (y < 0) // Negative y => West of Greenwich meridian
915  *lon = - *lon;
916 
917  // Compute an initial value for geodetic latitude
918  new_lat = std::atan2(z, xy_dist);
919 
920  do
921  {
922  *geodetic_lat = new_lat;
923 
924  // Compute the sin and cos of the approx geodetic lat
925 
926  // Compute the distance N along the height
927  N = A / std::sqrt(1 - ee * std::sin(*geodetic_lat) * std::sin(*geodetic_lat));
928 
929  // new_lat is already in the range -PI to PI
930  new_lat = std::atan2( (z + N * ee * std::sin(*geodetic_lat)), xy_dist );
931  }
932  while (std::fabs(new_lat - *geodetic_lat) > EPSILON);
933 
934  *geodetic_lat = new_lat;
935  *el = xy_dist/std::cos(new_lat) - N;
936 }
937 
938 
939 //**********************************************************************
940 // latlong_to_GRS (lat, lon, el, x, y, z, A, B):
941 // Computes GRS coordinates of a point from its the lat, long,
942 // and elevation in a GEODETIC coordinate system.
943 //
944 // Nov. 28, 2000 - JLM
945 // The input lat and lon values are assumed to be in radians
946 // The input elevation and A, B values are assumed to be in meters
947 // The output x, y, z, values are in meters
948 //
949 //
950 // An alternate way of doing this computation would be:
951 //
952 // N = A/std::sqrt(1 - e*e*std::sin(lat)*std::sin(lat))
953 // x = (N + el)std::cos(lat)std::cos(lon)
954 // y = (N + el)std::cos(lat)std::sin(lon)
955 // z = (N(1-e*e) + el) std::sin(lat)
956 //**********************************************************************
957 
958 void latlong_to_GRS
959  (double geodetic_lat,
960  double lon,
961  double el, //!< Input coordinates of point
962  double *x,
963  double *y,
964  double *z, //!< Output GRS coords
965  double A,
966  double B) //!< Major and minor axes of earth
967 {
968  double geocentric_lat;
969  double local_radius; //!< Local distance to earth center
970  double c, s;
971 
972  geocentric_lat = geo_detic2centric(geodetic_lat, A, B);
973 
974  // Compute the sin and cos of the latitude
975  s = std::sin(geocentric_lat);
976  c = std::cos(geocentric_lat);
977 
978 
979  // Compute the distance to the centre of the earth
980  local_radius = (A*B) / std::sqrt(B*B*c*c + A*A*s*s);
981 
982  *x = (local_radius*c + el*std::cos(geodetic_lat))*std::cos(lon);
983  *y = (local_radius*c + el*std::cos(geodetic_lat))*std::sin(lon);
984  *z = (local_radius*s + el*std::sin(geodetic_lat));
985 
986 #ifdef DEBUG_DATUM
987  double tlat, tlon, tel;
988  GRS_to_latlong(*x, *y, *z, &tlat, &tlon, &tel, A, B);
989  std::cout << "Error in computation: (" << (tlat-geodetic_lat) << ", " << (tlon-lon) << ", " << (tel-el) << ")\n";
990 #endif
991 }
992 
993 //*************************************************************************
994 // FUNCTION: nad27n_to_wgs72()
995 //*************************************************************************
996 
997 void nad27n_to_wgs72
998  (
999  double phi, //!< input lat, lon, elev coord (degrees)
1000  double lamda,
1001  double height,
1002  double *wgs72_phi, //!< lat in wgs72 (degrees)
1003  double *wgs72_lamda, //!< lon in wgs72 (degrees)
1004  double *wgs72_hgt) //!< elev in wgs72 (meters)
1005 {
1006  double wgs84_phi, wgs84_lamda, wgs84_hgt;
1007 
1008  nad27n_to_wgs84(phi, lamda, height,
1009  &wgs84_phi, &wgs84_lamda, &wgs84_hgt);
1010 
1011  wgs84_to_wgs72(wgs84_phi, wgs84_lamda, wgs84_hgt,
1012  wgs72_phi, wgs72_lamda, wgs72_hgt);
1013 }
1014 
1015 
1016 //*************************************************************************
1017 // FUNCTION: wgs72_to_nad27n()
1018 //*************************************************************************
1019 
1020 void wgs72_to_nad27n
1021  (
1022  double phi, //!< input lat, lon, elev coord (degrees)
1023  double lamda,
1024  double height,
1025  double *nad27n_phi, //!< lat in nad27n (degrees)
1026  double *nad27n_lamda, //!< lon in nad27n (degrees)
1027  double *nad27n_hgt) //!< elev in nad27n (meters)
1028 {
1029  double wgs84_phi, wgs84_lamda, wgs84_hgt;
1030 
1031  wgs72_to_wgs84(phi, lamda, height,
1032  &wgs84_phi, &wgs84_lamda, &wgs84_hgt);
1033 
1034  wgs84_to_nad27n(wgs84_phi, wgs84_lamda, wgs84_hgt,
1035  nad27n_phi, nad27n_lamda, nad27n_hgt);
1036 }
geo_centric2detic
double geo_centric2detic(double geocentric_lat, double A, double B)
Definition: vpgl_datum_conversion.cxx:867
degree_to_rad
#define degree_to_rad
Definition: vpgl_datum_conversion.cxx:33
dsin
#define dsin(x)
Definition: vpgl_datum_conversion.cxx:35
vpgl_datum_conversion.h
wgs72_to_wgs84_deltas
void wgs72_to_wgs84_deltas(double phi, double *delta_phi, double *delta_lamda, double *delta_hgt)
Definition: vpgl_datum_conversion.cxx:55
GRS_to_latlong
void GRS_to_latlong(double x, double y, double z, double *geodetic_lat, double *lon, double *el, double A, double B)
Major and minor axes of earth.
Definition: vpgl_datum_conversion.cxx:894
latlong_to_GRS
void latlong_to_GRS(double geodetic_lat, double lon, double el, double *x, double *y, double *z, double A, double B)
Major and minor axes of earth.
Definition: vpgl_datum_conversion.cxx:959
nad27n_to_wgs72
void nad27n_to_wgs72(double phi, double lamda, double height, double *wgs72_phi, double *wgs72_lamda, double *wgs72_hgt)
elev in wgs72 (meters)
Definition: vpgl_datum_conversion.cxx:998
nad27m_to_wgs84
void nad27m_to_wgs84(double phi, double lamda, double height, double *wgs84_phi, double *wgs84_lamda, double *wgs84_hgt)
elev new (meters)
Definition: vpgl_datum_conversion.cxx:196
nad27m_to_wgs84_deltas
void nad27m_to_wgs84_deltas(double phi, double lamda, double, double *delta_phi, double *delta_lamda, double *delta_hgt)
Definition: vpgl_datum_conversion.cxx:146
vnl_math.h
EPSILON
#define EPSILON
Definition: vpgl_datum_conversion.cxx:36
wgs84_to_nad27n
void wgs84_to_nad27n(double phi, double lamda, double height, double *nad27n_phi, double *nad27n_lamda, double *nad27n_hgt)
elev new (meters)
Definition: vpgl_datum_conversion.cxx:312
v
#define v
geo_detic2centric
double geo_detic2centric(double geodetic_lat, double A, double B)
Major and minor axes of earth.
Definition: vpgl_datum_conversion.cxx:856
nad27n_to_wgs84_alternate
void nad27n_to_wgs84_alternate(double nad27_lat, double nad27_lon, double nad27_el, double *wgs84_lat, double *wgs84_lon, double *wgs84_el)
Definition: vpgl_datum_conversion.cxx:347
nad27n_to_wgs84_deltas
void nad27n_to_wgs84_deltas(double phi, double lamda, double, double *delta_phi, double *delta_lamda, double *delta_hgt)
Definition: vpgl_datum_conversion.cxx:242
wgs72_to_wgs84
void wgs72_to_wgs84(double phi, double lamda, double height, double *wgs84_phi, double *wgs84_lamda, double *wgs84_hgt)
elev new (meters)
Definition: vpgl_datum_conversion.cxx:101
nad27n_to_wgs84
void nad27n_to_wgs84(double phi, double lamda, double height, double *wgs84_phi, double *wgs84_lamda, double *wgs84_hgt)
elev new (meters)
Definition: vpgl_datum_conversion.cxx:290
wgs84_to_nad27n_alternate
void wgs84_to_nad27n_alternate(double wgs84_lat, double wgs84_lon, double wgs84_el, double *nad27n_lat, double *nad27n_lon, double *nad27n_el)
Definition: vpgl_datum_conversion.cxx:603
dcos
#define dcos(x)
Definition: vpgl_datum_conversion.cxx:34
ipow
double ipow(double x, int i)
Definition: vpgl_datum_conversion.cxx:39
wgs72_to_nad27n
void wgs72_to_nad27n(double phi, double lamda, double height, double *nad27n_phi, double *nad27n_lamda, double *nad27n_hgt)
elev in nad27n (meters)
Definition: vpgl_datum_conversion.cxx:1021
wgs84_to_wgs72
void wgs84_to_wgs72(double phi, double lamda, double height, double *wgs72_phi, double *wgs72_lamda, double *wgs72_hgt)
elev new (meters)
Definition: vpgl_datum_conversion.cxx:123
wgs84_to_nad27m
void wgs84_to_nad27m(double phi, double lamda, double height, double *nad27m_phi, double *nad27m_lamda, double *nad27m_hgt)
elev new (meters)
Definition: vpgl_datum_conversion.cxx:219
Generated by   doxygen 1.8.15