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 Bonne projection code.
11 //*******************************************************************
12 // $Id: ossimBonneProjection.cpp 17815 2010-08-03 13:23:14Z dburken $
13 #include <ossim/projection/ossimBonneProjection.h>
14 #include <ossim/projection/ossimSinusoid.h>
15 #include <ossim/base/ossimKeywordNames.h>
16
17 RTTI_DEF1(ossimBonneProjection, "ossimBonneProjection", ossimMapProjection)
18
19
20 /***************************************************************************/
21 /*
22 * DEFINES
23 */
24 #define BONN_NO_ERROR 0x0000
25 #define BONN_LAT_ERROR 0x0001
26 #define BONN_LON_ERROR 0x0002
27 #define BONN_EASTING_ERROR 0x0004
28 #define BONN_NORTHING_ERROR 0x0008
29 #define BONN_ORIGIN_LAT_ERROR 0x0010
30 #define BONN_CENT_MER_ERROR 0x0020
31 #define BONN_A_ERROR 0x0040
32 #define BONN_B_ERROR 0x0080
33 #define BONN_A_LESS_B_ERROR 0x0100
34
35 #define PI_OVER_2 (M_PI / 2.0)
36 #define BONN_m(coslat,sinlat) (coslat/sqrt(1.0 - es2*sinlat*sinlat))
37 #define BONN_M(c0lat,c1s2lat,c2s4lat,c3s6lat) (Bonn_a*(c0lat-c1s2lat+c2s4lat-c3s6lat))
38 #define COEFF_TIMES_BONN_SIN(coeff,x,latit) (coeff*(sin(x * latit)))
39 #define FLOAT_EQ(x,v,epsilon) (((v - epsilon) < x) && (x < (v + epsilon)))
40
ossimBonneProjection(const ossimEllipsoid & ellipsoid,const ossimGpt & origin)41 ossimBonneProjection::ossimBonneProjection(const ossimEllipsoid& ellipsoid,
42 const ossimGpt& origin)
43 :ossimMapProjection(ellipsoid, origin)
44 {
45 setDefaults();
46 update();
47 }
48
ossimBonneProjection(const ossimEllipsoid & ellipsoid,const ossimGpt & origin,const double falseEasting,const double falseNorthing)49 ossimBonneProjection::ossimBonneProjection(const ossimEllipsoid& ellipsoid,
50 const ossimGpt& origin,
51 const double falseEasting,
52 const double falseNorthing)
53 :ossimMapProjection(ellipsoid, origin)
54 {
55 Bonn_False_Easting = falseEasting;
56 Bonn_False_Northing = falseNorthing;
57 Bonn_Delta_Northing = 25000000.0;
58
59 update();
60 }
61
update()62 void ossimBonneProjection::update()
63 {
64 Set_Bonne_Parameters(theEllipsoid.getA(),
65 theEllipsoid.getFlattening(),
66 theOrigin.latr(),
67 theOrigin.lonr(),
68 Bonn_False_Easting,
69 Bonn_False_Northing);
70
71 theFalseEastingNorthing.x = Bonn_False_Easting;
72 theFalseEastingNorthing.y = Bonn_False_Northing;
73
74 ossimMapProjection::update();
75 }
76
setFalseEasting(double falseEasting)77 void ossimBonneProjection::setFalseEasting(double falseEasting)
78 {
79 Bonn_False_Easting = falseEasting;
80 update();
81 }
82
setFalseNorthing(double falseNorthing)83 void ossimBonneProjection::setFalseNorthing(double falseNorthing)
84 {
85 Bonn_False_Northing = falseNorthing;
86 update();
87 }
88
setFalseEastingNorthing(double falseEasting,double falseNorthing)89 void ossimBonneProjection::setFalseEastingNorthing(double falseEasting,
90 double falseNorthing)
91 {
92 Bonn_False_Easting = falseEasting;
93 Bonn_False_Northing = falseNorthing;
94 update();
95 }
96
setDefaults()97 void ossimBonneProjection::setDefaults()
98 {
99
100 Bonn_False_Easting = 0.0;
101 Bonn_False_Northing = 0.0;
102 Bonn_Delta_Northing = 25000000.0;
103 if(theOrigin.latd() == 0.0)
104 {
105 // we can't have the origin of lat 0 for Bonne
106 // so bump it up an arc second.
107 //
108 theOrigin.latd(1.0/3600.0);
109 }
110 }
111
inverse(const ossimDpt & eastingNorthing) const112 ossimGpt ossimBonneProjection::inverse(const ossimDpt &eastingNorthing)const
113 {
114 double lat, lon;
115
116 Convert_Bonne_To_Geodetic(eastingNorthing.x,
117 eastingNorthing.y,
118 &lat,
119 &lon);
120 return ossimGpt(lat*DEG_PER_RAD, lon*DEG_PER_RAD, 0, theDatum);
121 }
122
forward(const ossimGpt & latLon) const123 ossimDpt ossimBonneProjection::forward(const ossimGpt &latLon)const
124 {
125 double easting = 0.0;
126 double northing = 0.0;
127 ossimGpt gpt = latLon;
128
129 if (theDatum)
130 {
131 if (theDatum->code() != latLon.datum()->code())
132 {
133 gpt.changeDatum(theDatum); // Shift to our datum.
134 }
135 }
136
137 Convert_Geodetic_To_Bonne(gpt.latr(),
138 gpt.lonr(),
139 &easting,
140 &northing);
141 return ossimDpt(easting, northing);
142 }
143
saveState(ossimKeywordlist & kwl,const char * prefix) const144 bool ossimBonneProjection::saveState(ossimKeywordlist& kwl, const char* prefix) const
145 {
146 return ossimMapProjection::saveState(kwl, prefix);
147 }
148
loadState(const ossimKeywordlist & kwl,const char * prefix)149 bool ossimBonneProjection::loadState(const ossimKeywordlist& kwl, const char* prefix)
150 {
151 // Must do this first.
152 bool flag = ossimMapProjection::loadState(kwl, prefix);
153
154 const char* type = kwl.find(prefix, ossimKeywordNames::TYPE_KW);
155
156 setDefaults();
157
158 if(ossimString(type) == STATIC_TYPE_NAME(ossimBonneProjection))
159 {
160 Bonn_False_Easting = theFalseEastingNorthing.x;
161 Bonn_False_Northing = theFalseEastingNorthing.y;
162 }
163
164 update();
165 return flag;
166 }
167
168
169 /*
170 * These state variables are for optimization purposes. The only function
171 * that should modify them is Set_Bonne_Parameters.
172 */
173
174
175 /***************************************************************************/
176 /*
177 * FUNCTIONS
178 */
179
180
Set_Bonne_Parameters(double a,double f,double Origin_Latitude,double Central_Meridian,double False_Easting,double False_Northing)181 long ossimBonneProjection::Set_Bonne_Parameters(double a,
182 double f,
183 double Origin_Latitude,
184 double Central_Meridian,
185 double False_Easting,
186 double False_Northing)
187 { /* Begin Set_Bonne_Parameters */
188 /*
189 * The function Set_Bonne_Parameters receives the ellipsoid parameters and
190 * Bonne projection parameters as inputs, and sets the corresponding state
191 * variables. If any errors occur, the error code(s) are returned by the
192 * function, otherwise BONN_NO_ERROR is returned.
193 *
194 * a : Semi-major axis of ellipsoid, in meters (input)
195 * f : Flattening of ellipsoid (input)
196 * Origin_Latitude : Latitude in radians at which the (input)
197 * point scale factor is 1.0
198 * Central_Meridian : Longitude in radians at the center of (input)
199 * the projection
200 * False_Easting : A coordinate value in meters assigned to the
201 * central meridian of the projection. (input)
202 * False_Northing : A coordinate value in meters assigned to the
203 * origin latitude of the projection (input)
204 */
205
206 double j, three_es4;
207 double x,e1,e2,e3,e4;
208 double clat;
209 double sin2lat, sin4lat, sin6lat, lat;
210 // double inv_f = 1 / f;
211 long Error_Code = BONN_NO_ERROR;
212
213 // if (a <= 0.0)
214 // { /* Semi-major axis must be greater than zero */
215 // Error_Code |= BONN_A_ERROR;
216 // }
217 // if ((inv_f < 250) || (inv_f > 350))
218 // { /* Inverse flattening must be between 250 and 350 */
219 // Error_Code |= BONN_INV_F_ERROR;
220 // }
221 // if ((Origin_Latitude < -PI_OVER_2) || (Origin_Latitude > PI_OVER_2))
222 // { /* origin latitude out of range */
223 // Error_Code |= BONN_ORIGIN_LAT_ERROR;
224 // }
225 // if ((Central_Meridian < -PI) || (Central_Meridian > TWO_PI))
226 // { /* origin longitude out of range */
227 // Error_Code |= BONN_CENT_MER_ERROR;
228 // }
229 if (!Error_Code)
230 { /* no errors */
231 Bonn_a = a;
232 Bonn_f = f;
233 Bonn_Origin_Lat = Origin_Latitude;
234 // if (Central_Meridian > PI)
235 // Central_Meridian -= TWO_PI;
236 Bonn_Origin_Long = Central_Meridian;
237 Bonn_False_Northing = False_Northing;
238 Bonn_False_Easting = False_Easting;
239 if (Bonn_Origin_Lat == 0.0)
240 {
241 if (Bonn_Origin_Long > 0)
242 {
243 Bonn_Max_Easting = 19926189.0;
244 Bonn_Min_Easting = -20037509.0;
245 }
246 else if (Bonn_Origin_Long < 0)
247 {
248 Bonn_Max_Easting = 20037509.0;
249 Bonn_Min_Easting = -19926189.0;
250 }
251 else
252 {
253 Bonn_Max_Easting = 20037509.0;
254 Bonn_Min_Easting = -20037509.0;
255 }
256 Bonn_Delta_Northing = 10001966.0;
257 Set_Sinusoidal_Parameters(Bonn_a, Bonn_f, Bonn_Origin_Long, Bonn_False_Easting, Bonn_False_Northing);
258 }
259 else
260 {
261 Sin_Bonn_Origin_Lat = sin(Bonn_Origin_Lat);
262
263 es2 = 2 * Bonn_f - Bonn_f * Bonn_f;
264 es4 = es2 * es2;
265 es6 = es4 * es2;
266 j = 45.0 * es6 / 1024.0;
267 three_es4 = 3.0 * es4;
268 c0 = 1 - es2 / 4.0 - three_es4 / 64.0 - 5.0 * es6 / 256.0;
269 c1 = 3.0 * es2 / 8.0 + three_es4 / 32.0 + j;
270 c2 = 15.0 * es4 / 256.0 + j;
271 c3 = 35.0 * es6 / 3072.0;
272
273 clat = cos(Bonn_Origin_Lat);
274 m1 = BONN_m(clat, Sin_Bonn_Origin_Lat);
275
276 lat = c0 * Bonn_Origin_Lat;
277 sin2lat = COEFF_TIMES_BONN_SIN(c1, 2.0, Bonn_Origin_Lat);
278 sin4lat = COEFF_TIMES_BONN_SIN(c2, 4.0, Bonn_Origin_Lat);
279 sin6lat = COEFF_TIMES_BONN_SIN(c3, 6.0, Bonn_Origin_Lat);
280 M1 = BONN_M(lat, sin2lat, sin4lat, sin6lat);
281
282 x = sqrt (1.0 - es2);
283 e1 = (1.0 - x) / (1.0 + x);
284 e2 = e1 * e1;
285 e3 = e2 * e1;
286 e4 = e3 * e1;
287 a0 = 3.0 * e1 / 2.0 - 27.0 * e3 / 32.0;
288 a1 = 21.0 * e2 / 16.0 - 55.0 * e4 / 32.0;
289 a2 = 151.0 * e3 / 96.0;
290 a3 = 1097.0 * e4 / 512.0;
291 if (Sin_Bonn_Origin_Lat == 0.0)
292 Bonn_am1sin = 0.0;
293 else
294 Bonn_am1sin = Bonn_a * m1 / Sin_Bonn_Origin_Lat;
295
296 Bonn_Max_Easting = 20027474.0;
297 Bonn_Min_Easting = -20027474.0;
298 Bonn_Delta_Northing = 20003932.0;
299
300 }
301
302 } /* End if(!Error_Code) */
303 return (Error_Code);
304 } /* End Set_Bonne_Parameters */
305
306
Get_Bonne_Parameters(double * a,double * f,double * Origin_Latitude,double * Central_Meridian,double * False_Easting,double * False_Northing) const307 void ossimBonneProjection::Get_Bonne_Parameters(double *a,
308 double *f,
309 double *Origin_Latitude,
310 double *Central_Meridian,
311 double *False_Easting,
312 double *False_Northing)const
313 { /* Begin Get_Bonne_Parameters */
314 /*
315 * The function Get_Bonne_Parameters returns the current ellipsoid
316 * parameters and Bonne projection parameters.
317 *
318 * a : Semi-major axis of ellipsoid, in meters (output)
319 * f : Flattening of ellipsoid (output)
320 * Origin_Latitude : Latitude in radians at which the (output)
321 * point scale factor is 1.0
322 * Central_Meridian : Longitude in radians at the center of (output)
323 * the projection
324 * False_Easting : A coordinate value in meters assigned to the
325 * central meridian of the projection. (output)
326 * False_Northing : A coordinate value in meters assigned to the
327 * origin latitude of the projection (output)
328 */
329
330 *a = Bonn_a;
331 *f = Bonn_f;
332 *Origin_Latitude = Bonn_Origin_Lat;
333 *Central_Meridian = Bonn_Origin_Long;
334 *False_Easting = Bonn_False_Easting;
335 *False_Northing = Bonn_False_Northing;
336 return;
337 } /* End Get_Bonne_Parameters */
338
339
Convert_Geodetic_To_Bonne(double Latitude,double Longitude,double * Easting,double * Northing) const340 long ossimBonneProjection::Convert_Geodetic_To_Bonne (double Latitude,
341 double Longitude,
342 double *Easting,
343 double *Northing)const
344 { /* Begin Convert_Geodetic_To_Bonne */
345 /*
346 * The function Convert_Geodetic_To_Bonne converts geodetic (latitude and
347 * longitude) coordinates to Bonne projection (easting and northing)
348 * coordinates, according to the current ellipsoid and Bonne projection
349 * parameters. If any errors occur, the error code(s) are returned by the
350 * function, otherwise BONN_NO_ERROR is returned.
351 *
352 * Latitude : Latitude (phi) in radians (input)
353 * Longitude : Longitude (lambda) in radians (input)
354 * Easting : Easting (X) in meters (output)
355 * Northing : Northing (Y) in meters (output)
356 */
357
358 double dlam; /* Longitude - Central Meridan */
359 double mm;
360 double MM;
361 double rho;
362 double EE;
363 double clat = cos(Latitude);
364 double slat = sin(Latitude);
365 double lat, sin2lat, sin4lat, sin6lat;
366 long Error_Code = BONN_NO_ERROR;
367
368 // if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
369 // { /* Latitude out of range */
370 // Error_Code |= BONN_LAT_ERROR;
371 // }
372 // if ((Longitude < -PI) || (Longitude > TWO_PI))
373 // { /* Longitude out of range */
374 // Error_Code |= BONN_LON_ERROR;
375 // }
376 if (!Error_Code)
377 { /* no errors */
378 if (Bonn_Origin_Lat == 0.0)
379 Convert_Geodetic_To_Sinusoidal(Latitude, Longitude, Easting, Northing);
380 else
381 {
382 dlam = Longitude - Bonn_Origin_Long;
383 if (dlam > M_PI)
384 {
385 dlam -= TWO_PI;
386 }
387 if (dlam < -M_PI)
388 {
389 dlam += TWO_PI;
390 }
391 if ((Latitude - Bonn_Origin_Lat) == 0.0 && FLOAT_EQ(fabs(Latitude),PI_OVER_2,.00001))
392 {
393 *Easting = 0.0;
394 *Northing = 0.0;
395 }
396 else
397 {
398 mm = BONN_m(clat, slat);
399
400 lat = c0 * Latitude;
401 sin2lat = COEFF_TIMES_BONN_SIN(c1, 2.0, Latitude);
402 sin4lat = COEFF_TIMES_BONN_SIN(c2, 4.0, Latitude);
403 sin6lat = COEFF_TIMES_BONN_SIN(c3, 6.0, Latitude);
404 MM = BONN_M(lat, sin2lat, sin4lat, sin6lat);
405
406 rho = Bonn_am1sin + M1 - MM;
407 if (rho == 0)
408 EE = 0;
409 else
410 EE = Bonn_a * mm * dlam / rho;
411 *Easting = rho * sin(EE) + Bonn_False_Easting;
412 *Northing = Bonn_am1sin - rho * cos(EE) + Bonn_False_Northing;
413 }
414 }
415 }
416 return (Error_Code);
417 } /* End Convert_Geodetic_To_Bonne */
418
419
Convert_Bonne_To_Geodetic(double Easting,double Northing,double * Latitude,double * Longitude) const420 long ossimBonneProjection::Convert_Bonne_To_Geodetic(double Easting,
421 double Northing,
422 double *Latitude,
423 double *Longitude)const
424 { /* Begin Convert_Bonne_To_Geodetic */
425 /*
426 * The function Convert_Bonne_To_Geodetic converts Bonne projection
427 * (easting and northing) coordinates to geodetic (latitude and longitude)
428 * coordinates, according to the current ellipsoid and Bonne projection
429 * coordinates. If any errors occur, the error code(s) are returned by the
430 * function, otherwise BONN_NO_ERROR is returned.
431 *
432 * Easting : Easting (X) in meters (input)
433 * Northing : Northing (Y) in meters (input)
434 * Latitude : Latitude (phi) in radians (output)
435 * Longitude : Longitude (lambda) in radians (output)
436 */
437
438 double dx; /* Delta easting - Difference in easting (easting-FE) */
439 double dy; /* Delta northing - Difference in northing (northing-FN) */
440 double mu;
441 double MM;
442 double mm;
443 double am1sin_dy;
444 double rho;
445 double sin2mu, sin4mu, sin6mu, sin8mu;
446 double clat, slat;
447 long Error_Code = BONN_NO_ERROR;
448
449 // if ((Easting < (Bonn_False_Easting + Bonn_Min_Easting))
450 // || (Easting > (Bonn_False_Easting + Bonn_Max_Easting)))
451 // { /* Easting out of range */
452 // Error_Code |= BONN_EASTING_ERROR;
453 // }
454 // if ((Northing < (Bonn_False_Northing - Bonn_Delta_Northing))
455 // || (Northing > (Bonn_False_Northing + Bonn_Delta_Northing)))
456 // { /* Northing out of range */
457 // Error_Code |= BONN_NORTHING_ERROR;
458 // }
459 if (!Error_Code)
460 { /* no errors */
461 if (Bonn_Origin_Lat == 0.0)
462 Convert_Sinusoidal_To_Geodetic(Easting, Northing, Latitude, Longitude);
463 else
464 {
465 dy = Northing - Bonn_False_Northing;
466 dx = Easting - Bonn_False_Easting;
467 am1sin_dy = Bonn_am1sin - dy;
468 rho = sqrt(dx * dx + am1sin_dy * am1sin_dy);
469 if (Bonn_Origin_Lat < 0.0)
470 rho = -rho;
471 MM = Bonn_am1sin + M1 - rho;
472
473 mu = MM / (Bonn_a * c0);
474 sin2mu = COEFF_TIMES_BONN_SIN(a0, 2.0, mu);
475 sin4mu = COEFF_TIMES_BONN_SIN(a1, 4.0, mu);
476 sin6mu = COEFF_TIMES_BONN_SIN(a2, 6.0, mu);
477 sin8mu = COEFF_TIMES_BONN_SIN(a3, 8.0, mu);
478 *Latitude = mu + sin2mu + sin4mu + sin6mu + sin8mu;
479
480 if (FLOAT_EQ(fabs(*Latitude),PI_OVER_2,.00001))
481 {
482 *Longitude = Bonn_Origin_Long;
483 }
484 else
485 {
486 clat = cos(*Latitude);
487 slat = sin(*Latitude);
488 mm = BONN_m(clat, slat);
489
490 if (Bonn_Origin_Lat < 0.0)
491 {
492 dx = -dx;
493 am1sin_dy = -am1sin_dy;
494 }
495 *Longitude = Bonn_Origin_Long + rho * (atan2(dx, am1sin_dy)) /
496 (Bonn_a * mm);
497 }
498
499 if (*Latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
500 *Latitude = PI_OVER_2;
501 else if (*Latitude < -PI_OVER_2)
502 *Latitude = -PI_OVER_2;
503
504 // if (*Longitude > PI)
505 // *Longitude -= TWO_PI;
506 // if (*Longitude < -PI)
507 // *Longitude += TWO_PI;
508
509 if (*Longitude > M_PI)/* force distorted values to 180, -180 degrees */
510 *Longitude = M_PI;
511 else if (*Longitude < -M_PI)
512 *Longitude = -M_PI;
513 }
514 }
515 return (Error_Code);
516 } /* End Convert_Bonne_To_Geodetic */
517
518 //*************************************************************************************************
519 //! Returns TRUE if principal parameters are within epsilon tolerance.
520 //*************************************************************************************************
operator ==(const ossimProjection & proj) const521 bool ossimBonneProjection::operator==(const ossimProjection& proj) const
522 {
523 if (!ossimMapProjection::operator==(proj))
524 return false;
525
526 const ossimBonneProjection* p = dynamic_cast<const ossimBonneProjection*>(&proj);
527 if (!p)
528 return false;
529
530 if (!ossim::almostEqual(Bonn_False_Easting, p->Bonn_False_Easting) ) return false;
531 if (!ossim::almostEqual(Bonn_False_Northing,p->Bonn_False_Northing)) return false;
532 if (!ossim::almostEqual(Bonn_Delta_Northing,p->Bonn_Delta_Northing)) return false;
533
534 return true;
535 }
536