USGS

Isis 3.0 Application Source Code Reference

Home

MdisCalUtils.h
Go to the documentation of this file.
1 #ifndef MdisCalUtils_h
2 #define MdisCalUtils_h
3 /**
4  * @file
5  * $Revision: 6715 $
6  * $Date: 2016-04-28 10:58:43 -0700 (Thu, 28 Apr 2016) $
7  *
8  * Unless noted otherwise, the portions of Isis written by the USGS are
9  * public domain. See individual third-party library and package descriptions
10  * for intellectual property information, user agreements, and related
11  * information.
12  *
13  * Although Isis has been used by the USGS, no warranty, expressed or
14  * implied, is made by the USGS as to the accuracy and functioning of such
15  * software and related material nor shall the fact of distribution
16  * constitute any such warranty, and no responsibility is assumed by the
17  * USGS in connection therewith.
18  *
19  * For additional information, launch
20  * $ISISROOT/doc//documents/Disclaimers/Disclaimers.html
21  * in a browser or see the Privacy & Disclaimers page on the Isis website,
22  * http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
23  * http://www.usgs.gov/privacy.html.
24  */
25 #include <cmath>
26 #include <string>
27 #include <vector>
28 
29 #include "CSVReader.h"
30 #include "FileName.h"
31 #include "IString.h"
32 #include "Spice.h"
33 /**
34  * @author ????-??-?? Kris Becker
35  *
36  * @internal
37  * @history 2015-08-31 Jeannie Backer - Changed method name from loadContaminationEvent() to
38  * loadEmpiricalCorrection(). Brought code closer to coding standards.
39  */
40 namespace Isis {
41  /**
42  * @brief Helper function to convert values to doubles
43  *
44  * @param T Type of value to convert
45  * @param value Value to convert
46  *
47  * @return double Converted value
48  */
49  template <typename T> double toDouble(const T &value) {
50  return toDouble(QString(value).trimmed());
51  }
52 
53  template <typename T> int toInteger(const T &value) {
54  return toInt(QString(value).trimmed());
55  }
56 
57  template <typename T> inline T MIN(const T &A, const T &B) {
58  return (((A) < (B)) ? (A) : (B));
59  }
60 
61  template <typename T> inline T MAX(const T &A, const T &B) {
62  return (((A) > (B)) ? (A) : (B));
63  }
64 
65  inline QString quote(const QString &value) {
66  if (value.isEmpty()) return (value);
67  if (value[0] == '"') return (value);
68  return (QString('"' + value + '"'));
69  }
70 
71  /**
72  * @brief Load required NAIF kernels required for timing needs
73  *
74  * This method maintains the loading of kernels for MESSENGER timing and
75  * planetary body ephemerides to support time and relative positions of planet
76  * bodies.
77  */
78  static void loadNaifTiming() {
79  static bool naifLoaded = false;
80  if (!naifLoaded) {
81 // Load the NAIF kernels to determine timing data
82  Isis::FileName leapseconds("$base/kernels/lsk/naif????.tls");
83  leapseconds = leapseconds.highestVersion();
84 
85  Isis::FileName sclk("$messenger/kernels/sclk/messenger_????.tsc");
86  sclk = sclk.highestVersion();
87 
88  Isis::FileName pck("$base/kernels/spk/de???.bsp");
89  pck = pck.highestVersion();
90 
91 // Load the kernels
92  QString leapsecondsName(leapseconds.expanded());
93  QString sclkName(sclk.expanded());
94  QString pckName(pck.expanded());
95  furnsh_c(leapsecondsName.toLatin1().data());
96  furnsh_c(sclkName.toLatin1().data());
97  furnsh_c(pckName.toLatin1().data());
98 
99 // Ensure it is loaded only once
100  naifLoaded = true;
101  }
102  return;
103  }
104 
105  /**
106  * @brief Computes the distance from the Sun to the observed body
107  *
108  * This method requires the appropriate NAIK kernels to be loaded that
109  * provides instrument time support, leap seconds and planet body ephemeris.
110  *
111  * @return double Distance in AU between Sun and observed body
112  */
113  static bool sunDistanceAU(const QString &scStartTime,
114  const QString &target,
115  double &sunDist) {
116 
117  // Ensure NAIF kernels are loaded
118  loadNaifTiming();
119  sunDist = 1.0;
120 
121  // Determine if the target is a valid NAIF target
122  SpiceInt tcode;
123  SpiceBoolean found;
124  bodn2c_c(target.toLatin1().data(), &tcode, &found);
125  if (!found) return (false);
126 
127  // Convert starttime to et
128  double obsStartTime;
129  scs2e_c(-236, scStartTime.toLatin1().data(), &obsStartTime);
130 
131  // Get the vector from target to sun and determine its length
132  double sunv[3];
133  double lt;
134  spkpos_c(target.toLatin1().data(), obsStartTime, "J2000", "LT+S", "sun",
135  sunv, &lt);
136  double sunkm = vnorm_c(sunv);
137 
138  // Return in AU units
139  sunDist = sunkm / 1.49597870691E8;
140  return (true);
141  }
142 
143  std::vector<double> loadWACCSV(const QString &fname, int filter,
144  int nvalues, bool header=true, int skip=0) {
145  // Open the CSV file
146  FileName csvfile(fname);
147  CSVReader csv(csvfile.expanded(), header, skip);
148  for (int i = 0; i < csv.rows(); i++) {
149  CSVReader::CSVAxis row = csv.getRow(i);
150  if (toInteger(row[0]) == filter) {
151  if ((row.dim1() - 1) < nvalues) {
152  QString mess = "Number values (" + QString(row.dim1() - 1) +
153  ") in file " + fname +
154  " less than number requested (" +
155  QString(nvalues) + ")!";
156  throw IException(IException::User, mess, _FILEINFO_);
157  }
158 
159  std::vector<double> rsp;
160  for (int c = 0; c < nvalues; c++) {
161  rsp.push_back(toDouble(row[1+c]));
162  }
163  return (rsp);
164  }
165  }
166 
167  // If it reaches here, the filter was not found
168  std::ostringstream mess;
169  mess << "CSV Vector MDIS filter " << filter << ", not found in file "
170  << fname << "!";
171  throw IException(IException::User, mess.str(), _FILEINFO_);
172  }
173 
174 
175  std::vector<double> loadNACCSV(const QString &fname, int nvalues,
176  bool header=true, int skip=0) {
177  // Open the CSV file
178  FileName csvfile(fname);
179  CSVReader csv(csvfile.expanded(), header, skip);
180  CSVReader::CSVAxis row = csv.getRow(0);
181  if (row.dim1() < nvalues) {
182  QString mess = "Number values (" + QString(row.dim1()) +
183  ") in file " + fname + " less than number requested (" +
184  QString(nvalues) + ")!";
185  throw IException(IException::User, mess, _FILEINFO_);
186 
187  }
188  std::vector<double> rsp;
189  for (int i = 0; i < nvalues; i++) {
190  rsp.push_back(toDouble(row[i]));
191  }
192  return (rsp);
193  }
194 
195 
196  std::vector<double> loadResponsivity(bool isNAC, bool binned, int filter,
197  QString &fname) {
198 
199  FileName resfile(fname);
200  if (fname.isEmpty()) {
201  QString camstr = (isNAC) ? "NAC" : "WAC";
202  QString binstr = (binned) ? "_BINNED" : "_NOTBIN";
203  QString base = "$messenger/calibration/RESPONSIVITY/";
204  resfile = base + "MDIS" + camstr + binstr + "_RESP_?.TAB";
205  resfile = resfile.highestVersion();
206  fname = resfile.originalPath() + "/" + resfile.name();
207  }
208 
209  // Unfortunately NAC has a slightly different format so must do it
210  // explicitly
211  if (isNAC) {
212  return (loadNACCSV(fname, 4, false, 0));
213  }
214  else {
215  // Load the WAC parameters
216  return (loadWACCSV(fname, filter, 4, false, 0));
217  }
218  }
219 
220 
221  std::vector<double> loadSolarIrr(bool isNAC, bool binned, int filter,
222  QString &fname) {
223 
224  FileName solfile(fname);
225  if (fname.isEmpty()) {
226  QString camstr = (isNAC) ? "NAC" : "WAC";
227  QString base = "$messenger/calibration/SOLAR/";
228  solfile = base + "MDIS" + camstr + "_SOLAR_?.TAB";
229  solfile = solfile.highestVersion();
230  fname = solfile.originalPath() + "/" + solfile.name();
231  }
232 
233  if (isNAC) {
234  return (loadNACCSV(fname, 3, false, 0));
235  }
236  else {
237  return (loadWACCSV(fname, filter, 3, false, 0));
238  }
239  }
240 
241  double loadSmearComponent(bool isNAC, int filter, QString &fname) {
242 
243  FileName smearfile(fname);
244  if (fname.isEmpty()) {
245  QString camstr = (isNAC) ? "NAC" : "WAC";
246  QString base = "$messenger/calibration/smear/";
247  smearfile = base + "MDIS" + camstr + "_FRAME_TRANSFER_??.TAB";
248  smearfile = smearfile.highestVersion();
249  fname = smearfile.originalPath() + "/" + smearfile.name();
250  }
251 
252  std::vector<double> smear;
253  if (isNAC) {
254  smear = loadNACCSV(fname, 1, false, 0);
255  }
256  else {
257  smear = loadWACCSV(fname, filter, 1, false, 0);
258  }
259  return (smear[0]);
260  }
261 
262 /**
263  * @brief Load and retrieve empirical correction factor
264  *
265  * This function determines the empirical correction factor for changes that
266  * that occured on the spacecraft after Mercury orbit insertion. The
267  * affected dates are May 24, 2011 to January 3, 2012.
268  *
269  * The table of correction factors is expected to be stored in
270  * $messenger/calibration/events/event_table_ratioed_v?.txt. However, the
271  * caller may provide a table that conforms to the expected format. The
272  * expected format for the empirical correction file is a comma separated value
273  * (CSV) table that contains 13 columns of data per row. The first column is
274  * the UTC time during the event. The next 12 columns contain multiplicative
275  * correction factors for each WAC filter (NAC correction factors are not
276  * provided). These factors are expected to be around 1.0 (the default) as it
277  * is expected to directly scale DN values.
278  *
279  * Below is the expected mapping of column indexes to filter numbers as
280  * specfied in the BandBin/Number keyword from MDIS ISIS cube labels. Index is
281  * the column index from each row for a given filter, Number is the value of
282  * the BandBin/Number keyword from the label designating the filter number
283  * (corresponding to the filter parameter passed to this routine) and Letter is
284  * the filter letter designation used in the last alpha numeric character in
285  * MDIS filenames:
286  *
287  * Index Number Letter Wavelength
288  * 1 6 F 430 nm
289  * 2 3 C 480 nm
290  * 3 4 D 560 nm
291  * 4 5 E 630 nm
292  * 5 7 G 750 nm
293  * 6 12 L 830 nm
294  * 7 10 J 900 nm
295  * 8 9 I 1000 nm
296  * 9 1 A Filter 1 (700 nm)
297  * 10 2 B Filter 2 (clear)
298  * 11 8 H Filter 8 (950 nm)
299  * 12 11 K Filter 11 (1010 nm)
300  *
301  * The UTC dates in the first column are assumed to be strictly increasing in
302  * time. The initial table (*_v2) contains dates that span the complete
303  * expected timeframe of the mission (launch at 2004-08-04T10:00:00.000000,
304  * termination at 2015-01-03T09:00:00.000000).
305  *
306  * The spacecraft clock time is provided as input (scStartTime) to this
307  * function. This value is converted to ET (SCET) and used to determine the
308  * corresponding event time in the first column of the table. The first table
309  * column time is represented in UTC time. This time is converted to ET and
310  * then compared with the start time in ET.
311  *
312  * The algorithm searches linearly through the table essentially storing the
313  * time slot prior to the SCET and the next occuring one. Ultimately, the
314  * factor returned by the algorithm is the one whose event time is closest to
315  * the SCET.
316  *
317  * The empirical correction model and algorithm was developed by
318  * Mary Ruth Keller of JHA/APL.
319  *
320  * @author Kris Becker - 10/23/2012
321  *
322  * @param scStartTime - Start time of the image in SCLK format
323  * @param filter - WAC filter number to return event correction factor for
324  * @param ename - Returns the name of the event table file if not
325  * provided by caller. If non-empty string is passed by
326  * caller, it is assumed to be a fully qualified filename
327  * of the event table.
328  * @param eDate - Returns the UTC date entry of the selected correction
329  * event factor
330  *
331  * @return double - Event correction factor at the selected time to apply
332  * to WAC filter data.
333  */
334  double loadEmpiricalCorrection(const QString &scStartTime, const int filter,
335  QString &ename, QString &eDate) {
336 
337  // This table maps the filter number extracted from BandBin/Number keyword
338  // to the columns (index) in the empirical correction table
339  const int filterMap[12] = { 6, 3, 4, 5, 7, 12, 10, 9, 1, 2, 8, 11 };
340 
341  // Find the WAC filter column index
342  int ncols = sizeof(filterMap)/sizeof(filterMap[0]);
343  int column = -1;
344  for (int c = 0; c < ncols; c++) {
345  if (filterMap[c] == filter) {
346  column = c + 1; // indexes start after 1st (time) column
347  break;
348  }
349  }
350 
351  // Ensure we have a valid filter number
352  if (column <= 0) {
353  std::ostringstream mess;
354  mess << "Invalid MDIS WAC filter number (" << filter
355  << " - range:1-12) for determining index into empirical correction table.";
356  throw IException(IException::User, mess.str(), _FILEINFO_);
357  }
358 
359  // File name not provided by caller. Determine the event table name
360  if ( ename.isEmpty() ) {
361  FileName eventfile("$messenger/calibration/events/event_table_ratioed_v?.txt");
362  eventfile = eventfile.highestVersion();
363  ename = eventfile.originalPath() + "/" + eventfile.name();
364  }
365 
366  // Open/read the CSV empirical correction file
367  FileName csvfile(ename);
368  const bool header = false; // No header in file
369  const int skip = 0; // No lines to skip to data
370  const int nvalues = 13; // Expected columns in table
371  CSVReader csv(csvfile.expanded(), header, skip);
372  if (csv.columns() < nvalues) { // All rows should have same # columns
373  QString mess = "Number values (" + QString(csv.columns()) +
374  ") in file " + ename + " less than number requested (" +
375  QString(nvalues) + ")!";
376  throw IException(IException::User, mess, _FILEINFO_);
377  }
378 
379  // Ensure NAIF kernels are loaded for NAIF time computations
380  loadNaifTiming();
381 
382  // Convert s/c clock start time to et
383  double obsStartTime;
384  scs2e_c(-236, scStartTime.toLatin1().data(), &obsStartTime);
385 
386  // Set initial conditions and loop through all rows in the event table
387  double evalue = 1.0;
388  eDate = "N/A"; // Will attain a valid time on guaranteed first pass
389  double preEventTime = 0.0;
390  for (int i = 0; i < csv.rows(); i++) {
391  CSVReader::CSVAxis eRow = csv.getRow(i);
392  QString utcTime = eRow[0];
393  double eTime;
394  utc2et_c(utcTime.toLatin1().data(), &eTime);
395 
396  // If current time is greater than start time this is the post event case
397  if (eTime > obsStartTime) {
398  // Get closest pre or post event correction factor
399  if ( fabs(obsStartTime-preEventTime) > fabs(eTime-obsStartTime) ) {
400  // Post-event time closer to SCLK than Pre-event time
401  eDate = utcTime;
402  evalue = toDouble(eRow[column]);
403  }
404 
405  break; // Terminate loop and return
406  }
407 
408  // Record pre-event time slot - Sets return variables as well
409  eDate = utcTime;
410  preEventTime = eTime;
411  evalue = toDouble(eRow[column]);
412  }
413 
414  // Return the factor
415  return (evalue);
416  }
417 
418 
419 };
420 #endif
double base
Definition: tonematch.cpp:18
double sunDist
Definition: mar10cal.cpp:25
T MIN(const T &A, const T &B)
Implement templatized MIN fumnction.
Definition: Hillier.h:38
std::vector< double > loadWACCSV(const QString &fname, int filter, int nvalues, bool header=true, int skip=0)
Definition: MdisCalUtils.h:143
QString utcTime
Definition: apollo2isis.cpp:36
T MAX(const T &A, const T &B)
Implement templatized MAX function.
Definition: Hillier.h:48
double loadSmearComponent(bool isNAC, int filter, QString &fname)
Definition: MdisCalUtils.h:241
double loadEmpiricalCorrection(const QString &scStartTime, const int filter, QString &ename, QString &eDate)
Load and retrieve empirical correction factor.
Definition: MdisCalUtils.h:334
int toInteger(const T &value)
Definition: MdisCalUtils.h:53
std::vector< double > loadResponsivity(bool isNAC, bool binned, int filter, QString &fname)
Definition: MdisCalUtils.h:196
std::vector< double > loadSolarIrr(bool isNAC, bool binned, int filter, QString &fname)
Definition: MdisCalUtils.h:221
double toDouble(const T &value)
Helper function to convert values to doubles.
Definition: MdisCalUtils.h:49
QString quote(const QString &value)
Definition: MdisCalUtils.h:65
Isis::CubeCalculator c
Definition: fx.cpp:35
void filter(Buffer &in, double &v)
The filter depends on the user input kernel.
Definition: kernfilter.cpp:71
std::vector< double > loadNACCSV(const QString &fname, int nvalues, bool header=true, int skip=0)
Definition: MdisCalUtils.h:175
bool trimmed
Definition: trimfilter.cpp:10