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

Mission Quality Monitor

HELP

How to Use Earthdata Search

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

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster