OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
l1agen_oci.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <string.h>
5 #include <unistd.h>
6 
7 #include <sstream>
8 #include <fstream>
9 #include <iomanip>
10 #include <getopt.h>
11 #include <libgen.h>
12 
13 #include "nc4utils.h"
14 #include "l1agen_oci.h"
15 
16 using namespace std;
17 using namespace netCDF;
18 using namespace netCDF::exceptions;
19 
20 int ccsds_sec_to_yds( uint8_t *cctime,
21  int32_t *iyear, int32_t *iday, double *sec) {
22 
23  uint32_t ui32;
24  uint32_t ccsec;
25 
26  memcpy( &ui32, cctime, 4);
27  ccsec = SWAP_4(ui32);
28 
29  double dsec = (double) ccsec;
30  int leap = leapseconds_since_1993( dsec);
31  ccsec -= (leap + 27);
32 
33  *iday = ccsec / 86400;
34  int32_t jday = *iday + 2436205; // Jan. 1, 1958 is Julian day 2436205
35  jdate(jday, iyear, iday);
36 
37  // Get milliseconds
38  int32_t msec = cctime[4]*256 + cctime[5];
39  int32_t isec = ccsec % 86400;
40  *sec = isec + msec/65536.0;
41 
42  return 0;
43 
44 }
45 
46 #define VERSION "1.05.00_2022-06-24"
47 
48 // Modification history:
49 // Programmer Organization Date Ver Description of change
50 // ---------- ------------ ---- --- ---------------------
51 // Joel Gales FutureTech 09/20/18 0.10 Original development
52 // based on IDL routines
53 // developed by F. Patt
54 // Joel Gales SAIC 10/26/18 0.20 Complete initial alpha version
55 // Joel Gales SAIC 12/20/18 0.30 Complete initial version
56 // of SWIR bands
57 // Joel Gales SAIC 04/26/19 0.40 Implement code changes from
58 // F. Patt
59 // Joel Gales SAIC 05/20/19 0.50 Implement code changes from
60 // F. Patt
61 // Joel Gales SAIC 06/04/19 0.60 Implement code changes from
62 // F. Patt
63 // Joel Gales SAIC 07/01/19 0.70 Add support for outlist
64 // Joel Gales SAIC 07/23/19 0.75 Implement code changes from
65 // F. Patt
66 // Joel Gales SAIC 08/02/19 0.80 Flag and ignore packets
67 // greater than 1200 bytes.
68 // Remove common routines and
69 // Place in common.cpp
70 // Joel Gales SAIC 08/09/19 0.81 Fix memory overwrite bug in
71 // unpack_ccd_packet() in the
72 // ossdata array. Initialize
73 // "lines" and "bands" arrays.
74 // Add support for granules with
75 // missing blue/red bands.
76 // Joel Gales SAIC 10/25/19 0.82 Exit if EOF before finding
77 // good packet
78 // Joel Gales SAIC 11/20/19 0.85 Implement code changes from
79 // F. Patt
80 // Joel Gales SAIC 11/25/19 0.86 Change 60 to maxsc fpr dspn
81 // comparison
82 // Joel Gales SAIC 11/27/19 0.87 Change L1A name to standard
83 // Trap granules with no ancillary
84 // granules
85 // Joel Gales SAIC 12/05/19 0.88 Add nametage command line
86 // option
87 // Joel Gales SAIC 01/03/20 0.90 Update and add additional
88 // telemetry fields
89 // Joel Gales SAIC 01/22/20 0.91 Add scomp/maxgap
90 // Joel Gales SAIC 01/23/20 0.92 Check for EOF when checking
91 // for zero science pixels
92 // Zero out sci/dark fields
93 // Joel Gales SAIC 01/28/20 0.93 Bug fixes and updates
94 // Joel Gales SAIC 02/07/20 0.94 Implement changes from F.Patt
95 // (020720)
96 // Joel Gales SAIC 02/20/20 0.95 Fixed bugs with bagerr, ragerr,
97 // digerr and dautemp
98 // Joel Gales SAIC 02/25/20 0.96 Fixed bug in ccsds..
99 // Writing 6 bytes into uint32_t
100 // Joel Gales SAIC 03/04/20 0.97 Implement changes from F.Patt
101 // (022820)
102 // Joel Gales SAIC 03/24/20 0.98 Fix HAM_side issue
103 // Joel Gales SAIC 04/02/20 0.99 Implement SWIR mode changes
104 // Liang Hong SAIC 04/28/20 0.9901 Fixed last scan_line_attributes
105 // records in last output granule
106 // Liang Hong SAIC 05/20/20 0.9902 Handle dark zone science data
107 // Liang Hong SAIC 08/11/20 0.9903 Handle 0 pix 0 band sci/cal data in l1a
108 // "no data" and order variation in input
109 // Liang Hong SAIC 08/24/20 0.9904 fixed a bug to close file with 0-scan
110 // ; demand non-negative indices
111 // Liang Hong SAIC 09/02/20 0.9905 HKT and science data overlap check; order
112 // SWIR bands in ascending wavelength
113 // Liang Hong SAIC 10/28/20 0.9906 handle fill values in SWIR
114 // Liang Hong SAIC 10/29/20 0.9907 APID for the MCE HK packet changed to 713
115 // Liang Hong SAIC 11/23/20 0.9908 fixed rare start/end time error; SWIR band
116 // -specific pixel shifts as input option;
117 // run with optional granule start time;
118 // fixed science packet sequence error flag
119 // fixed HAM_side value after index nmce
120 // Liang Hong SAIC 12/01/20 0.99.00 fixed number_of_filled_scans in metadata
121 // Liang Hong SAIC 04/22/21 0.99.10 fixed no ancillary data exit conditions
122 // Liang Hong SAIC 06/17/21 0.99.20 generate 1 telemetry entry when no data
123 // bug fix in duplicated reading of last tlm
124 // Liang Hong SAIC 07/27/21 0.99.21 return 110 for No ancillary packets
125 // Liang Hong SAIC 01/07/22 1.00.00 temperature fields update in HKT packets
126 // SWAP_4 in common.h updated
127 // Liang Hong SAIC 01/11/22 1.00.01 OCI SWIR Raw mode reading correction
128 // Liang Hong SAIC 01/25/22 1.00.11 blue and red spectral mode order correction
129 // Liang Hong SAIC 03/11/22 1.01.00 added telemetry for the solar calibrator;
130 // references in metadata; write ancillary_tlm
131 // Liang Hong SAIC 03/30/22 1.02.00 SPCA and lunar stare data types added
132 // Liang Hong SAIC 04/14/22 1.03.00 update error checking and handling
133 // Liang Hong SAIC 04/21/22 1.03.01 updated packet # threshold; mode table check
134 // Liang Hong SAIC 04/24/22 1.03.02 exit read packets when endfile
135 // Liang Hong SAIC 05/11/22 1.04.00 update of APID list and check all for spin #
136 // added navigation data from hkt input
137 // Liang Hong SAIC 05/27/22 1.04.01 limit packet reading to maxpkts
138 // Liang Hong SAIC 06/24/22 1.05.00 added CCD masks; fixed bugs in nav data
139 
140 ofstream tempOut;
141 
142 int main (int argc, char* argv[])
143 {
144 
145  cout << "l1agen_oci " << VERSION << " ("
146  << __DATE__ << " " << __TIME__ << ")" << endl;
147 
148  if ( argc == 1) {
149  cout << endl <<
150  "l1agen_oci OCI_packet_file granule_len " <<
151  "[-k SC_HKT_input_list] " << // file containing list of Spacecraft housekeeping telemetry file names
152  "[--hktlist SC_HKT_input_list] "
153  "[-s SWIR_LOFF_config] " <<
154  "[--swir_loff_set] " <<
155  "[-t YYYYmmddTHH:MM:SS] " <<
156  "[--start_time YYYYmmddTHH:MM:SS] " <<
157  "[-o output_list_file] " <<
158  "[--outlist output_list_file]" <<
159  endl;
160  return 0;
161  }
162 
163  int c;
164  string hktlist = "";
165  string swir_loff_set = "";
166  string time_start = "";
167  string outlist = "";
168  string nametag = "";
169  while (1) {
170  static struct option long_options[] = {
171  {"hktlist", required_argument, 0, 'k'},
172  {"swir_loff_set", required_argument,0,'s'},
173  {"start_time", required_argument, 0, 't'},
174  {"outlist", required_argument, 0, 'o'},
175  {"nametag", required_argument, 0, 'p'},
176  {0, 0, 0, 0}
177  };
178 
179  /* getopt_long stores the option index here. */
180  int option_index = 0;
181 
182  c = getopt_long( argc, argv, "k:s:t:o:p:", long_options, &option_index);
183 
184  /* Detect the end of the options. */
185  if (c == -1)
186  break;
187 
188  switch (c)
189  {
190  case 'k':
191  hktlist.assign(optarg);
192  break;
193 
194  case 's':
195  swir_loff_set.assign(optarg);
196  break;
197 
198  case 't':
199  time_start.assign(optarg);
200  break;
201 
202  case 'o':
203  // printf ("option -o with value `%s'\n", optarg);
204  outlist.assign( optarg);
205  break;
206 
207  case 'p':
208  nametag.assign( optarg);
209  break;
210 
211  default:
212  abort ();
213  }
214  }
215 
216  ofstream fout;
217  if ( outlist.compare("") != 0)
218  fout.open( outlist.c_str());
219 
220  fstream tfileStream;
221 
222  //Read S/C telemetry data
223  //read_pace_telemetry(hktlist,iyrs,idays,otime,pos,vel,atime,quat,arate,ttime,tilt);
224 
225  // OCI packet file
226  tfileStream.open( argv[optind+0], fstream::in | fstream::binary);
227  if ( tfileStream.fail()) {
228  cout << argv[optind+0] << " not found" << endl;
229  exit(1);
230  }
231 
232  uint32_t apid = 0;
233  uint32_t len = 0;
234  int32_t endfile = 0;
235  uint8_t fpacket[PKTSIZE];
236  uint8_t apacket[ANCSIZE];
237  uint8_t apacket0[ANCSIZE];
238  uint32_t maxpkts = 30000; // LH, 4/21/2022, v1.03.01
239  int acomp;
240  int maxgap = 10;
241 
242  uint8_t **pbuffer0 = new uint8_t *[maxpkts];
243  pbuffer0[0] = new uint8_t[PKTSIZE*maxpkts];
244  for (size_t i=1; i<maxpkts; i++) pbuffer0[i] = pbuffer0[i-1] + PKTSIZE;
245 
246  uint8_t **pbuffer1 = new uint8_t *[maxpkts];
247  pbuffer1[0] = new uint8_t[PKTSIZE*maxpkts];
248  for (size_t i=1; i<maxpkts; i++) pbuffer1[i] = pbuffer1[i-1] + PKTSIZE;
249 
250  uint32_t npkts0;
251  int32_t ancind0;
252  // uint32_t spn0, spn;
253  int32_t spn0, spn, spnp;
254  vector<int32_t> tlmind0;
255 
256 
257  // Get first science or ancillary packet
258  read_packet( &tfileStream, NULL, len, apid, endfile);
259  if (len > PKTSIZE) {
260  cout << "Packet too big (" << len << ") for buffer (" << PKTSIZE << ")"
261  << endl;
262  exit(1);
263  }
264  read_packet( &tfileStream, fpacket, len, apid, endfile);
265  apid = (fpacket[0] % 8)*256 + fpacket[1];
266  int nsk=0;
267  while (apid != 636 && apid != 700 && apid != 720 && !endfile) {
268  nsk++;
269  read_packet( &tfileStream, NULL, len, apid, endfile);
270  if (len > PKTSIZE) {
271  cout << "Packet too big (" << len << ") for buffer (" << PKTSIZE << ")"
272  << endl;
273  exit(1);
274  }
275  read_packet( &tfileStream, fpacket, len, apid, endfile);
276  apid = (fpacket[0] % 8)*256 + fpacket[1];
277  }
278  if (endfile) {
279  cout << "No science packets found in file" << endl;
280  exit(1);
281  }
282  if (nsk > 0) cout << nsk << " packets skipped" << endl;
283 
284 
285  // Read first scan and check for ancillary packet
286  ancind0 = -1;
287  tlmind0.clear();
288 
289  int32_t iyear, iday;
290  double stime;
291 
292  // Get granule period in minutes
293  int32_t mper;
294  string str = argv[optind+1];
295  istringstream( str) >> mper;
296 
297  itab itable[10];
298  string dtypes[] = {"", "", "_DARK", "_SOL-D", "_SOL-M", "_LIN", "_LUN",
299  "_DIAG", "_STAT", "_SPEC", "","_SNAP-X", "_SNAP-I",
300  "SPCA","LUN-ST"};
301  string smodes[] = {"", "_SDIAG", "_SRAW", "_STEST"};
302 
303  uint8_t **pbuffer = new uint8_t *[maxpkts];
304  pbuffer[0] = new uint8_t[PKTSIZE*maxpkts];
305  for (size_t i=1; i<maxpkts; i++) pbuffer[i] = pbuffer[i-1] + PKTSIZE;
306 
307  uint16_t ncps, nbbs, nrbs, nsps, ndcs, ndss, btaps[16], rtaps[16], msps;
308 
309  uint8_t *linerr = new uint8_t[maxpkts];
310  uint8_t *seqerr = new uint8_t[maxpkts];
311  for (size_t i=1; i<maxpkts; i++) {
312  linerr[i] = 0;
313  seqerr[i] = 0;
314  }
315  uint8_t noseq=255;
316  uint16_t smode = 0;
317  uint16_t sstop = 0;
318 
319  while ((ancind0 == -1) && !endfile) {
320  read_oci_scan_packets( &tfileStream, fpacket,
321  (uint8_t (*)[PKTSIZE]) &pbuffer0[0][0],
322  npkts0, spn0, ancind0, tlmind0, noseq, endfile);
323  }
324  if ( endfile) {
325  cout << "No ancillary packets found in file" << endl;
326  exit(110); // ver 0.99.21
327  }
328 
330  while (!endfile) {
331 
332  // ncps - number of CCD band pixels (ccd_pixels)
333  // nsps - number of SWIR band pixels
334  // ndcs - number of dark collect pixels
335  // ndss - number of dark SWIR pixels
336  // nbbs - number of blue bands
337  // nrbs - number of red bands
338  // btaps - Blue CCD tap enable flags (1 = enabled)
339  // rtaps - Red CCD tap enable flags (1 = enabled)
340 
341  // Check for zero science pixels
342  if (ancind0 != -1) {
343  memcpy( apacket0, &pbuffer0[ancind0][0], ANCSIZE);
344  get_band_dims( apacket0, ncps, nbbs, nrbs, nsps, ndcs, ndss,
345  btaps, rtaps, itable);
346  }
347 
348  while ((ncps == 1 || ancind0 == -1) && !endfile) {
349  if (ancind0 != -1) cout<<"Ancillary packet has zero science pixels at spin "<<spn0<<endl;
350  read_oci_scan_packets( &tfileStream, fpacket,
351  (uint8_t (*)[PKTSIZE]) &pbuffer0[0][0],
352  npkts0, spn0, ancind0, tlmind0, noseq, endfile);
353  if (ancind0 != -1) {
354  memcpy( apacket0, &pbuffer0[ancind0][0], ANCSIZE);
355  get_band_dims( apacket0, ncps, nbbs, nrbs, nsps, ndcs, ndss,
356  btaps, rtaps, itable);
357  }
358  }
359 
360  if (ancind0 == -1) {
361  cout << "No ancillary packets for last granule" << endl;
362  break;
363  }
364 
365  cout << endl;
366  cout << "ncps - number of CCD band pixels (ccd_pixels): " << ncps << endl;
367  cout << "nsps - number of SWIR band pixels: " << nsps << endl;
368  cout << "ndcs - number of dark collect pixels: " << ndcs << endl;
369  cout << "ndss - number of dark SWIR pixels: " << ndss << endl;
370  cout << "nbbs - number of blue bands: " << nbbs << endl;
371  cout << "nrbs - number of red bands: " << nrbs << endl;
372 
373  get_anc_packet_time( apacket0, iyear, iday, stime);
374 
375  double scanp = 0.1755; // changed from 1.0 / 5.737 to 0.1755 in v1.04.00
376  double stimp = stime - scanp;
377 
378  time_struct starttime, endtime;
379  starttime.iyear = iyear;
380  starttime.iday = iday;
381  starttime.sec = stime;
382 
383  uint32_t ltime, mtime;
384  uint16_t maxsc;
385 
386  // optional input start time for first L1A granule
387  // LH, 11/20/2020
388  struct tm tm_start;
389  if (!time_start.empty()) {
390  strptime(time_start.c_str(), "%Y-%m-%dT%H:%M:%S", &tm_start);
391  ltime = (int32_t)(tm_start.tm_hour*3600+tm_start.tm_min*60+tm_start.tm_sec);
392  mtime = mper>0?ltime + mper*60:ltime + 600;
393  while (((stime <ltime) || (ncps==1)) && !endfile) {
394  read_oci_scan_packets( &tfileStream, fpacket,
395  (uint8_t (*)[PKTSIZE]) &pbuffer0[0][0],
396  npkts0, spn0, ancind0, tlmind0, noseq, endfile);
397  if (ancind0!=-1) {
398  memcpy( apacket0, &pbuffer0[ancind0][0], ANCSIZE);
399  get_anc_packet_time( apacket0, iyear, iday, stime);
400 
401  get_band_dims( apacket0, ncps, nbbs, nrbs, nsps, ndcs, ndss,
402  btaps, rtaps, itable);
403  if (ncps==1) cout<<"Ancillary packet has zero science pixels at spin "<<spn0<<endl;
404  }
405  }
406 
407  if (endfile || (stime>mtime)) {
408  cout << "No science data in time range" << endl;
409  exit(1);
410  }
411 
412  starttime.sec = stime;
413 
414  sstop = 1;
415  }
416 
417  if (mper > 0) {
418  ltime = (((int32_t) (stime)) / 60 / mper) * (mper*60);
419  mtime = ltime + mper*60;
420  maxsc = (uint16_t) ((mper*60/scanp) + 2);
421  } else {
422  cout << "Processing with single file option" << endl;
423  ltime = (int32_t) stime;
424  mtime = ltime + 600;
425  maxsc = 3600;
426  }
427 
428  acomp = 1;
429 
430  // number of SWIR bands
431  nsps = ((nsps+7)/8)*8; //Round up SWIR number of pixels to a multiple of 8, LH, 4/15/2022, v1.03.00
432  unsigned short nswb = 9;
433 
434  int16_t cindex[32768];
435  int16_t sindex[32768*nswb];
436  int16_t cdindex[32768];
437  int16_t sdindex[32768];
438  int16_t swir_loff[9];
439 
440  for (size_t i=0; i<32768; i++) {
441  cindex[i] = -1;
442  for (size_t j=0; j<nswb; j++) sindex[i*nswb+j] = -1;
443  cdindex[i] = -1;
444  sdindex[i] = -1;
445  }
446 
447  for (size_t i=0; i<9; i++) swir_loff[i] = 0;
448  if (swir_loff_set.compare("ETU")==0){
449  memcpy(swir_loff, SWIR_LOFF_ETU, 9*sizeof(uint16_t));
450  }
451  //int32_t ldark = 32768;
452 
453  //if (stime < 15299.458) {
454  // cout << "make_oci_line_index" << endl;
455  // make_oci_line_index( itable, cindex, sindex, ldark);
456  //}
457  // make_oci_line_index( itable, cindex, sindex, ldark);
458  make_oci_line_index( itable, cindex, sindex, cdindex, sdindex,swir_loff);
459 
460 
461 
462  // Get SWIR band data mode
463  // uint16_t smode = 0; // LH, 09/10/2020
464  get_swir_mode( (uint8_t (*)[PKTSIZE]) &pbuffer0[0][0], npkts0, smode);
465  uint16_t smodep = smode;
466  int scomp = 1;
467 
468  uint16_t cdsmode;
469 
470  // Determine start and end time of granule
471  uint32_t jd0 = jday(iyear, 1, iday);
472  int16_t yr16 = (int16_t) iyear;
473  int16_t doy = (int16_t) iday;
474  int16_t month, dom;
475  yd2md( yr16, doy, &month, &dom);
476 
477  int32_t ih = (int32_t) (stime / 3600);
478  int32_t mn = (int32_t) ((stime - ih*3600) / 60);
479  int32_t isec = (int32_t) (stime - ih*3600 - mn*60);
480 
481  stringstream timestr, datestr;
482  timestr << setfill('0') << setw(2) << ih
483  << setfill('0') << setw(2) << mn
484  << setfill('0') << setw(2) << isec;
485 
486  datestr << setfill('0') << setw(4) << iyear
487  << setfill('0') << setw(2) << month
488  << setfill('0') << setw(2) << dom;
489 
490  static l1aFile outfile;
491  string basenme;
492  basenme.assign(basename(argv[optind+0]));
493 
494  if (nametag == "") nametag.assign("PACE_OCI");
495  short dtype = itable[1].dtype; // ; Get data type for file name and metadata
496  short maxdtype = 2;
497  for (size_t i=0; i<10; i++) {
498  if (itable[i].dtype>maxdtype) maxdtype = itable[i].dtype;
499  }
500  // if (maxdtype > 2) dtype = maxdtype;
501  if ((maxdtype != 2) && (maxdtype != 10)) dtype = maxdtype; // LH, 8/24/2020
502 
503  // data type mod for ETU before June 2020
504  if ((jd0 < 2459000) && dtype == 11) dtype = 9;
505 
506  string l1a_name = nametag + dtypes[dtype] + smodes[smode] +
507  "_" + basenme.substr(0,4) + basenme.substr(5,3) + basenme.substr(9,3);
508  l1a_name += "." + datestr.str() + "T" + timestr.str() + ".L1A.nc";
509 
510  // Initialize data arrays
511 
512  // blue band dark collect data for granule
513  uint16_t **dark_b = new uint16_t *[maxsc];
514  dark_b[0] = new uint16_t[ndcs*(nbbs>0?nbbs:1)*maxsc];
515  for (size_t i=1; i<maxsc; i++) dark_b[i] = dark_b[i-1] + ndcs*(nbbs>0?nbbs:1);
516 
517  // red band dark collect data for granule
518  uint16_t **dark_r = new uint16_t *[maxsc];
519  dark_r[0] = new uint16_t[ndcs*(nrbs>0?nrbs:1)*maxsc];
520  for (size_t i=1; i<maxsc; i++) dark_r[i] = dark_r[i-1] + ndcs*(nrbs>0?nrbs:1);
521 
522  // swir band dark collect data for granule
523  uint32_t **dark_s = new uint32_t *[maxsc];
524  dark_s[0] = new uint32_t[ndss*(nswb>0?nswb:1)*maxsc];
525  for (size_t i=1; i<maxsc; i++) dark_s[i] = dark_s[i-1] + ndss*(nswb>0?nswb:1);
526 
527  uint16_t **bdark = new uint16_t *[nbbs];
528  uint16_t **rdark = new uint16_t *[nrbs];
529  uint32_t **sdark = new uint32_t *[nswb];
530 
531  uint16_t **bsci = new uint16_t *[nbbs];
532  bsci[0] = new uint16_t[ncps*(nbbs>0?nbbs:1)];
533  for (size_t i=1; i<nbbs; i++) bsci[i] = bsci[i-1] + ncps;
534 
535  uint16_t **rsci = new uint16_t *[nrbs];
536  rsci[0] = new uint16_t[ncps*(nrbs>0?nrbs:1)];
537  for (size_t i=1; i<nrbs; i++) rsci[i] = rsci[i-1] + ncps;
538 
539  uint32_t **ssci = new uint32_t *[nswb];
540  ssci[0] = new uint32_t[nsps*(nswb>0?nswb:1)];
541  for (size_t i=1; i<nswb; i++) ssci[i] = ssci[i-1] + nsps;
542 
543 
544  uint8_t **ancdata = new uint8_t *[maxsc+1];
545  ancdata[0] = new uint8_t[ANCSIZE*(maxsc+1)];
546  for (size_t i=1; i< (size_t) (maxsc+1); i++)
547  ancdata[i] = ancdata[i-1] + ANCSIZE;
548 
549  uint32_t maxtlm = maxsc*10;
550  uint8_t **tlmdata = new uint8_t *[maxtlm];
551  tlmdata[0] = new uint8_t[TLMSIZE*(maxtlm)];
552  for (size_t i=1; i< (size_t) (maxtlm); i++)
553  tlmdata[i] = tlmdata[i-1] + TLMSIZE;
554 
555  // uint16_t ncpt = ncps + ndcs;
556  // uint16_t nspt0 = nsps + ndss;
557  uint16_t ncpt = ncps; // Dark pixels now included in science pixels
558  uint16_t nspt0 = nsps; // Dark pixels now included in science pixels
559 
560  // Round up SWIR number of pixels to a multiple of 8
561  uint16_t nspt = ((nspt0 + 7) / 8) * 8;
562 
563  // Note: "bands" arrays here are the reverse order as IDL versions
564  uint16_t **bbands = new uint16_t*[ncpt];
565  uint16_t **rbands = new uint16_t*[ncpt];
566  for (size_t i=0; i<ncpt; i++) {
567  bbands[i] = NULL;
568  rbands[i] = NULL;
569  }
570 
571  int16_t *blines = new int16_t[ncpt];
572  int16_t *rlines = new int16_t[ncpt];
573 
574  for (size_t i=0; i<ncpt; i++) blines[i] = -1;
575  for (size_t i=0; i<ncpt; i++) rlines[i] = -1;
576 
577 
578  msps = 4096; // May be full-scan SWIR pixels in ETU data
579  uint32_t **sbands = new uint32_t*[msps];
580  for (size_t i=0; i<msps; i++) {
581  sbands[i] = new uint32_t[nswb];
582  }
583 
584  int16_t *slines = new int16_t[msps];
585  // for (size_t i=0; i<msps; i++) slines[i] = 0;
586 
587  // swir frame type for dark collect
588  int8_t **sdfrms = new int8_t *[maxsc];
589  sdfrms[0] = new int8_t[ndss*maxsc];
590  for (size_t i=1; i<maxsc; i++) sdfrms[i] = sdfrms[i-1] + ndss;
591  memset(sdfrms[0], -1, sizeof(int8_t)*ndss*maxsc);
592 
593  int8_t *sfrm = new int8_t[msps];
594  int8_t *sfrms = new int8_t[msps];
595 
596  // Initialize output file and get object IDs for EV data
597  cout << "Creating: " << l1a_name.c_str() << endl;
598  cout << "Starting at spin number "<<spn0<<endl;
599  outfile.createl1( (char *) l1a_name.c_str(),
600  maxsc, ncps, nbbs, nrbs, nsps, ndcs);
601 
602  // Write output filename to outlist if needed.
603  long outlistpos;
604  if ( outlist.compare("") != 0) {
605  outlistpos = fout.tellp();
606  fout << l1a_name.c_str() << " ";
607  }
608 
609  // Read and process OCI scans
610  int iret = 0;
611  int icheck = 0;
612  uint32_t isc = 0;
613  uint32_t npkts,npkt1;
614  int32_t ancind;
615  vector<int32_t> tlmind;
616  spnp = spn0;
617  int32_t enddata = 0;
618  int dspn = 1;
619  uint32_t ntlm = 0;
620 
621  uint16_t btype, rtype;
622  uint16_t bagg[16], ragg[16];
623  uint16_t nbands;
624 
625  // tempOut.open ("bsci.bin", ios::out | ios::trunc | ios::binary);
626 
627  // ---------------------------------------------------------------------
628  // ---------------------------------------------------------------------
629  // ---------------------------------------------------------------------
630  while ( stime < mtime && acomp && !enddata && scomp && (dspn <= maxgap) && isc<maxsc) {
631  // sometimes, data contain more scans than possible, set isc<maxsc to avoid overflow, LH 8/26/2020
632  if ((isc % 100) == 0) cout << "Processing scan " << isc << endl;
633  // cout << "Processing scan " << isc << endl;
634  //cout << tfileStream.tellg() << endl;
635 
636  memcpy( &pbuffer1[0][0], &pbuffer0[0][0], PKTSIZE*maxpkts);
637  ancind = ancind0;
638  tlmind = tlmind0;
639  npkt1 = npkts0;
640  spn = spn0;
641  enddata = endfile;
642 
643  // Read next scan
644  read_oci_scan_packets( &tfileStream, fpacket,
645  (uint8_t (*)[PKTSIZE]) &pbuffer0[0][0],
646  npkts0, spn0, ancind0, tlmind0, seqerr[isc],
647  endfile);
648 
649  // Disabled maxsc check for I&T
650  if (dspn >= 1 && npkt1 > 1) { // ver 1.03.00
651 
652  // Save ancillary packet in array
653  if (ancind != -1) {
654  memcpy(ancdata[isc], pbuffer1[ancind], ANCSIZE);
655  } else {
656  // insert spin # if missing anc packet
657  uint32_t ui32;
658  ui32 = SWAP_4(spn);
659  memcpy(&ancdata[isc][24], &ui32, 4);
660 
661  // set apid to 0 where there is missing anc packet, LH 4/28/2020
662  ui32 = SWAP_4(0);
663  memcpy(&ancdata[isc][0], &ui32, 4);
664  }
665 
666  // Save telemetry packet in array
667  int ntind = tlmind.size();
668  if ( ntind > 0 && ntlm < maxtlm) {
669  int itt = 0;
670  while (itt < ntind && ntlm < maxtlm) {
671  memcpy(tlmdata[ntlm], pbuffer1[tlmind[itt]], TLMSIZE);
672  ntlm++;
673  itt++;
674  }
675  if (ntlm >= maxtlm)
676  cout << "Maximum number of telemetry packets exceeded at spin: "
677  << spn << endl;
678  }
679 
680  // Zero out slines, sfrm, and sbands arrays
681  for (size_t i=0; i<msps; i++) {
682  slines[i] = -1;
683  sfrm[i] = 0;
684  for (size_t j=0; j<nswb; j++) sbands[i][j] = 0;
685  }
686 
687  // Unpack science data from packets
688  npkts = npkt1 + npkts0;
689  if (npkts>maxpkts) { // LH, 5/27/2022, ver 1.04.01
690  cout<<"Packets exceed max ["<<maxpkts<<"]."<<endl;
691  npkts = maxpkts;
692  npkts0 = npkts - npkt1;
693  }
694  memcpy( &pbuffer[0][0], &pbuffer1[0][0], PKTSIZE*npkt1);
695  if (npkts0 > 0) memcpy( &pbuffer[npkt1][0], &pbuffer0[0][0], PKTSIZE*npkts0); // pbuffer[*,npkt1:npkts-1] = pbuffer2[*,0:npkt2-1]
696 
697  unpack_oci_sci( npkts, spn, ncpt, nspt, msps, nbands, btaps, rtaps,
698  (uint8_t (*)[PKTSIZE]) &pbuffer[0][0],
699  bbands, rbands, sbands, blines, rlines, slines,
700  btype, bagg, rtype, ragg, sfrm, icheck);
701  if (icheck != 0) cout << "Science data unpacking error in spin: "
702  << spn << endl;
703  iret = iret | icheck;
704 
705  bdark[0] = dark_b[isc];
706  for (size_t i=1; i<nbbs; i++) bdark[i] = bdark[i-1] + ndcs;
707 
708  rdark[0] = dark_r[isc];
709  for (size_t i=1; i<nrbs; i++) rdark[i] = rdark[i-1] + ndcs;
710 
711  sdark[0] = dark_s[isc];
712  for (size_t i=1; i<nswb; i++) sdark[i] = sdark[i-1] + ndss;
713 
714  // Initialize science/dark array with fill value
715  std::fill(bsci[0], bsci[0]+ncps*((nbbs>0)?nbbs:1), 65535);
716  std::fill(rsci[0], rsci[0]+ncps*((nrbs>0)?nrbs:1), 65535);
717  std::fill(ssci[0],ssci[0]+nsps*((nswb>0)?nswb:1), 1048575);
718  memset(sfrms, -1, sizeof(int8_t)*msps);
719 
720  // Initialize with fill value, LH 8/6/2020
721  std::fill(bdark[0],bdark[0]+ndcs*((nbbs>0)?nbbs:1), 65535);
722  std::fill(rdark[0],rdark[0]+ndcs*((nrbs>0)?nrbs:1), 65535);
723  std::fill(sdark[0],sdark[0]+ndss*((nswb>0)?nswb:1), 1048575);
724 
725  // If science data in scan, check data for gaps or inconsistencies and determine data types
726  if ((blines[0] != -1) || (rlines[0] != -1) || (slines[0] != -1)) {
727  check_load_oci_data( dtype, ncpt, nsps, ndcs, ndss, nbbs, nrbs, nswb, // msps -> nsps, LH, ver 1.03.00
728  cindex, sindex, cdindex,sdindex, bbands, rbands, sbands,
729  blines, rlines, slines, bsci, rsci, ssci,
730  bdark, rdark, sdark, linerr[isc], icheck);
731  if (icheck >= 2) cout << "Science data checking error in spin: "
732  << spn << endl;
733  iret = iret | icheck;
734 
735  // Store SWIR band frame types
736  for (size_t i=0; i<msps; i++) {
737  if (sindex[slines[i]*nswb] > -1) { // LH, 8/24/2020; use SWIR band0, LH 11/20/2020
738  sfrms[sindex[slines[i]*nswb]] = sfrm[i];
739  }
740  }
741 
742  for (size_t i=0; i<msps; i++) {
743  if (sdindex[slines[i]] > -1) { // LH, 8/24/2020
744  sdfrms[isc][sdindex[slines[i]]] = sfrm[i];
745  }
746  }
747 
748  // Write science data to file
749  outfile.write_oci_science_data( isc, nbbs, nrbs, nswb,
750  ncps, nsps, bsci, rsci, ssci, sfrms);
751 
752  // Check for instrument HKT packets in scan (placeholder for now)
753 
754  isc++;
755  stimp = stime;
756  endtime.iyear = iyear;
757  endtime.iday = iday;
758  endtime.sec = stime;
759  spnp = spn;
760  }
761  } else {
762  ih = (int32_t) (stime / 3600);
763  mn = (int32_t) ((stime - ih*3600) / 60);
764  isec = (int32_t) (stime - ih*3600 - mn*60);
765 
766  if(stime <= stimp && ancind != -1)
767  cout << "Scan " << spn << " out of order at" <<
768  ih << " " << mn << " " << isec << endl;
769  } // if (dspn >= 1 && npkts >
770 
771  // Get scan time and band dimensions for next scan
772  if (!enddata && ancind0 != -1) {
773  memcpy( apacket, &pbuffer0[ancind0][0], ANCSIZE);
774  get_anc_packet_time( apacket, iyear, iday, stime);
775  uint32_t jd = jday(iyear, 1, iday);
776  stime += (jd - jd0) * 86400;
777 
778  get_swir_mode( (uint8_t (*)[PKTSIZE]) &pbuffer0[0][0], npkts0, smode);
779 
780  scomp = (smode == smodep);
781 
782  if ( !scomp) {
783  ih = (int32_t) (stime / 3600);
784  mn = (int32_t) ((stime - ih*3600) / 60);
785  isec = (int32_t) (stime - ih*3600 - mn*60);
786  cout << "SWIR data mode change at: "
787  << ih << " " << mn << " " << isec << endl;
788  }
789  acomp = anc_compare( apacket0, apacket);
790  // cout << "acomp: " << acomp << endl;
791  }
792 
793  dspn = spn0 - spnp;
794  if (dspn > maxgap)
795  cout << "Spin number gap at spin " << spnp << endl;
796 
797  } // while ( stime < mtime && acomp && !enddata)
798  // ---------------------------------------------------------------------
799  // ---------------------------------------------------------------------
800  // ---------------------------------------------------------------------
801 
802  if (mper > 0 && isc < (size_t) (maxsc-10) && acomp) {
803  cout << "(mper > 0 && isc < (size_t) (maxsc-10) && acomp)" << endl;
804  iret = iret | 1;
805  }
806  cout << "Scans in file: " << isc << endl;
807  cout << "Complete flag: " << iret << endl;
808 
809  if (isc > 0) {
810 
811  // Need to include ancillary packet from next scan
812  if (ancind0 != -1) memcpy(ancdata[isc], apacket, ANCSIZE);
813 
814  // Write calibration data to file
815  outfile.write_oci_cal_data( isc, nbbs, nrbs, nswb, ndcs, ndss,
816  dark_b[0], dark_r[0], dark_s[0],
817  sdfrms[0]);
818 
819  // Get scan metadata and write to file
820  // Save spid id computed in this function for later use
821  int32_t *spinID = new int32_t[isc+100];
822  for (size_t i=0; i<(isc+100); i++) spinID[i] = -999;
823  outfile.write_oci_scan_metadata( isc, ancdata[0], seqerr, linerr, spinID);
824 
825  // Write ancillary data to file
826  outfile.write_oci_ancil_data( isc, ancdata[0]);
827 
828  // Unpack OCI telemetry data and write to file
829  cdsmode = 0;
830  int ntind = tlmind0.size();
831  if (tlmind0.size() != 0 && dspn <= maxgap && !endfile) { // 0.99.20
832  int itt = 0;
833  while (itt < ntind && ntlm < maxtlm) {
834  memcpy(tlmdata[ntlm], pbuffer0[tlmind0[itt]], TLMSIZE);
835  ntlm++;
836  itt++;
837  }
838  }
839 
840  if (ntlm > 0) {
841  cout << ntlm << " HKT packets" << endl;
842  outfile.write_oci_tlm_data( itable, ntlm,
843  (uint8_t (*)[TLMSIZE]) &tlmdata[0][0],
844  spinID, cdsmode,isc);
845  }
846  delete[] spinID;
847 
848  // Locate navigation data and write to file
849  if ( hktlist.compare("") != 0)
850  outfile.write_navigation(hktlist,starttime,endtime);
851 
852  // if ((endfile && mper <= 0) || !acomp)
853  // mtime = (uint32_t) floor(stimp);
854 
855  // Generate granule metadata and write to file
856  string sdir = "Ascending"; // Hard-code until we have spacecraft data
857  string edir = "Ascending"; // Hard-code until we have spacecraft data
858  outfile.write_oci_global_metadata( starttime, endtime, l1a_name,
859  sdir, edir, isc, dtype, //itable[1].dtype changed to dtype
860  smode, cdsmode, fout);
861 
862  // Write complete flag to outlist if needed
863  if ( outlist.compare("") != 0)
864  fout << iret << "\n";
865 
866  // Close file
867  outfile.close();
868 
869  } else {
870  // Remove 0-scan file
871  outfile.close(); // LH, 08/19/2020
872  int status = remove( l1a_name.c_str());
873  if( status == 0) {
874  cout << "Removing 0-scan file: " << l1a_name.c_str() << endl;
875  if ( outlist.compare("") != 0) fout.seekp( outlistpos);
876  } else {
877  cout << "Error removing " << l1a_name.c_str() << endl;
878  exit(1);
879  }
880  }
881 
882  delete[] dark_b[0];
883  delete[] dark_b;
884  delete[] dark_r[0];
885  delete[] dark_r;
886  delete[] dark_s[0];
887  delete[] dark_s;
888 
889  delete[] bdark;
890  delete[] rdark;
891  delete[] sdark;
892 
893  delete[] bsci[0];
894  delete[] bsci;
895  delete[] rsci[0];
896  delete[] rsci;
897  delete[] ssci[0];
898  delete[] ssci;
899 
900  for (size_t i=0; i<ncps; i++) if ( bbands[i] != NULL) delete[] bbands[i];
901  delete[] bbands;
902  for (size_t i=0; i<ncps; i++) if ( rbands[i] != NULL) delete[] rbands[i];
903  delete[] rbands;
904  for (size_t i=0; i<msps; i++) delete[] sbands[i];
905  delete[] sbands;
906 
907  delete[] blines;
908  delete[] rlines;
909  delete[] slines;
910 
911  delete[] sdfrms[0];
912  delete[] sdfrms;
913 
914  delete[] sfrm;
915  delete[] sfrms;
916 
917  delete[] ancdata[0];
918  delete[] ancdata;
919 
920  if (stime>mtime && sstop==1) break;
921  } // while (!endfile)
922 
923  // tempOut.close();
924 
925 
927 
928  if ( outlist.compare("") != 0) fout.close();
929 
930  tfileStream.close();
931 
932  //Deallocate
933  delete[] pbuffer0[0];
934  delete[] pbuffer0;
935  delete[] pbuffer[0];
936  delete[] pbuffer;
937 
938  delete[] linerr;
939  delete[] seqerr;
940 
941  return 0;
942 }
943 
944 int make_oci_line_index( itab *itable, int16_t *cindex, int16_t *sindex,
945  // int32_t &ldark) {
946  // cdindex: CCD dark line index array
947  // sdindex: SWIR dark line index array
948  int16_t *cdindex, int16_t *sdindex, int16_t *swir_loff) {
949 
950  uint16_t maxlines = 32768;
951  uint16_t nagg[4] = {1,2,4,8};
952  unsigned short nswb = 9;
953 
954  // Loop through data zones
955  uint16_t loff = 0;
956  uint16_t cpix = 0;
957  uint16_t spix = 0;
958  uint16_t cdpix = 0;
959  uint16_t sdpix = 0;
960  uint16_t sloff = 0;
961 
962  for (size_t i=0; i<10; i++) {
963 
964  // If not "no data" type, add indices to array
965  if ((itable[i].dtype != 0) && (itable[i].dtype != 10)) {
966  if (itable[i].lines > 0) {
967 
968  // CCD pixel index
969  uint16_t iagg = nagg[itable[i].iagg];
970  uint16_t lines = itable[i].lines;
971  // Check for total lines within limit
972  if ((loff+lines) > maxlines) {
973  cout << "Mode table entry " << i << " exceeds max lines" << endl;
974  lines = maxlines - loff - iagg;
975  }
976 
977  uint16_t cp = lines / iagg;
978  for (size_t j=0; j<cp; j++) {
979  uint16_t cind = loff + j*iagg;
980  cindex[cind] = cpix + j;
981  }
982  cpix += cp;
983 
984  // SWIR pixel index (fixed aggregation of 8)
985  uint16_t sp = lines / 8;
986  for (size_t j=0; j<sp; j++) {
987  uint16_t sind = (loff/8 + j)*8;
988  // sindex[sind] = spix + j;
989  for (size_t k=0;k<nswb;k++) {
990  //If Earth view, use line offsets, LH 11/20/2020
991  if (itable[i].dtype == 1) {
992  sloff = swir_loff[k];
993  } else {
994  sloff = 0;
995  }
996  sindex[(sind-sloff)*nswb+k] = spix + j;
997  }
998  }
999 
1000  spix += sp;
1001 
1002  // Dark collect
1003  if (itable[i].dtype == 2) {
1004  for (size_t j=0; j<cp; j++) {
1005  uint16_t cind = loff + j*iagg;
1006  cdindex[cind] = cdpix + j;
1007  }
1008  cdpix += cp;
1009 
1010  for (size_t j=0; j<sp; j++) {
1011  uint16_t sind = loff + j*8;
1012  sdindex[sind] = sdpix + j;
1013  }
1014  sdpix += sp;
1015  }
1016  // Check for dark collect
1017  // if (itable[i].dtype == 2 && ldark == 32768) ldark = loff;
1018  } else {
1019  // if (itable[i].dtype != 0 && itable[i].lines == 0)
1020  cout << "Data zone " << i << " type " << itable[i].dtype <<
1021  " has zero lines" << endl;
1022  }
1023  }
1024  loff += itable[i].lines;
1025  if ( loff >= maxlines) break;
1026  //cout << "itable[i].lines: " << i << " " << itable[i].lines << endl;
1027  }
1028 
1029  return 0;
1030 }
1031 
1032 int unpack_oci_sci( uint32_t npkts, int32_t spin, uint16_t ncps, uint16_t nsps,
1033  uint16_t msps,
1034  uint16_t &nbands,
1035  uint16_t btaps[16], uint16_t rtaps[16],
1036  uint8_t (*pbuffer)[PKTSIZE],
1037  uint16_t **bbands, uint16_t **rbands, uint32_t **sbands,
1038  int16_t *blines, int16_t *rlines, int16_t *slines,
1039  uint16_t &btype, uint16_t bagg[16],
1040  uint16_t &rtype, uint16_t ragg[16],
1041  int8_t *sfrm, int &iret) {
1042 // Unpack all of the science data for an OCI scan.
1043 // The SWIR bands are re-ordered in ascending wavelength order.
1044 // Order of dual-gain bands is (SG, HG).
1045  for (size_t i=0;i<ncps;i++) {
1046  blines[i] = -1;
1047  rlines[i] = -1;
1048  }
1049  for (size_t i=0;i<msps;i++) slines[i] = -1;
1050 
1051  std::vector<int> iswav = {0,1,2,3,4,5,6,7,8};
1052  int ibpix = 0;
1053  int irpix = 0;
1054  int ispix = 0;
1055  iret = 0;
1056 
1057  uint16_t ccdid;
1058  uint32_t line;
1059  uint16_t iagg;
1060  uint16_t jagg[16];
1061 
1062  uint16_t **ccddata = new uint16_t*;
1063  uint16_t ossdata[16];
1064 
1065  uint16_t nbbnds;
1066  uint16_t nrbnds;
1067  uint16_t nb;
1068 
1069  uint32_t apid;
1070  uint32_t ui32;
1071  int32_t spn;
1072 
1073  //tempOut.open ("ccddata.bin", ios::out | ios::trunc | ios::binary);
1074 
1075  // SWIR band sorting indices for science mode
1076  uint16_t smode = 0;
1077  get_swir_mode( (uint8_t (*)[PKTSIZE]) &pbuffer[0][0], npkts, smode);
1078 
1079  if (smode == 0) {
1080  iswav = {3,0,1,2,8,6,7,5,4};
1081  }
1082 
1083  for (size_t ipkt=0; ipkt<npkts; ipkt++) {
1084  apid = (pbuffer[ipkt][0] % 8)*256 + pbuffer[ipkt][1];
1085  uint16_t dtype = (pbuffer[ipkt][12] % 64) / 4;
1086  memcpy( &ui32, &pbuffer[ipkt][6], 4);
1087  spn = SWAP_4( ui32);
1088 
1089  // If CCD science APID
1090  if (apid == 700 && dtype > 0 && spn == spin) {
1091  unpack_ccd_packet( pbuffer[ipkt], btaps, rtaps, ccdid, line, dtype,
1092  iagg, jagg, nbands, ccddata, ossdata);
1093 
1094  // tempOut.write((char *) &nbands, sizeof(uint16_t));
1095  //tempOut.write((char *) &(*ccddata)[0], nbands*sizeof(uint16_t));
1096 
1097  // If blue
1098  if (ccdid) {
1099  if (ibpix == 0) {
1100  nbbnds = nbands;
1101  btype = dtype;
1102  memcpy( bagg, jagg, 16*sizeof(uint16_t));
1103  for (size_t i=0; i<ncps; i++) {
1104  bbands[i] = new uint16_t[nbands];
1105  for (size_t j=0; j<nbands; j++) bbands[i][j] = 65535; // LH, 11/23/2020
1106  }
1107 
1108  } else {
1109  // Check for inconsistent non-dark data type or spectral aggregation
1110  uint8_t agg = 0;
1111  for (size_t i=0; i<16; i++)
1112  if (jagg[i] != bagg[i]) agg++;
1113 
1114  if ((dtype != btype && dtype != 2 && dtype !=5) || (agg > 0)) { // Ignore dark and linearity
1115  cout << "Data type or spectral aggregation error, CCDID: " <<
1116  ccdid << " Line: " << line << endl;
1117  iret = 4;
1118  }
1119  }
1120  if (nbands <= nbbnds) nb = nbands; else nb = nbbnds;
1121  if (ibpix < ncps) {
1122  memcpy( bbands[ibpix], &(*ccddata)[0], nb*sizeof(uint16_t));
1123  blines[ibpix] = (line/iagg) * iagg;
1124  } else {
1125  //memcpy( &ui32, &pbuffer[ipkt][6], 4);
1126  //int32_t spn = SWAP_4( ui32);
1127  cout << "Number of blue pixels exceeded in spin: " << spin
1128  << " in packet (0-based): " << ipkt << endl;
1129  }
1130  ibpix++;
1131 
1132  } else {
1133 
1134  if (irpix == 0) {
1135  nrbnds = nbands;
1136  rtype = dtype;
1137  memcpy( ragg, jagg, 16*sizeof(uint16_t));
1138  for (size_t i=0; i<ncps; i++) {
1139  rbands[i] = new uint16_t[nbands];
1140  for (size_t j=0; j<nbands; j++) rbands[i][j] = 65535; // LH, 11/23/2020
1141  }
1142  } else {
1143  // Check for inconsistent non-dark data type or spectral aggregation
1144  uint8_t agg = 0;
1145  for (size_t i=0; i<16; i++)
1146  if (jagg[i] != ragg[i]) agg++;
1147 
1148  if ((dtype != rtype && dtype != 2 && dtype !=5) || (agg > 0)) { // Ignore dark and linearity
1149  cout << "Data type or spectral aggregation error, CCDID: " <<
1150  ccdid << " Line: " << line << endl;
1151  iret = 4;
1152  }
1153  }
1154 
1155  if (nbands <= nrbnds) nb = nbands; else nb = nrbnds;
1156  if (irpix < ncps) {
1157  memcpy( rbands[irpix], &(*ccddata)[0], nb*sizeof(uint16_t));
1158  rlines[irpix] = (line/iagg) * iagg;
1159  } else {
1160  //memcpy( &ui32, &pbuffer[ipkt][6], 4);
1161  //int32_t spn = SWAP_4( ui32);
1162  cout << "Number of red pixels exceeded in spin: " << spin
1163  << " in packet (0-based): " << ipkt << endl;
1164  // cout << "irpix="<<irpix<<", ncps="<<ncps<<endl;
1165  // exit(1); // ver 1.03.02
1166  }
1167  irpix++;
1168  } // if (ccdid)
1169 
1170  // cout << "delete ccddata" << endl;
1171  delete[] *ccddata;
1172 
1173  } // if (apid == 700 && dtype > 0)
1174 
1175  // If SWIR science APID
1176  if (apid == 720 && (ispix+7) < msps && spn == spin) {
1177  int16_t lines[8];
1178  uint8_t swirfrm[8];
1179  uint32_t swirdata[8*9];
1180 
1181  // swirdata: 8 rows of 9 columns
1182  unpack_swir_packet( pbuffer[ipkt], lines, swirfrm, swirdata);
1183 
1184  for (size_t i=0; i<8; i++) {
1185  slines[ispix+i] = (lines[i]/8) * 8;
1186  for (size_t j=0; j<9; j++) sbands[ispix+i][j] = swirdata[iswav[j]+9*i]; // bands in ascending order
1187  }
1188  memcpy(&sfrm[ispix], swirfrm, 8*sizeof(uint8_t));
1189 
1190  ispix += 8;
1191  } // if (apid == 720 && (ispix+7) < msps)
1192 
1193  if ((ibpix <= 0) || (ibpix > ncps)) blines[0] = -1;
1194  if ((irpix <= 0) || (irpix > ncps)) rlines[0] = -1;
1195  if ((ispix <= 0) || (ispix > msps)) slines[0] = -1;
1196 
1197  } // ipkt loop
1198 
1199  delete ccddata;
1200  // tempOut.close();
1201 
1202  return 0;
1203 }
1204 
1205 
1206 int unpack_ccd_packet( uint8_t *packet, uint16_t btaps[16], uint16_t rtaps[16],
1207  uint16_t &ccdid, uint32_t &line, uint16_t &dtype,
1208  uint16_t &iagg, uint16_t jagg[16], uint16_t &nbands,
1209  uint16_t **ccddata, uint16_t ossdata[16]) {
1210 
1211  // Get CCD ID, line number, data type and spatial aggregation
1212  ccdid = (packet[12] & 64) / 64;
1213  line = packet[10]*256 + packet[11];
1214  dtype = (packet[12] % 64) / 4;
1215  uint16_t oss = packet[17] % 16;
1216  iagg = packet[12] % 4;
1217  uint16_t agg[4] = {1,2,4,8};
1218  iagg = agg[iagg];
1219 
1220  // Get tap enable flags for focal plane
1221  uint16_t *ftaps;
1222  if (ccdid) ftaps = btaps; else ftaps = rtaps;
1223 
1224  // Get spectral aggregation factors for all taps
1225  uint16_t its[4] = {0,4,8,12};
1226  uint32_t taps[16];
1227  for (size_t i=0; i<4; i++) {
1228  jagg[its[i]] = agg[(packet[13+i] & 192) / 64];
1229  jagg[its[i]+1] = agg[(packet[13+i] & 48) / 16];
1230  jagg[its[i]+2] = agg[(packet[13+i] & 12) / 4];
1231  jagg[its[i]+3] = agg[(packet[13+i] & 3)];
1232  }
1233  for (size_t i=0; i<16; i++) taps[i] = 32*ftaps[i] / jagg[i];
1234 
1235  // Allocate output buffer
1236  nbands = 0;
1237  for (size_t i=0; i<16; i++) nbands += taps[i];
1238 
1239  //cout << "new ccddata" << endl;
1240  *ccddata = new uint16_t[nbands];
1241 
1242  int ioff = 18;
1243  uint16_t ibnd = nbands - 1;
1244 
1245  // For each tap (1-16):
1246  for (size_t j=0; j<16; j++) {
1247 
1248  // Copy band data from packet to output buffer
1249  // Packet data are in reverse spectral order, so swap order here
1250  if (ftaps[j]) {
1251  for (size_t i=0; i<taps[j]; i++) {
1252  uint16_t ui16;
1253  memcpy( &ui16, &packet[ioff+2*i], 2);
1254  ui16 = SWAP_2(ui16);
1255  memcpy( &(*ccddata)[ibnd-i], &ui16, 2);
1256  }
1257  ibnd -= taps[j];
1258  ioff += 2*taps[j];
1259  }
1260  }
1261 
1262  if (oss > 0) {
1263  for (size_t j=0; j<16; j++) {
1264  uint16_t ui16;
1265  memcpy( &ui16, &packet[ioff+2*j], 2);
1266  ui16 = SWAP_2(ui16);
1267  memcpy(&ossdata[j], &ui16, 2);
1268  }
1269  }
1270 
1271  return 0;
1272 }
1273 
1274 
1275 int unpack_swir_packet( uint8_t *packet, int16_t *slines, uint8_t *swirfrm,
1276  uint32_t *swirdata) {
1277 
1278  // Program to unpack one OCI SWIR data packet
1279  // Each packet contains data for eight science pixels
1280  // Reference: OCI-ELEC-SPEC-0028
1281 
1282  // Allocate output buffers (9 bands for 8 pixels)
1283 
1284  int ioff = 10;
1285  uint8_t smeta;;
1286 
1287  for (size_t i=0; i<8; i++) {
1288 
1289  // Get line number and metadata
1290  if ((unsigned(packet[ioff])==255) && (unsigned(packet[ioff+1])==255)) {
1291  // handle fill values in SWIR, LH 10/28/2020
1292  slines[i]=0;
1293  } else {
1294  slines[i] = packet[ioff]*256 + packet[ioff+1];
1295  }
1296 
1297  smeta = packet[ioff+2];
1298  swirfrm[i] = smeta % 8;
1299 
1300  // Extract SWIR band data
1301  eight20( &packet[ioff+3], &swirdata[9*i]);
1302 
1303  ioff += 26;
1304  }
1305 
1306  return 0;
1307 }
1308 
1309 
1310 int check_load_oci_data( short dtype, uint16_t ncpt, uint16_t nspt0,
1311  uint16_t ndcs, uint16_t ndss,
1312  uint16_t nbbs, uint16_t nrbs, uint16_t nswb,
1313  int16_t *cindex, int16_t *sindex,
1314  int16_t *cdindex, int16_t *sdindex,
1315  uint16_t **bbands, uint16_t **rbands,
1316  uint32_t **sbands,
1317  int16_t *blines, int16_t *rlines, int16_t *slines,
1318  uint16_t **bsci, uint16_t **rsci, uint32_t **ssci,
1319  uint16_t **bdark, uint16_t **rdark, uint32_t **sdark,
1320  uint8_t &linerr, int &icheck) {
1321 
1322  icheck = 0;
1323 
1324  // CCD line check except in linearity mode
1325  bool linchk = (dtype != 5);
1326 
1327  // Check for presence of blue and red focal plane data
1328  //dbb = size(bbands)
1329  bool bb = ((bbands[0] != NULL) && (blines[0] != -1));
1330  // drb = size(rbands)
1331  bool rb = ((rbands[0] != NULL) && (rlines[0] != -1));
1332  //drs = size(sbands)
1333  bool sb = ((sbands[0] != NULL) && (slines[0] != -1));
1334 
1336  uint16_t nbb;
1337  if (bbands[0] == NULL) nbb = 0; else nbb = nbbs;
1338  int mb=-1;
1339 
1340  // Check for invalid line numbers
1341  if ( bb && linchk) {
1342  for (size_t i=0; i<ncpt; i++) {
1343  if (blines[i] == -1) continue;
1344  if (cindex[blines[i]] == -1) mb++;
1345  }
1346  for (size_t i=0; i<ncpt; i++)
1347  if (cindex[blines[i]] != (cindex[blines[0]] + (int) i)) {
1348  cout << "Blue CCD line sequence error" << endl;
1349  linerr = 1;
1350  break;
1351  }
1352  }
1353 
1354  // Identify pixels and load into output arrays
1355  if (bb) {
1356  uint16_t nbs=0;
1357  for (size_t i=0; i<ncpt; i++) {
1358  if (blines[i] == -1) continue;
1359  // if (blines[i] < ldark && cindex[blines[i]] != -1) nbs++;
1360  if (cindex[blines[i]] > -1) nbs++; // LH, 8/24/2020, replaced != with >
1361  }
1362  if ((nbs < ncpt-ndcs) && linchk) {
1363  cout << "Missing blue band science pixels" << endl;
1364  icheck = 2;
1365  }
1366 
1367  for (size_t i=0; i<ncpt; i++) {
1368  if (blines[i] == -1) continue;
1369  // if (blines[i] < ldark && cindex[blines[i]] != -1) {
1370  if (cindex[blines[i]] > -1) { // LH, 8/24/2020
1371  for (size_t j=0; j<nbb; j++) {
1372  bsci[j][cindex[blines[i]]] = bbands[i][j];
1373  }
1374  }
1375  }
1376 
1377  // Identify dark count pixels and load into output arrays
1378  uint16_t nbd=0;
1379  for (size_t i=0; i<ncpt; i++) {
1380  if (blines[i] == -1) continue;
1381  // if (blines[i] >= ldark && cindex[blines[i]] != -1) nbd++;
1382  if (cdindex[blines[i]] > -1) nbd++; // LH, 8/24/2020
1383  }
1384  if (nbd < ndcs) {
1385  cout << "Missing blue band dark pixels" << endl;
1386  icheck = 2;
1387  }
1388 
1389  for (size_t i=0; i<ncpt; i++) {
1390  if (blines[i] == -1) continue;
1391  //if (blines[i] >= ldark && cindex[blines[i]] != -1) {
1392  if (cdindex[blines[i]] > -1) { // LH, 8/24/2020
1393  // int k = cindex[blines[i]]-cindex[ldark];
1394  int k = cdindex[blines[i]];
1395  if ( k >= ndcs || k < 0) {
1396  //cout << "bdark indice out of bounds: " << k << endl;
1397  continue;
1398  // exit(1);
1399  }
1400  for (size_t j=0; j<nbb; j++) {
1401  bdark[j][k] = bbands[i][j];
1402  }
1403  }
1404  }
1405  }
1406 
1408  uint16_t nrb;
1409  if (rbands[0] == NULL) nrb = 0; else nrb = nrbs;
1410  int mr=-1;
1411 
1412  // Check for invalid line numbers
1413  if ( rb && linchk) {
1414  for (size_t i=0; i<ncpt; i++) {
1415  if (rlines[i] == -1) continue;
1416  if (cindex[rlines[i]] == -1) mr++;
1417  }
1418  for (size_t i=0; i<ncpt; i++) {
1419  if (rlines[i] == -1) continue;
1420  if (cindex[rlines[i]] != (cindex[rlines[0]] + (int) i)) {
1421  cout << "Red CCD line sequence error" << endl;
1422  linerr = 1;
1423  break;
1424  }
1425  }
1426  }
1427 
1428  // Identify pixels and load into output arrays
1429  if (rb) {
1430  uint16_t nrs=0;
1431  for (size_t i=0; i<ncpt; i++) {
1432  if (rlines[i] == -1) continue;
1433  // if (rlines[i] < ldark && cindex[rlines[i]] != -1) nrs++;
1434  if (cindex[rlines[i]] > -1) nrs++; // LH, 8/24/2020
1435  }
1436 
1437  if ((nrs < ncpt-ndcs) && linchk) {
1438  cout << "Missing red band science pixels" << endl;
1439  icheck = 2;
1440  }
1441 
1442  for (size_t i=0; i<ncpt; i++) {
1443  // if (rlines[i] < ldark && cindex[rlines[i]] != -1) {
1444  if (rlines[i] == -1) continue;
1445  if (rlines[i] > -1 && cindex[rlines[i]] > -1) { // LH, 8/24/2020
1446  for (size_t j=0; j<nrb; j++) {
1447  rsci[j][cindex[rlines[i]]] = rbands[i][j];
1448  }
1449  }
1450  }
1451 
1452  // Identify dark count pixels and load into output arrays
1453  uint16_t nrd=0;
1454  for (size_t i=0; i<ncpt; i++) {
1455  if (rlines[i] == -1) continue;
1456  // if (rlines[i] >= ldark && cindex[rlines[i]] != -1) nrd++;
1457  if (cdindex[rlines[i]] > -1) nrd++; // LH, 8/24/2020
1458  }
1459 
1460  if (nrd < ndcs) {
1461  cout << "Missing red band dark pixels" << endl;
1462  icheck = 2;
1463  }
1464 
1465  for (size_t i=0; i<ncpt; i++) {
1466  if (rlines[i] == -1) continue;
1467  // if (rlines[i] >= ldark && cindex[rlines[i]] != -1) {
1468  if (cdindex[rlines[i]] > -1) { // LH, 8/24/2020
1469  // int k = cindex[rlines[i]]-cindex[ldark];
1470  int k = cdindex[rlines[i]];
1471  if ( k >= ndcs || k < 0) {
1472  //cout << "rdark indice out of bounds: " << k << endl;
1473  continue;
1474  // exit(1);
1475  }
1476  for (size_t j=0; j<nrb; j++) {
1477  rdark[j][k] = rbands[i][j];
1478  }
1479  }
1480  }
1481  }
1482 
1483  if (mb != -1 || mr !=-1) {
1484  // Disable SWIR line check for ETU
1485  cout << "Invalid line numbers" << endl;
1486  icheck = 4;
1487  }
1488 
1490 
1491  // Identify pixels and load into output arrays
1492  // SWIR science data line index has band-by-band offsets, LH, 11/20/2020
1493  if (sb) {
1494  uint16_t nss=0;
1495  for (size_t i=0; i<nspt0; i++) {
1496  if (slines[i] > -1) {
1497  nss++;
1498  for (size_t j=0; j<nswb; j++) {
1499  if (sindex[slines[i]*nswb+j] > -1) { // LH, 8/24/2020; sindex expanded by * nswb bands, LH 11/20/2020
1500  ssci[j][sindex[slines[i]*nswb+j]] = sbands[i][j];
1501  }
1502  }
1503  }
1504  }
1505 
1506  if (nss < nspt0-ndss) {
1507  cout << "Missing SWIR band science pixels" << endl;
1508  icheck = 2;
1509  }
1510 
1511  // Identify dark count pixels and load into output arrays
1512  uint16_t nsd=0;
1513  for (size_t i=0; i<nspt0; i++) {
1514  if (slines[i] == -1) continue;
1515  // if (slines[i] >= ldark && sindex[slines[i]] != -1) nsd++;
1516  if (sdindex[slines[i]] > -1) nsd++; // LH, 8/24/2020
1517  }
1518  //if (nsd < ndss) {
1519  // cout << "Missing SWIR band dark pixels" << endl;
1520  // icheck = 1;
1521  //}
1522 
1523  for (size_t i=0; i<nspt0; i++) {
1524  // if (slines[i] >= ldark && sindex[slines[i]] != -1) {
1525  if (slines[i] == -1) continue;
1526  if (sdindex[slines[i]] > -1) { // LH, 8/24/2020
1527  // int k = sindex[slines[i]]-sindex[ldark];
1528  int k = sdindex[slines[i]];
1529  if ( k >= ndss || k < 0) {
1530  cout << "sdark indice out of bounds: " << k << endl;
1531  exit(1);
1532  }
1533 
1534  for (size_t j=0; j<nswb; j++) {
1535  sdark[j][k] = sbands[i][j];
1536  }
1537  }
1538  }
1539  }
1540  return 0;
1541 }
1542 
1543 uint8_t check_sum( int32_t nc, uint8_t *dat, uint8_t *chk) {
1544 
1545  // Function to check data against checksum
1546  // Checksum is 4 bytes computed by XOR
1547 
1548  uint8_t chks[4], tmp[4];
1549  memcpy( chks, &dat[0], 4);
1550 
1551  for ( int i=1; i<nc; i++) {
1552  memcpy( &tmp, &dat[4*i], 4);
1553  for ( int j=0; j<4; j++) chks[j] = chks[j] ^ tmp[j];
1554  }
1555 
1556  uint8_t check[4];
1557  for ( int i=0; i<4; i++) check[i] = (chk[i] == chks[i]);
1558 
1559  uint8_t check_sum = (check[0] & check[1] & check[2] & check[3]);
1560 
1561  return check_sum;
1562 }
1563 
1564 
1565 /*----------------------------------------------------------------- */
1566 /* Create an Generic NETCDF4 level1 file */
1567 /* ---------------------------------------------------------------- */
1568 int l1aFile::createl1( char* l1_filename, uint16_t maxsc,
1569  uint16_t ncps, uint16_t nbbs, uint16_t nrbs,
1570  uint16_t nsps, uint16_t ndcs) {
1571 
1572  try {
1573  l1afile = new NcFile( l1_filename, NcFile::replace);
1574  }
1575  catch ( NcException& e) {
1576  e.what();
1577  cerr << "Failure creating OCI L1A file: " << l1_filename << endl;
1578  exit(1);
1579  }
1580 
1581  fileName.assign( l1_filename);
1582 
1583  ifstream oci_l1a_data_structure;
1584  string line;
1585  string dataStructureFile;
1586  dataStructureFile.assign("$OCDATAROOT/oci/OCI_Level-1A_Data_Structure.cdl");
1587  expandEnvVar( &dataStructureFile);
1588 
1589  oci_l1a_data_structure.open( dataStructureFile.c_str(), ifstream::in);
1590  if ( oci_l1a_data_structure.fail() == true) {
1591  cout << "\"" << dataStructureFile.c_str() << "\" not found" << endl;
1592  exit(1);
1593  }
1594 
1595  // Find "dimensions" section of CDL file
1596  while(1) {
1597  getline( oci_l1a_data_structure, line);
1598  size_t pos = line.find("dimensions:");
1599  if ( pos == 0) break;
1600  }
1601 
1602  // Define dimensions from "dimensions" section of CDL file
1603  ndims = 0;
1604 
1605  while(1) {
1606  getline( oci_l1a_data_structure, line);
1607  if (line.substr(0,2) == "//") continue;
1608 
1609  size_t pos = line.find(" = ");
1610  size_t semi = line.find(" ;");
1611  if ( pos == string::npos) break;
1612 
1613  uint32_t dimSize;
1614  istringstream iss;
1615  // istringstream iss(line.substr(pos+2, string::npos));
1616  string dimString = line.substr(pos+2, semi-(pos+2));
1617  if (dimString.find("UNLIMITED") == string::npos) {
1618  iss.str(dimString);
1619  iss >> dimSize;
1620  } else {
1621  dimSize = NC_UNLIMITED;
1622  }
1623 
1624  iss.clear();
1625  iss.str( line);
1626  iss >> skipws >> line;
1627 
1628  // cout << "Dimension Name: " << line.c_str() << " Dimension Size: "
1629  // << dimSize << endl;
1630 
1631  if (line.compare("ccd_pixels") == 0) {
1632  dimSize = ncps;
1633  }
1634 
1635  if (line.compare("SWIR_pixels") == 0) {
1636  dimSize = nsps;
1637  }
1638 
1639  if (line.compare("DC_pixels") == 0) {
1640  dimSize = ndcs;
1641  }
1642 
1643  if (line.compare("blue_bands") == 0) {
1644  dimSize = nbbs;
1645  }
1646 
1647  if (line.compare("red_bands") == 0) {
1648  dimSize = nrbs;
1649  }
1650 
1651  try {
1652  ncDims[ndims++] = l1afile->addDim( line, dimSize);
1653  }
1654  catch ( NcException& e) {
1655  e.what();
1656  cerr << "Failure creating dimension: " << line.c_str() << endl;
1657  exit(1);
1658  }
1659 
1660  } // while loop
1661 
1662  // Read global attributes (string attributes only)
1663  while(1) {
1664  getline( oci_l1a_data_structure, line);
1665  size_t pos = line.find("// global attributes");
1666  if ( pos == 0) break;
1667  }
1668 
1669  while(1) {
1670  getline( oci_l1a_data_structure, line);
1671  size_t pos = line.find(" = ");
1672  if ( pos == string::npos) break;
1673 
1674  string attValue;
1675 
1676  // Remove leading and trailing quotes
1677  attValue.assign(line.substr(pos+4));
1678  size_t posQuote = attValue.find("\"");
1679  attValue.assign(attValue.substr(0, posQuote));
1680 
1681  istringstream iss(line.substr(pos+2));
1682  iss.clear();
1683  iss.str( line);
1684  iss >> skipws >> line;
1685 
1686  // Skip commented out attributes
1687  if (line.substr(0,2) == "//") continue;
1688 
1689  string attName;
1690  attName.assign(line.substr(1).c_str());
1691 
1692  // Skip non-string attributes
1693  if (attName.compare("orbit_number") == 0) continue;
1694  if (attName.compare("history") == 0) continue;
1695  if (attName.compare("format_version") == 0) continue;
1696  if (attName.compare("instrument_number") == 0) continue;
1697  if (attName.compare("pixel_offset") == 0) continue;
1698  if (attName.compare("number_of_filled_scans") == 0) continue;
1699 
1700  // cout << attName.c_str() << " " << attValue.c_str() << endl;
1701  try {
1702  l1afile->putAtt(attName, attValue);
1703  }
1704  catch ( NcException& e) {
1705  e.what();
1706  cerr << "Failure creating attribute: " + attName << endl;
1707  exit(1);
1708  }
1709 
1710  } // while(1)
1711 
1712 
1713  ngrps = 0;
1714  // Loop through groups
1715  while(1) {
1716  getline( oci_l1a_data_structure, line);
1717 
1718  // Check if end of CDL file
1719  // If so then close CDL file and return
1720  if (line.substr(0,1).compare("}") == 0) {
1721  oci_l1a_data_structure.close();
1722  return 0;
1723  }
1724 
1725  // Check for beginning of new group
1726  size_t pos = line.find("group:");
1727 
1728  // If found then create new group and variables
1729  if ( pos == 0) {
1730 
1731  // Parse group name
1732  istringstream iss(line.substr(6, string::npos));
1733  iss >> skipws >> line;
1734 
1735  ncGrps[ngrps++] = l1afile->addGroup(line);
1736 
1737  int numDims=0;
1738  string sname;
1739  string lname;
1740  string standard_name;
1741  string units;
1742  string flag_values;
1743  string flag_meanings;
1744  string reference;
1745  double valid_min=0.0;
1746  double valid_max=0.0;
1747  double fill_value=0.0;
1748 
1749  vector<NcDim> varVec;
1750 
1751  int ntype=0;
1752  NcType ncType;
1753 
1754  // Loop through datasets in group
1755  getline( oci_l1a_data_structure, line); // skip "variables:"
1756  while(1) {
1757  getline( oci_l1a_data_structure, line);
1758 
1759  if (line.substr(0,2) == "//") continue;
1760  if (line.length() == 0) continue;
1761  if (line.substr(0,1).compare("\r") == 0) continue;
1762  if (line.substr(0,1).compare("\n") == 0) continue;
1763 
1764  size_t pos = line.find(":");
1765 
1766  // No ":" found, new dataset or empty line or end-of-group
1767  if ( pos == string::npos) {
1768 
1769  if ( numDims > 0) {
1770  // Create previous dataset
1771 
1772  createNCDF( ncGrps[ngrps-1], sname.c_str(), lname.c_str(),
1773  standard_name.c_str(), units.c_str(),
1774  (void *) &fill_value,
1775  flag_values.c_str(), flag_meanings.c_str(),
1776  reference.c_str(),
1777  valid_min, valid_max, ntype, varVec);
1778 
1779  flag_values.assign("");
1780  flag_meanings.assign("");
1781  reference.assign("");
1782  units.assign("");
1783  varVec.clear();
1784  }
1785 
1786  valid_min=0.0;
1787  valid_max=0.0;
1788  fill_value=0.0;
1789 
1790  if (line.substr(0,10).compare("} // group") == 0) break;
1791 
1792  // Parse variable type
1793  string varType;
1794  istringstream iss(line);
1795  iss >> skipws >> varType;
1796 
1797  // Get corresponding NC variable type
1798  if ( varType.compare("char") == 0) ntype = NC_CHAR;
1799  else if ( varType.compare("byte") == 0) ntype = NC_BYTE;
1800  else if ( varType.compare("short") == 0) ntype = NC_SHORT;
1801  else if ( varType.compare("int") == 0) ntype = NC_INT;
1802  else if ( varType.compare("long") == 0) ntype = NC_INT;
1803  else if ( varType.compare("float") == 0) ntype = NC_FLOAT;
1804  else if ( varType.compare("real") == 0) ntype = NC_FLOAT;
1805  else if ( varType.compare("double") == 0) ntype = NC_DOUBLE;
1806  else if ( varType.compare("ubyte") == 0) ntype = NC_UBYTE;
1807  else if ( varType.compare("ushort") == 0) ntype = NC_USHORT;
1808  else if ( varType.compare("uint") == 0) ntype = NC_UINT;
1809  else if ( varType.compare("ulong") == 0) ntype = NC_UINT;
1810  else if ( varType.compare("int64") == 0) ntype = NC_INT64;
1811  else if ( varType.compare("uint64") == 0) ntype = NC_UINT64;
1812 
1813  // Parse short name (sname)
1814  pos = line.find("(");
1815  size_t posSname = line.substr(0, pos).rfind(" ");
1816  sname.assign(line.substr(posSname+1, pos-posSname-1));
1817  // cout << "sname: " << sname.c_str() << endl;
1818 
1819  // Parse variable dimension info
1820  this->parseDims( line.substr(pos+1, string::npos), varVec);
1821  numDims = varVec.size();
1822 
1823  } else {
1824  // Parse variable attributes
1825  size_t posEql = line.find("=");
1826  size_t pos1qte = line.find("\"");
1827  size_t pos2qte = line.substr(pos1qte+1, string::npos).find("\"");
1828  // cout << line.substr(pos+1, posEql-pos-2).c_str() << endl;
1829 
1830  string attrName = line.substr(pos+1, posEql-pos-2);
1831 
1832  // Get long_name
1833  if ( attrName.compare("long_name") == 0) {
1834  lname.assign(line.substr(pos1qte+1, pos2qte));
1835  // cout << "lname: " << lname.c_str() << endl;
1836  }
1837 
1838  // Get units
1839  else if ( attrName.compare("units") == 0) {
1840  units.assign(line.substr(pos1qte+1, pos2qte));
1841  // cout << "units: " << units.c_str() << endl;
1842  }
1843 
1844  // Get _FillValue
1845  else if ( attrName.compare("_FillValue") == 0) {
1846  iss.clear();
1847  iss.str( line.substr(posEql+1, string::npos));
1848  iss >> fill_value;
1849  }
1850 
1851  // Get flag_values
1852  else if ( attrName.compare("flag_values") == 0) {
1853  flag_values.assign(line.substr(pos1qte+1, pos2qte));
1854  }
1855 
1856  // Get flag_meanings
1857  else if ( attrName.compare("flag_meanings") == 0) {
1858  flag_meanings.assign(line.substr(pos1qte+1, pos2qte));
1859  }
1860 
1861  // Get reference
1862  else if ( attrName.compare("reference") == 0) {
1863  reference.assign(line.substr(pos1qte+1, pos2qte));
1864  }
1865 
1866  // Get valid_min
1867  else if ( attrName.compare("valid_min") == 0) {
1868  iss.clear();
1869  iss.str( line.substr(posEql+1, string::npos));
1870  iss >> valid_min;
1871  // cout << "valid_min: " << valid_min << endl;
1872  }
1873 
1874  // Get valid_max
1875  else if ( attrName.compare("valid_max") == 0) {
1876  iss.clear();
1877  iss.str( line.substr(posEql+1, string::npos));
1878  iss >> valid_max;
1879  // cout << "valid_max: " << valid_max << endl;
1880  }
1881 
1882  } // if ( pos == string::npos)
1883  } // datasets in group loop
1884  } // New Group loop
1885  } // Main Group loop
1886 
1887  return 0;
1888 }
1889 
1890 
1891 int l1aFile::parseDims( string dimString, vector<NcDim>& varDims) {
1892 
1893  size_t curPos=0;
1894  // char dimName[NC_MAX_NAME+1];
1895  string dimName;
1896 
1897  while(1) {
1898  size_t pos = dimString.find(",", curPos);
1899  if ( pos == string::npos)
1900  pos = dimString.find(")");
1901 
1902  string varDimName;
1903  istringstream iss(dimString.substr(curPos, pos-curPos));
1904  iss >> skipws >> varDimName;
1905 
1906  for (int i=0; i<ndims; i++) {
1907 
1908  try {
1909  dimName = ncDims[i].getName();
1910  }
1911  catch ( NcException& e) {
1912  e.what();
1913  cerr << "Failure accessing dimension: " + dimName << endl;
1914  exit(1);
1915  }
1916 
1917  if ( varDimName.compare(dimName) == 0) {
1918  // cout << " " << dimName << " " << ncDims[i].getSize() << endl;
1919  varDims.push_back(ncDims[i]);
1920  break;
1921  }
1922  }
1923  if ( dimString.substr(pos, 1).compare(")") == 0) break;
1924 
1925  curPos = pos + 1;
1926  }
1927 
1928  return 0;
1929 }
1930 
1931 
1932 int l1aFile::write_oci_science_data( uint32_t isc,
1933  uint16_t nbbs, uint16_t nrbs,
1934  uint16_t nswb,
1935  uint16_t ncps, uint16_t nsps,
1936  uint16_t **bsci, uint16_t **rsci,
1937  uint32_t **ssci, int8_t *sfrms) {
1938 
1939  // Writes one scan at a time
1940 
1941  NcVar blu_bands;
1942  NcVar red_bands;
1943  NcVar swir_bands;
1944  NcVar swir_frms;
1945  vector<size_t> start;
1946  vector<size_t> count_b;
1947  vector<size_t> count_r;
1948  vector<size_t> count_s;
1949  vector<size_t> count_0; // for 1-D fill value variables
1950 
1951  vector<size_t> start_frms;
1952  vector<size_t> count_frms;
1953 
1954  NcGroup gid = l1afile->getGroup( "science_data");
1955  blu_bands = gid.getVar("sci_blue");
1956  red_bands = gid.getVar("sci_red");
1957  swir_bands = gid.getVar("sci_SWIR");
1958  swir_frms = gid.getVar("frm_type_SWIR");
1959 
1960  start.push_back(0);
1961  start.push_back(0);
1962  start.push_back(0);
1963 
1964  count_b.push_back(1);
1965  count_b.push_back(nbbs);
1966  count_b.push_back(ncps);
1967 
1968  count_r.push_back(1);
1969  count_r.push_back(nrbs);
1970  count_r.push_back(ncps);
1971 
1972  count_s.push_back(1);
1973  count_s.push_back(nswb);
1974  count_s.push_back(nsps);
1975 
1976  count_0.push_back(1);
1977  count_0.push_back(1);
1978  count_0.push_back(1);
1979 
1980  //NcDim dim = blu_bands.getDim(0);
1981  //cout << dim.getName() << endl;
1982  //cout << dim.getSize() << endl;
1983  //cout << dim.isUnlimited() << endl;
1984 
1985  start[0] = isc;
1986 
1987  // create 1-D sci_blue/sci_red/sci_SWIR with fill value if number of blue bands is 0
1988  (count_b[1] > 0) ? blu_bands.putVar(start, count_b, &bsci[0][0]) : blu_bands.putVar(start, count_0, &bsci[0][0]);
1989  (count_r[1] > 0) ? red_bands.putVar(start, count_r, &rsci[0][0]) : red_bands.putVar(start, count_0, &rsci[0][0]);
1990  (count_s[1] > 0) ? swir_bands.putVar(start, count_s, &ssci[0][0]) : swir_bands.putVar(start, count_0, &ssci[0][0]);
1991 
1992  start_frms.push_back(isc);
1993  start_frms.push_back(0);
1994 
1995  count_frms.push_back(1);
1996  count_frms.push_back(nsps);
1997 
1998  if ( count_frms[1] > 0)
1999  swir_frms.putVar(start_frms, count_frms, sfrms);
2000 
2001  return 0;
2002 }
2003 
2004 
2005 int l1aFile::write_oci_cal_data( uint32_t isc,
2006  uint16_t nbbs, uint16_t nrbs, uint16_t nswb,
2007  uint16_t ndcs, uint16_t ndss,
2008  uint16_t *dark_b, uint16_t *dark_r,
2009  uint32_t *dark_s,
2010  int8_t *sdfrms) {
2011 
2012  NcVar blu_dark;
2013  NcVar red_dark;
2014  NcVar swr_dark;
2015  vector<size_t> start;
2016  vector<size_t> count_b;
2017  vector<size_t> count_r;
2018  vector<size_t> count_s;
2019  vector<size_t> count_0; // for 1-D fill value variables
2020 
2021  start.push_back(0);
2022  start.push_back(0);
2023  start.push_back(0);
2024 
2025  count_b.push_back(isc);
2026  count_b.push_back(nbbs);
2027  count_b.push_back(ndcs);
2028 
2029  count_r.push_back(isc);
2030  count_r.push_back(nrbs);
2031  count_r.push_back(ndcs);
2032 
2033  count_s.push_back(isc);
2034  count_s.push_back(nswb);
2035  count_s.push_back(ndss);
2036 
2037  count_0.push_back(1);
2038  count_0.push_back(1);
2039  count_0.push_back(1);
2040 
2041  NcGroup gid = l1afile->getGroup( "onboard_calibration_data");
2042  blu_dark = gid.getVar("DC_blue");
2043  red_dark = gid.getVar("DC_red");
2044  swr_dark = gid.getVar("DC_SWIR");
2045 
2046  // create 1-D DC_blue/DC_red/DC_SWIR with fill value if number of blue bands is 0
2047  (count_b[1] > 0) ? blu_dark.putVar(start, count_b, dark_b) : blu_dark.putVar(start, count_0, dark_b);
2048  (count_r[1] > 0) ? red_dark.putVar(start, count_r, dark_r) : red_dark.putVar(start, count_0, dark_r);
2049  (count_s[1] > 0) ? swr_dark.putVar(start, count_s, dark_s) : swr_dark.putVar(start, count_0, dark_s);
2050 
2051  NcVar sdfrm_dark;
2052  vector<size_t> start_sdfrm;
2053  vector<size_t> count_sdfrm;
2054 
2055  start_sdfrm.push_back(0);
2056  start_sdfrm.push_back(0);
2057  count_sdfrm.push_back(isc);
2058  count_sdfrm.push_back(ndss);
2059 
2060  sdfrm_dark = gid.getVar("frm_type_DC_SWIR");
2061  // if ( count_sdfrm[1] > 0 && sdfrm_dark.isNull())
2062  if ( count_sdfrm[1] > 0)
2063  sdfrm_dark.putVar(start_sdfrm, count_sdfrm, sdfrms);
2064 
2065  return 0;
2066 }
2067 
2068 
2069 int l1aFile::write_oci_scan_metadata( uint32_t isc, uint8_t *ancdata,
2070  uint8_t *seqerr, uint8_t *linerr,
2071  int32_t *spinID) {
2072 
2073  vector<size_t> start;
2074  vector<size_t> count;
2075  start.push_back(0);
2076  count.push_back(isc);
2077 
2078  uint32_t ancap = 636;
2079 
2080  int32_t iyear, iday;
2081  double stime;
2082 
2083  // Extract and convert times to seconds of day
2084  // Extract CCSDS scan start times
2085  double *stimes = new double[isc];
2086  short int toff = 24;
2087  uint32_t *scss = new uint32_t[isc];
2088  int32_t *scsu = new int32_t[isc];
2089 
2090  for (size_t i=0; i<isc; i++) {
2091  uint32_t apid = (ancdata[i*ANCSIZE] % 8)*256 + ancdata[i*ANCSIZE+1];
2092  if (apid == ancap) {
2093  get_anc_packet_time( &ancdata[i*ANCSIZE], iyear, iday, stime);
2094  stimes[i] = stime;
2095 
2096  uint32_t ui32;
2097  memcpy( &ui32, &ancdata[i*ANCSIZE+toff+4], 4);
2098  scss[i] = SWAP_4( ui32);
2099  memcpy( &ui32, &ancdata[i*ANCSIZE+toff+8], 4);
2100  scsu[i] = SWAP_4( ui32) / 4096;
2101  } else {
2102  stimes[i] = -999;
2103  scss[i] = 0;
2104  scsu[i] = -999;
2105  }
2106  }
2107 
2108  NcGroup gid = l1afile->getGroup( "scan_line_attributes");
2109 
2110  NcVar var;
2111 
2112  // Scan start time (seconds of day)
2113  var = gid.getVar("scan_start_time");
2114  var.putVar(start, count, stimes);
2115 
2116  // Scan start time (CCSDS)
2117 
2118  // Seconds since 1958
2119  var = gid.getVar("scan_start_CCSDS_sec");
2120  var.putVar(start, count, scss);
2121 
2122  // Microseconds
2123  var = gid.getVar("scan_start_CCSDS_usec");
2124  var.putVar(start, count, scsu);
2125 
2126 
2127  // Extract and write HAM side
2128  uint8_t *hamSide = new uint8_t[isc];
2129 
2130  for (size_t i=0; i<isc; i++) {
2131  // This part is modified to fix bogus values for the last two records in last chunk of L1A
2132  // Liang Hong, 4/24/2020
2133  //int ham =
2134  // (ancdata[(i+1)*ANCSIZE+toff+86] % 128) * 256 +
2135  // ancdata[(i+1)*ANCSIZE+toff+87];
2136  int ham = 255;
2137  //ham /= 16384; // Need to revisit this when better known
2138  hamSide[i] = (uint8_t) ham;
2139  }
2140  var = gid.getVar("HAM_side");
2141  var.putVar(start, count, hamSide);
2142 
2143 
2144  // Extract and write instrument spin ID
2145  // int32_t *spinID = new int32_t[isc];
2146 
2147  for (size_t i=0; i<isc; i++) {
2148  uint32_t ui32;
2149  memcpy( &ui32, &ancdata[i*ANCSIZE+toff], 4);
2150  spinID[i] = SWAP_4( ui32);
2151  }
2152  var = gid.getVar("spin_ID");
2153  var.putVar(start, count, spinID);
2154 
2155 
2156  // Packet and line number sequence error flags
2157  var = gid.getVar("pseq_flag");
2158  var.putVar(start, count, seqerr);
2159  var = gid.getVar("line_flag");
2160  var.putVar(start, count, linerr);
2161 
2162  delete[] stimes;
2163  delete[] scss;
2164  delete[] scsu;
2165  delete[] hamSide;
2166  //delete[] spinID;
2167 
2168  return 0;
2169 }
2170 
2171 
2172 int l1aFile::write_oci_ancil_data( uint32_t isc, uint8_t *ancdata) {
2173 
2174  vector<size_t> start;
2175  vector<size_t> count;
2176  start.push_back(0);
2177 
2178  uint16_t ui16;
2179 
2180  // Extract ancillary telemetry status flags and error counts
2181  short ioff = 12;
2182  int16_t *anctlm = new int16_t[6*isc];
2183  for (size_t i=0;i<6*isc;i++) anctlm[i]=-999;
2184  for (size_t i=0; i<isc; i++) {
2185  for (size_t j=0; j<6; j++) {
2186  memcpy( &ui16, &ancdata[i*ANCSIZE+ioff+j*2], 2);
2187  anctlm[i*6+j] = SWAP_2( ui16);
2188  }
2189  }
2190 
2191  // Extract spatial aggregation / data type table
2192  short dtype[10];
2193  short iagg[10];
2194  short lines[10];
2195  short nagg[4] = {1,2,4,8};
2196  ioff = 36;
2197 
2198  int32_t iyear, iday;
2199  double stime;
2200  get_anc_packet_time( &ancdata[0], iyear, iday, stime);
2201  uint32_t jd = jday(iyear, 1, iday);
2202 
2203  for (size_t i=0; i<10; i++) {
2204  dtype[i] = ancdata[ioff+3] % 16;
2205  if ((jd < 2459000) && (dtype[i] == 11)) dtype[i] = 9; // data type mod for ETU before June 2020
2206  iagg[i] = ancdata[ioff+2] % 4;
2207  if (dtype[i] != 0) iagg[i] = nagg[iagg[i]];
2208  lines[i] = ancdata[ioff] * 256 + ancdata[ioff+1];
2209  ioff += 4;
2210  }
2211  ioff += 4;
2212 
2213  // Extract spectral aggregation and compute numbers of bands
2214 
2215  // Tap enable flags
2216  uint16_t btap = ancdata[ioff+2] * 256 + ancdata[ioff+3];
2217  uint16_t rtap = ancdata[ioff+0] * 256 + ancdata[ioff+1];
2218 
2219  // Tap aggregation factors
2220  uint32_t ui32;
2221  int32_t i32;
2222 
2223  memcpy(&ui32, &ancdata[ioff+8], 4);
2224  uint32_t bagg = SWAP_4(ui32);
2225 
2226  memcpy(&ui32, &ancdata[ioff+4], 4);
2227  uint32_t ragg = SWAP_4(ui32);
2228 
2229  int16_t baggs[16]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
2230  int16_t raggs[16]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
2231 
2232  // Compute number of bands for enabled taps
2233  uint16_t ken = 1;
2234  uint32_t kag = 3;
2235  uint32_t lag = 1;
2236 
2237  for (size_t i=0; i<16; i++) { // Put tap information in ascending spectral order
2238  uint16_t btaps = (btap & ken) / ken;
2239  if (btaps) baggs[15-i] = nagg[(bagg & kag) / lag];
2240  uint16_t rtaps = (rtap & ken) / ken;
2241  if (rtaps) raggs[15-i] = nagg[(ragg & kag) / lag];
2242 
2243  ken *= 2;
2244  kag *= 4;
2245  lag *= 4;
2246  }
2247 
2248 
2249  NcGroup gid;
2250  NcVar var;
2251 
2252  // Write to file
2253 
2254  gid = l1afile->getGroup( "spatial_spectral_modes");
2255 
2256  count.push_back(10);
2257 
2258  // Data type
2259  var = gid.getVar("spatial_zone_data_type");
2260  var.putVar(start, count, dtype);
2261 
2262  // Spatial aggregation
2263  var = gid.getVar("spatial_aggregation");
2264  var.putVar(start, count, iagg);
2265 
2266  // Number of lines
2267  var = gid.getVar("spatial_zone_lines");
2268  var.putVar(start, count, lines);
2269 
2270  count.pop_back();
2271  count.push_back(16);
2272 
2273  // Blue spectral aggregation
2274  var = gid.getVar("blue_spectral_mode");
2275  var.putVar(start, count, baggs);
2276 
2277  // Red spectral aggregation
2278  var = gid.getVar("red_spectral_mode");
2279  var.putVar(start, count, raggs);
2280 
2281 
2282  // Extract MCE data and write to file
2283  // Use ancillary packets AFTER science data
2284 
2285  gid = l1afile->getGroup( "engineering_data");
2286 
2287  count.pop_back();
2288  count.push_back(isc);
2289 
2290  // write ancillary telemetry status flags and error counts
2291  start.push_back(0);
2292  count.push_back(6);
2293  var = gid.getVar("ancillary_tlm");
2294  var.putVar(start, count, anctlm);
2295  delete[] anctlm;
2296  start.pop_back();
2297  count.pop_back();
2298 
2299  // Aggregation control
2300  int32_t *aggcon = new int32_t[isc];
2301  for (size_t i=0; i<isc; i++) {
2302  memcpy( &i32, &ancdata[i*ANCSIZE+ioff-4], 4);
2303  aggcon[i] = SWAP_4( i32);
2304  }
2305  var = gid.getVar("agg_control");
2306  var.putVar(start, count, aggcon);
2307  delete[] aggcon;
2308 
2309 
2310  // Aggregation errors
2311  uint16_t *bagerr = new uint16_t[isc];
2312  uint16_t *ragerr = new uint16_t[isc];
2313  for (size_t i=0; i<isc; i++) {
2314  memcpy( &ui16, &ancdata[(i+1)*ANCSIZE+ioff+12], 2);
2315  bagerr[i] = SWAP_2( ui16);
2316 
2317  memcpy( &ui16, &ancdata[(i+1)*ANCSIZE+ioff+14], 2);
2318  ragerr[i] = SWAP_2( ui16);
2319  }
2320  var = gid.getVar("blue_agg_error");
2321  var.putVar(start, count, bagerr);
2322  var = gid.getVar("red_agg_error");
2323  var.putVar(start, count, ragerr);
2324  delete[] bagerr;
2325  delete[] ragerr;
2326 
2327 
2328  // Digital card error status
2329  int32_t *digerr = new int32_t[isc];
2330  for (size_t i=0; i<isc; i++) {
2331  memcpy( &i32, &ancdata[(i+1)*ANCSIZE+ioff+20], 4);
2332  digerr[i] = SWAP_4( i32);
2333  }
2334  var = gid.getVar("dig_card_error");
2335  var.putVar(start, count, digerr);
2336  delete[] digerr;
2337 
2338  start.push_back(0);
2339  count.push_back(9);
2340 
2341  return 0;
2342 }
2343 
2344 
2345 int l1aFile::write_oci_tlm_data( itab *itable, uint32_t ntlm, uint8_t (*tlmdata)[TLMSIZE],
2346  int32_t *spinID, uint16_t &cdsmode, uint32_t isc) {
2347 // ntlm: Number of telemetry packets
2348 // tlmdata: OCI telemetry data packets
2349 // cdsmode: CDS mode (0 = enabled, 1 = reset, 2 = video)
2350 // tlmzs(2): Telmetry zone start times (msec)
2351 // tlmzd(2): Telmetry zone durations (msec)
2352 // tditime: CCD data line TDI time (clock cycles)
2353 
2354  // Set up output array and extract/convert data from DAU packets
2355  const uint16_t dauapid = 723;
2356  const uint16_t ddcapid = 701;
2357  const uint16_t mceapid = 713; // 711->713, changed by Liang Hong, 10/29/2020
2358  const uint16_t encapid = 712;
2359  const uint16_t scaapid = 717;
2360  const uint16_t senapid = 716;
2361  const uint16_t tcapid = 656;
2362  const uint16_t ddctapid = 703;
2363  const uint16_t dauctapid = 744;
2364  const uint16_t icdumcetapid = 745;
2365 
2366  //uint16_t nfptmps = 14;
2367  //uint16_t nsdtmps = 16;
2368  //uint16_t nabtmps = 9;
2369  //uint16_t ndtmps = 14;
2370  uint16_t nicthrm = 74;
2371  uint16_t ndctmps = 69;
2372  uint16_t nimtmps = 16;
2373  uint16_t nadclat = 4;
2374  uint16_t ndau = 0;
2375  uint16_t nddc = 0;
2376  uint16_t nmce = 0;
2377  uint16_t nenc = 0;
2378  uint16_t nsca = 0;
2379  uint16_t nsen = 0;
2380  uint16_t ntc = 0;
2381  uint16_t ndct = 0;
2382  uint16_t nddt = 0;
2383  uint16_t nimt = 0;
2384  int16_t *tlmzs = new int16_t[2];
2385  int16_t *tlmzd = new int16_t[2];
2386  uint16_t tditime = 0;
2387 
2388  double *dausec = new double[ntlm];
2389  for (size_t i=0;i<ntlm;i++) dausec[i]=-999;
2390  int32_t *dauspin = new int32_t[ntlm];
2391  for (size_t i=0;i<ntlm;i++) dauspin[i]=-999;
2392  uint8_t *dautlm = new uint8_t[620*ntlm];
2393  for (size_t i=0;i<620*ntlm;i++) dautlm[i]=0;
2394  double *ddcsec = new double[ntlm];
2395  for (size_t i=0;i<ntlm;i++) ddcsec[i]=-999;
2396  uint8_t *ddctlm = new uint8_t[524*ntlm];
2397  for (size_t i=0;i<524*ntlm;i++) ddctlm[i]=0;
2398  int32_t *mcespin = new int32_t[ntlm];
2399  for (size_t i=0;i<ntlm;i++) mcespin[i]=-999;
2400  uint8_t *mcetlm = new uint8_t[480*ntlm];
2401  for (size_t i=0;i<480*ntlm;i++) mcetlm[i]=0;
2402  int32_t *encspin = new int32_t[ntlm];
2403  for (size_t i=0;i<ntlm;i++) encspin[i]=-999;
2404  int16_t *encoder = new int16_t[4*200*ntlm];
2405  for (size_t i=0;i<4*200*ntlm;i++) encoder[i]=-999;
2406  int16_t *auxparms = new int16_t[33];
2407  for (size_t i=0; i<33;i++) auxparms[i] = -999;
2408  int32_t *scaspin = new int32_t[ntlm];
2409  for (size_t i=0;i<ntlm;i++) scaspin[i]=-999;
2410  uint8_t *scatlm = new uint8_t[480*ntlm];
2411  for (size_t i=0;i<480*ntlm;i++) scatlm[i]=0;
2412  int32_t *senspin = new int32_t[ntlm];
2413  for (size_t i=0;i<ntlm;i++) senspin[i]=-999;
2414  int16_t *sencoder = new int16_t[4*200*ntlm];
2415  for (size_t i=0;i<4*200*ntlm;i++) sencoder[i]=-999;
2416  double *tcsec = new double[ntlm];
2417  for (size_t i=0;i<ntlm;i++) tcsec[i]=-999;
2418  uint8_t *tctlm = new uint8_t[1216*ntlm];
2419  for (size_t i=0;i<1216*ntlm;i++) tctlm[i]=0;
2420  double *dauctsec = new double[ntlm];
2421  for (size_t i=0;i<ntlm;i++) dauctsec[i]=-999;
2422  double *icdumcetsec = new double[ntlm];
2423  for (size_t i=0;i<ntlm;i++) icdumcetsec[i]=-999;
2424  uint8_t *icdumcettlm = new uint8_t[76*ntlm];
2425  for (size_t i=0;i<76*ntlm;i++) icdumcettlm[i]=0;
2426  //uint16_t *redtemp = new uint16_t[nfptmps*ntlm];
2427  //for (size_t i=0;i<nfptmps*ntlm;i++) redtemp[i]=0;
2428  //uint16_t *blutemp = new uint16_t[nfptmps*ntlm];
2429  //for (size_t i=0;i<nfptmps*ntlm;i++) blutemp[i]=0;
2430  //uint16_t *sdstemp = new uint16_t[nsdtmps*ntlm];
2431  //for (size_t i=0;i<nsdtmps*ntlm;i++) sdstemp[i]=0;
2432  //uint16_t *aobtemp = new uint16_t[nabtmps*ntlm];
2433  //for (size_t i=0;i<nabtmps*ntlm;i++) aobtemp[i]=0;
2434  //uint16_t *dautemp = new uint16_t[ndtmps*ntlm];
2435  //for (size_t i=0;i<ndtmps*ntlm;i++) dautemp[i]=0;
2436  float *icdutherm = new float[nicthrm*ntlm];
2437  for (size_t i=0;i<nicthrm*ntlm;i++) icdutherm[i]=-999;
2438  float *dauctemp = new float[ndctmps*ntlm];
2439  for (size_t i=0;i<ndctmps*ntlm;i++) dauctemp[i]=-999;
2440  float *icdumcetemp = new float[nimtmps*ntlm];
2441  for (size_t i=0;i<nimtmps*ntlm;i++) icdumcetemp[i]=-999;
2442  uint8_t *adclat = new uint8_t[nadclat*ntlm];
2443  for (size_t i=0;i<nadclat*ntlm;i++) adclat[i]=0;
2444  uint8_t *cdsdis = new uint8_t[ntlm];
2445  for (size_t i=0;i<ntlm;i++) cdsdis[i]=0;
2446  uint8_t *hamside = new uint8_t[ntlm];
2447  uint32_t *redmask = new uint32_t[16];
2448  uint32_t *bluemask = new uint32_t[16];
2449 
2450  for (size_t i=0; i<ntlm; i++) {
2451  uint32_t apid = ((uint8_t) tlmdata[i][0] % 8)*256 + (uint8_t) tlmdata[i][1];
2452  int32_t iy, idy;
2453  double sc;
2454  uint8_t cctime[8]={0,0,0,0,0,0,0,0};
2455  uint32_t ui32;
2456  // uint16_t ui16;
2457  int16_t i16;
2458  float f32;
2459  cctime[6] = 0;
2460  cctime[7] = 0;
2461 
2462  switch (apid) {
2463 
2464  case dauapid:
2465  memcpy(cctime, &tlmdata[i][6], 6);
2466  ccsds_sec_to_yds(cctime, &iy, &idy, &sc);
2467  dausec[ndau] = sc;
2468  memcpy( &ui32, &tlmdata[i][12], 4);
2469  dauspin[ndau] = SWAP_4(ui32);
2470  memcpy( &dautlm[ndau*620], &tlmdata[i][16], 620);
2471 
2472  // Extract and convert temperatures
2473  /*
2474  for (size_t j=0; j<4; j++) {
2475  memcpy( &ui16, &tlmdata[i][312+2*j], 2);
2476  redtemp[ndau*nfptmps+j] = SWAP_2(ui16);
2477 
2478  memcpy( &ui16, &tlmdata[i][376+2*j], 2);
2479  blutemp[ndau*nfptmps+j] = SWAP_2(ui16);
2480  }
2481  for (size_t j=0; j<3; j++) {
2482  memcpy( &ui16, &tlmdata[i][434+2*j], 2);
2483  redtemp[ndau*nfptmps+4+j] = SWAP_2(ui16);
2484 
2485  memcpy( &ui16, &tlmdata[i][466+2*j], 2);
2486  blutemp[ndau*nfptmps+4+j] = SWAP_2(ui16);
2487  }
2488  for (size_t j=0; j<7; j++) {
2489  memcpy( &ui16, &tlmdata[i][482+2*j], 2);
2490  redtemp[ndau*nfptmps+7+j] = SWAP_2(ui16);
2491 
2492  memcpy( &ui16, &tlmdata[i][546+2*j], 2);
2493  blutemp[ndau*nfptmps+7+j] = SWAP_2(ui16);
2494  }
2495 
2496  for (size_t j=0; j<nsdtmps; j++) {
2497  memcpy( &ui16, &tlmdata[i][162+2*j], 2);
2498  sdstemp[ndau*nsdtmps+j] = SWAP_2(ui16);
2499  }
2500 
2501  for (size_t j=0; j<nabtmps; j++) {
2502  memcpy( &ui16, &tlmdata[i][194+2*j], 2);
2503  aobtemp[ndau*nabtmps+j] = SWAP_2(ui16);
2504  }
2505 
2506  for (size_t j=0; j<4; j++) {
2507  memcpy( &ui16, &tlmdata[i][154+2*j], 2);
2508  dautemp[ndau*ndtmps+j] = SWAP_2(ui16);
2509  }
2510  for (size_t j=0; j<4; j++) {
2511  memcpy( &ui16, &tlmdata[i][212+2*j], 2);
2512  dautemp[ndau*ndtmps+4+j] = SWAP_2(ui16);
2513  }
2514  for (size_t j=0; j<3; j++) {
2515  memcpy( &ui16, &tlmdata[i][402+2*j], 2);
2516  dautemp[ndau*ndtmps+8+j] = SWAP_2(ui16);
2517  }
2518  for (size_t j=0; j<3; j++) {
2519  memcpy( &ui16, &tlmdata[i][608+2*j], 2);
2520  dautemp[ndau*ndtmps+11+j] = SWAP_2(ui16);
2521  }
2522  */
2523 
2524  ndau++;
2525  break;
2526 
2527  case ddcapid:
2528  // DDC telemetry
2529  memcpy(cctime, &tlmdata[i][6], 6);
2530  ccsds_sec_to_yds(cctime, &iy, &idy, &sc);
2531  ddcsec[nddc] = sc;
2532  memcpy( &ddctlm[nddc*524], &tlmdata[i][12], 524);
2533  cdsdis[nddc] = tlmdata[i][29];
2534  memcpy( &adclat[nddc*nadclat], &tlmdata[i][176], nadclat);
2535  if (nddc==0) {
2536  for (size_t j=0; j<16; j++) {
2537  memcpy( &ui32, &tlmdata[i][196+j*4], 4);
2538  redmask[j] = (uint32_t) SWAP_4(ui32);
2539 
2540  memcpy( &ui32, &tlmdata[i][260+j*4], 4);
2541  bluemask[j] = (uint32_t) SWAP_4(ui32);
2542  }
2543  }
2544 
2545  nddc++;
2546  break;
2547 
2548  case mceapid:
2549  // RTA/HAM MCE telemetry
2550  memcpy( &ui32, &tlmdata[i][12], 4);
2551  mcespin[nmce] = SWAP_4(ui32);
2552  memcpy( &mcetlm[nmce*480], &tlmdata[i][16], 480);
2553  hamside[nmce] = (tlmdata[i][49] & 8) / 8;
2554 
2555  nmce++;
2556  break;
2557 
2558  case encapid:
2559  // RTA/HAM MCE encoder data
2560  memcpy( &ui32, &tlmdata[i][12], 4);
2561  encspin[nenc] = SWAP_4(ui32);
2562 
2563  for (size_t j=0; j<4*200; j++) {
2564  memcpy( &i16, &tlmdata[i][16+2*j], 2);
2565  encoder[nenc*4*200+j] = SWAP_2(i16);
2566  }
2567 
2568  nenc++;
2569  break;
2570 
2571  case scaapid:
2572  // SCA MCE telemetry
2573  memcpy( &ui32, &tlmdata[i][12], 4);
2574  scaspin[nsca] = SWAP_4(ui32);
2575  memcpy( &scatlm[nsca*480], &tlmdata[i][16], 480);
2576 
2577  nsca++;
2578  break;
2579 
2580  case senapid:
2581  // SCA MCE encoder data
2582  memcpy( &ui32, &tlmdata[i][12], 4);
2583  senspin[nsen] = SWAP_4(ui32);
2584  for (size_t j=0; j<4*200; j++) {
2585  memcpy( &i16, &tlmdata[i][16+2*j], 2);
2586  sencoder[nsen*4*200+j] = SWAP_2(i16);
2587  }
2588 
2589  nsen++;
2590  break;
2591 
2592  case tcapid:
2593  // ICDU TC telemetry and thermistors
2594  memcpy(cctime, &tlmdata[i][6], 6);
2595  ccsds_sec_to_yds(cctime, &iy, &idy, &sc);
2596  tcsec[ntc] = sc;
2597  memcpy( &tctlm[ntc*1216], &tlmdata[i][12], 1216);
2598 
2599  for (size_t j=0; j<nicthrm; j++) {
2600  memcpy( &f32, &tlmdata[i][184+4*j], 4);
2601  icdutherm[ntc*nicthrm+j] = ReverseFloat(f32);
2602  }
2603 
2604  ntc++;
2605  break;
2606 
2607  case ddctapid:
2608  // DDC table telemetry
2609  // Only need one set of these
2610  if (nddt == 0) {
2611  for (size_t j=0; j<33; j++) {
2612  memcpy( &ui32, &tlmdata[i][64+j*4], 4);
2613  auxparms[j] = (int16_t) SWAP_4(ui32);
2614  }
2615  nddt = 1;
2616  }
2617  break;
2618 
2619  case dauctapid:
2620  // FSW DAUC temperatures
2621  memcpy(cctime, &tlmdata[i][6], 6);
2622  ccsds_sec_to_yds(cctime, &iy, &idy, &sc);
2623  dauctsec[ndct] = sc;
2624 
2625  for (size_t j=0; j<ndctmps; j++) {
2626  memcpy( &f32, &tlmdata[i][12+4*j], 4);
2627  dauctemp[ndct*ndctmps+j] = ReverseFloat(f32);
2628  }
2629  ndct++;
2630  break;
2631 
2632  case icdumcetapid:
2633  // FSW ICCU/MCE telemetry and temperatures
2634  memcpy(cctime, &tlmdata[i][6], 6);
2635  ccsds_sec_to_yds(cctime, &iy, &idy, &sc);
2636  icdumcetsec[nimt] = sc;
2637  memcpy( &icdumcettlm[nimt*76], &tlmdata[i][12], 76);
2638  for (size_t j=0; j<nimtmps; j++) {
2639  memcpy( &f32, &tlmdata[i][24+4*j], 4);
2640  icdumcetemp[nimt*nimtmps+j] = ReverseFloat(f32);
2641  }
2642  nimt++;
2643  break;
2644 
2645  }
2646  } // tlm loop
2647 
2648  NcGroup gid;
2649  NcVar var;
2650 
2651  vector<size_t> start;
2652  vector<size_t> count;
2653 
2654  // Write to file
2655 
2656  gid = l1afile->getGroup( "engineering_data");
2657 
2658  if (ndau==0) ndau = 1; // 0.99.20
2659 
2660  //if (ndau > 0) {
2661  start.push_back(0);
2662  count.push_back(ndau);
2663 
2664  var = gid.getVar("DAU_tlm_time");
2665  var.putVar(start, count, dausec);
2666 
2667  var = gid.getVar("DAU_spin_ID");
2668  var.putVar(start, count, dauspin);
2669 
2670  start.push_back(0);
2671  count.push_back(620);
2672 
2673  var = gid.getVar("DAU_telemetry");
2674  var.putVar(start, count, dautlm);
2675 
2676  //count.pop_back();
2677  //count.push_back(nfptmps);
2678 
2679  //var = gid.getVar("blue_FPA_temperatures");
2680  //var.putVar(start, count, blutemp);
2681  //var = gid.getVar("red_FPA_temperatures");
2682  //var.putVar(start, count, redtemp);
2683 
2684  //count.pop_back();
2685  //count.push_back(nsdtmps);
2686 
2687  //var = gid.getVar("SDS_det_temperatures");
2688  //var.putVar(start, count, sdstemp);
2689 
2690  //count.pop_back();
2691  //count.push_back(nabtmps);
2692 
2693  //var = gid.getVar("AOB_temperatures");
2694  //var.putVar(start, count, aobtemp);
2695 
2696  //count.pop_back();
2697  //count.push_back(ndtmps);
2698 
2699  //var = gid.getVar("DAU_temperatures");
2700  //var.putVar(start, count, dautemp);
2701 
2702  // Get telemetry zone start times and durations
2703  tlmzs[0] = dautlm[123];
2704  tlmzs[1] = dautlm[125];
2705  tlmzd[0] = dautlm[122];
2706  tlmzd[1] = dautlm[124];
2707  //}
2708 
2709  if (nddc==0) {
2710  for (size_t i=0;i<16;i++) {
2711  redmask[i] = 4294967295;
2712  bluemask[i] = 4294967295;
2713  }
2714  nddc = 1; // 0.99.20
2715  }
2716  //if (nddc > 0) {
2717  start.clear();
2718  count.clear();
2719  start.push_back(0);
2720  count.push_back(nddc);
2721 
2722  // DDC tlm time
2723  var = gid.getVar("DDC_tlm_time");
2724  var.putVar(start, count, ddcsec);
2725 
2726  // CDS disable flag
2727  var = gid.getVar("CDS_disable");
2728  var.putVar(start, count, cdsdis);
2729 
2730  start.push_back(0);
2731  count.push_back(524);
2732 
2733  // DDC telemetry
2734  var = gid.getVar("DDC_telemetry");
2735  var.putVar(start, count, ddctlm);
2736 
2737  count.pop_back();
2738  count.push_back(nadclat);
2739 
2740  // ADC latency
2741  var = gid.getVar("ADC_latency");
2742  var.putVar(start, count, adclat);
2743 
2744  // Red and blue channel masks
2745  start.clear();
2746  count.clear();
2747  start.push_back(0);
2748  count.push_back(16);
2749  var = gid.getVar("blue_channel_mask");
2750  var.putVar(start, count, bluemask);
2751  var = gid.getVar("red_channel_mask");
2752  var.putVar(start, count, redmask);
2753 
2754  // Get CDS mode for metadata
2755  cdsmode = cdsdis[0] * (adclat[0]-14);
2756 
2757  // Get TDI time
2758  uint16_t ui16;
2759  memcpy( &ui16, &ddctlm[346], 2);
2760  tditime = SWAP_2(ui16);
2761  //}
2762 
2763  if (nmce==0) nmce=1; // 0.99.20
2764  //if (nmce > 0) {
2765  start.clear();
2766  count.clear();
2767  start.push_back(0);
2768  count.push_back(nmce);
2769 
2770  // RTA/HAM MCE spin ID
2771  var = gid.getVar("MCE_spin_ID");
2772  var.putVar(start, count, mcespin);
2773 
2774  // HAM side
2775  //if (nmce > 1) {
2776  uint8_t *mside = new uint8_t[isc];
2777  for (size_t i=0; i<isc; i++) mside[i] = 255;
2778  for (size_t i=0; i<isc; i++) {
2779  for (size_t j=0; j<ntlm; j++) {
2780  if ( (int) mcespin[j] == spinID[i]) {
2781  mside[i] = hamside[j];
2782  break;
2783  }
2784  }
2785  }
2786  count.pop_back();
2787  count.push_back(isc); // LH, 11/20/2020
2788  NcGroup scgid = l1afile->getGroup( "scan_line_attributes");
2789 
2790  var = scgid.getVar("HAM_side");
2791  // var.putVar(start, count, &hamside[1]);
2792  var.putVar(start, count, mside);
2793  delete[] mside;
2794  //}
2795 
2796  start.push_back(0);
2797  count.clear(); // LH, 11/20/2020
2798  count.push_back(nmce); // 1.00.00
2799  count.push_back(480);
2800 
2801  // RTA/HAM MCE telemetry
2802  var = gid.getVar("MCE_telemetry");
2803  var.putVar(start, count, mcetlm);
2804  //}
2805 
2806  if (nenc==0) nenc = 1; // 0.99.20
2807  //if (nenc > 0) {
2808  start.clear();
2809  count.clear();
2810  start.push_back(0);
2811  count.push_back(nenc);
2812 
2813  // RTA/HAM encoder spin ID
2814  var = gid.getVar("encoder_spin_ID");
2815  var.putVar(start, count, encspin);
2816 
2817  start.push_back(0);
2818  count.push_back(200);
2819  start.push_back(0);
2820  count.push_back(4);
2821 
2822  // RTA/HAM encoder data
2823  var = gid.getVar("MCE_encoder_data");
2824  var.putVar(start, count, encoder);
2825  //}
2826 
2827 
2828  if (nsca==0) nsca = 1; // 1.01.00
2829  start.clear();
2830  count.clear();
2831  start.push_back(0);
2832  count.push_back(nsca);
2833  // SCA MCE spin ID
2834  var = gid.getVar("SCA_spin_ID");
2835  var.putVar(start, count, scaspin);
2836 
2837  // SCA MCE telemetry
2838  start.clear();
2839  start.push_back(0);
2840  count.clear();
2841  count.push_back(nsca);
2842  count.push_back(480);
2843 
2844  var = gid.getVar("SCA_telemetry");
2845  var.putVar(start, count, scatlm);
2846 
2847  if (nsen==0) nsen = 1;
2848  start.clear();
2849  count.clear();
2850  start.push_back(0);
2851  count.push_back(nsen);
2852 
2853  // RTA/HAM encoder spin ID
2854  var = gid.getVar("SCA_encoder_spin_ID");
2855  var.putVar(start, count, senspin);
2856 
2857  // SCA encoder data
2858  start.push_back(0);
2859  count.push_back(200);
2860  start.push_back(0);
2861  count.push_back(4);
2862 
2863  // RTA/HAM encoder data
2864  var = gid.getVar("SCA_encoder_data");
2865  var.putVar(start, count, sencoder);
2866 
2867  if (ntc==0) ntc = 1; // 0.99.20
2868  //if (ntc > 0) {
2869  start.clear();
2870  count.clear();
2871  start.push_back(0);
2872  count.push_back(ntc);
2873 
2874  // TC tlm time
2875  var = gid.getVar("TC_tlm_time");
2876  var.putVar(start, count, tcsec);
2877 
2878  start.push_back(0);
2879  count.push_back(1216);
2880 
2881  // TC telemetry
2882  var = gid.getVar("TC_telemetry");
2883  var.putVar(start, count, tctlm);
2884 
2885  count.pop_back();
2886  count.push_back(nicthrm);
2887 
2888  // ICDU thermisters
2889  var = gid.getVar("ICDU_thermisters");
2890  var.putVar(start, count, icdutherm);
2891  //}
2892 
2893  // Auxilary parameter table
2894  gid = l1afile->getGroup( "spatial_spectral_modes");
2895 
2896  //if (nddt > 0) {
2897  start.clear();
2898  count.clear();
2899  start.push_back(0);
2900  count.push_back(33);
2901 
2902  var = gid.getVar("aux_param_table");
2903  var.putVar(start, count, auxparms);
2904  //}
2905 
2906  if (ndct==0) ndct=1;
2907  start.clear();
2908  count.clear();
2909  start.push_back(0);
2910  count.push_back(ndct);
2911 
2912  gid = l1afile->getGroup( "engineering_data");
2913 
2914  // FSW DAUC temperature time
2915  var = gid.getVar("DAUC_temp_time");
2916  var.putVar(start, count, dauctsec);
2917 
2918  start.push_back(0);
2919  count.push_back(ndctmps);
2920 
2921  // FSW DAUC temperatures
2922  var = gid.getVar("DAUC_temperatures");
2923  var.putVar(start, count, dauctemp);
2924 
2925  if (nimt==0) nimt=1;
2926  start.clear();
2927  count.clear();
2928  start.push_back(0);
2929  count.push_back(nimt);
2930 
2931  // FSW ICDU/MCE temperature time
2932  var = gid.getVar("ICDU_MCE_temp_time");
2933  var.putVar(start, count, icdumcetsec);
2934 
2935  start.push_back(0);
2936  count.push_back(76);
2937 
2938  // FSW ICDU/MCE temperature telemetry
2939  var = gid.getVar("ICDU_MCE_temp_tlm");
2940  var.putVar(start, count, icdumcettlm);
2941 
2942  count.pop_back();
2943  count.push_back(nimtmps);
2944 
2945  // FSW ICDU/MCE temperatures
2946  var = gid.getVar("ICDU_MCE_temperatures");
2947  var.putVar(start, count, icdumcetemp);
2948 
2949 
2950  // check the science data and telemetry collection zones for conflicts and write the results
2951  // as an attribute to the engineering data group. There are two telemetry collection zones
2952  string conflict = "No";
2953  gid = l1afile->getGroup("engineering_data");
2954 
2955  if ((tlmzd[0] == 0) && (tlmzd[1] == 0)) {
2956  cout<<"No telemetry zone fields in file"<<endl;
2957  gid.putAtt("science_telemetry_zone_conflict", conflict);
2958  } else {
2959 
2960  float clock = 136000; // Master clock frequency in msec
2961  float secpline = (tditime + 1)/clock; // msec per line
2962 
2963  // Find no-data zones
2964  float *znd = new float[2*10];
2965  float *zndcp = new float[2*10];
2966  int ndz = 0;
2967  uint16_t line = 0;
2968  for (size_t i=0; i<10; i++) {
2969  if (((itable[i].dtype == 0) || (itable[i].dtype == 10)) && (itable[i].lines > 0)) {
2970  znd[ndz] = line*secpline;
2971  znd[1*10+ndz] = (line + itable[i].lines)*secpline;
2972  //Check for consecutive no-data zones
2973  if ((ndz > 0) && (znd[ndz] == znd[1*10+ndz-1])) {
2974  znd[1*10+ndz-1] = znd[1*10+ndz];
2975  } else {
2976  ndz++;
2977  }
2978  }
2979  line += itable[i].lines;
2980  }
2981  memcpy(zndcp,znd,20);
2982  // znd new dimension: 2*ndz
2983  for (int i=0;i<ndz;i++) znd[ndz+i]=znd[10+i]; // znd = znd(*,0:ndz-1)
2984 
2985  // If first no-data zone is at start of spin, extend last zone by that amount.
2986  if ((itable[0].dtype == 0) || (itable[0].dtype == 10)) znd[1*ndz+ndz-1] += znd[1*ndz+0];
2987 
2988  if (ndz > 0) {
2989  // Loop through telemetry zones
2990  for (size_t i=0; i<2; i++) {
2991  // Find overlapping no-data zone
2992  int16_t tlmze = tlmzs[i] + tlmzd[i] + 3;
2993  int k;
2994  for (k=ndz-1; k>=0; k--) {
2995  // k = max(where(tlmze gt znd(0,*)))
2996  if (tlmze>znd[k]) break;
2997  }
2998  if ((k == -1) || (tlmzs[i] < znd[k]) || (tlmze > znd[1*ndz+k])) {
2999  conflict = "Yes";
3000  cout<<"Telemetry zone "<< tlmzs[i]<<", " << tlmze<<endl;
3001  cout<<"No-data zone "<<znd[k]<<", "<<znd[1*ndz+k]<<endl;
3002  }
3003  }
3004  } else conflict = "Yes";
3005 
3006  // Write attribute to engineering_data_group
3007  // Open group
3008  NcGroup gid = l1afile->getGroup("engineering_data");
3009  gid.putAtt("science_telemetry_zone_conflict", conflict);
3010 
3011  delete[] znd;
3012  delete[] zndcp;
3013  } //if (tlmzd[0] == 0) && (tlmzd[1] == 0) else
3014 
3015 
3016  delete[] dausec;
3017  delete[] dauspin;
3018  delete[] dautlm;
3019  delete[] ddcsec;
3020  delete[] ddctlm;
3021  delete[] mcespin;
3022  delete[] mcetlm;
3023  delete[] encspin;
3024  delete[] encoder;
3025  delete[] auxparms;
3026  delete[] tcsec;
3027  delete[] tctlm;
3028  delete[] dauctsec;
3029  delete[] icdumcetsec;
3030  delete[] icdumcettlm;
3031  //delete[] redtemp;
3032  //delete[] blutemp;
3033  //delete[] sdstemp;
3034  //delete[] aobtemp;
3035  //delete[] dautemp;
3036  delete[] icdutherm;
3037  delete[] adclat;
3038  delete[] cdsdis;
3039  delete[] hamside;
3040  delete[] tlmzs;
3041  delete[] tlmzd;
3042 
3043  return 0;
3044 }
3045 
3046 int l1aFile::write_navigation(std::string hktlist,
3047  time_struct &starttime,
3048  time_struct &endtime) {
3049  uint16_t nnavmax = 1000;
3050  size_t nnav = 0;
3051 
3052  double usec_l1a_start, usec_l1a_end, usec_hkt;
3053  usec_l1a_start = yds2unix(starttime.iyear,starttime.iday,starttime.sec);
3054  usec_l1a_end = yds2unix(endtime.iyear,endtime.iday,endtime.sec);
3055 
3056  // read HKT file list, loop through HKT files
3057  ifstream file(hktlist);
3058  string strHKTfile;
3059 
3060  double *atime = new double[nnavmax];
3061  for (size_t i=0;i<nnavmax;i++) atime[i]=-999;
3062  double *otime = new double[nnavmax];
3063  for (size_t i=0;i<nnavmax;i++) otime[i]=-999;
3064  double *ttime = new double[nnavmax];
3065  for (size_t i=0;i<nnavmax;i++) ttime[i]=-999;
3066 
3067  float *arate = new float[3*nnavmax];
3068  for (size_t i=0;i<3*nnavmax;i++) arate[i]=-999;
3069  float *quat = new float[4*nnavmax];
3070  for (size_t i=0;i<4*nnavmax;i++) quat[i]=-999;
3071  float *pos = new float[3*nnavmax];
3072  for (size_t i=0;i<3*nnavmax;i++) pos[i]=-9999999;
3073  float *vel = new float[3*nnavmax];
3074  for (size_t i=0;i<3*nnavmax;i++) vel[i]=-9999999;
3075  float *tlt = new float[nnavmax];
3076  for (size_t i=0;i<nnavmax;i++) tlt[i]=-999;
3077 
3078  NcVar var;
3079 
3080  while (getline(file, strHKTfile)) {
3081  try {
3082  // read HKT netcdf file
3083  int ncid, gid;
3084  int atttid, attqid,attrid;
3085  int orbpid, orbtid, orbvid;
3086  int tid, ttid;
3087  size_t nHKTlen=0;
3088  char *hkt_t_start = (char*)malloc(100 * sizeof(char));
3089 
3090  nc_open(strHKTfile.c_str(), NC_NOWRITE, &ncid);
3091  cout<<"Reading PACE HKT file "<<strHKTfile<<"..."<<endl;
3092 
3093  // find HKT data year, month, date from global attributes
3094  nc_get_att_text(ncid, NC_GLOBAL, "time_coverage_start", hkt_t_start); // e.g. 2022-04-20T17:26:00.000000Z
3095 
3096  int hkt_yr, hkt_mon, hkt_day;
3097  sscanf(hkt_t_start, "%4d-%2d-%2d", &hkt_yr, &hkt_mon, &hkt_day);
3098 
3099  nc_inq_grp_ncid(ncid,"navigation_data",&gid);
3100 
3101  nc_type attt_type;
3102  int attt_ndims;
3103  int attt_dimids[NC_MAX_VAR_DIMS];
3104  int attt_natts;
3105  nc_inq_varid (gid, "att_time", &atttid);
3106  nc_inq_var (gid, atttid, 0, &attt_type, &attt_ndims, attt_dimids,&attt_natts);
3107  nc_inq_dimlen(gid, attt_dimids[0], &nHKTlen);
3108  if (nHKTlen<1) {
3109  nc_close(ncid);
3110  continue;
3111  }
3112 
3113  // read data from HKT file
3114  double *at = new double[nHKTlen];
3115  double *ot = new double[nHKTlen];
3116  double *tt = new double[nHKTlen];
3117 
3118  float *r = new float[3*nHKTlen];
3119  float *q = new float[4*nHKTlen];
3120  float *p = new float[3*nHKTlen];
3121  float *v = new float[3*nHKTlen];
3122  float *t = new float[nHKTlen];
3123 
3124  nc_get_var(gid,atttid,at);
3125 
3126  // find values within L1A file time range, i.e. [starttime-10 sec, endtime+10 sec]
3127  int ind_start=1e6, ind_end=-1;
3128  for (size_t i=0;i<nHKTlen;i++) {
3129  usec_hkt = ymds2unix(hkt_yr,hkt_mon,hkt_day,at[i]);
3130  if ((usec_hkt>usec_l1a_start-10) && (ind_start==1e6)) {
3131  ind_start =i;
3132  break;
3133  }
3134  }
3135  for (size_t i=nHKTlen-1;i>=0;i--) {
3136  usec_hkt = ymds2unix(hkt_yr,hkt_mon,hkt_day,at[i]);
3137  if ((usec_hkt<usec_l1a_end+10) && (ind_end==-1)) {
3138  ind_end =i;
3139  break;
3140  }
3141  }
3142 
3143  if ((ind_end - ind_start)>0) {
3144  nc_inq_varid (gid, "att_quat", &attqid);
3145  nc_get_var(gid,attqid,q);
3146  nc_inq_varid (gid, "att_rate", &attrid);
3147  nc_get_var(gid,attrid,r);
3148 
3149  nc_inq_varid (gid, "orb_pos", &orbpid);
3150  nc_get_var(gid,orbpid,p);
3151  nc_inq_varid (gid, "orb_time", &orbtid);
3152  nc_get_var(gid,orbtid,ot);
3153  nc_inq_varid (gid, "orb_vel", &orbvid);
3154  nc_get_var(gid,orbvid,v);
3155 
3156  nc_inq_varid (gid, "tilt", &tid);
3157  nc_get_var(gid,tid,t);
3158  nc_inq_varid (gid, "tilt_time", &ttid);
3159  nc_get_var(gid,ttid,tt);
3160 
3161  // concatenate arrays
3162  size_t nnav0 = ind_end - ind_start + 1;
3163  memcpy(atime+nnav,at+ind_start,nnav0*sizeof(double));
3164  memcpy(otime+nnav,ot+ind_start,nnav0*sizeof(double));
3165  memcpy(ttime+nnav,tt+ind_start,nnav0*sizeof(double));
3166  memcpy(tlt+nnav,t+ind_start,nnav0*sizeof(float));
3167 
3168  memcpy(quat+nnav*4,q+ind_start*4,4*nnav0*sizeof(float));
3169  memcpy(arate+nnav*3,r+ind_start*3,3*nnav0*sizeof(float));
3170  memcpy(pos+nnav*3,p+ind_start*3,3*nnav0*sizeof(float));
3171  memcpy(vel+nnav*3,v+ind_start*3,3*nnav0*sizeof(float));
3172 
3173  nnav += nnav0;
3174 
3175  } // if ((ind_end - ind_start)>0)
3176 
3177  nc_close(ncid);
3178 
3179  delete[] hkt_t_start;
3180  delete[] at;
3181  delete[] ot;
3182  delete[] tt;
3183  delete[] r;
3184  delete[] q;
3185  delete[] p;
3186  delete[] v;
3187  delete[] t;
3188  }
3189  catch (NcException& e){
3190  //e.what();
3191  cout<<"Error reading "<<strHKTfile<<"!"<<endl;
3192  return NC2_ERR;
3193  }
3194 
3195  } // while (getline(file, strHKTfile))
3196 
3197  // write to L1A file
3198  vector<size_t> start;
3199  vector<size_t> count;
3200 
3201  NcGroup l1a_gid = l1afile->getGroup("navigation_data");
3202 
3203  if (nnav==0) nnav = 1;
3204  start.push_back(0);
3205  count.push_back(nnav);
3206  var = l1a_gid.getVar("att_time");
3207  var.putVar(start, count, atime);
3208 
3209  var = l1a_gid.getVar("orb_time");
3210  var.putVar(start, count, otime);
3211 
3212  var = l1a_gid.getVar("tilt_time");
3213  var.putVar(start, count, ttime);
3214 
3215  var = l1a_gid.getVar("tilt");
3216  var.putVar(start, count, tlt);
3217 
3218  start.clear();
3219  count.clear();
3220 
3221  start.push_back(0);
3222  start.push_back(0);
3223  count.push_back(nnav);
3224  count.push_back(4);
3225  var = l1a_gid.getVar("att_quat");
3226  var.putVar(start, count, quat);
3227 
3228  count.pop_back();
3229  count.push_back(3);
3230  var = l1a_gid.getVar("att_rate");
3231  var.putVar(start, count, arate);
3232 
3233  var = l1a_gid.getVar("orb_pos");
3234  var.putVar(start, count, pos);
3235 
3236  var = l1a_gid.getVar("orb_vel");
3237  var.putVar(start, count, vel);
3238 
3239  delete[] atime;
3240  delete[] otime;
3241  delete[] ttime;
3242  delete[] arate;
3243  delete[] quat;
3244  delete[] pos;
3245  delete[] vel;
3246  delete[] tlt;
3247 
3248  return 0;
3249 }
3250 
3251 
3252 
3253 int l1aFile::write_oci_global_metadata( time_struct &starttime,
3254  time_struct &endtime,
3255  std::string l1a_name,
3256  std::string sdir, std::string edir,
3257  uint32_t isc, short dtype,
3258  uint16_t smode, uint16_t cdsmode,
3259  std::ofstream &fout) {
3260 
3261  string dtypes[] = {"", "Earth Collect", "Dark Collect", "Solar Cal Daily",
3262  "Solar Cal Monthly", "Response Curve", "Lunar Cal",
3263  "Diagnostic","Static", "Earth Spectral","",
3264  "External Snapshot Trigger","Internal Snapshot Trigger",
3265  "SPCA","Lunar Stare"};
3266 
3267  string smodes[] = {"Science", "Diagnostic", "Single-image raw",
3268  "Test pattern"};
3269 
3270  string cdsmodes[] = {"CDS", "Reset", "Video"};
3271 
3272  // Date created
3273  char buf[32];
3274  strcpy( buf, unix2isodate(now(), 'G'));
3275  l1afile->putAtt("date_created", buf);
3276 
3277  // Write start, end, create time attributes
3278  int16_t mon, idm;
3279  int32_t ih, mn;
3280  stringstream ss;
3281 
3282  yd2md((int16_t) starttime.iyear, (int16_t) starttime.iday, &mon, &idm);
3283 
3284  starttime.sec = floor(starttime.sec*1000)/1000; // LH, 11/18/2020
3285  ih = (int) (starttime.sec/3600);
3286  starttime.sec -= ih * 3600;
3287  mn = (int) (starttime.sec/60);
3288  starttime.sec -= mn * 60;
3289 
3290  // yyyy-mn-dyThr:mn:ss.sss
3291  ss = stringstream();
3292  ss << setw(4) << to_string(starttime.iyear) << "-";
3293  ss << setw(2) << setfill('0') << mon << "-";
3294  ss << setw(2) << setfill('0') << idm << "T";
3295 
3296  ss << setw(2) << setfill('0') << ih << ":";
3297  ss << setw(2) << setfill('0') << mn << ":";
3298  ss << fixed << setw(6) << setprecision(3) <<
3299  setfill('0') << starttime.sec;
3300 
3301  cout << ss.str() << endl;
3302  l1afile->putAtt("time_coverage_start", ss.str().c_str());
3303 
3304  // Write to outlist if needed
3305  if ( fout.is_open())
3306  fout << ss.str().c_str() << " ";
3307 
3308  yd2md((int16_t) endtime.iyear, (int16_t) endtime.iday, &mon, &idm);
3309 
3310  endtime.sec = floor(endtime.sec*1000)/1000; // LH, 11/18/2020, to handle 2020-11-05T20:06:59.999980000
3311  ih = (int) (endtime.sec/3600);
3312  endtime.sec -= ih * 3600;
3313  mn = (int) (endtime.sec/60);
3314  endtime.sec -= mn * 60;
3315 
3316  // yyyy-mn-dyThr:mn:ss.sss
3317  ss = stringstream();
3318  ss << setw(4) << to_string(endtime.iyear) << "-";
3319  ss << setw(2) << setfill('0') << mon << "-";
3320  ss << setw(2) << setfill('0') << idm << "T";
3321 
3322  ss << setw(2) << setfill('0') << ih << ":";
3323  ss << setw(2) << setfill('0') << mn << ":";
3324  ss << fixed << setw(6) << setprecision(3) <<
3325  setfill('0') << endtime.sec;
3326 
3327  cout << ss.str() << endl;
3328  l1afile->putAtt("time_coverage_end", ss.str().c_str());
3329 
3330  // Write product file name
3331  l1afile->putAtt("product_name", l1a_name.c_str());
3332 
3333  // Write orbit direction
3334  l1afile->putAtt("startDirection", sdir.c_str());
3335  l1afile->putAtt("endDirection", edir.c_str());
3336 
3337  // Write number of scans
3338  l1afile->putAtt("number_of_filled_scans", to_string(isc));
3339 
3340  // Write data collect type and SWIR mode
3341  l1afile->putAtt("data_collect_mode", dtypes[dtype].c_str());
3342  l1afile->putAtt("SWIR_data_mode", smodes[smode].c_str());
3343  l1afile->putAtt("CDS_mode", cdsmodes[cdsmode].c_str());
3344 
3345  // Write to outlist if needed
3346  if ( fout.is_open())
3347  fout << ss.str().c_str() << " ";
3348 
3349  return 0;
3350 }
3351 
3352 int l1aFile::close() {
3353 
3354  try {
3355  l1afile->close();
3356  }
3357  catch ( NcException& e) {
3358  cout << e.what() << endl;
3359  cerr << "Failure closing: " + fileName << endl;
3360  exit(1);
3361  }
3362 
3363  return 0;
3364 }
3365 
3366 
3367 int createNCDF( NcGroup &ncGrp, const char *sname, const char *lname,
3368  const char *standard_name, const char *units,
3369  void *fill_value, const char *flag_values,
3370  const char *flag_meanings, const char *reference,
3371  double low, double high, int nt, vector<NcDim>& varVec) {
3372 
3373  size_t dimlength;
3374 
3375 
3376  /* Create the NCDF dataset */
3377  NcVar ncVar;
3378  try {
3379  ncVar = ncGrp.addVar(sname, nt, varVec);
3380  }
3381  catch ( NcException& e) {
3382  cout << e.what() << endl;
3383  cerr << "Failure creating variable: " << sname << endl;
3384  exit(1);
3385  }
3386 
3387  // Set fill value
3388  double fill_value_dbl;
3389  memcpy( &fill_value_dbl, fill_value, sizeof(double));
3390 
3391  int8_t i8;
3392  uint8_t ui8;
3393  int16_t i16;
3394  uint16_t ui16;
3395  int32_t i32;
3396  uint32_t ui32;
3397  float f32;
3398 
3399  if ( low != fill_value_dbl) {
3400  if ( nt == NC_BYTE) {
3401  i8 = fill_value_dbl;
3402  ncVar.setFill(true, (void *) &i8);
3403  } else if ( nt == NC_UBYTE) {
3404  ui8 = fill_value_dbl;
3405  ncVar.setFill(true, (void *) &ui8);
3406  } else if ( nt == NC_SHORT) {
3407  i16 = fill_value_dbl;
3408  ncVar.setFill(true, (void *) &i16);
3409  } else if ( nt == NC_USHORT) {
3410  ui16 = fill_value_dbl;
3411  ncVar.setFill(true, (void *) &ui16);
3412  } else if ( nt == NC_INT) {
3413  i32 = fill_value_dbl;
3414  ncVar.setFill(true, (void *) &i32);
3415  } else if ( nt == NC_UINT) {
3416  ui32 = fill_value_dbl;
3417  ncVar.setFill(true, (void *) &ui32);
3418  } else if ( nt == NC_FLOAT) {
3419  f32 = fill_value_dbl;
3420  ncVar.setFill(true, (void *) &f32);
3421  } else {
3422  ncVar.setFill(true, (void *) &fill_value_dbl);
3423  }
3424  }
3425 
3426  /* vary chunck size based on dimensions */
3427  int do_deflate = 0;
3428  vector<size_t> chunkVec;
3429  if ( varVec.size() == 3 && (strncmp(sname, "EV_", 3) == 0)) {
3430  dimlength = varVec[2].getSize();
3431 
3432  chunkVec.push_back(1);
3433  chunkVec.push_back(16);
3434  chunkVec.push_back(dimlength/10);
3435 
3436  do_deflate = 1;
3437  }
3438 
3439  /* Set compression */
3440  if ( do_deflate) {
3441  /* First set chunking */
3442  try {
3443  ncVar.setChunking(ncVar.nc_CHUNKED, chunkVec);
3444  }
3445  catch ( NcException& e) {
3446  e.what();
3447  cerr << "Failure setting chunking: " << sname << endl;
3448  exit(1);
3449  }
3450 
3451  try {
3452  ncVar.setCompression(true, true, 5);
3453  }
3454  catch ( NcException& e) {
3455  e.what();
3456  cerr << "Failure setting compression: " << sname << endl;
3457  exit(1);
3458  }
3459  }
3460 
3461 
3462  /* Add a "long_name" attribute */
3463  try {
3464  ncVar.putAtt("long_name", lname);
3465  }
3466  catch ( NcException& e) {
3467  e.what();
3468  cerr << "Failure creating 'long_name' attribute: " << lname << endl;
3469  exit(1);
3470  }
3471 
3472  if ( strcmp( flag_values, "") != 0) {
3473 
3474  size_t curPos=0;
3475 
3476  string fv;
3477  fv.assign( flag_values);
3478  size_t pos = fv.find("=", curPos);
3479  fv = fv.substr(pos+1);
3480 
3481  size_t semicln = fv.find(";");
3482  pos = 0;
3483 
3484  int8_t vec[1024];
3485  int n = 0;
3486  while(pos != semicln) {
3487  pos = fv.find(",", curPos);
3488  if ( pos == string::npos)
3489  pos = semicln;
3490 
3491  string flag_value;
3492  istringstream iss(fv.substr(curPos, pos-curPos));
3493  iss >> skipws >> flag_value;
3494  vec[n++] = atoi( flag_value.c_str());
3495  curPos = pos + 1;
3496  }
3497 
3498  try {
3499  ncVar.putAtt("flag_values", NC_BYTE, n, vec);
3500  }
3501  catch ( NcException& e) {
3502  e.what();
3503  cerr << "Failure creating 'flag_values' attribute: " << lname << endl;
3504  exit(1);
3505  }
3506  }
3507 
3508  /* Add a "flag_meanings" attribute if specified*/
3509  if ( strcmp( flag_meanings, "") != 0) {
3510 
3511  try {
3512  ncVar.putAtt("flag_meanings", flag_meanings);
3513  }
3514  catch ( NcException& e) {
3515  e.what();
3516  cerr << "Failure creating 'flag_meanings' attribute: "
3517  << flag_meanings << endl;
3518  exit(1);
3519  }
3520  }
3521 
3522  /* Add a "reference" attribute if specified*/
3523  if ( strcmp( reference, "") != 0) {
3524 
3525  try {
3526  ncVar.putAtt("reference", reference);
3527  }
3528  catch ( NcException& e) {
3529  e.what();
3530  cerr << "Failure creating 'reference' attribute: "
3531  << reference << endl;
3532  exit(1);
3533  }
3534  }
3535 
3536  /* Add "valid_min/max" attributes if specified */
3537  if (low < high) {
3538  switch(nt) { /* Use the appropriate number type */
3539  case NC_BYTE:
3540  {
3541  uint8_t vr[2];
3542  vr[0] = (uint8_t)low;
3543  vr[1] = (uint8_t)high;
3544 
3545  try {
3546  ncVar.putAtt("valid_min", NC_BYTE, 1, &vr[0]);
3547  }
3548  catch ( NcException& e) {
3549  e.what();
3550  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3551  exit(1);
3552  }
3553 
3554  try {
3555  ncVar.putAtt("valid_max", NC_BYTE, 1, &vr[1]);
3556  }
3557  catch ( NcException& e) {
3558  e.what();
3559  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3560  exit(1);
3561  }
3562  }
3563  break;
3564  case NC_UBYTE:
3565  {
3566  uint8_t vr[2];
3567  vr[0] = (uint8_t)low;
3568  vr[1] = (uint8_t)high;
3569 
3570  try {
3571  ncVar.putAtt("valid_min", NC_UBYTE, 1, &vr[0]);
3572  }
3573  catch ( NcException& e) {
3574  e.what();
3575  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3576  exit(1);
3577  }
3578 
3579  try {
3580  ncVar.putAtt("valid_max", NC_UBYTE, 1, &vr[1]);
3581  }
3582  catch ( NcException& e) {
3583  e.what();
3584  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3585  exit(1);
3586  }
3587  }
3588  break;
3589  case NC_SHORT:
3590  {
3591  int16_t vr[2];
3592  vr[0] = (int16_t)low;
3593  vr[1] = (int16_t)high;
3594 
3595  try {
3596  ncVar.putAtt("valid_min", NC_SHORT, 1, &vr[0]);
3597  }
3598  catch ( NcException& e) {
3599  e.what();
3600  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3601  exit(1);
3602  }
3603 
3604  try {
3605  ncVar.putAtt("valid_max", NC_SHORT, 1, &vr[1]);
3606  }
3607  catch ( NcException& e) {
3608  e.what();
3609  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3610  exit(1);
3611  }
3612  }
3613  break;
3614  case NC_USHORT:
3615  {
3616  uint16_t vr[2];
3617  vr[0] = (uint16_t)low;
3618  vr[1] = (uint16_t)high;
3619 
3620  try {
3621  ncVar.putAtt("valid_min", NC_USHORT, 1, &vr[0]);
3622  }
3623  catch ( NcException& e) {
3624  e.what();
3625  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3626  exit(1);
3627  }
3628 
3629  try {
3630  ncVar.putAtt("valid_max", NC_USHORT, 1, &vr[1]);
3631  }
3632  catch ( NcException& e) {
3633  e.what();
3634  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3635  exit(1);
3636  }
3637  }
3638  break;
3639  case NC_INT:
3640  {
3641  int32_t vr[2];
3642  vr[0] = (int32_t)low;
3643  vr[1] = (int32_t)high;
3644 
3645  try {
3646  ncVar.putAtt("valid_min", NC_INT, 1, &vr[0]);
3647  }
3648  catch ( NcException& e) {
3649  e.what();
3650  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3651  exit(1);
3652  }
3653 
3654  try {
3655  ncVar.putAtt("valid_max", NC_INT, 1, &vr[1]);
3656  }
3657  catch ( NcException& e) {
3658  e.what();
3659  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3660  exit(1);
3661  }
3662 
3663  }
3664  break;
3665  case NC_UINT:
3666  {
3667  uint32_t vr[2];
3668  vr[0] = (uint32_t)low;
3669  vr[1] = (uint32_t)high;
3670 
3671  try {
3672  ncVar.putAtt("valid_min", NC_UINT, 1, &vr[0]);
3673  }
3674  catch ( NcException& e) {
3675  e.what();
3676  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3677  exit(1);
3678  }
3679 
3680  try {
3681  ncVar.putAtt("valid_max", NC_UINT, 1, &vr[1]);
3682  }
3683  catch ( NcException& e) {
3684  e.what();
3685  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3686  exit(1);
3687  }
3688 
3689  }
3690  break;
3691  case NC_FLOAT:
3692  {
3693  float vr[2];
3694  vr[0] = (float)low;
3695  vr[1] = (float)high;
3696 
3697  try {
3698  ncVar.putAtt("valid_min", NC_FLOAT, 1, &vr[0]);
3699  }
3700  catch ( NcException& e) {
3701  e.what();
3702  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3703  exit(1);
3704  }
3705 
3706  try {
3707  ncVar.putAtt("valid_max", NC_FLOAT, 1, &vr[1]);
3708  }
3709  catch ( NcException& e) {
3710  e.what();
3711  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3712  exit(1);
3713  }
3714  }
3715  break;
3716  case NC_DOUBLE:
3717  {
3718  double vr[2];
3719  vr[0] = low;
3720  vr[1] = high;
3721 
3722  try {
3723  ncVar.putAtt("valid_min", NC_DOUBLE, 1, &vr[0]);
3724  }
3725  catch ( NcException& e) {
3726  e.what();
3727  cerr << "Failure creating 'valid_min' attribute: " << vr[0] << endl;
3728  exit(1);
3729  }
3730 
3731  try {
3732  ncVar.putAtt("valid_max", NC_DOUBLE, 1, &vr[1]);
3733  }
3734  catch ( NcException& e) {
3735  e.what();
3736  cerr << "Failure creating 'valid_max' attribute: " << vr[1] << endl;
3737  exit(1);
3738  }
3739  }
3740  break;
3741  default:
3742  fprintf(stderr,"-E- %s line %d: ",__FILE__,__LINE__);
3743  fprintf(stderr,"Got unsupported number type (%d) ",nt);
3744  fprintf(stderr,"while trying to create NCDF variable, \"%s\", ",sname);
3745  return(1);
3746  }
3747  }
3748 
3749  /* Add a "units" attribute if one is specified */
3750  if(units != NULL && *units != 0) {
3751 
3752  try {
3753  ncVar.putAtt("units", units);
3754  }
3755  catch ( NcException& e) {
3756  e.what();
3757  cerr << "Failure creating 'units' attribute: " << units << endl;
3758  exit(1);
3759  }
3760  }
3761 
3762  /* Add a "standard_name" attribute if one is specified */
3763  if(standard_name != NULL && *standard_name != 0) {
3764  try {
3765  ncVar.putAtt("standard_name", units);
3766  }
3767  catch ( NcException& e) {
3768  e.what();
3769  cerr << "Failure creating 'standard_name' attribute: "
3770  << standard_name << endl;
3771  exit(1);
3772  }
3773  }
3774 
3775  return 0;
3776 }
3777 
3778 
3779 
3780 int eight20( uint8_t *inbytes, uint32_t *outsamples) {
3781 
3782  // This routine takes groups of 9 20-bit samples that are packed into
3783  // 23 bytes and unpacks them into a 4-byte integer array
3784 
3785  for (size_t i=0; i<5; i++) {
3786  outsamples[i*2] = 4096*inbytes[i*5] + 16*inbytes[i*5+1] + inbytes[i*5+2]/16;
3787  }
3788 
3789  for (size_t i=0; i<4; i++) {
3790  outsamples[i*2+1] = 65536*(inbytes[i*5+2] % 16) + 256*inbytes[i*5+3] +
3791  inbytes[i*5+4];
3792  }
3793 
3794  return 0;
3795 }
3796 
3797 
unpack_ccd_packet
int unpack_ccd_packet(uint8_t *packet, uint16_t btaps[16], uint16_t rtaps[16], uint16_t &ccdid, uint32_t &line, uint16_t &dtype, uint16_t &iagg, uint16_t jagg[16], uint16_t &nbands, uint16_t **ccddata, uint16_t ossdata[16])
Definition: l1agen_oci.cpp:1206
q
data_t q
Definition: decode_rs.h:74
check_load_oci_data
int check_load_oci_data(short dtype, uint16_t ncpt, uint16_t nspt0, uint16_t ndcs, uint16_t ndss, uint16_t nbbs, uint16_t nrbs, uint16_t nswb, int16_t *cindex, int16_t *sindex, int16_t *cdindex, int16_t *sdindex, uint16_t **bbands, uint16_t **rbands, uint32_t **sbands, int16_t *blines, int16_t *rlines, int16_t *slines, uint16_t **bsci, uint16_t **rsci, uint32_t **ssci, uint16_t **bdark, uint16_t **rdark, uint32_t **sdark, uint8_t &linerr, int &icheck)
Definition: l1agen_oci.cpp:1310
r
int r
Definition: decode_rs.h:73
t
data_t t[NROOTS+1]
Definition: decode_rs.h:77
get_swir_mode
int get_swir_mode(uint8_t(*pbuffer)[PKTSIZE], uint32_t npkts, uint16_t &smode)
Definition: common.cpp:370
SWIR_LOFF_ETU
const uint16_t SWIR_LOFF_ETU[9]
Definition: l1agen_oci.h:12
itab::lines
short lines
Definition: common.h:12
j
int j
Definition: decode_rs.h:73
status
int status
Definition: l1_czcs_hdf.c:32
yds2unix
double yds2unix(int16_t year, int16_t day, double secs)
Definition: yds2unix.c:7
l1aFile::write_navigation
int write_navigation(std::string hktlist, time_struct &starttime, time_struct &endtime)
Definition: l1agen_oci.cpp:3046
ReverseFloat
float ReverseFloat(const float inFloat)
Definition: l1agen_oci.h:16
l1aFile::createl1
int createl1(char *l1_filename, uint32_t nSC, uint32_t imgWidth, uint32_t imgHeight, uint32_t fndWidth, uint32_t fndHeight)
Definition: l1agen_hawkeye.cpp:1207
NULL
#define NULL
Definition: decode_rs.h:63
spin
void spin(double st, double *pos1, double *pos2)
check_sum
uint8_t check_sum(int32_t nc, uint8_t *dat, uint8_t *chk)
Definition: l1agen_oci.cpp:1543
viirs_ler_luts::sb
real, dimension(260) sb
Definition: viirs_ler_luts.f95:12
color_dtdb.v
v
Definition: color_dtdb.py:430
itab::dtype
short dtype
Definition: common.h:10
main
int main(int argc, char *argv[])
Definition: l1agen_oci.cpp:142
TLMSIZE
#define TLMSIZE
Definition: l1agen_oci.h:10
hico.proc_hico.valid_min
valid_min
Definition: proc_hico.py:246
pos
float32 * pos
Definition: l1_czcs_hdf.c:35
msec
int32 * msec
Definition: l1_czcs_hdf.c:31
SWAP_4
#define SWAP_4(x)
Definition: common.h:18
tm
float tm[MODELMAX]
Definition: atrem_corl1v2.h:248
create_images.time_start
time_start
Definition: create_images.py:87
jdate
int jdate(int32_t julian, int32_t *year, int32_t *doy)
Definition: jdate.c:5
jday
int32_t jday(int16_t i, int16_t j, int16_t k)
Definition: jday.c:4
hico.proc_hico.valid_max
valid_max
Definition: proc_hico.py:247
string
@ string
Definition: GEO_read_param_file.c:64
seadasutils.ProcUtils.remove
def remove(file_to_delete)
Definition: ProcUtils.py:319
file
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 file
Definition: HISTORY.txt:413
time_struct::iday
int32_t iday
Definition: l0info_oci.h:7
seadasutils.ProcUtils.mtime
def mtime(the_file)
Definition: ProcUtils.py:344
l1aFile::write_oci_cal_data
int write_oci_cal_data(uint32_t isc, uint16_t nbbs, uint16_t nrbs, uint16_t nswb, uint16_t ndcs, uint16_t ndss, uint16_t *dark_b, uint16_t *dark_r, uint32_t *dark_s, int8_t *sdfrms)
Definition: l1agen_oci.cpp:2005
eight20
int eight20(uint8_t *inbytes, uint32_t *outsamples)
Definition: l1agen_oci.cpp:3780
yd2md
void yd2md(int16_t year, int16_t doy, int16_t *month, int16_t *dom)
Definition: yd2md.c:6
tempOut
ofstream tempOut
Definition: l1agen_oci.cpp:140
mside
int32_t mside
Definition: l1_hmodis_hdf.c:752
tmp
data_t tmp
Definition: decode_rs.h:74
leapseconds_since_1993
int leapseconds_since_1993(double tai93time)
Definition: leapsecond.c:58
l1agen_oci.h
l1aFile::write_oci_tlm_data
int write_oci_tlm_data(itab *itable, uint32_t ntlm, uint8_t(*tlmdata)[TLMSIZE], int32_t *spinID, uint16_t &cdsmode, uint32_t isc)
Definition: l1agen_oci.cpp:2345
get_anc_packet_time
int get_anc_packet_time(uint8_t *apacket, int32_t &iyear, int32_t &iday, double &stime)
Definition: common.cpp:97
ANCSIZE
#define ANCSIZE
Definition: l0info_oci.h:3
l1aFile::parseDims
int parseDims(string dimString, int *numDims, int *varDims)
Definition: l1agen_hawkeye.cpp:1510
read_packet
int read_packet(FILE *infile, uint8_t packet[], int *len, long int *endfile)
Definition: read_packet.c:18
expandEnvVar
int expandEnvVar(string *sValue)
Definition: hawkeyeUtil.h:41
jd
Definition: jd.py:1
tt
subroutine tt(NMAX, NCHECK)
Definition: ampld.lp.f:1852
ioff
int32 ioff
Definition: get_dem_height.c:17
time_struct::iyear
int32_t iyear
Definition: l0info_oci.h:6
parseDims
int parseDims(int ncid, int ndims, string dimString, int *numDims, int *dimid, int *varDims)
Definition: hawkeyeUtil.cpp:339
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
fix_mac_rpath.line
line
Definition: fix_mac_rpath.py:55
ih
int32_t ih
Definition: atrem_corl1.h:161
spacetrack_tles.doy
dictionary doy
Definition: spacetrack_tles.py:171
SWAP_2
#define SWAP_2(x)
Definition: l0gen_hawkeye.cpp:54
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
anc_compare
int anc_compare(uint8_t *apacket0, uint8_t *apacket)
Definition: common.cpp:285
basename
#define basename(s)
Definition: l0chunk_modis.c:29
hico.proc_hico.units
units
Definition: proc_hico.py:245
str
const char * str
Definition: l1c_msi.cpp:35
fix_mac_rpath.lines
lines
Definition: fix_mac_rpath.py:34
unpack_oci_sci
int unpack_oci_sci(uint32_t npkts, int32_t spin, uint16_t ncps, uint16_t nsps, uint16_t msps, uint16_t &nbands, uint16_t btaps[16], uint16_t rtaps[16], uint8_t(*pbuffer)[PKTSIZE], uint16_t **bbands, uint16_t **rbands, uint32_t **sbands, int16_t *blines, int16_t *rlines, int16_t *slines, uint16_t &btype, uint16_t bagg[16], uint16_t &rtype, uint16_t ragg[16], int8_t *sfrm, int &iret)
Definition: l1agen_oci.cpp:1032
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
dtype
dtype
Definition: DDataset.hpp:31
read_oci_scan_packets
int read_oci_scan_packets(fstream *tfileStream, uint8_t *apacket, uint8_t(*pbuffer)[PKTSIZE], uint32_t &npkts, int32_t &spnum, int32_t &ancind, vector< int32_t > &tlmind, uint8_t &seqerr, int32_t &endfile)
Definition: common.cpp:132
itab
Definition: common.h:9
unpack_swir_packet
int unpack_swir_packet(uint8_t *packet, int16_t *slines, uint8_t *swirfrm, uint32_t *swirdata)
Definition: l1agen_oci.cpp:1275
time_struct
Definition: l0info_oci.h:5
time_struct::sec
double sec
Definition: l0info_oci.h:8
spix
int32 spix
Definition: l1_czcs_hdf.c:21
now
u5 which has been done in the LOCALGRANULEID metadata should have an extension NRT It is requested to identify the NRT production Changes from v6 which may affect scientific the sector rotation may actually occur during one of the scans earlier than the one where it is first reported As a the b1 values are about the LOCALGRANULEID metadata should have an extension NRT It is requested to identify the NRT to fill pixels affected by dead subframes with a special value Output the metadata of noisy and dead subframe Dead Subframe EV and Detector Quality Flag2 Removed the function call of Fill_Dead_Detector_SI to stop interpolating SI values for dead but also for all downstream products for science test only Changes from v5 which will affect scientific to conform to MODIS requirements Removed the Mixed option from the ScanType in the code because the L1A Scan Type is never Mixed Changed for ANSI C compliance and comments to better document the fact that when the HDF_EOS metadata is stricly the and products are off by and in the track respectively Corrected some misspelling of RCS swir_oob_sending_detector to the Reflective LUTs to enable the SWIR OOB correction detector so that if any of the sending detectors becomes noisy or non near by good detectors from the same sending band can be specified as the substitute in the new look up table Code change for adding an additional dimension of mirror side to the Band_21_b1 LUT to separate the coefficient of the two mirror sides for just like other thermal emissive so that the L1B code can calibrate Band scan to scan with mirror side dependency which leads better calibration result Changes which do not affect scientific when the EV data are not provided in this Crosstalk Correction will not be performed to the Band calibration data Changes which do not affect scientific and BB_500m in L1A Logic was added to turn off the or to spatial aggregation processes and the EV_250m_Aggr1km_RefSB and EV_500m_Aggr1km_RefSB fields were set to fill values when SDSs EV_250m and EV_500m are absent in L1A file Logic was added to skip the processing and turn off the output of the L1B QKM and HKM EV data when EV_250m and EV_500m are absent from L1A In this the new process avoids accessing and reading the and L1A EV skips and writing to the L1B and EV omits reading and subsampling SDSs from geolocation file and writing them to the L1B and omits writing metadata to L1B and EV and skips closing the L1A and L1B EV and SDSs Logic was added to turn off the L1B OBC output when the high resolution OBC SDSs are absent from L1A This is accomplished by skipping the openning the writing of metadata and the closing of the L1B OBC hdf which is Bit in the scan by scan bit QA has been changed Until now
Definition: HISTORY.txt:361
unix2isodate
char * unix2isodate(double dtime, char zone)
Definition: unix2isodate.c:10
l1aFile::write_oci_science_data
int write_oci_science_data(uint32_t isc, uint16_t nbbs, uint16_t nrbs, uint16_t nswb, uint16_t ncps, uint16_t nsps, uint16_t **bsci, uint16_t **rsci, uint32_t **ssci, int8_t *sfrms)
Definition: l1agen_oci.cpp:1932
l1aFile::write_oci_global_metadata
int write_oci_global_metadata(time_struct &starttime, time_struct &endtime, std::string l1a_name, std::string sdir, std::string edir, uint32_t isc, short dtype, uint16_t smode, uint16_t cdsmode, std::ofstream &fout)
Definition: l1agen_oci.cpp:3253
itab::iagg
short iagg
Definition: common.h:11
ymds2unix
double ymds2unix(short year, short month, short day, double secs)
leap
logical function leap(YEAR)
Definition: leap.f:10
c
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
idy
int32_t idy
Definition: atrem_corl1.h:161
standard_name
These two strings are used for the product XML output If product_id is not set then prefix is used If the last char of the name_prefix is _ then it is removed If algorithm_id is not set then name_suffix is used If the first char is _ then it is removed l2prod standard_name[0]
Definition: HOWTO_Add_a_product.txt:109
i
int i
Definition: decode_rs.h:71
PKTSIZE
#define PKTSIZE
Definition: common.h:7
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
createNCDF
int createNCDF(NcGroup &ncGrp, const char *sname, const char *lname, const char *standard_name, const char *units, void *fill_value, const char *flag_values, const char *flag_meanings, const char *reference, double low, double high, int nt, vector< NcDim > &varVec)
Definition: l1agen_oci.cpp:3367
ccsds_sec_to_yds
int ccsds_sec_to_yds(uint8_t *cctime, int32_t *iyear, int32_t *iday, double *sec)
Definition: l1agen_oci.cpp:20
VERSION
#define VERSION
Definition: l1agen_oci.cpp:46
get_band_dims
int get_band_dims(uint8_t *apacket, uint16_t &ncp, uint16_t &nbb, uint16_t &nrb, uint16_t &nsp, uint16_t &ndc, uint16_t &nds, uint16_t *btaps, uint16_t *rtaps, itab *itable)
Definition: common.cpp:17
k
int k
Definition: decode_rs.h:73
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
l1aFile::write_oci_scan_metadata
int write_oci_scan_metadata(uint32_t isc, uint8_t *ancdata, uint8_t *seqerr, uint8_t *linerr, int32_t *spinID)
Definition: l1agen_oci.cpp:2069
f32
float32 f32
Definition: l2bin.cpp:104
nb
int32_t nb
Definition: atrem_corl1.h:132
l1aFile::write_oci_ancil_data
int write_oci_ancil_data(uint32_t isc, uint8_t *ancdata)
Definition: l1agen_oci.cpp:2172
make_oci_line_index
int make_oci_line_index(itab *itable, int16_t *cindex, int16_t *sindex, int16_t *cdindex, int16_t *sdindex, int16_t *swir_loff)
Definition: l1agen_oci.cpp:944
count
int count
Definition: decode_rs.h:79


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster