Isis 3 Programmer Reference
Robinson.cpp
Go to the documentation of this file.
1 
23 #include "Robinson.h"
24 
25 #include <cmath>
26 #include <cfloat>
27 
28 #include "Constants.h"
29 #include "IException.h"
30 #include "Projection.h"
31 #include "Pvl.h"
32 #include "PvlGroup.h"
33 #include "PvlKeyword.h"
34 
35 const double EPSILON = 1.0e-10;
36 
37 using namespace std;
38 namespace Isis {
39 
56  Robinson::Robinson(Pvl &label, bool allowDefaults) :
57  TProjection::TProjection(label) {
58  try {
59 
60  // Initialize coefficients
61  // PR Add xtra element at beginning to mimic 1-based arrays to match Snyder's code for
62  // easier debugging.
63  m_pr << 0 << -.062 <<
64  0 <<
65  0.062 <<
66  0.124 <<
67  0.186 <<
68  0.248 <<
69  0.310 <<
70  0.372 <<
71  0.434 <<
72  0.4958 <<
73  0.5571 <<
74  0.6176 <<
75  0.6769 <<
76  0.7346 <<
77  0.7903 <<
78  0.8435 <<
79  0.8936 <<
80  0.9394 <<
81  0.9761 <<
82  1.;
83  // XLR
84  m_xlr << 0 << 0.9986 <<
85  1.0 <<
86  0.9986 <<
87  0.9954 <<
88  0.99 <<
89  0.9822 <<
90  0.973 <<
91  0.96 <<
92  0.9427 <<
93  0.9216 <<
94  0.8962 <<
95  0.8679 <<
96  0.835 <<
97  0.7986 <<
98  0.7597 <<
99  0.7186 <<
100  0.6732 <<
101  0.6213 <<
102  0.5722 <<
103  0.5322;
104 
105  // Try to read the mapping group
106  PvlGroup &mapGroup = label.findGroup("Mapping", Pvl::Traverse);
107 
108  // Compute and write the default center longitude if allowed and
109  // necessary
110  if ((allowDefaults) && (!mapGroup.hasKeyword("CenterLongitude"))) {
111  double lon = (m_minimumLongitude + m_maximumLongitude) / 2.0;
112  mapGroup += PvlKeyword("CenterLongitude", toString(lon));
113  }
114 
115  // Get the center longitude
116  m_centerLongitude = mapGroup["CenterLongitude"];
117 
118  // convert to radians, adjust for longitude direction
121  }
122  catch(IException &e) {
123  QString message = "Invalid label group [Mapping]";
124  throw IException(e, IException::Io, message, _FILEINFO_);
125  }
126  }
127 
130  }
131 
140  bool Robinson::operator== (const Projection &proj) {
141  if (!Projection::operator==(proj)) return false;
142  // dont do the below it is a recusive plunge
143  // if (Projection::operator!=(proj)) return false;
144  Robinson *robin = (Robinson *) &proj;
145  if (robin->m_centerLongitude != m_centerLongitude) return false;
146  return true;
147  }
148 
154  QString Robinson::Name() const {
155  return "Robinson";
156  }
157 
164  QString Robinson::Version() const {
165  return "1.0";
166  }
167 
180  bool Robinson::SetGround(const double lat, const double lon) {
181  // Convert to radians
182  m_latitude = lat;
183  m_longitude = lon;
184  double latRadians = lat * DEG2RAD;
185  double lonRadians = lon * DEG2RAD;
186  if (m_longitudeDirection == PositiveWest) lonRadians *= -1.0;
187 
188  // Compute the coordinate
189  double deltaLon = (lonRadians - m_centerLongitude);
190  double p2 = fabs(latRadians / 5.0 / DEG2RAD);
191  long ip1 = (long) (p2 - EPSILON);
192  if (ip1 > 17) {
193  return false;
194  }
195 
196  // Stirling's interpolation formula (using 2nd Diff.)
197  //
198  p2 -= (double) ip1;
199  double x = 0.8487 * m_equatorialRadius * (m_xlr[ip1 + 2] + p2 * (m_xlr[ip1 + 3] -
200  m_xlr[ip1 + 1]) / 2.0 +
201  p2 * p2 * (m_xlr[ip1 + 3] - 2.0 * m_xlr[ip1 + 2] +
202  m_xlr[ip1 + 1])/2.0) * deltaLon;
203 
204  double y = 1.3523 * m_equatorialRadius * (m_pr[ip1 + 2] + p2 * (m_pr[ip1 + 3] -
205  m_pr[ip1 +1]) / 2.0 + p2 * p2 * (m_pr[ip1 + 3] -
206  2.0 * m_pr[ip1 + 2] + m_pr[ip1 + 1]) / 2.0);
207  if (lat < 0) y *= -1.;
208 
209  SetComputedXY(x, y);
210  m_good = true;
211  return m_good;
212  }
213 
214 
228  bool Robinson::SetCoordinate(const double x, const double y) {
229  m_latitude = 0;
230  m_longitude = 0;
231  // Save the coordinate
232  SetXY(x, y);
233 
234  double yy = y / m_equatorialRadius / 1.3523;
235  double phid = yy * 90.0;
236  double p2 = fabs(phid / 5.0);
237  long ip1 = (long) (p2 - EPSILON);
238  if (ip1 == 0) ip1 = 1;
239  if (ip1 > 17) {
240  return false;
241  }
242 
243  // Stirling's interpolation formula as used in forward transformation is
244  // reversed for first estimation of LAT. from rectangular coordinates. LAT.
245  // is then adjusted by iteration until use of forward series provides correct
246  // value of Y within tolerance.
247  //
248  double u, v, t, c;
249  double y1;
250 
251  for (int i = 0;;) {
252  u = m_pr[ip1 + 3] - m_pr[ip1 + 1];
253  v = m_pr[ip1 + 3] - 2.0 * m_pr[ip1 + 2] + m_pr[ip1 + 1];
254  t = 2.0 * (fabs(yy) - m_pr[ip1 + 2]) / u;
255  c = v / u;
256  p2 = t * (1.0 - c * t * (1.0 - 2.0 * c * t));
257 
258  if ((p2 >= 0.0) || (ip1 == 1)) {
259  phid = (p2 + (double) ip1 ) * 5.0;
260  if (y < 0) phid *= -1;
261 
262  do {
263  p2 = fabs(phid / 5.0);
264  ip1 = (long) (p2 - EPSILON);
265  if (ip1 > 17) {
266  return false;
267  }
268  p2 -= (double) ip1;
269 
270  y1 = 1.3523 * m_equatorialRadius * (m_pr[ip1 +2] + p2 *(m_pr[ip1 + 3] -
271  m_pr[ip1 +1]) / 2.0 + p2 * p2 * (m_pr[ip1 + 3] -
272  2.0 * m_pr[ip1 + 2] + m_pr[ip1 + 1])/2.0);
273  if (y < 0) y1 *= -1.;
274 
275  phid -= (90.0 * (y1 - y) / m_equatorialRadius / 1.3523);
276  i++;
277  if (i > 75) {
278  return false;
279  }
280  } while (fabs(y1 - y) > .00001);
281  break;
282  }
283  else {
284  ip1--;
285  if (ip1 < 0) {
286  return false;
287  }
288  }
289  }
290 
291  // Compute latitude and make sure it is not above 90
292  m_latitude = phid;
293 
294  // Calculate longitude. using final latitude. with transposed forward Stirling's
295  // interpolation formula.
296  m_longitude = m_centerLongitude + x / m_equatorialRadius / 0.8487 / (m_xlr[ip1 + 2] +
297  p2 * (m_xlr[ip1 + 3] - m_xlr[ip1 + 1]) / 2.0 +
298  p2 * p2 * (m_xlr[ip1 + 3] - 2.0 * m_xlr[ip1 + 2] +
299  m_xlr[ip1 + 1]) / 2.0);
300  m_longitude *= RAD2DEG;
302 
303  // Our double precision is not good once we pass a certain magnitude of
304  // longitude. Prevent failures down the road by failing now.
305  m_good = (fabs(m_longitude) < 1E10);
306  return m_good;
307  }
308 
309 
333  bool Robinson::XYRange(double &minX, double &maxX,
334  double &minY, double &maxY) {
335 
336  // Check the corners of the lat/lon range
341 
342  // If the latitude crosses the equator check there
343  if ((m_minimumLatitude < 0.0) && (m_maximumLatitude > 0.0)) {
346  }
347 
348  // Make sure everything is ordered
349  if (m_minimumX >= m_maximumX) return false;
350  if (m_minimumY >= m_maximumY) return false;
351 
352  // Return X/Y min/maxs
353  minX = m_minimumX;
354  maxX = m_maximumX;
355  minY = m_minimumY;
356  maxY = m_maximumY;
357  return true;
358  }
359 
360 
367  PvlGroup mapping = TProjection::Mapping();
368 
369  mapping += m_mappingGrp["CenterLongitude"];
370 
371  return mapping;
372  }
373 
381 
382  return mapping;
383  }
384 
392 
393  mapping += m_mappingGrp["CenterLongitude"];
394 
395  return mapping;
396  }
397 
398 } // end namespace isis
399 
413  bool allowDefaults) {
414  return new Isis::Robinson(lab, allowDefaults);
415 }
bool SetCoordinate(const double x, const double y)
This method is used to set the projection x/y.
Definition: Robinson.cpp:228
bool hasKeyword(const QString &name) const
Check to see if a keyword exists.
PvlGroup MappingLatitudes()
This function returns the latitude keywords that this projection uses.
Definition: Robinson.cpp:379
Isis::Projection * RobinsonPlugin(Isis::Pvl &lab, bool allowDefaults)
This is the function that is called in order to instantiate a Robinson object.
Definition: Robinson.cpp:412
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition: PvlObject.h:141
Longitude values increase in the westerly direction.
Definition: TProjection.h:241
double m_centerLongitude
The center longitude for the map projection.
Definition: Robinson.h:89
Robinson Map Projection.
Definition: Robinson.h:68
Base class for Map TProjections.
Definition: TProjection.h:182
bool XYRange(double &minX, double &maxX, double &minY, double &maxY)
This method is used to determine the x/y range which completely covers the area of interest specified...
Definition: Robinson.cpp:333
Namespace for the standard library.
PvlGroup MappingLongitudes()
This function returns the longitude keywords that this projection uses.
Definition: Robinson.cpp:390
Search child objects.
Definition: PvlObject.h:170
double m_minimumX
The data elements m_minimumX, m_minimumY, m_maximumX, and m_maximumY are convience data elements when...
Definition: Projection.h:333
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition: IString.cpp:226
A type of error that occurred when performing an actual I/O operation.
Definition: IException.h:171
void SetComputedXY(double x, double y)
This protected method is a helper for derived classes.
Definition: Projection.cpp:795
double m_latitude
This contains the currently set latitude value.
Definition: TProjection.h:332
virtual PvlGroup MappingLatitudes()
This function returns the latitude keywords that this projection uses.
QString Version() const
Returns the version of the map projection.
Definition: Robinson.cpp:164
double m_maximumY
See minimumX description.
Definition: Projection.h:344
double m_longitude
This contains the currently set longitude value.
Definition: TProjection.h:334
double m_maximumLongitude
Contains the maximum longitude for the entire ground range.
Definition: TProjection.h:376
double m_minimumLatitude
Contains the minimum latitude for the entire ground range.
Definition: TProjection.h:370
double m_maximumLatitude
Contains the maximum latitude for the entire ground range.
Definition: TProjection.h:372
Base class for Map Projections.
Definition: Projection.h:171
double m_minimumY
See minimumX description.
Definition: Projection.h:343
double m_equatorialRadius
Polar radius of the target.
Definition: TProjection.h:351
Contains multiple PvlContainers.
Definition: PvlGroup.h:57
#define _FILEINFO_
Macro for the filename and line number.
Definition: IException.h:40
A single keyword-value pair.
Definition: PvlKeyword.h:98
bool SetGround(const double lat, const double lon)
This method is used to set the latitude/longitude (assumed to be of the correct LatitudeType, LongitudeDirection, and LongitudeDomain.
Definition: Robinson.cpp:180
Container for cube-like labels.
Definition: Pvl.h:135
~Robinson()
Destroys the Robinson object.
Definition: Robinson.cpp:129
virtual PvlGroup Mapping()
This function returns the keywords that this projection uses.
const double DEG2RAD
Multiplier for converting from degrees to radians.
Definition: Constants.h:59
bool m_good
Indicates if the contents of m_x, m_y, m_latitude, and m_longitude are valid.
Definition: Projection.h:316
LongitudeDirection m_longitudeDirection
An enumerated type indicating the LongitudeDirection read from the labels.
Definition: TProjection.h:340
Isis exception class.
Definition: IException.h:107
PvlGroup Mapping()
This function returns the keywords that this projection uses.
Definition: Robinson.cpp:366
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31
QString Name() const
Returns the name of the map projection, "Robinson".
Definition: Robinson.cpp:154
const double RAD2DEG
Multiplier for converting from radians to degrees.
Definition: Constants.h:60
virtual PvlGroup MappingLongitudes()
This function returns the longitude keywords that this projection uses.
void SetXY(double x, double y)
This protected method is a helper for derived classes.
Definition: Projection.cpp:819
void XYRangeCheck(const double latitude, const double longitude)
This convience function is established to assist in the development of the XYRange virtual method...
double m_minimumLongitude
Contains the minimum longitude for the entire ground range.
Definition: TProjection.h:374
double m_maximumX
See minimumX description.
Definition: Projection.h:342
bool operator==(const Projection &proj)
Compares two Projection objects to see if they are equal.
Definition: Robinson.cpp:140
PvlGroup m_mappingGrp
Mapping group that created this projection.
Definition: Projection.h:345