Isis 3 Programmer Reference
TProjection.cpp
1
6/* SPDX-License-Identifier: CC0-1.0 */
7#include "TProjection.h"
8
9#include <QObject>
10
11#include <cfloat>
12#include <cmath>
13#include <iomanip>
14#include <sstream>
15#include <vector>
16
17#include <SpiceUsr.h>
18
19#include "Constants.h"
20#include "Displacement.h"
21#include "FileName.h"
22#include "IException.h"
23#include "IString.h"
24#include "Longitude.h"
25#include "NaifStatus.h"
26#include "Pvl.h"
27#include "PvlGroup.h"
28#include "PvlKeyword.h"
29#include "SpecialPixel.h"
30#include "Target.h"
31#include "WorldMapper.h"
32
33using namespace std;
34namespace Isis {
91 try {
92 // Mapping group is read by the parent
93 // Get the radii from the EquatorialRadius and PolarRadius keywords
94 if ((m_mappingGrp.hasKeyword("EquatorialRadius")) &&
95 (m_mappingGrp.hasKeyword("PolarRadius"))) {
96 m_equatorialRadius = m_mappingGrp["EquatorialRadius"];
97 m_polarRadius = m_mappingGrp["PolarRadius"];
98 }
99 // Get the radii
100 try {
102 m_equatorialRadius = radii["EquatorialRadius"];
103 m_polarRadius = radii["PolarRadius"];
104 }
105 catch (IException &e) {
106 QString msg = "Projection failed. No target radii are available "
107 "through keywords [EquatorialRadius and PolarRadius] "
108 "or [TargetName].";
109 throw IException(e, IException::Unknown, msg, _FILEINFO_);
110 }
111
112 // Check the radii for validity
113 if (m_equatorialRadius <= 0.0) {
114 QString msg = "Projection failed. Invalid value for keyword "
115 "[EquatorialRadius]. It must be greater than zero";
116 throw IException(IException::Unknown, msg, _FILEINFO_);
117 }
118 if (m_polarRadius <= 0.0) {
119 QString msg = "Projection failed. Invalid value for keyword "
120 "[PolarRadius]. It must be greater than zero";
121 throw IException(IException::Unknown, msg, _FILEINFO_);
122 }
123
124 // Get the LatitudeType
125 if ((QString) m_mappingGrp["LatitudeType"] == "Planetographic") {
127 }
128 else if ((QString) m_mappingGrp["LatitudeType"] == "Planetocentric") {
130 }
131 else {
132 QString msg = "Projection failed. Invalid value for keyword "
133 "[LatitudeType] must be "
134 "[Planetographic or Planetocentric]";
135 throw IException(IException::Unknown, msg, _FILEINFO_);
136 }
137
138 // Get the LongitudeDirection
139 if ((QString) m_mappingGrp["LongitudeDirection"] == "PositiveWest") {
141 }
142 else if ((QString) m_mappingGrp["LongitudeDirection"] == "PositiveEast") {
144 }
145 else {
146 QString msg = "Projection failed. Invalid value for keyword "
147 "[LongitudeDirection] must be "
148 "[PositiveWest or PositiveEast]";
149 throw IException(IException::Unknown, msg, _FILEINFO_);
150 }
151
152 // Get the LongitudeDomain
153 if ((QString) m_mappingGrp["LongitudeDomain"] == "360") {
154 m_longitudeDomain = 360;
155 }
156 else if ((QString) m_mappingGrp["LongitudeDomain"] == "180") {
157 m_longitudeDomain = 180;
158 }
159 else {
160 QString msg = "Projection failed. Invalid value for keyword "
161 "[LongitudeDomain] must be [180 or 360]";
162 throw IException(IException::Unknown, msg, _FILEINFO_);
163 }
164
165 // Get the ground range if it exists
166 m_groundRangeGood = false;
167 if ((m_mappingGrp.hasKeyword("MinimumLatitude")) &&
168 (m_mappingGrp.hasKeyword("MaximumLatitude")) &&
169 (m_mappingGrp.hasKeyword("MinimumLongitude")) &&
170 (m_mappingGrp.hasKeyword("MaximumLongitude"))) {
171 m_minimumLatitude = m_mappingGrp["MinimumLatitude"];
172 m_maximumLatitude = m_mappingGrp["MaximumLatitude"];
173 m_minimumLongitude = m_mappingGrp["MinimumLongitude"];
174 m_maximumLongitude = m_mappingGrp["MaximumLongitude"];
175
176 if ((m_minimumLatitude < -90.0) || (m_minimumLatitude > 90.0)) {
177 QString msg = "Projection failed. "
178 "[MinimumLatitude] of [" + toString(m_minimumLatitude)
179 + "] is outside the range of [-90:90]";
180 throw IException(IException::Unknown, msg, _FILEINFO_);
181 }
182
183 if ((m_maximumLatitude < -90.0) || (m_maximumLatitude > 90.0)) {
184 QString msg = "Projection failed. "
185 "[MaximumLatitude] of [" + toString(m_maximumLatitude)
186 + "] is outside the range of [-90:90]";
187 throw IException(IException::Unknown, msg, _FILEINFO_);
188 }
189
191 QString msg = "Projection failed. "
192 "[MinimumLatitude,MaximumLatitude] of ["
194 + toString(m_maximumLatitude) + "] are not "
195 + "properly ordered";
196 throw IException(IException::Unknown, msg, _FILEINFO_);
197 }
198
200 QString msg = "Projection failed. "
201 "[MinimumLongitude,MaximumLongitude] of ["
203 + toString(m_maximumLongitude) + "] are not "
204 + "properly ordered";
205 throw IException(IException::Unknown, msg, _FILEINFO_);
206 }
207
208 m_groundRangeGood = true;
209 }
210 else {
211 // if no ground range is given, initialize the min/max lat/lon to 0
212 m_minimumLatitude = 0.0;
213 m_maximumLatitude = 0.0;
214 m_minimumLongitude = 0.0;
215 m_maximumLongitude = 0.0;
216 }
217
218 // Initialize miscellaneous protected data elements
220 QString msg = "Projection failed. Invalid keyword value(s). "
221 "[EquatorialRadius] = " + toString(m_equatorialRadius)
222 + " must be greater than or equal to [PolarRadius] = "
224 throw IException(IException::Unknown, msg, _FILEINFO_);
225 }
226 else {
227 m_eccentricity = 1.0 -
231 }
232
233 // initialize the rest of the x,y,lat,lon member variables
236
237 // If we made it to this point, we have what we need for a triaxial projection
239 }
240 catch (IException &e) {
241 QString msg = "Projection failed. Invalid label group [Mapping]";
242 throw IException(e, IException::Unknown, msg, _FILEINFO_);
243 }
244 }
245
249
261 if (!Projection::operator==(proj)) return false;
262 TProjection *tproj = (TProjection *) &proj;
263 if (EquatorialRadius() != tproj->EquatorialRadius()) return false;
264 if (PolarRadius() != tproj->PolarRadius()) return false;
265 if (IsPlanetocentric() != tproj->IsPlanetocentric()) return false;
266 if (IsPositiveWest() != tproj->IsPositiveWest()) return false;
267 return true;
268 }
269
270
278 return m_equatorialRadius;
279 }
280
288 return m_polarRadius;
289 }
290
305 return m_eccentricity;
306 }
307
326 double TProjection::LocalRadius(double latitude) const {
327 if (latitude == Null) {
329 "Unable to calculate local radius. The given latitude value ["
330 + toString(latitude) + "] is invalid.",
331 _FILEINFO_);
332 }
333 double a = m_equatorialRadius;
334 double c = m_polarRadius;
335 // to save calculations, if the target is spherical, return the eq. rad
336 if (a - c < DBL_EPSILON) {
337 return a;
338 }
339 else {
340 double lat = latitude * PI / 180.0;
341 return a * c / sqrt(pow(c * cos(lat), 2) + pow(a * sin(lat), 2));
342 }
343 }
344
354 return LocalRadius(m_latitude);
355 }
356
357
369 return 0.0;
370 }
371
382 return false;
383 }
384
395
406
418 double TProjection::ToPlanetocentric(const double lat) const {
420 }
421
434 double TProjection::ToPlanetocentric(const double lat,
435 double eRadius, double pRadius) {
436 if (lat == Null || abs(lat) > 90.0) {
438 "Unable to convert to Planetocentric. The given latitude value ["
439 + toString(lat) + "] is invalid.",
440 _FILEINFO_);
441 }
442 double mylat = lat;
443 if (abs(mylat) < 90.0) { // So tan doesn't fail
444 mylat *= PI / 180.0;
445 mylat = atan(tan(mylat) * (pRadius / eRadius) *
446 (pRadius / eRadius));
447 mylat *= 180.0 / PI;
448 }
449 return mylat;
450 }
451
463 double TProjection::ToPlanetographic(const double lat) const {
465 }
466
481 double eRadius, double pRadius) {
482 //Account for double rounding error.
483 if (qFuzzyCompare(fabs(lat), 90.0)) {
484 lat = qRound(lat);
485 }
486 if (lat == Null || fabs(lat) > 90.0) {
488 "Unable to convert to Planetographic. The given latitude value ["
489 + toString(lat) + "] is invalid.",
490 _FILEINFO_);
491 }
492 double mylat = lat;
493 if (fabs(mylat) < 90.0) { // So tan doesn't fail
494 mylat *= PI / 180.0;
495 mylat = atan(tan(mylat) * (eRadius / pRadius) *
496 (eRadius / pRadius));
497 mylat *= 180.0 / PI;
498 }
499 return mylat;
500 }
501
509 if (m_latitudeType == Planetographic) return "Planetographic";
510 return "Planetocentric";
511 }
512
523
534
548 double TProjection::ToPositiveEast(const double lon, const int domain) {
549 if (lon == Null) {
551 "Unable to convert to PositiveEast. The given longitude value ["
552 + toString(lon) + "] is invalid.",
553 _FILEINFO_);
554 }
555 double mylon = lon;
556
557 mylon *= -1;
558
559 if (domain == 360) {
560 mylon = To360Domain(mylon);
561 }
562 else if (domain == 180) {
563 mylon = To180Domain(mylon);
564 }
565 else {
566 QString msg = "Unable to convert longitude. Domain [" + toString(domain)
567 + "] is not 180 or 360.";
568 throw IException(IException::Unknown, msg, _FILEINFO_);
569 }
570
571 return mylon;
572 }
573
587 double TProjection::ToPositiveWest(const double lon, const int domain) {
588 if (lon == Null) {
590 "Unable to convert to PositiveWest. The given longitude value ["
591 + toString(lon) + "] is invalid.",
592 _FILEINFO_);
593 }
594 double mylon = lon;
595
596 mylon *= -1;
597
598 if (domain == 360) {
599 mylon = To360Domain(mylon);
600 }
601 else if (domain == 180) {
602 mylon = To180Domain(mylon);
603 }
604 else {
605 QString msg = "Unable to convert longitude. Domain [" + toString(domain)
606 + "] is not 180 or 360.";
607 throw IException(IException::Unknown, msg, _FILEINFO_);
608 }
609
610 return mylon;
611 }
612
621 if (m_longitudeDirection == PositiveEast) return "PositiveEast";
622 return "PositiveWest";
623 }
624
633 return m_longitudeDomain == 180;
634 }
635
644 return m_longitudeDomain == 360;
645 }
646
657 double TProjection::To180Domain(const double lon) {
658 if (lon == Null) {
660 "Unable to convert to 180 degree domain. The given longitude value ["
661 + toString(lon) + "] is invalid.",
662 _FILEINFO_);
663 }
665 }
666
675 double TProjection::To360Domain(const double lon) {
676 if (lon == Null) {
678 "Unable to convert to 360 degree domain. The given longitude value ["
679 + toString(lon) + "] is invalid.",
680 _FILEINFO_);
681 }
682 double result = lon;
683
684 if ( (lon < 0.0 || lon > 360.0) &&
685 !qFuzzyCompare(lon, 0.0) && !qFuzzyCompare(lon, 360.0)) {
687 }
688
689 return result;
690 }
691
699 if (m_longitudeDomain == 360) return "360";
700 return "180";
701 }
702
711 return m_minimumLatitude;
712 }
713
722 return m_maximumLatitude;
723 }
724
733 return m_minimumLongitude;
734 }
735
744 return m_maximumLongitude;
745 }
746
760 bool TProjection::SetGround(const double lat, const double lon) {
761 if (lat == Null || lon == Null) {
762 m_good = false;
763 return m_good;
764 }
765 else {
766 m_latitude = lat;
767 m_longitude = lon;
768 m_good = true;
769 SetComputedXY(lon, lat);
770 }
771 return m_good;
772 }
773
789 bool TProjection::SetCoordinate(const double x, const double y) {
790 if (x == Null || y == Null) {
791 m_good = false;
792 }
793 else {
794 m_good = true;
795 SetXY(x, y);
796 m_latitude = XCoord();
798 }
799 return m_good;
800 }
801
802
811 double TProjection::Latitude() const {
812 return m_latitude;
813 }
814
823 double TProjection::Longitude() const {
824 return m_longitude;
825 }
826
827
839 bool TProjection::SetUniversalGround(const double lat, const double lon) {
840 if (lat == Null || lon == Null) {
841 m_good = false;
842 return m_good;
843 }
844 // Deal with the longitude first
845 m_longitude = lon;
847 if (m_longitudeDomain == 180) {
849 }
850 else {
851 // Do this because longitudeDirection could cause (-360,0)
853 }
854
855 // Deal with the latitude
858 }
859 else {
860 m_latitude = lat;
861 }
862
863 // Now the lat/lon are in user defined coordinates so set them
865 }
866
867
879 bool TProjection::SetUnboundUniversalGround(const double lat, const double lon) {
880 if (lat == Null || lon == Null) {
881 m_good = false;
882 return m_good;
883 }
884 // Deal with the longitude first
885 m_longitude = lon;
887
888 // Deal with the latitude
891 }
892 else {
893 m_latitude = lat;
894 }
895
896 // Now the lat/lon are in user defined coordinates so set them
898 }
899
900
909 double lat = m_latitude;
911 return lat;
912 }
913
923 double lon = m_longitude;
924 if (m_longitudeDirection == PositiveWest) lon = -lon;
925 lon = To360Domain(lon);
926 return lon;
927 }
928
929
940 double TProjection::Scale() const {
941 if (m_mapper != NULL) {
942 double lat = TrueScaleLatitude() * PI / 180.0;
943 double a = m_polarRadius * cos(lat);
944 double b = m_equatorialRadius * sin(lat);
945 double localRadius = m_equatorialRadius * m_polarRadius /
946 sqrt(a * a + b * b);
947
948 return localRadius / m_mapper->Resolution();
949 }
950 else {
951 return 1.0;
952 }
953 }
954
955
993 bool TProjection::XYRange(double &minX, double &maxX,
994 double &minY, double &maxY) {
995 if (minX == Null || maxX == Null || minY == Null || maxY == Null) {
996 return false;
997 }
998 if (m_groundRangeGood) {
999 minX = m_minimumLongitude;
1000 maxX = m_maximumLongitude;
1001 minY = m_minimumLatitude;
1002 maxY = m_maximumLatitude;
1003 return true;
1004 }
1005 return false;
1006 }
1007
1062 void TProjection::XYRangeCheck(const double latitude, const double longitude) {
1063
1064 if (latitude == Null || longitude == Null) {
1065 m_good = false;
1066 return;
1067 }
1068
1069 SetGround(latitude, longitude);
1070 if (!IsGood()) return;
1071
1072 if (XCoord() < m_minimumX) m_minimumX = XCoord();
1073 if (XCoord() > m_maximumX) m_maximumX = XCoord();
1074 if (YCoord() < m_minimumY) m_minimumY = YCoord();
1075 if (YCoord() > m_maximumY) m_maximumY = YCoord();
1076 return;
1077 }
1078
1079
1098 double maxLon,
1099 double longitude) {
1100 // get the min/max range closest to 0.0 lon
1101 double adjustedLon = To360Domain(longitude);
1102 double adjustedMinLon = To360Domain(minLon);
1103 double adjustedMaxLon = To360Domain(maxLon);
1104
1105 if (adjustedMinLon > adjustedMaxLon) {
1106 if (adjustedLon > adjustedMinLon) {
1107 adjustedLon -= 360;
1108 }
1109 adjustedMinLon -= 360;
1110 }
1111
1112 // if this range covers all longitudes, then the given longitude is clearly in range
1113 if (qFuzzyCompare(maxLon - minLon, 360.0)) {
1114 return true;
1115 }
1116 else if (adjustedMinLon <= adjustedLon && adjustedLon <= adjustedMaxLon) {
1117 return true;
1118 }
1119 else {
1120 return false;
1121 }
1122 }
1123
1124
1140 bool TProjection::inLongitudeRange(double longitude) {
1141 return inLongitudeRange(MinimumLongitude(), MaximumLongitude(), longitude);
1142 }
1143
1144
1160 bool TProjection::inLatitudeRange(double latitude) {
1161 if (MaximumLatitude() - MinimumLatitude() == 180) {
1162 return true;
1163 }
1164 else if (MinimumLatitude() <= latitude && latitude <= MaximumLatitude()) {
1165 return true;
1166 }
1167 else {
1168 return false;
1169 }
1170 }
1171
1172
1195 bool TProjection::xyRangeOblique(double &minX, double &maxX,
1196 double &minY, double &maxY) {
1197 if (minX == Null || maxX == Null || minY == Null || maxY == Null) {
1198 return false;
1199 }
1200 //For oblique, we'll have to walk all 4 sides to find out min/max x/y values
1201 if (!HasGroundRange()) return false; // Don't have min/max lat/lon,
1202 //can't continue
1203
1204 m_specialLatCases.clear();
1205 m_specialLonCases.clear();
1206
1207 // First, search longitude for min X/Y
1208 double minFoundX1, minFoundX2;
1209 double minFoundY1, minFoundY2;
1210
1211 // Search for minX between minlat and maxlat along minlon
1213 minFoundX1, MinimumLongitude(), true, true, true);
1214 // Search for minX between minlat and maxlat along maxlon
1216 minFoundX2, MaximumLongitude(), true, true, true);
1217 // Search for minY between minlat and maxlat along minlon
1219 minFoundY1, MinimumLongitude(), false, true, true);
1220 // Search for minY between minlat and maxlat along maxlon
1222 minFoundY2, MaximumLongitude(), false, true, true);
1223
1224 // Second, search latitude for min X/Y
1225 double minFoundX3, minFoundX4;
1226 double minFoundY3, minFoundY4;
1227
1228 // Search for minX between minlon and maxlon along minlat
1230 minFoundX3, MinimumLatitude(), true, false, true);
1231 // Search for minX between minlon and maxlon along maxlat
1233 minFoundX4, MaximumLatitude(), true, false, true);
1234 // Search for minY between minlon and maxlon along minlat
1236 minFoundY3, MinimumLatitude(), false, false, true);
1237 // Search for minY between minlon and maxlon along maxlat
1239 minFoundY4, MaximumLatitude(), false, false, true);
1240
1241 // We've searched all possible minimums, go ahead and store the lowest
1242 double minFoundX5 = min(minFoundX1, minFoundX2);
1243 double minFoundX6 = min(minFoundX3, minFoundX4);
1244 m_minimumX = min(minFoundX5, minFoundX6);
1245
1246 double minFoundY5 = min(minFoundY1, minFoundY2);
1247 double minFoundY6 = min(minFoundY3, minFoundY4);
1248 m_minimumY = min(minFoundY5, minFoundY6);
1249
1250 // Search longitude for max X/Y
1251 double maxFoundX1, maxFoundX2;
1252 double maxFoundY1, maxFoundY2;
1253
1254 // Search for maxX between minlat and maxlat along minlon
1256 maxFoundX1, MinimumLongitude(), true, true, false);
1257 // Search for maxX between minlat and maxlat along maxlon
1259 maxFoundX2, MaximumLongitude(), true, true, false);
1260 // Search for maxY between minlat and maxlat along minlon
1262 maxFoundY1, MinimumLongitude(), false, true, false);
1263 // Search for maxY between minlat and maxlat along maxlon
1265 maxFoundY2, MaximumLongitude(), false, true, false);
1266
1267 // Search latitude for max X/Y
1268 double maxFoundX3, maxFoundX4;
1269 double maxFoundY3, maxFoundY4;
1270
1271 // Search for maxX between minlon and maxlon along minlat
1273 maxFoundX3, MinimumLatitude(), true, false, false);
1274 // Search for maxX between minlon and maxlon along maxlat
1276 maxFoundX4, MaximumLatitude(), true, false, false);
1277 // Search for maxY between minlon and maxlon along minlat
1279 maxFoundY3, MinimumLatitude(), false, false, false);
1280 // Search for maxY between minlon and maxlon along maxlat
1282 maxFoundY4, MaximumLatitude(), false, false, false);
1283
1284 // We've searched all possible maximums, go ahead and store the highest
1285 double maxFoundX5 = max(maxFoundX1, maxFoundX2);
1286 double maxFoundX6 = max(maxFoundX3, maxFoundX4);
1287 m_maximumX = max(maxFoundX5, maxFoundX6);
1288
1289 double maxFoundY5 = max(maxFoundY1, maxFoundY2);
1290 double maxFoundY6 = max(maxFoundY3, maxFoundY4);
1291 m_maximumY = max(maxFoundY5, maxFoundY6);
1292
1293 // Look along discontinuities for more extremes
1294 vector<double> specialLatCases = m_specialLatCases;
1295 for (unsigned int specialLatCase = 0;
1296 specialLatCase < specialLatCases.size();
1297 specialLatCase ++) {
1298 double minX, maxX, minY, maxY;
1299
1300 // Search for minX between minlon and maxlon along latitude discontinuities
1302 minX, specialLatCases[specialLatCase], true, false, true);
1303 // Search for minY between minlon and maxlon along latitude discontinuities
1305 minY, specialLatCases[specialLatCase], false, false, true);
1306 // Search for maxX between minlon and maxlon along latitude discontinuities
1308 maxX, specialLatCases[specialLatCase], true, false, false);
1309 // Search for maxX between minlon and maxlon along latitude discontinuities
1311 maxY, specialLatCases[specialLatCase], false, false, false);
1312
1313 m_minimumX = min(minX, m_minimumX);
1314 m_maximumX = max(maxX, m_maximumX);
1315 m_minimumY = min(minY, m_minimumY);
1316 m_maximumY = max(maxY, m_maximumY);
1317 }
1318
1319 vector<double> specialLonCases = m_specialLonCases;
1320 for (unsigned int specialLonCase = 0;
1321 specialLonCase < specialLonCases.size();
1322 specialLonCase ++) {
1323 double minX, maxX, minY, maxY;
1324
1325 // Search for minX between minlat and maxlat along longitude discontinuities
1327 minX, specialLonCases[specialLonCase], true, true, true);
1328 // Search for minY between minlat and maxlat along longitude discontinuities
1330 minY, specialLonCases[specialLonCase], false, true, true);
1331 // Search for maxX between minlat and maxlat along longitude discontinuities
1333 maxX, specialLonCases[specialLonCase], true, true, false);
1334 // Search for maxY between minlat and maxlat along longitude discontinuities
1336 maxY, specialLonCases[specialLonCase], false, true, false);
1337
1338 m_minimumX = min(minX, m_minimumX);
1339 m_maximumX = max(maxX, m_maximumX);
1340 m_minimumY = min(minY, m_minimumY);
1341 m_maximumY = max(maxY, m_maximumY);
1342 }
1343
1344 m_specialLatCases.clear();
1345 m_specialLonCases.clear();
1346
1347 // Make sure everything is ordered
1348 if (m_minimumX >= m_maximumX) return false;
1349 if (m_minimumY >= m_maximumY) return false;
1350
1351 // Return X/Y min/maxs
1352 minX = m_minimumX;
1353 maxX = m_maximumX;
1354 minY = m_minimumY;
1355 maxY = m_maximumY;
1356
1357 return true;
1358 }
1359
1397 void TProjection::doSearch(double minBorder, double maxBorder,
1398 double &extremeVal, const double constBorder,
1399 bool searchX, bool searchLongitude, bool findMin) {
1400 if (minBorder == Null || maxBorder == Null || constBorder == Null) {
1401 return;
1402 }
1403 const double TOLERANCE = PixelResolution()/2;
1404 const int NUM_ATTEMPTS = (unsigned int)DBL_DIG; // It's unsafe to go past
1405 // this precision
1406
1407 double minBorderX, minBorderY, maxBorderX, maxBorderY;
1408 int attempts = 0;
1409
1410 do {
1411 findExtreme(minBorder, maxBorder, minBorderX, minBorderY, maxBorderX,
1412 maxBorderY, constBorder, searchX, searchLongitude, findMin);
1413 if (minBorderX == Null && maxBorderX == Null
1414 && minBorderY == Null && maxBorderY == Null ) {
1415 attempts = NUM_ATTEMPTS;
1416 continue;
1417 }
1418 attempts ++;
1419 }
1420 while ((fabs(minBorderX - maxBorderX) > TOLERANCE
1421 || fabs(minBorderY - maxBorderY) > TOLERANCE)
1422 && (attempts < NUM_ATTEMPTS));
1423 // check both x and y distance in case symmetry of map
1424 // For example, if minBorderX = maxBorderX but minBorderY = -maxBorderY,
1425 // these points may not be close enough.
1426
1427 if (attempts >= NUM_ATTEMPTS) {
1428 // We zoomed in on a discontinuity because our range never shrank, this
1429 // will need to be rechecked later.
1430 // *min and max border should be nearly identical, so it doesn't matter
1431 // which is used here
1432 if (searchLongitude) {
1433 m_specialLatCases.push_back(minBorder);
1434 }
1435 else {
1436 m_specialLonCases.push_back(minBorder);
1437 }
1438 }
1439
1440 // These values will always be accurate, even over a discontinuity
1441 if (findMin) {
1442 if (searchX) extremeVal = min(minBorderX, maxBorderX);
1443 else extremeVal = min(minBorderY, maxBorderY);
1444 }
1445 else {
1446 if (searchX) extremeVal = max(minBorderX, maxBorderX);
1447 else extremeVal = max(minBorderY, maxBorderY);
1448 }
1449 return;
1450 }
1451
1511 void TProjection::findExtreme(double &minBorder, double &maxBorder,
1512 double &minBorderX, double &minBorderY,
1513 double &maxBorderX, double &maxBorderY,
1514 const double constBorder, bool searchX,
1515 bool searchLongitude, bool findMin) {
1516 if (minBorder == Null || maxBorder == Null || constBorder == Null) {
1517 minBorderX = Null;
1518 minBorderY = minBorderX;
1519 minBorderY = minBorderX;
1520 return;
1521 }
1522 if (!searchLongitude && (fabs(fabs(constBorder) - 90.0) < DBL_EPSILON)) {
1523 // it is impossible to search "along" a pole
1524 setSearchGround(minBorder, constBorder, searchLongitude);
1525 minBorderX = XCoord();
1526 minBorderY = YCoord();
1527 maxBorderY = minBorderY;
1528 return;
1529 }
1530 // Always do 10 steps
1531 const double STEP_SIZE = (maxBorder - minBorder) / 10.0;
1532 const double LOOP_END = maxBorder + (STEP_SIZE / 2.0); // This ensures we do
1533 // all of the steps
1534 // properly
1535 double currBorderVal = minBorder;
1536 setSearchGround(minBorder, constBorder, searchLongitude);
1537
1538 // this makes sure that the initial currBorderVal is valid before entering
1539 // the loop below
1540 if (!m_good){
1541 // minBorder = currBorderVal+STEP_SIZE < LOOP_END until setGround is good?
1542 // then, if still not good return?
1543 while (!m_good && currBorderVal <= LOOP_END) {
1544 currBorderVal+=STEP_SIZE;
1545 if (searchLongitude && (currBorderVal - 90.0 > DBL_EPSILON)) {
1546 currBorderVal = 90.0;
1547 }
1548 setSearchGround(currBorderVal, constBorder, searchLongitude);
1549 }
1550 if (!m_good) {
1551 minBorderX = Null;
1552 minBorderY = Null;
1553 return;
1554 }
1555 }
1556
1557 // save the values of three consecutive steps from the minBorder towards
1558 // the maxBorder along the constBorder. initialize these three border
1559 // values (the non-constant lat or lon)
1560 double border1 = currBorderVal;
1561 double border2 = currBorderVal;
1562 double border3 = currBorderVal;
1563
1564 // save the coordinate (x or y) values that correspond to the first
1565 // two borders that are being saved.
1566 // initialize these two coordinate values (x or y)
1567 double value1 = (searchX) ? XCoord() : YCoord();
1568 double value2 = value1;
1569
1570 // initialize the extreme coordinate value
1571 // -- this is the largest coordinate value found so far
1572 double extremeVal2 = value2;
1573
1574 // initialize the extreme border values
1575 // -- these are the borders on either side of the extreme coordinate value
1576 double extremeBorder1 = minBorder;
1577 double extremeBorder3 = minBorder;
1578
1579 while (currBorderVal <= LOOP_END) {
1580
1581 // this conditional was added to prevent trying to SetGround with an
1582 // invalid latitude greater than 90 degrees. There is no need check for
1583 // latitude less than -90 since we start at the minBorder (already
1584 // assumed to be valid) and step forward toward (and possibly past)
1585 // maxBorder
1586 if (searchLongitude && (currBorderVal - 90.0 > DBL_EPSILON)) {
1587 currBorderVal = 90.0;
1588 }
1589
1590 // update the current border value along constBorder
1591 currBorderVal += STEP_SIZE;
1592 setSearchGround(currBorderVal, constBorder, searchLongitude);
1593 if (!m_good){
1594 continue;
1595 }
1596
1597 // update the border and coordinate values
1598 border3 = border2;
1599 border2 = border1;
1600 border1 = currBorderVal;
1601 value2 = value1;
1602 value1 = (searchX) ? XCoord() : YCoord();
1603
1604 if ((findMin && value2 < extremeVal2)
1605 || (!findMin && value2 > extremeVal2)) {
1606 // Compare the coordinate value associated with the center border with
1607 // the current extreme. If the updated coordinate value is more extreme
1608 // (smaller or larger, depending on findMin), then we update the
1609 // extremeVal and it's borders.
1610 extremeVal2 = value2;
1611
1612 extremeBorder3 = border3;
1613 extremeBorder1 = border1;
1614 }
1615 }
1616
1617 // update min/max border values to the values on either side of the most
1618 // extreme coordinate found in this call to this method
1619
1620 minBorder = extremeBorder3; // Border 3 is lagging and thus smaller
1621
1622 // since the loop steps past the original maxBorder, we want to retain
1623 // the original maxBorder value so we don't go outside of the original
1624 // min/max range given
1625 if (extremeBorder1 <= maxBorder ) {
1626 maxBorder = extremeBorder1; // Border 1 is leading and thus larger
1627 }
1628
1629 // update minBorder coordinate values
1630 setSearchGround(minBorder, constBorder, searchLongitude);
1631 // if (!m_good){
1632 // this should not happen since minBorder has already been verified in
1633 // the while loop above
1634 // }
1635
1636 minBorderX = XCoord();
1637 minBorderY = YCoord();
1638
1639 // update maxBorder coordinate values
1640 setSearchGround(maxBorder, constBorder, searchLongitude);
1641 // if (!m_good){
1642 // this should not happen since maxBorder has already been verified in
1643 // the while loop above
1644 // }
1645
1646 maxBorderX = XCoord();
1647 maxBorderY = YCoord();
1648 return;
1649 }
1650
1673 void TProjection::setSearchGround(const double variableBorder,
1674 const double constBorder,
1675 bool variableIsLat) {
1676 if (variableBorder == Null || constBorder == Null) {
1677 return;
1678 }
1679 double lat, lon;
1680 if (variableIsLat) {
1681 lat = variableBorder;
1682 lon = constBorder;
1683 }
1684 else {
1685 lat = constBorder;
1686 lon = variableBorder;
1687 }
1688 SetGround(lat, lon);
1689 return;
1690 }
1691
1692
1699 PvlGroup mapping("Mapping");
1700
1701 QStringList keyNames;
1702 keyNames << "TargetName" << "ProjectionName" << "EquatorialRadius" << "PolarRadius"
1703 << "LatitudeType" << "LongitudeDirection" << "LongitudeDomain"
1704 << "PixelResolution" << "Scale" << "UpperLeftCornerX" << "UpperLeftCornerY"
1705 << "MinimumLatitude" << "MaximumLatitude" << "MinimumLongitude" << "MaximumLongitude"
1706 << "Rotation";
1707
1708 foreach (QString keyName, keyNames) {
1709 if (m_mappingGrp.hasKeyword(keyName)) {
1710 mapping += m_mappingGrp[keyName];
1711 }
1712 }
1713
1714 return mapping;
1715 }
1716
1717
1724 PvlGroup mapping("Mapping");
1725
1726 if (HasGroundRange()) {
1727 mapping += m_mappingGrp["MinimumLatitude"];
1728 mapping += m_mappingGrp["MaximumLatitude"];
1729 }
1730
1731 return mapping;
1732 }
1733
1740 PvlGroup mapping("Mapping");
1741
1742 if (HasGroundRange()) {
1743 mapping += m_mappingGrp["MinimumLongitude"];
1744 mapping += m_mappingGrp["MaximumLongitude"];
1745 }
1746
1747 return mapping;
1748 }
1749
1770 double TProjection::qCompute(const double sinPhi) const {
1771 if (m_eccentricity < DBL_EPSILON) {
1772 QString msg = "Snyder's q variable should only be computed for "
1773 "ellipsoidal projections.";
1774 throw IException(IException::Unknown, msg, _FILEINFO_);
1775 }
1776 double eSinPhi = m_eccentricity * sinPhi;
1777 return (1 - m_eccentricity * m_eccentricity)
1778 * (sinPhi / (1 - eSinPhi * eSinPhi)
1779 - 1 / (2 * m_eccentricity) * log( (1 - eSinPhi) / (1 + eSinPhi) ));
1780 // Note: We know that q is well defined since
1781 // 0 < e < 1 and -1 <= sin(phi) <= 1
1782 // implies that -1 < e*sin(phi) < 1
1783 // Thus, there are no 0 denominators and the log domain is
1784 // satisfied, (1-e*sin(phi))/(1+e*sin(phi)) > 0
1785 }
1786
1803 double TProjection::phi2Compute(const double t) const {
1804 double localPhi = HALFPI - 2.0 * atan(t);
1805 double halfEcc = 0.5 * Eccentricity();
1806 double difference = DBL_MAX;
1807 int iteration = 0;
1808 // a failure in this loop means an exception will be thrown, which is
1809 // an expensive operation. So letting let the loop iterate quite a bit
1810 // is not a big deal. Also, the user will be unable to project at all
1811 // if this fails so more reason to let this loop iterate: better to
1812 // function slow than not at all.
1813 const int MAX_ITERATIONS = 45;
1814
1815 while ((iteration < MAX_ITERATIONS) && (difference > 0.0000000001)) {
1816 double eccTimesSinphi = Eccentricity() * sin(localPhi);
1817 double newPhi = HALFPI -
1818 2.0 * atan(t * pow((1.0 - eccTimesSinphi) /
1819 (1.0 + eccTimesSinphi), halfEcc));
1820 difference = fabs(newPhi - localPhi);
1821 localPhi = newPhi;
1822 iteration++;
1823 }
1824
1825 if (iteration >= MAX_ITERATIONS) {
1826 QString msg = "Failed to converge in TProjection::phi2Compute()";
1827 throw IException(IException::Unknown, msg, _FILEINFO_);
1828 }
1829
1830 return localPhi;
1831 }
1832
1847 double TProjection::mCompute(const double sinphi, const double cosphi) const {
1848 double eccTimesSinphi = Eccentricity() * sinphi;
1849 double denominator = sqrt(1.0 - eccTimesSinphi * eccTimesSinphi);
1850 return cosphi / denominator;
1851 }
1852
1870 double TProjection::tCompute(const double phi, const double sinphi) const {
1871 if ((HALFPI) - fabs(phi) < DBL_EPSILON) return 0.0;
1872
1873 double eccTimesSinphi = Eccentricity() * sinphi;
1874 double denominator = pow((1.0 - eccTimesSinphi) /
1875 (1.0 + eccTimesSinphi),
1876 0.5 * Eccentricity());
1877 return tan(0.5 * (HALFPI - phi)) / denominator;
1878 }
1879
1891 double TProjection::e4Compute() const {
1892 double onePlusEcc = 1.0 + Eccentricity();
1893 double oneMinusEcc = 1.0 - Eccentricity();
1894
1895 return sqrt(pow(onePlusEcc, onePlusEcc) *
1896 pow(oneMinusEcc, oneMinusEcc));
1897 }
1898
1899} //end namespace isis
1900
1901
double degrees() const
Get the angle in units of Degrees.
Definition Angle.h:232
@ Degrees
Degrees are generally considered more human readable, 0-360 is one circle, however most math does not...
Definition Angle.h:56
Isis exception class.
Definition IException.h:91
@ Unknown
A type of error that cannot be classified as any of the other error types.
Definition IException.h:118
This class is designed to encapsulate the concept of a Longitude.
Definition Longitude.h:40
Longitude force180Domain() const
This returns a longitude that is constricted to -180 to 180 degrees.
Longitude force360Domain() const
This returns a longitude that is constricted to 0-360 degrees.
Base class for Map Projections.
Definition Projection.h:155
@ Triaxial
These projections are used to map triaxial and irregular-shaped bodies.
Definition Projection.h:166
WorldMapper * m_mapper
This points to a mapper passed into the SetWorldMapper method.
Definition Projection.h:292
double XCoord() const
This returns the projection X provided SetGround, SetCoordinate, SetUniversalGround,...
double m_maximumX
See minimumX description.
Definition Projection.h:326
virtual bool HasGroundRange() const
This indicates if the longitude direction type is positive west (as opposed to postive east).
bool m_groundRangeGood
Indicates if the ground range (min/max lat/lons) were read from the labels.
Definition Projection.h:313
double PixelResolution() const
Returns the pixel resolution value from the PVL mapping group in meters/pixel.
double YCoord() const
This returns the projection Y provided SetGround, SetCoordinate, SetUniversalGround,...
bool m_good
Indicates if the contents of m_x, m_y, m_latitude, and m_longitude are valid.
Definition Projection.h:300
double m_minimumX
The data elements m_minimumX, m_minimumY, m_maximumX, and m_maximumY are convience data elements when...
Definition Projection.h:317
PvlGroup m_mappingGrp
Mapping group that created this projection.
Definition Projection.h:329
double m_minimumY
See minimumX description.
Definition Projection.h:327
void SetXY(double x, double y)
This protected method is a helper for derived classes.
bool IsGood() const
This indicates if the last invocation of SetGround, SetCoordinate, SetUniversalGround,...
void setProjectionType(const ProjectionType ptype)
Sets the projection subclass type.
double m_maximumY
See minimumX description.
Definition Projection.h:328
void SetComputedXY(double x, double y)
This protected method is a helper for derived classes.
bool hasKeyword(const QString &name) const
Check to see if a keyword exists.
Contains multiple PvlContainers.
Definition PvlGroup.h:41
Container for cube-like labels.
Definition Pvl.h:119
Base class for Map TProjections.
double m_longitude
This contains the currently set longitude value.
virtual bool SetGround(const double lat, const double lon)
This method is used to set the latitude/longitude (assumed to be of the correct LatitudeType,...
bool xyRangeOblique(double &minX, double &maxX, double &minY, double &maxY)
This method is used to find the XY range for oblique aspect projections (non-polar projections) by "w...
bool SetUnboundUniversalGround(const double coord1, const double coord2)
This method is used to set the latitude/longitude.
bool IsPlanetocentric() const
This indicates if the latitude type is planetocentric (as opposed to planetographic).
double m_minimumLatitude
Contains the minimum latitude for the entire ground range.
virtual 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...
double m_polarRadius
Polar radius of the target.
std::vector< double > m_specialLatCases
Constant Latitudes that intersect a discontinuity.
double m_maximumLongitude
Contains the maximum longitude for the entire ground range.
static double To180Domain(const double lon)
This method converts a longitude into the -180 to 180 domain.
double m_equatorialRadius
Polar radius of the target.
LongitudeDirection m_longitudeDirection
An enumerated type indicating the LongitudeDirection read from the labels.
virtual PvlGroup MappingLongitudes()
This function returns the longitude keywords that this projection uses.
virtual double MaximumLatitude() const
This returns the maximum latitude of the area of interest.
virtual ~TProjection()
Destroys the TProjection object.
void findExtreme(double &minBorder, double &maxBorder, double &minBorderX, double &minBorderY, double &maxBorderX, double &maxBorderY, const double constBorder, bool searchX, bool searchLongitude, bool findMin)
Searches for extreme (min/max/discontinuity) coordinate values across latitudes/longitudes.
double LocalRadius() const
This method returns the local radius in meters at the current latitude position.
virtual bool SetCoordinate(const double x, const double y)
This method is used to set the projection x/y.
int m_longitudeDomain
This integer is either 180 or 360 and is read from the labels.
double mCompute(const double sinphi, const double cosphi) const
A convience method to compute Snyder's m equation (14-15) for a given latitude, .
double ToPlanetocentric(const double lat) const
This method converts a planetographic latitude to a planetocentric latitude.
bool Has180Domain() const
This indicates if the longitude domain is -180 to 180 (as opposed to 0 to 360).
virtual bool SetUniversalGround(const double lat, const double lon)
This method is used to set the latitude/longitude which must be Planetocentric (latitude) and Positiv...
std::vector< double > m_specialLonCases
Constant Longitudes that intersect a discontinuity.
void setSearchGround(const double variableBorder, const double constBorder, bool variableIsLat)
This function sets the ground for the given border values.
virtual double MinimumLongitude() const
This returns the minimum longitude of the area of interest.
virtual double UniversalLongitude()
This returns a universal longitude (positive east in 0 to 360 domain).
double PolarRadius() const
This returns the polar radius of the target.
bool inLongitudeRange(double longitude)
Determine whether the given longitude is within the range of the MinimumLongitude and MaximumLongitud...
void XYRangeCheck(const double latitude, const double longitude)
This convience function is established to assist in the development of the XYRange virtual method.
double Eccentricity() const
This returns the eccentricity of the target,.
bool IsPositiveWest() const
This indicates if the longitude direction type is positive east (as opposed to postive west).
QString LongitudeDirectionString() const
This method returns the longitude direction as a string.
bool inLatitudeRange(double latitude)
Determine whether the given latitude is within the range of the MinimumLatitude and MaximumLatitude r...
virtual PvlGroup MappingLatitudes()
This function returns the latitude keywords that this projection uses.
double m_eccentricity
The eccentricity of the target body.
double qCompute(const double sinPhi) const
A convience method to compute Snyder's q equation (3-12) for a given latitude, .
virtual bool operator==(const Projection &proj)
This method determines whether two map projection objects are equal by comparing the equatorial radiu...
double phi2Compute(const double t) const
A convience method to compute latitude angle phi2 given small t, from Syder's recursive equation (7-9...
double m_minimumLongitude
Contains the minimum longitude for the entire ground range.
void doSearch(double minBorder, double maxBorder, double &extremeVal, const double constBorder, bool searchX, bool searchLongitude, bool findMin)
This method searches for extreme (min/max/discontinuity) coordinate values along the constBorder line...
double e4Compute() const
A convience method to compute.
double tCompute(const double phi, const double sinphi) const
A convience method to compute Snyder's t equation (15-9) for a given latitude, .
@ PositiveWest
Longitude values increase in the westerly direction.
@ PositiveEast
Longitude values increase in the easterly direction.
double m_maximumLatitude
Contains the maximum latitude for the entire ground range.
virtual double Latitude() const
This returns a latitude with correct latitude type as specified in the label object.
static double ToPositiveWest(const double lon, const int domain)
This method converts a longitude into the positive west direction.
static double To360Domain(const double lon)
This method converts a longitude into the 0 to 360 domain.
TProjection(Pvl &label)
Constructs an empty TProjection object.
QString LongitudeDomainString() const
This method returns the longitude domain as a string.
bool IsPlanetographic() const
This indicates if the latitude type is planetographic (as opposed to planetocentric).
virtual bool IsEquatorialCylindrical()
This method returns true if the projection is equatorial cylindrical.
QString LatitudeTypeString() const
This method returns the latitude type as a string.
bool Has360Domain() const
This indicates if the longitude domain is 0 to 360 (as opposed to -180 to 180).
@ Planetocentric
Latitudes are measured as the angle from the equatorial plane to the plane through the center of the ...
@ Planetographic
Latitudes are measured as the angle from the equatorial plane to the normal to the surface of the pla...
virtual double MinimumLatitude() const
This returns the minimum latitude of the area of interest.
LatitudeType m_latitudeType
An enumerated type indicating the LatitudeType read from the labels.
virtual double MaximumLongitude() const
This returns the maximum longitude of the area of interest.
virtual PvlGroup Mapping()
This function returns the keywords that this projection uses.
static double ToPositiveEast(const double lon, const int domain)
This method converts a longitude into the positive east direction.
virtual double UniversalLatitude()
This returns a universal latitude (planetocentric).
double EquatorialRadius() const
This returns the equatorial radius of the target.
bool IsPositiveEast() const
This indicates if the longitude direction type is positive west (as opposed to postive east).
double Scale() const
This method returns the scale for mapping world coordinates into projection coordinates.
virtual double Longitude() const
This returns a longitude with correct longitude direction and domain as specified in the label object...
double ToPlanetographic(const double lat) const
This method converts a planetocentric latitude to a planetographic latitude.
double m_latitude
This contains the currently set latitude value.
virtual double TrueScaleLatitude() const
This method returns the latitude of true scale.
static PvlGroup radiiGroup(QString target)
Creates a Pvl Group with keywords TargetName, EquitorialRadius, and PolarRadius.
Definition Target.cpp:428
virtual double Resolution() const
This virtual method will the resolution of the world system relative to one unit in the projection sy...
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition IString.cpp:211
const double Null
Value for an Isis Null pixel.
const double HALFPI
The mathematical constant PI/2.
Definition Constants.h:41
const double PI
The mathematical constant PI.
Definition Constants.h:40
Namespace for the standard library.