OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
scpad_viirs.cpp
Go to the documentation of this file.
1 #include <math.h>
2 #include <stdio.h>
3 #include <libgen.h>
4 #include <stdlib.h>
5 #include <getopt.h>
6 #include <stdint.h>
7 #include <string.h>
8 
9 #include <sstream>
10 #include <iostream>
11 #include <fstream>
12 #include "nc4utils.h"
13 #include "l1agen_viirs.h"
14 
15 #define SWAP_2(x) ( (((x) & 0xff) << 8) | ((unsigned short)(x) >> 8) )
16 
17 #define VERSION "1.15"
18 
19 // Modification history:
20 // Programmer Organization Date Ver Description of change
21 // ---------- ------------ ---- --- ---------------------
22 // Joel Gales FutureTech 09/25/15 0.10 Original development
23 // Joel Gales FutureTech 09/30/15 0.11 Return 110 if time gaps
24 // in output file
25 // Joel Gales FutureTech 10/01/15 0.12 Return 120 if no overlap
26 // with prev/following files
27 // Joel Gales FutureTech 10/01/15 0.13 Set default nPad to 10"
28 // Joel Gales FutureTech 10/01/15 0.14 Add support in case files
29 // abut with no overlap
30 // Joel Gales FutureTech 10/02/15 0.15 Assume prev/foll files abut
31 // Joel Gales FutureTech 10/05/15 0.16 Remove path from filename
32 // before determining filetype
33 // Joel Gales FutureTech 10/06/15 0.17 Change toff to 6 for S/C recs
34 // Joel Gales FutureTech 10/06/15 0.18 Change get_pkt_sec() to
35 // support both input and output
36 // of julian date
37 // Joel Gales FutureTech 10/09/15 0.19 Change time comparison to
38 // (secEnd+nPad + 0.5) for the
39 // following file to avoid
40 // roundoff error
41 // Joel Gales FutureTech 10/15/15 0.19 Change time comparison to
42 // (secEnd+nPad - 0.5) for the
43 // preceeding file
44 // Joel Gales FutureTech 12/07/16 1.00 Add support for JPSS1
45 // Joel Gales FutureTech 02/27/18 1.10 Exit with 1 if all times in
46 // following file less than
47 // endtime (plus pad) of file
48 // to be padded
49 // Joel Gales SAIC 03/11/20 1.11 Exit with 110 if all times in
50 // following file less than
51 // endtime (plus pad) of file
52 // to be padded
53 // Sean Bailey NASA 03/18/20 1.12 Changed to exit with 120 for
54 // above
55 // Joel Gales SAIC 08/16/21 1.13 Add support for JPSS2
56 // Joel Gales SAIC 09/29/21 1.14 Add support for JPSS2 GPS
57 // Joel Gales SAIC 10/04/21 1.15 Fix problems with 10 Hz files
58 
59 using namespace std;
60 
61 extern "C" int32_t jday( int16_t i, int16_t j, int16_t k);
62 extern "C" int jdate( int32_t julian, int32_t *year, int32_t *doy);
63 extern "C" int ccsds_to_yds( uint8_t *cctime, int32_t *iyear, int32_t *iday,
64  double *sec);
65 
66 double get_pkt_sec( uint8_t *pkt_time, int32_t in_jd, int32_t *out_jd) {
67 
68  uint16_t ui16a;
69  double sec;
70  int16_t iy16, idy16;
71  int32_t iy, idy, jd, iyr, iday, jd0;
72 
73  if ( in_jd == -1) {
74  memcpy( &ui16a, (uint16_t *) pkt_time, 2);
75  jd0 = SWAP_2( ui16a) + 2436205;
76  } else {
77  jd0 = in_jd;
78  }
79 
80  if ( out_jd != NULL) *out_jd = jd0;
81 
82  ccsds_to_yds( pkt_time, &iy, &idy, &sec);
83  iy16 = iy;
84  idy16 = idy;
85  jd = jday( iy16, 1, idy16);
86 
87  return sec + (jd - jd0)*86400;
88 }
89 
90 
91 int main (int argc, char* argv[])
92 {
93 
94  cout << "scpad_viirs " << VERSION << " ("
95  << __DATE__ << " " << __TIME__ << ")" << endl;
96 
97  if ( argc == 1) {
98  cout << endl <<
99  "scpad_viirs [-n secs-to-pad]" << endl <<
100  " { [-p previous-file]" << endl <<
101  " [-f following-file] }" << endl <<
102  " file-to-pad" << endl <<
103  " [output-file]" << endl;
104  return 0;
105  }
106 
107  int c;
108  int index;
109  int nSec2Pad = 10;
110  char *pfile = NULL;
111  char *ffile = NULL;
112  char *file_to_pad = NULL;
113  while ((c = getopt (argc, argv, "n:p:f:")) != -1)
114  switch (c)
115  {
116  case 'n':
117  nSec2Pad = atoi( optarg);
118  break;
119 
120  case 'p':
121  pfile = optarg;
122  break;
123 
124  case 'f':
125  ffile = optarg;
126  break;
127 
128  default:
129  abort ();
130  }
131 
132  // cout << nSec2Pad << endl;
133 
134  if ( pfile == 0x0 && ffile == 0x0) {
135  cout << "Neither previous nor following file specified" << endl;
136  return 120;
137  }
138 
139  fstream xfileStream, pfileStream, ffileStream, ofileStream;
140 
141  xfileStream.open( argv[optind], fstream::in | fstream::binary);
142  if ( xfileStream.fail()) {
143  cout << argv[optind] << " not found" << endl;
144  exit(1);
145  }
146 
147 
148  xfileStream.seekg (0, ios::end);
149  int xFileSize = xfileStream.tellg();
150  xfileStream.seekg (0, ios::beg);
151 
152  char fileType = 0;
153  if (strncmp(basename(argv[optind]), "P1570011", 8) == 0) fileType = 'S';
154  if (strncmp(basename(argv[optind]), "P1570008", 8) == 0) fileType = 'A';
155  if (strncmp(basename(argv[optind]), "P1570000", 8) == 0) fileType = 'B';
156 
157  if (strncmp(basename(argv[optind]), "P1590011", 8) == 0) fileType = 'S';
158  if (strncmp(basename(argv[optind]), "P1590008", 8) == 0) fileType = 'A';
159  if (strncmp(basename(argv[optind]), "P1590000", 8) == 0) fileType = 'B';
160 
161  if (strncmp(basename(argv[optind]), "P1770011", 8) == 0) fileType = 'S';
162  if (strncmp(basename(argv[optind]), "P1770030", 8) == 0) fileType = 'A';
163  if (strncmp(basename(argv[optind]), "P1770037", 8) == 0) fileType = 'B';
164  if (strncmp(basename(argv[optind]), "P1770034", 8) == 0) fileType = 'G';
165 
166  if (fileType == 0) {
167  cout << "File type cannot be determined from filename" << endl;
168  exit(1);
169  }
170 
171  int nPad = nSec2Pad;
172  if (strncmp(basename(argv[optind]), "P1770011", 8) == 0) nPad = 10*nSec2Pad;
173  if (strncmp(basename(argv[optind]), "P1770030", 8) == 0) nPad = 10*nSec2Pad;
174 
175  float timeGap = 0.5;
176  if (strncmp(basename(argv[optind]), "P1770011", 8) == 0) timeGap = 0.05;
177  if (strncmp(basename(argv[optind]), "P1770030", 8) == 0) timeGap = 0.05;
178 
179  float gapCheck = 1.5;
180  if (strncmp(basename(argv[optind]), "P1770011", 8) == 0) gapCheck = 0.15;
181  if (strncmp(basename(argv[optind]), "P1770030", 8) == 0) gapCheck = 0.15;
182 
183  int32_t recSize, toff;
184  toff = 6;
185 
186  if (strncmp(basename(argv[optind]), "P157", 4) == 0) {
187  if (fileType == 'S') {
188  recSize = 71;
189  } else if (fileType == 'A') {
190  recSize = 355;
191  } else if (fileType == 'B') {
192  recSize = 207;
193  }
194  }
195 
196  if (strncmp(basename(argv[optind]), "P159", 4) == 0) {
197  if (fileType == 'S') {
198  recSize = 71;
199  } else if (fileType == 'A') {
200  recSize = 393;
201  } else if (fileType == 'B') {
202  recSize = 183;
203  }
204  }
205 
206  if (strncmp(basename(argv[optind]), "P177", 4) == 0) {
207  if (fileType == 'S') {
208  recSize = 71;
209  } else if (fileType == 'A') {
210  recSize = 493;
211  } else if (fileType == 'B') {
212  recSize = 212;
213  } else if (fileType == 'G') {
214  recSize = 374;
215  }
216  }
217 
218  int32_t nRec = xFileSize / recSize;
219 
220  // Allocate continuous 2-dim array
221  uint8_t **pkts = new uint8_t*[nRec];
222  pkts[0] = new uint8_t[recSize*nRec];
223  for (int i=1; i<nRec; i++) pkts[i] = pkts[i-1] + recSize;
224 
225  // Read entire file into buffer
226  xfileStream.read( (char *) &pkts[0][0], xFileSize);
227  xfileStream.close();
228 
229  uint16_t ui16a;
230  double sec;
231  int16_t iy16, idy16;
232  int32_t iy, idy, jd, iyr, iday, jdtest;
233 
234  // Get start time of file to pad
235 
236  double secStart = get_pkt_sec( &pkts[0][toff], -1, &jdtest);
237  cout << "Start of file to pad: " << secStart << endl;
238 
239  double secEnd = get_pkt_sec( &pkts[nRec-1][toff], jdtest, NULL);
240  cout << "Stop of file to pad: " << secEnd << endl;
241 
242  //
243  // Read previous file
244  //
245  int iret = 0;
246  int iPad = 0;
247  uint8_t **pPkts;
248  if (pfile != 0x0) {
249  pfileStream.open( pfile, fstream::in | fstream::binary);
250  if ( pfileStream.fail()) {
251  cout << pfile << " not found" << endl;
252  exit(1);
253  }
254 
255  pPkts = new uint8_t*[nPad];
256  pPkts[0] = new uint8_t[recSize*nPad];
257  for (int i=1; i<nPad; i++) pPkts[i] = pPkts[i-1] + recSize;
258 
259  // pfileStream.seekg (0, ios::beg);
260  int i = 0;
261  // Find record just greater than start time of file to pad
262  while (1) {
263  pfileStream.read( (char *) &pPkts[0][0], recSize);
264  if (pfileStream.eof()) {
265  cout << "No padding from previous file within " << nSec2Pad <<
266  " seconds" << endl;
267  iPad = iPad | 1;
268  iret = 110;
269  break;
270  }
271 
272  double sec0 = get_pkt_sec( &pPkts[0][toff], jdtest, NULL);
273 
274  // cout << "prev: " << i << " " << sec0 << endl;
275  if (sec0 > (secStart-nSec2Pad - timeGap) && (sec0 < secStart)) {
276  pfileStream.seekg (-recSize*1, ios::cur);
277  pfileStream.read( (char *) &pPkts[0][0], recSize*nPad);
278  break;
279  }
280  i++;
281  }
282  pfileStream.close();
283  } // if (pfile != 0x0)
284 
285 
286  //
287  // Read following file
288  //
289  uint8_t **fPkts;
290  if (ffile != 0x0) {
291  ffileStream.open( ffile, fstream::in | fstream::binary);
292  if ( ffileStream.fail()) {
293  cout << ffile << " not found" << endl;
294  exit(1);
295  }
296 
297  fPkts = new uint8_t*[nPad];
298  fPkts[0] = new uint8_t[recSize*nPad];
299  for (int i=1; i<nPad; i++) fPkts[i] = fPkts[i-1] + recSize;
300 
301  ffileStream.seekg (0, ios::beg);
302  int i = 0;
303  // Find record just greater than start time of file to pad
304  while (1) {
305  ffileStream.read( (char *) &fPkts[0][0], recSize);
306 
307  double sec0 = get_pkt_sec( &fPkts[0][toff], jdtest, NULL);
308  // cout << "foll: " << i << " " << sec0 << endl;
309  if (sec0 > (secEnd+nSec2Pad + timeGap)) {
310  ffileStream.seekg (-recSize*(nPad+1), ios::cur);
311  ffileStream.read( (char *) &fPkts[0][0], recSize*nPad);
312 
313  break;
314  }
315  i++;
316  if (i == nRec) {
317  cout << "All times in following file less than endtime (plus pad) of file to be padded." << endl;
318  exit(120); // status 120 to be interpreted as - nothing to do - moving along with no padding done
319  }
320  }
321 
322  // Read nPad records of following file so that first record
323  // is right after end time
324  if (i == 0) {
325  iPad = iPad | 2;
326  iret = 110;
327  cout << "No padding from following file within " << nSec2Pad << " seconds"
328  << endl;
329  }
330 
331  ffileStream.close();
332  } // if (ffile != 0x0)
333 
334  if (ffile == 0x0 && (iPad & 1) == 1) {
335  delete[] pkts[0];
336  delete[] pkts;
337 
338  delete[] pPkts[0];
339  delete[] pPkts;
340 
341  return 120;
342  }
343 
344 
345  if (pfile == 0x0 && (iPad & 2) == 2) {
346  delete[] pkts[0];
347  delete[] pkts;
348 
349  delete[] fPkts[0];
350  delete[] fPkts;
351 
352  return 120;
353  }
354 
355 
356  //
357  // Open and write to output file
358  //
359  ofileStream.open( argv[optind+1], fstream::out | fstream::binary);
360 
361  if (pfile != 0x0 && (iPad & 1) == 0) {
362  cout << "Writing " << nPad << " records from previous file" << endl;
363  ofileStream.write( (char *) &pPkts[0][0], recSize*nPad);
364  }
365 
366  cout << "Writing " << nRec << " records from original file" << endl;
367  ofileStream.write( (char *) &pkts[0][0], recSize*nRec);
368 
369  if (ffile != 0x0 && (iPad & 2) == 0) {
370  cout << "Writing " << nPad << " records from following file" << endl;
371  ofileStream.write( (char *) &fPkts[0][0], recSize*nPad);
372  }
373 
374  ofileStream.close();
375 
376  delete[] pkts[0];
377  delete[] pkts;
378 
379  if (pfile != 0x0) {
380  delete[] pPkts[0];
381  delete[] pPkts;
382  }
383 
384  if (ffile != 0x0) {
385  delete[] fPkts[0];
386  delete[] fPkts;
387  }
388 
389 
390  // Check
391  uint8_t **oPkts;
392  int32_t oRec = nRec;
393  if (pfile != 0x0 && (iPad & 1) == 0) oRec += nPad;
394  if (ffile != 0x0 && (iPad & 2) == 0) oRec += nPad;
395  oPkts = new uint8_t*[oRec];
396  oPkts[0] = new uint8_t[recSize*oRec];
397  for (int i=1; i<oRec; i++) oPkts[i] = oPkts[i-1] + recSize;
398 
399  ofileStream.open( argv[optind+1], fstream::in | fstream::binary);
400  ofileStream.read( (char *) &oPkts[0][0], recSize*oRec);
401  ofileStream.close();
402 
403  //Days since 1/1/1958
404  uint16_t ui16;
405  uint32_t ui32;
406  memcpy( &ui16, (uint16_t *) &oPkts[0][toff], 2);
407  int32_t jd0 = SWAP_2( ui16) + 2436205;
408 
409  jdate( jd0, &iyr, &iday);
410 
411  // Loop through packets
412  ccsds_to_yds( &oPkts[0][toff], &iy, &idy, &sec);
413  iy16 = iy;
414  idy16 = idy;
415  jd = jday( iy16, 1, idy16);
416  double prevsec = sec + (jd - jd0)*86400;
417 
418  for (int i=0; i<oRec; i++) {
419  ccsds_to_yds( &oPkts[i][toff], &iy, &idy, &sec);
420  iy16 = iy;
421  idy16 = idy;
422  jd = jday( iy16, 1, idy16);
423  //cout << i << " " << sec + (jd - jd0)*86400 << endl;
424  if (fabs(sec + (jd - jd0)*86400 - prevsec) > gapCheck) {
425  cout << "Time gap in output file at record: " << i <<
426  " " << prevsec << " " << sec + (jd - jd0)*86400 << endl;
427  iret = 110;
428  }
429  prevsec = sec + (jd - jd0)*86400;
430  }
431 
432  delete[] oPkts[0];
433  delete[] oPkts;
434 
435  return iret;
436 }
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
int j
Definition: decode_rs.h:73
int32_t jday(int16_t i, int16_t j, int16_t k)
Definition: jday.c:4
function jd(i, j, k)
Definition: jd.f:2
#define NULL
Definition: decode_rs.h:63
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in out
Definition: HISTORY.txt:422
int ccsds_to_yds(uint8_t *cctime, int32_t *iyear, int32_t *iday, double *sec)
Definition: ccsds_to_yds.c:5
int main(int argc, char *argv[])
Definition: scpad_viirs.cpp:91
Definition: jd.py:1
double get_pkt_sec(uint8_t *pkt_time, int32_t in_jd, int32_t *out_jd)
Definition: scpad_viirs.cpp:66
#define basename(s)
Definition: l0chunk_modis.c:29
#define SWAP_2(x)
Definition: scpad_viirs.cpp:15
#define fabs(a)
Definition: misc.h:93
int jdate(int32_t julian, int32_t *year, int32_t *doy)
Definition: jdate.c:5
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
integer function julian(DAY, MONTH, YEAR)
Definition: julian.f:2
int32_t idy
Definition: atrem_corl1.h:161
int i
Definition: decode_rs.h:71
int32_t iyr
Definition: atrem_corl1.h:161
#define VERSION
Definition: scpad_viirs.cpp:17
int k
Definition: decode_rs.h:73