OSSIM - Open Source Software Image Map  Version 1.9.0 (20180803)
ossimTransMercatorProjection.cpp
Go to the documentation of this file.
1 //*******************************************************************
2 // Copyright (C) 2000 ImageLinks Inc.
3 //
4 // License: See top LICENSE.txt file.
5 //
6 // Author: Garrett Potts
7 //
8 // Description:
9 //
10 // Calls Geotrans Transverse Mercator projection code.
11 //*******************************************************************
12 // $Id: ossimTransMercatorProjection.cpp 23002 2014-11-24 17:11:17Z dburken $
13 #include <cmath>
14 using namespace std;
15 
16 #include <ossim/projection/ossimTransMercatorProjection.h>
17 #include <ossim/base/ossimKeywordNames.h>
18 
19 RTTI_DEF1(ossimTransMercatorProjection, "ossimTransMercatorProjection", ossimMapProjection)
20 
21 /******************************* DEFINES *********************************/
22 
23 #define TRANMERC_NO_ERROR 0x0000
24 #define TRANMERC_LAT_ERROR 0x0001
25 #define TRANMERC_LON_ERROR 0x0002
26 #define TRANMERC_EASTING_ERROR 0x0004
27 #define TRANMERC_NORTHING_ERROR 0x0008
28 #define TRANMERC_ORIGIN_LAT_ERROR 0x0010
29 #define TRANMERC_CENT_MER_ERROR 0x0020
30 #define TRANMERC_A_ERROR 0x0040
31 #define TRANMERC_B_ERROR 0x0080
32 #define TRANMERC_A_LESS_B_ERROR 0x0100
33 #define TRANMERC_SCALE_FACTOR_ERROR 0x0200
34 #define TRANMERC_LON_WARNING 0x0400
35 
36 #define MAX_LAT ((M_PI * 90.0)/180.0) /* 90 degrees in radians */
37 #define MAX_DELTA_LONG ((M_PI * 90.0)/180.0) /* 90 degrees in radians */
38 #define MIN_SCALE_FACTOR 0.3
39 #define MAX_SCALE_FACTOR 3.0
40 
41 #define SPHTMD(Latitude) ((double) (TranMerc_ap * Latitude \
42  - TranMerc_bp * sin(2.e0 * Latitude) + TranMerc_cp * sin(4.e0 * Latitude) \
43  - TranMerc_dp * sin(6.e0 * Latitude) + TranMerc_ep * sin(8.e0 * Latitude) ) )
44 
45 #define SPHSN(Latitude) ((double) (getA() / sqrt( 1.e0 - TranMerc_es * \
46  pow(sin(Latitude), 2))))
47 
48 #define SPHSR(Latitude) ((double) (getA() * (1.e0 - TranMerc_es) / \
49  pow(DENOM(Latitude), 3)))
50 
51 #define DENOM(Latitude) ((double) (sqrt(1.e0 - TranMerc_es * pow(sin(Latitude),2))))
52 
53 ossimTransMercatorProjection::ossimTransMercatorProjection(const ossimEllipsoid& ellipsoid,
54  const ossimGpt& origin)
55  :
56  ossimMapProjection(ellipsoid, origin),
57  TranMerc_a(6378137.0),
58  TranMerc_f(1.0/298.257223563),
59  TranMerc_es(0.0066943799901413800),
60  TranMerc_ebs(0.0067394967565869),
61  TranMerc_Origin_Lat(origin.latr()),
62  TranMerc_Origin_Long(origin.lonr()),
63  TranMerc_False_Northing(0.0),
64  TranMerc_False_Easting(0.0),
65  TranMerc_Scale_Factor(1.0),
66  TranMerc_ap(6367449.1458008),
67  TranMerc_bp(16038.508696861),
68  TranMerc_cp(16.832613334334),
69  TranMerc_dp(0.021984404273757),
70  TranMerc_Delta_Easting(40000000.0),
71  TranMerc_Delta_Northing(40000000.0)
72 {
73  update();
74 }
75 
76 ossimTransMercatorProjection::ossimTransMercatorProjection(const ossimEllipsoid& ellipsoid,
77  const ossimGpt& origin,
78  double falseEasting,
79  double falseNorthing,
80  double scaleFactor)
81  :
82  ossimMapProjection(ellipsoid, origin),
83  TranMerc_a(6378137.0),
84  TranMerc_f(1.0/298.257223563),
85  TranMerc_es(0.0066943799901413800),
86  TranMerc_ebs(0.0067394967565869),
87  TranMerc_Origin_Lat(origin.latr()),
88  TranMerc_Origin_Long(origin.lonr()),
89  TranMerc_False_Northing(falseNorthing),
90  TranMerc_False_Easting(falseEasting),
91  TranMerc_Scale_Factor(scaleFactor),
92  TranMerc_ap(6367449.1458008),
93  TranMerc_bp(16038.508696861),
94  TranMerc_cp(16.832613334334),
95  TranMerc_dp(0.021984404273757),
96  TranMerc_Delta_Easting(40000000.0),
97  TranMerc_Delta_Northing(40000000.0)
98 {
99  update();
100 }
101 
102 void ossimTransMercatorProjection::update()
103 {
104  Set_Transverse_Mercator_Parameters(theEllipsoid.getA(),
105  theEllipsoid.getFlattening(),
106  theOrigin.latr(),
107  theOrigin.lonr(),
108  TranMerc_False_Easting,
109  TranMerc_False_Northing,
110  TranMerc_Scale_Factor);
111 
112  theFalseEastingNorthing.x = TranMerc_False_Easting;
113  theFalseEastingNorthing.y = TranMerc_False_Northing;
114 
115  ossimMapProjection::update();
116 }
117 
118 void ossimTransMercatorProjection::setFalseEasting(double falseEasting)
119 {
120  TranMerc_False_Easting = falseEasting;
121 
122  update();
123 }
124 
125 
126 void ossimTransMercatorProjection::setFalseNorthing(double falseNorthing)
127 {
128  TranMerc_False_Northing = falseNorthing;
129 
130  update();
131 }
132 
133 void ossimTransMercatorProjection::setFalseEastingNorthing(double falseEasting,
134  double falseNorthing)
135 {
136  TranMerc_False_Easting = falseEasting;
137  TranMerc_False_Northing = falseNorthing;
138 
139  update();
140 }
141 
142 void ossimTransMercatorProjection::setScaleFactor(double scaleFactor)
143 {
144  TranMerc_Scale_Factor = scaleFactor;
145 
146  update();
147 }
148 
149 void ossimTransMercatorProjection::setParameters(double falseEasting,
150  double falseNorthing,
151  double scaleFactor)
152 {
153  TranMerc_False_Easting = falseEasting;
154  TranMerc_False_Northing = falseNorthing;
155  TranMerc_Scale_Factor = scaleFactor;
156 
157  update();
158 }
159 
160 void ossimTransMercatorProjection::setDefaults()
161 {
162  TranMerc_False_Easting = 0.0;
163  TranMerc_False_Northing = 0.0;
164  TranMerc_Scale_Factor = 1.0;
165  TranMerc_Delta_Easting = 40000000.0;
166  TranMerc_Delta_Northing = 40000000.0;
167 }
168 
169 ossimGpt ossimTransMercatorProjection::inverse(const ossimDpt &eastingNorthing)const
170 {
171  double lat = 0.0;
172  double lon = 0.0;
173 
174  Convert_Transverse_Mercator_To_Geodetic(eastingNorthing.x,
175  eastingNorthing.y,
176  &lat,
177  &lon);
178 
179  return ossimGpt(lat*DEG_PER_RAD, lon*DEG_PER_RAD, 0.0, theDatum);
180 }
181 
182 ossimDpt ossimTransMercatorProjection::forward(const ossimGpt &latLon)const
183 {
184  double easting = 0.0;
185  double northing = 0.0;
186  ossimGpt gpt = latLon;
187 
188  if (theDatum)
189  {
190  if (theDatum->code() != latLon.datum()->code())
191  {
192  gpt.changeDatum(theDatum); // Shift to our datum.
193  }
194  }
195 
196  Convert_Geodetic_To_Transverse_Mercator(gpt.latr(),
197  gpt.lonr(),
198  &easting,
199  &northing);
200 
201  return ossimDpt(easting, northing);
202 }
203 
204 bool ossimTransMercatorProjection::saveState(ossimKeywordlist& kwl, const char* prefix) const
205 {
206  kwl.add(prefix,
207  ossimKeywordNames::SCALE_FACTOR_KW,
208  TranMerc_Scale_Factor,
209  true);
210 
211  return ossimMapProjection::saveState(kwl, prefix);
212 }
213 
214 bool ossimTransMercatorProjection::loadState(const ossimKeywordlist& kwl,
215  const char* prefix)
216 {
217  bool flag = ossimMapProjection::loadState(kwl, prefix);
218  const char* type = kwl.find(prefix, ossimKeywordNames::TYPE_KW);
219  const char* scaleFactor = kwl.find(prefix, ossimKeywordNames::SCALE_FACTOR_KW);
220 
221  setDefaults();
222  if(ossimString(type) == STATIC_TYPE_NAME(ossimTransMercatorProjection))
223  {
224  TranMerc_False_Easting = theFalseEastingNorthing.x;
225  TranMerc_False_Northing = theFalseEastingNorthing.y;
226 
227  if(scaleFactor)
228  {
229  double d = ossimString(scaleFactor).toDouble();
230  if (d > 0.0) // Check to avoid divide by zero.
231  {
232  TranMerc_Scale_Factor = d;
233  }
234  }
235  }
236 
237  update();
238 
239  return flag;
240 }
241 
242 /************************************************************************/
243 /* FUNCTIONS
244  *
245  */
246 
247 
248 long ossimTransMercatorProjection::Set_Transverse_Mercator_Parameters(double a,
249  double f,
250  double Origin_Latitude,
251  double Central_Meridian,
252  double False_Easting,
253  double False_Northing,
254  double Scale_Factor)
255 
256 { /* BEGIN Set_Tranverse_Mercator_Parameters */
257  /*
258  * The function Set_Tranverse_Mercator_Parameters receives the ellipsoid
259  * parameters and Tranverse Mercator projection parameters as inputs, and
260  * sets the corresponding state variables. If any errors occur, the error
261  * code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
262  * returned.
263  *
264  * a : Semi-major axis of ellipsoid, in meters (input)
265  * f : Flattening of ellipsoid (input)
266  * Origin_Latitude : Latitude in radians at the origin of the (input)
267  * projection
268  * Central_Meridian : Longitude in radians at the center of the (input)
269  * projection
270  * False_Easting : Easting/X at the center of the projection (input)
271  * False_Northing : Northing/Y at the center of the projection (input)
272  * Scale_Factor : Projection scale factor (input)
273  */
274 
275  double tn; /* True Meridianal distance constant */
276  double tn2;
277  double tn3;
278  double tn4;
279  double tn5;
280  double dummy_northing;
281  double TranMerc_b; /* Semi-minor axis of ellipsoid, in meters */
282 // double inv_f = 1 / f;
283  long Error_Code = TRANMERC_NO_ERROR;
284 
285 // if (a <= 0.0)
286 // { /* Semi-major axis must be greater than zero */
287 // Error_Code |= TRANMERC_A_ERROR;
288 // }
289 // if ((inv_f < 250) || (inv_f > 350))
290 // { /* Inverse flattening must be between 250 and 350 */
291 // Error_Code |= TRANMERC_INV_F_ERROR;
292 // }
293 // if ((Origin_Latitude < -MAX_LAT) || (Origin_Latitude > MAX_LAT))
294 // { /* origin latitude out of range */
295 // Error_Code |= TRANMERC_ORIGIN_LAT_ERROR;
296 // }
297 // if ((Central_Meridian < -M_PI) || (Central_Meridian > TWO_PI))
298 // { /* origin longitude out of range */
299 // Error_Code |= TRANMERC_CENT_MER_ERROR;
300 // }
301 // if ((Scale_Factor < MIN_SCALE_FACTOR) || (Scale_Factor > MAX_SCALE_FACTOR))
302 // {
303 // Error_Code |= TRANMERC_SCALE_FACTOR_ERROR;
304 // }
305  if (!Error_Code)
306  { /* no errors */
307  TranMerc_a = a;
308  TranMerc_f = f;
309  TranMerc_Origin_Lat = 0;
310  TranMerc_Origin_Long = 0;
311  TranMerc_False_Northing = 0;
312  TranMerc_False_Easting = 0;
313  TranMerc_Scale_Factor = 1;
314 
315  /* Eccentricity Squared */
316  TranMerc_es = 2 * TranMerc_f - TranMerc_f * TranMerc_f;
317  /* Second Eccentricity Squared */
318  TranMerc_ebs = (1 / (1 - TranMerc_es)) - 1;
319 
320  TranMerc_b = TranMerc_a * (1 - TranMerc_f);
321  /*True meridianal constants */
322  tn = (TranMerc_a - TranMerc_b) / (TranMerc_a + TranMerc_b);
323  tn2 = tn * tn;
324  tn3 = tn2 * tn;
325  tn4 = tn3 * tn;
326  tn5 = tn4 * tn;
327 
328  TranMerc_ap = TranMerc_a * (1.e0 - tn + 5.e0 * (tn2 - tn3)/4.e0
329  + 81.e0 * (tn4 - tn5)/64.e0 );
330  TranMerc_bp = 3.e0 * TranMerc_a * (tn - tn2 + 7.e0 * (tn3 - tn4)
331  /8.e0 + 55.e0 * tn5/64.e0 )/2.e0;
332  TranMerc_cp = 15.e0 * TranMerc_a * (tn2 - tn3 + 3.e0 * (tn4 - tn5 )/4.e0) /16.0;
333  TranMerc_dp = 35.e0 * TranMerc_a * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
334  TranMerc_ep = 315.e0 * TranMerc_a * (tn4 - tn5) / 512.e0;
335  Convert_Geodetic_To_Transverse_Mercator(MAX_LAT,
336  MAX_DELTA_LONG,
337  &TranMerc_Delta_Easting,
338  &TranMerc_Delta_Northing);
339  Convert_Geodetic_To_Transverse_Mercator(0,
340  MAX_DELTA_LONG,
341  &TranMerc_Delta_Easting,
342  &dummy_northing);
343  TranMerc_Origin_Lat = Origin_Latitude;
344  if (Central_Meridian > M_PI)
345  Central_Meridian -= TWO_PI;
346  TranMerc_Origin_Long = Central_Meridian;
347  TranMerc_False_Northing = False_Northing;
348  TranMerc_False_Easting = False_Easting;
349  TranMerc_Scale_Factor = Scale_Factor;
350  } /* END OF if(!Error_Code) */
351  return (Error_Code);
352 } /* END of Set_Transverse_Mercator_Parameters */
353 
354 
355 void ossimTransMercatorProjection::Get_Transverse_Mercator_Parameters(double *a,
356  double *f,
357  double *Origin_Latitude,
358  double *Central_Meridian,
359  double *False_Easting,
360  double *False_Northing,
361  double *Scale_Factor)const
362 
363 { /* BEGIN Get_Tranverse_Mercator_Parameters */
364  /*
365  * The function Get_Transverse_Mercator_Parameters returns the current
366  * ellipsoid and Transverse Mercator projection parameters.
367  *
368  * a : Semi-major axis of ellipsoid, in meters (output)
369  * f : Flattening of ellipsoid (output)
370  * Origin_Latitude : Latitude in radians at the origin of the (output)
371  * projection
372  * Central_Meridian : Longitude in radians at the center of the (output)
373  * projection
374  * False_Easting : Easting/X at the center of the projection (output)
375  * False_Northing : Northing/Y at the center of the projection (output)
376  * Scale_Factor : Projection scale factor (output)
377  */
378 
379  *a = TranMerc_a;
380  *f = TranMerc_f;
381  *Origin_Latitude = TranMerc_Origin_Lat;
382  *Central_Meridian = TranMerc_Origin_Long;
383  *False_Easting = TranMerc_False_Easting;
384  *False_Northing = TranMerc_False_Northing;
385  *Scale_Factor = TranMerc_Scale_Factor;
386 
387  return;
388 } /* END OF Get_Tranverse_Mercator_Parameters */
389 
390 
391 
392 long ossimTransMercatorProjection::Convert_Geodetic_To_Transverse_Mercator (double Latitude,
393  double Longitude,
394  double *Easting,
395  double *Northing)const
396 
397 { /* BEGIN Convert_Geodetic_To_Transverse_Mercator */
398 
399  /*
400  * The function Convert_Geodetic_To_Transverse_Mercator converts geodetic
401  * (latitude and longitude) coordinates to Transverse Mercator projection
402  * (easting and northing) coordinates, according to the current ellipsoid
403  * and Transverse Mercator projection coordinates. If any errors occur, the
404  * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
405  * returned.
406  *
407  * Latitude : Latitude in radians (input)
408  * Longitude : Longitude in radians (input)
409  * Easting : Easting/X in meters (output)
410  * Northing : Northing/Y in meters (output)
411  */
412 
413  double c; /* Cosine of latitude */
414  double c2;
415  double c3;
416  double c5;
417  double c7;
418  double dlam; /* Delta longitude - Difference in Longitude */
419  double eta; /* constant - TranMerc_ebs *c *c */
420  double eta2;
421  double eta3;
422  double eta4;
423  double s; /* Sine of latitude */
424  double sn; /* Radius of curvature in the prime vertical */
425  double t; /* Tangent of latitude */
426  double tan2;
427  double tan3;
428  double tan4;
429  double tan5;
430  double tan6;
431  double t1; /* Term in coordinate conversion formula - GP to Y */
432  double t2; /* Term in coordinate conversion formula - GP to Y */
433  double t3; /* Term in coordinate conversion formula - GP to Y */
434  double t4; /* Term in coordinate conversion formula - GP to Y */
435  double t5; /* Term in coordinate conversion formula - GP to Y */
436  double t6; /* Term in coordinate conversion formula - GP to Y */
437  double t7; /* Term in coordinate conversion formula - GP to Y */
438  double t8; /* Term in coordinate conversion formula - GP to Y */
439  double t9; /* Term in coordinate conversion formula - GP to Y */
440  double tmd; /* True Meridional distance */
441  double tmdo; /* True Meridional distance for latitude of origin */
442  long Error_Code = TRANMERC_NO_ERROR;
443 // double temp_Origin;
444 // double temp_Long;
445 
446 // if ((Latitude < -MAX_LAT) || (Latitude > MAX_LAT))
447 // { /* Latitude out of range */
448 // Error_Code|= TRANMERC_LAT_ERROR;
449 // }
450  if (Longitude > M_PI)
451  Longitude -= TWO_PI;
452 // if ((Longitude < (TranMerc_Origin_Long - MAX_DELTA_LONG))
453 // || (Longitude > (TranMerc_Origin_Long + MAX_DELTA_LONG)))
454 // {
455 // if (Longitude < 0)
456 // temp_Long = Longitude + TWO_PI;
457 // else
458 // temp_Long = Longitude;
459 // if (TranMerc_Origin_Long < 0)
460 // temp_Origin = TranMerc_Origin_Long + TWO_PI;
461 // else
462 // temp_Origin = TranMerc_Origin_Long;
463 // if ((temp_Long < (temp_Origin - MAX_DELTA_LONG))
464 // || (temp_Long > (temp_Origin + MAX_DELTA_LONG)))
465 // Error_Code|= TRANMERC_LON_ERROR;
466 // }
467  if (!Error_Code)
468  { /* no errors */
469 
470  /*
471  * Delta Longitude
472  */
473  dlam = Longitude - TranMerc_Origin_Long;
474 
475 // if (fabs(dlam) > (9.0 * M_PI / 180))
476 // { /* Distortion will result if Longitude is more than 9 degrees from the Central Meridian */
477 // Error_Code |= TRANMERC_LON_WARNING;
478 // }
479 
480  if (dlam > M_PI)
481  dlam -= TWO_PI;
482  if (dlam < -M_PI)
483  dlam += TWO_PI;
484  if (fabs(dlam) < 2.e-10)
485  dlam = 0.0;
486 
487  s = sin(Latitude);
488  c = cos(Latitude);
489  c2 = c * c;
490  c3 = c2 * c;
491  c5 = c3 * c2;
492  c7 = c5 * c2;
493  t = tan (Latitude);
494  tan2 = t * t;
495  tan3 = tan2 * t;
496  tan4 = tan3 * t;
497  tan5 = tan4 * t;
498  tan6 = tan5 * t;
499  eta = TranMerc_ebs * c2;
500  eta2 = eta * eta;
501  eta3 = eta2 * eta;
502  eta4 = eta3 * eta;
503 
504  /* radius of curvature in prime vertical */
505  sn = SPHSN(Latitude);
506 
507  /* True Meridianal Distances */
508  tmd = SPHTMD(Latitude);
509 
510  /* Origin */
511  tmdo = SPHTMD (TranMerc_Origin_Lat);
512 
513  /* northing */
514  t1 = (tmd - tmdo) * TranMerc_Scale_Factor;
515  t2 = sn * s * c * TranMerc_Scale_Factor/ 2.e0;
516  t3 = sn * s * c3 * TranMerc_Scale_Factor * (5.e0 - tan2 + 9.e0 * eta
517  + 4.e0 * eta2) /24.e0;
518 
519  t4 = sn * s * c5 * TranMerc_Scale_Factor * (61.e0 - 58.e0 * tan2
520  + tan4 + 270.e0 * eta - 330.e0 * tan2 * eta + 445.e0 * eta2
521  + 324.e0 * eta3 -680.e0 * tan2 * eta2 + 88.e0 * eta4
522  -600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4) / 720.e0;
523 
524  t5 = sn * s * c7 * TranMerc_Scale_Factor * (1385.e0 - 3111.e0 *
525  tan2 + 543.e0 * tan4 - tan6) / 40320.e0;
526 
527  *Northing = TranMerc_False_Northing + t1 + pow(dlam,2.e0) * t2
528  + pow(dlam,4.e0) * t3 + pow(dlam,6.e0) * t4
529  + pow(dlam,8.e0) * t5;
530 
531  /* Easting */
532  t6 = sn * c * TranMerc_Scale_Factor;
533  t7 = sn * c3 * TranMerc_Scale_Factor * (1.e0 - tan2 + eta ) /6.e0;
534  t8 = sn * c5 * TranMerc_Scale_Factor * (5.e0 - 18.e0 * tan2 + tan4
535  + 14.e0 * eta - 58.e0 * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3
536  - 64.e0 * tan2 * eta2 - 24.e0 * tan2 * eta3 )/ 120.e0;
537  t9 = sn * c7 * TranMerc_Scale_Factor * ( 61.e0 - 479.e0 * tan2
538  + 179.e0 * tan4 - tan6 ) /5040.e0;
539 
540  *Easting = TranMerc_False_Easting + dlam * t6 + pow(dlam,3.e0) * t7
541  + pow(dlam,5.e0) * t8 + pow(dlam,7.e0) * t9;
542  }
543  return (Error_Code);
544 } /* END OF Convert_Geodetic_To_Transverse_Mercator */
545 
546 
547 long ossimTransMercatorProjection::Convert_Transverse_Mercator_To_Geodetic (double Easting,
548  double Northing,
549  double *Latitude,
550  double *Longitude)const
551 { /* BEGIN Convert_Transverse_Mercator_To_Geodetic */
552 
553  /*
554  * The function Convert_Transverse_Mercator_To_Geodetic converts Transverse
555  * Mercator projection (easting and northing) coordinates to geodetic
556  * (latitude and longitude) coordinates, according to the current ellipsoid
557  * and Transverse Mercator projection parameters. If any errors occur, the
558  * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
559  * returned.
560  *
561  * Easting : Easting/X in meters (input)
562  * Northing : Northing/Y in meters (input)
563  * Latitude : Latitude in radians (output)
564  * Longitude : Longitude in radians (output)
565  */
566 
567  double c; /* Cosine of latitude */
568  double de; /* Delta easting - Difference in Easting (Easting-Fe) */
569  double dlam; /* Delta longitude - Difference in Longitude */
570  double eta; /* constant - TranMerc_ebs *c *c */
571  double eta2;
572  double eta3;
573  double eta4;
574  double ftphi; /* Footpoint latitude */
575  int i; /* Loop iterator */
576  // double s; /* Sine of latitude */
577  double sn; /* Radius of curvature in the prime vertical */
578  double sr; /* Radius of curvature in the meridian */
579  double t; /* Tangent of latitude */
580  double tan2;
581  double tan4;
582  double t10; /* Term in coordinate conversion formula - GP to Y */
583  double t11; /* Term in coordinate conversion formula - GP to Y */
584  double t12; /* Term in coordinate conversion formula - GP to Y */
585  double t13; /* Term in coordinate conversion formula - GP to Y */
586  double t14; /* Term in coordinate conversion formula - GP to Y */
587  double t15; /* Term in coordinate conversion formula - GP to Y */
588  double t16; /* Term in coordinate conversion formula - GP to Y */
589  double t17; /* Term in coordinate conversion formula - GP to Y */
590  double tmd; /* True Meridional distance */
591  double tmdo; /* True Meridional distance for latitude of origin */
592  long Error_Code = TRANMERC_NO_ERROR;
593 
594 // if ((Easting < (TranMerc_False_Easting - TranMerc_Delta_Easting))
595 // ||(Easting > (TranMerc_False_Easting + TranMerc_Delta_Easting)))
596 // { /* Easting out of range */
597 // Error_Code |= TRANMERC_EASTING_ERROR;
598 // }
599 // if ((Northing < (TranMerc_False_Northing - TranMerc_Delta_Northing))
600 // || (Northing > (TranMerc_False_Northing + TranMerc_Delta_Northing)))
601 // { /* Northing out of range */
602 // Error_Code |= TRANMERC_NORTHING_ERROR;
603 // }
604 
605  if (!Error_Code)
606  {
607  /* True Meridional Distances for latitude of origin */
608  tmdo = SPHTMD(TranMerc_Origin_Lat);
609 
610  /* Origin */
611  tmd = tmdo + (Northing - TranMerc_False_Northing) / TranMerc_Scale_Factor;
612 
613  /* First Estimate */
614  sr = SPHSR(0.e0);
615  ftphi = tmd/sr;
616 
617  for (i = 0; i < 5 ; i++)
618  {
619  t10 = SPHTMD (ftphi);
620  sr = SPHSR(ftphi);
621  ftphi = ftphi + (tmd - t10) / sr;
622  }
623 
624  /* Radius of Curvature in the meridian */
625  sr = SPHSR(ftphi);
626 
627  /* Radius of Curvature in the meridian */
628  sn = SPHSN(ftphi);
629 
630  /* Sine Cosine terms */
631  // s = sin(ftphi);
632  c = cos(ftphi);
633 
634  /* Tangent Value */
635  t = tan(ftphi);
636  tan2 = t * t;
637  tan4 = tan2 * tan2;
638  eta = TranMerc_ebs * pow(c,2);
639  eta2 = eta * eta;
640  eta3 = eta2 * eta;
641  eta4 = eta3 * eta;
642  de = Easting - TranMerc_False_Easting;
643  if (fabs(de) < 0.0001)
644  de = 0.0;
645 
646  /* Latitude */
647  t10 = t / (2.e0 * sr * sn * pow(TranMerc_Scale_Factor, 2));
648  t11 = t * (5.e0 + 3.e0 * tan2 + eta - 4.e0 * pow(eta,2)
649  - 9.e0 * tan2 * eta) / (24.e0 * sr * pow(sn,3)
650  * pow(TranMerc_Scale_Factor,4));
651  t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4
652  - 252.e0 * tan2 * eta - 3.e0 * eta2 + 100.e0
653  * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4
654  * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2
655  + 84.e0 * tan2* eta3 - 192.e0 * tan2 * eta4)
656  / ( 720.e0 * sr * pow(sn,5) * pow(TranMerc_Scale_Factor, 6) );
657  t13 = t * ( 1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0
658  * pow(t,6))/ (40320.e0 * sr * pow(sn,7) * pow(TranMerc_Scale_Factor,8));
659  *Latitude = ftphi - pow(de,2) * t10 + pow(de,4) * t11 - pow(de,6) * t12
660  + pow(de,8) * t13;
661 
662  t14 = 1.e0 / (sn * c * TranMerc_Scale_Factor);
663 
664  t15 = (1.e0 + 2.e0 * tan2 + eta) / (6.e0 * pow(sn,3) * c *
665  pow(TranMerc_Scale_Factor,3));
666 
667  t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2
668  + 8.e0 * tan2 * eta + 24.e0 * tan4 - 4.e0
669  * eta3 + 4.e0 * tan2 * eta2 + 24.e0
670  * tan2 * eta3) / (120.e0 * pow(sn,5) * c
671  * pow(TranMerc_Scale_Factor,5));
672 
673  t17 = (61.e0 + 662.e0 * tan2 + 1320.e0 * tan4 + 720.e0
674  * pow(t,6)) / (5040.e0 * pow(sn,7) * c
675  * pow(TranMerc_Scale_Factor,7));
676 
677  /* Difference in Longitude */
678  dlam = de * t14 - pow(de,3) * t15 + pow(de,5) * t16 - pow(de,7) * t17;
679 
680  /* Longitude */
681  (*Longitude) = TranMerc_Origin_Long + dlam;
682  while (*Latitude > (90.0 * RAD_PER_DEG))
683  {
684  *Latitude = M_PI - *Latitude;
685  *Longitude += M_PI;
686  if (*Longitude > M_PI)
687  *Longitude -= TWO_PI;
688  }
689 
690  while (*Latitude < (-90.0 * RAD_PER_DEG))
691  {
692  *Latitude = - (*Latitude + M_PI);
693  *Longitude += M_PI;
694  if (*Longitude > M_PI)
695  *Longitude -= TWO_PI;
696  }
697  if (*Longitude > TWO_PI)
698  *Longitude -= TWO_PI;
699  if (*Longitude < -M_PI)
700  *Longitude += TWO_PI;
701 
702 // if (fabs(dlam) > (9.0 * M_PI / 180))
703 // { /* Distortion will result if Longitude is more than 9 degrees from the Central Meridian */
704 // Error_Code |= TRANMERC_LON_WARNING;
705 // }
706  }
707  return (Error_Code);
708 } /* END OF Convert_Transverse_Mercator_To_Geodetic */
709 
710 std::ostream& ossimTransMercatorProjection::print(std::ostream& out) const
711 {
712  out << setiosflags(ios::fixed) << setprecision(15)
713  << "// ossimTransMercatorProjection::print\n"
714  << ossimKeywordNames::SCALE_FACTOR_KW << ": " << TranMerc_Scale_Factor
715  << endl;
716  return ossimMapProjection::print(out);
717 }
718 
719 //*************************************************************************************************
721 //*************************************************************************************************
722 bool ossimTransMercatorProjection::operator==(const ossimProjection& proj) const
723 {
724  if (!ossimMapProjection::operator==(proj))
725  return false;
726 
727  const ossimTransMercatorProjection* p =
728  dynamic_cast<const ossimTransMercatorProjection*>(&proj);
729  if (!p) return false;
730 
731  if (!ossim::almostEqual(TranMerc_Scale_Factor,p->TranMerc_Scale_Factor)) return false;
732 
733  return true;
734 }
ossimMapProjection::loadState
virtual bool loadState(const ossimKeywordlist &kwl, const char *prefix=0)
Method to the load (recreate) the state of an object from a keyword list.
Definition: ossimMapProjection.cpp:923
ossimTransMercatorProjection::Get_Transverse_Mercator_Parameters
void Get_Transverse_Mercator_Parameters(double *a, double *f, double *Origin_Latitude, double *Central_Meridian, double *False_Easting, double *False_Northing, double *Scale_Factor) const
Definition: ossimTransMercatorProjection.cpp:355
DEG_PER_RAD
#define DEG_PER_RAD
Definition: ossimConstants.h:452
ossimKeywordlist
Represents serializable keyword/value map.
Definition: ossimKeywordlist.h:50
ossimKeywordlist::find
const char * find(const char *key) const
Definition: ossimKeywordlist.cpp:599
ossim::almostEqual
bool almostEqual(T x, T y, T tolerance=FLT_EPSILON)
Definition: ossimCommon.h:53
ossimDpt::y
double y
Definition: ossimDpt.h:165
ossimDatum::code
virtual const ossimString & code() const
Definition: ossimDatum.h:57
ossimTransMercatorProjection::setFalseEastingNorthing
void setFalseEastingNorthing(double falseEasting, double falseNorthing)
Definition: ossimTransMercatorProjection.cpp:133
ossimTransMercatorProjection::ossimTransMercatorProjection
ossimTransMercatorProjection(const ossimEllipsoid &ellipsoid=ossimEllipsoid(6378137, 6356752.3142), const ossimGpt &origin=ossimGpt())
Definition: ossimTransMercatorProjection.cpp:53
ossimKeywordNames::TYPE_KW
static const char * TYPE_KW
Definition: ossimKeywordNames.h:190
ossimGpt::changeDatum
void changeDatum(const ossimDatum *datum)
This will actually perform a shift.
Definition: ossimGpt.cpp:316
ossimTransMercatorProjection::saveState
virtual bool saveState(ossimKeywordlist &kwl, const char *prefix=0) const
Definition: ossimTransMercatorProjection.cpp:204
STATIC_TYPE_NAME
#define STATIC_TYPE_NAME(T)
Definition: ossimRtti.h:325
ossimGpt::datum
const ossimDatum * datum() const
datum().
Definition: ossimGpt.h:196
SPHSN
#define SPHSN(Latitude)
Definition: ossimTransMercatorProjection.cpp:45
ossimKeywordlist::add
void add(const char *prefix, const ossimKeywordlist &kwl, bool overwrite=true)
Definition: ossimKeywordlist.cpp:177
SPHTMD
#define SPHTMD(Latitude)
Definition: ossimTransMercatorProjection.cpp:41
M_PI
#define M_PI
Definition: ossimConstants.h:446
ossimTransMercatorProjection::Convert_Transverse_Mercator_To_Geodetic
long Convert_Transverse_Mercator_To_Geodetic(double Easting, double Northing, double *Latitude, double *Longitude) const
Definition: ossimTransMercatorProjection.cpp:547
ossimEllipsoid::getA
const double & getA() const
Definition: ossimEllipsoid.h:66
ossimGpt::lonr
double lonr() const
Returns the longitude in radian measure.
Definition: ossimGpt.h:76
ossimString::toDouble
double toDouble() const
Definition: ossimString.cpp:870
ossimMapProjection::saveState
virtual bool saveState(ossimKeywordlist &kwl, const char *prefix=0) const
Method to save the state of an object to a keyword list.
Definition: ossimMapProjection.cpp:811
ossimTransMercatorProjection::setFalseNorthing
void setFalseNorthing(double falseNorthing)
Definition: ossimTransMercatorProjection.cpp:126
MAX_DELTA_LONG
#define MAX_DELTA_LONG
Definition: ossimTransMercatorProjection.cpp:37
ossimTransMercatorProjection::print
virtual std::ostream & print(std::ostream &out) const
Definition: ossimTransMercatorProjection.cpp:710
SPHSR
#define SPHSR(Latitude)
Definition: ossimTransMercatorProjection.cpp:48
TRANMERC_NO_ERROR
#define TRANMERC_NO_ERROR
Definition: ossimTransMercatorProjection.cpp:23
ossimTransMercatorProjection::operator==
virtual bool operator==(const ossimProjection &projection) const
Returns TRUE if principal parameters are within epsilon tolerance.
Definition: ossimTransMercatorProjection.cpp:722
ossimDpt::x
double x
Definition: ossimDpt.h:164
TWO_PI
#define TWO_PI
Definition: ossimConstants.h:449
ossimGpt::latr
double latr() const
latr().
Definition: ossimGpt.h:66
ossimTransMercatorProjection::Convert_Geodetic_To_Transverse_Mercator
long Convert_Geodetic_To_Transverse_Mercator(double Latitude, double Longitude, double *Easting, double *Northing) const
Definition: ossimTransMercatorProjection.cpp:392
ossimEllipsoid::getFlattening
const double & getFlattening() const
Definition: ossimEllipsoid.h:68
ossimTransMercatorProjection::setParameters
void setParameters(double falseEasting, double falseNorthing, double scaleFactor)
Definition: ossimTransMercatorProjection.cpp:149
ossimMapProjection::theEllipsoid
ossimEllipsoid theEllipsoid
This method verifies that the projection parameters match the current pcs code.
Definition: ossimMapProjection.h:312
ossimKeywordNames::SCALE_FACTOR_KW
static const char * SCALE_FACTOR_KW
Definition: ossimKeywordNames.h:165
RTTI_DEF1
#define RTTI_DEF1(cls, name, b1)
Definition: ossimRtti.h:485
ossimMapProjection::theFalseEastingNorthing
ossimDpt theFalseEastingNorthing
Hold the false easting northing.
Definition: ossimMapProjection.h:349
ossimTransMercatorProjection::Set_Transverse_Mercator_Parameters
long Set_Transverse_Mercator_Parameters(double a, double f, double Origin_Latitude, double Central_Meridian, double False_Easting, double False_Northing, double Scale_Factor)
Definition: ossimTransMercatorProjection.cpp:248
ossimTransMercatorProjection::inverse
virtual ossimGpt inverse(const ossimDpt &eastingNorthing) const
Will take a point in meters and convert it to ground.
Definition: ossimTransMercatorProjection.cpp:169
ossimMapProjection::print
virtual std::ostream & print(std::ostream &out) const
Prints data members to stream.
Definition: ossimMapProjection.cpp:1366
ossimTransMercatorProjection::forward
virtual ossimDpt forward(const ossimGpt &latLon) const
All map projections will convert the world coordinate to an easting northing (Meters).
Definition: ossimTransMercatorProjection.cpp:182
RAD_PER_DEG
#define RAD_PER_DEG
Definition: ossimConstants.h:458
ossimTransMercatorProjection::loadState
virtual bool loadState(const ossimKeywordlist &kwl, const char *prefix=0)
Definition: ossimTransMercatorProjection.cpp:214
ossim::ostream
std::basic_ostream< char > ostream
Base class for char output streams.
Definition: ossimIosFwd.h:23
ossimMapProjection::theDatum
const ossimDatum * theDatum
This is only set if we want to have built in datum shifting.
Definition: ossimMapProjection.h:318
Generated on Fri Aug 3 2018 08:46:54 for OSSIM - Open Source Software Image Map by   doxygen 1.8.14