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
common.cpp
Go to the documentation of this file.
1 #include <algorithm>
2 #include <iterator>
3 
4 #include "common.h"
5 
6 using namespace std;
7 
8 #define BAGGBASE 2863311530
9 #define RAGGBASE 2774100650
10 
11 // Modification history:
12 // Programmer Organization Date Ver Description of change
13 // ---------- ------------ ---- --- ---------------------
14 // Joel Gales FutureTech 08/02/19 0.10 Original development
15 // Liang Hong SAIC 04/14/22 0.20 l1agen_oci v1.03.00
16 // Liang Hong SAIC 04/24/22 0.21 l1agen_oci v1.03.02
17 // Liang Hong SAIC 05/05/22 0.30 l1agen_oci v1.04.00
18 // Liang Hong SAIC 11/18/22 0.40 l1agen_oci v1.06.00
19 // Jakob Lindo SSAI 03/16/24 0.50 Added aggregation indication
20 
21 int get_band_dims(uint8_t *apacket, uint16_t &ncp, uint16_t &nbb, uint16_t &nrb, uint16_t &nsp, uint16_t &ndc,
22  uint16_t &nds, uint16_t *btaps, uint16_t *rtaps, itab *itable) {
23  // Extract spatial aggregation table and compute numbers of pixels
24  short ioff = 36;
25  short nagg[4] = {1, 2, 4, 8};
26  ncp = 0;
27  nsp = 0;
28  ndc = 0;
29  nds = 0;
30 
31  for (size_t i = 0; i < 10; i++) {
32  itable[i].dtype = apacket[ioff + 3] % 16;
33  itable[i].iagg = apacket[ioff + 2] % 4;
34  if (itable[i].dtype == 5)
35  itable[i].iagg = 0;
36  itable[i].lines = apacket[ioff] * 256 + apacket[ioff + 1];
37  ioff += 4;
38 
39  if (itable[i].dtype > 0 && itable[i].dtype <= 14 &&
40  itable[i].dtype != 10) { // Changed dtype<=12 to 14, LH, 3/30/2022
41  if (itable[i].dtype == 2) {
42  ndc += itable[i].lines / nagg[itable[i].iagg];
43  nds += itable[i].lines / 8;
44  } // else { Changed to include dark pixels, Liang Hong, 5/12/2020
45  ncp += itable[i].lines / nagg[itable[i].iagg];
46  nsp += itable[i].lines / 8;
47  // }
48  }
49  }
50 
51  // to ensure that the science and dark count arrays will be created with at least one pixel
52  if (ncp == 0) {
53  ncp = 1;
54  nsp = 1;
55  }
56 
57  if (ndc == 0) {
58  ndc = 1;
59  nds = 1;
60  }
61 
62  ioff += 4;
63 
64  // Extract spectral aggregation and compute numbers of bands
65  // Tap enable flags
66  short btap = apacket[ioff + 2] * 256 + apacket[ioff + 3];
67  short rtap = apacket[ioff] * 256 + apacket[ioff + 1];
68 
69  // Tap aggregation factors
70  uint32_t bagg;
71  uint32_t ragg;
72  memcpy(&bagg, &apacket[ioff + 8], sizeof(uint32_t));
73  memcpy(&ragg, &apacket[ioff + 4], sizeof(uint32_t));
74  bagg = SWAP_4(bagg);
75  ragg = SWAP_4(ragg);
76 
77  if (itable[1].dtype == 1 && (bagg != BAGGBASE || ragg != RAGGBASE))
78  itable[1].dtype = 13;
79 
80  // Compute number of bands for enabled taps
81  // Bands are in reverse spectral order
82  nbb = 0;
83  nrb = 0;
84  uint16_t ken = 1;
85  uint32_t kag = 3;
86  uint32_t lag = 1;
87 
88  // for (size_t i=0; i<16; i++) {
89  for (int i = 15; i >= 0; i--) {
90  btaps[i] = (btap & ken) / ken;
91  if (btaps[i])
92  nbb += 32 / nagg[(bagg & kag) / lag];
93  rtaps[i] = (rtap & ken) / ken;
94  if (rtaps[i])
95  nrb += 32 / nagg[(ragg & kag) / lag];
96  ken *= 2;
97  kag *= 4;
98  lag *= 4;
99  }
100 
101  return 0;
102 }
103 
104 int get_anc_packet_time(uint8_t *apacket, int32_t &iyear, int32_t &iday, double &stime) {
105  // Unpack and convert the CCSDS segmented time code
106  // from the OCI ancillary packet
107 
108  // Get day count since Jan. 1, 1958 (Julian day 2436205)
109  // sec58 = seconds since 0 UT Jan. 1, 1958
110  short int toff = 28;
111  uint32_t sec58;
112  memcpy(&sec58, &apacket[toff], sizeof(uint32_t));
113  sec58 = SWAP_4(sec58);
114 
115  double dbl58 = (double)sec58;
116  int leap = leapseconds_since_1993(dbl58) + 27;
117  sec58 -= leap;
118 
119  // Convert to year and day
120  iday = sec58 / 86400;
121  int32_t jd = iday + 2436205; // Jan. 1, 1958 is Julian day 2436205
122  jdate(jd, &iyear, &iday);
123 
124  // Get microseconds
125  uint32_t usec;
126  memcpy(&usec, &apacket[toff + 4], sizeof(uint32_t));
127  usec = SWAP_4(usec);
128  usec = usec / 4096; // 20 MSBs used, last 12 are spares
129 
130  uint32_t isec = sec58 % 86400;
131  stime = isec + ((double)usec) * 1e-6;
132 
133  return 0;
134 }
135 
136 int read_oci_scan_packets(L0Stream *tfileStream, uint8_t *apacket, uint8_t (*pbuffer)[PKTSIZE],
137  uint32_t &npkts, int32_t &spnum, int32_t &ancind, vector<int32_t> &tlmind,
138  uint8_t &seqerr, int32_t &endfile, bool isSPW) {
139  // int pos = tfileStream->tellg();
140  if (endfile)
141  return 0; // LH, 4/24/2022, v1.03.02
142 
143  uint32_t apidmin = 636; // ver
144  uint32_t apidmax = 745; // ver 1.00.00, LH, 1/7/2022
145  uint32_t ui32;
146  uint32_t len;
147  uint32_t maxpkts = 30000; // ver 0.20, LH, 4/14/2022
148 
149  // Ancillary packet APID
150  static uint32_t apida = 636;
151 
152  // APIDs with time fields and spin numbers
153  static uint32_t apidt[8] = {
154  711, 712, 713, 715,
155  716, 717, 721, 723}; // telemetry APIDs with time fields and spin numbers, ver 0.30, LH, 5/5/2022
156 
157  // APIDs with spin numbers but no time fields
158  static uint32_t apidn[2] = {700, 720};
159 
160  // Get OCI spin number from first packet of next scan
161  uint32_t apid = (apacket[0] % 8) * 256 + apacket[1];
162 
163  npkts = 0;
164  tlmind.clear();
165  ancind = -1;
166  // int spne = 0;
167 
168  if (apid == apida) {
169  // Ancillary packet
170  memcpy(&ui32, &apacket[24], 4);
171  spnum = SWAP_4(ui32);
172  } else if (apid == apidn[0] || apid == apidn[1]) {
173  // Science Packet without time field, ignore SWIR packets for now
174  memcpy(&ui32, &apacket[6], 4);
175  spnum = SWAP_4(ui32);
176  // } else if (apid == apidt[0] || apid == apidt[1] || apid == apidt[2] ||apid == apidt[3]) {
177  } else if (std::find(std::begin(apidt), std::end(apidt), apid) != std::end(apidt)) {
178  // Packet with time field
179  memcpy(&ui32, &apacket[12], 4);
180  spnum = SWAP_4(ui32);
181  }
182 
183  // uint32_t spn = spnum;
184  // uint32_t spnt = spnum;
185  int32_t spn = spnum;
186  int32_t spnt = spnum;
187  uint8_t *packet = apacket;
188 
189  len = apacket[4] * 256 + apacket[5] + 7;
190 
191  int prev_pseq = 0;
192  seqerr = 0;
193 
194  // Read all of the packets with this spin number
195  // and store in the packet buffer
196  while (spn <= spnum && spnt <= (spnum + 1) && (npkts < maxpkts)) {
197  // Check for packet out of order
198  // Load packet into buffer
199  memcpy(&pbuffer[npkts][0], packet, len);
200  apid = (packet[0] % 8) * 256 + packet[1];
201 
202  // Check science packet sequence numbers, LH 11/20/2020
203  if (apid == apidn[0]) {
204  int pseq = (pbuffer[npkts][2] % 64) * 256 + pbuffer[npkts][3];
205  if ((prev_pseq > 0) && ((pseq - prev_pseq) != 1 && (pseq - prev_pseq) != -16383)) {
206  // cout<<"npkts="<<npkts<<"pseq,pseq-prev_pseq="<<pseq<<","<<pseq-prev_pseq<<endl;
207  seqerr = 1;
208  }
209  prev_pseq = pseq;
210  }
211 
212  if (apid == apida)
213  ancind = npkts;
214  // if (pos >= 527223898) {
215  // pos = tfileStream->tellg();
216  // cout << pos << " " << apid << " " << apida << " " << ancind
217  // << spnum << " " << spn << " " << spnt << endl;
218  //}
219 
220  if (apid >= apidmin & apid <= apidmax && apid != apidn[0] && apid != apidn[1]) {
221  tlmind.push_back(npkts);
222  }
223  npkts++;
224  if (npkts == maxpkts) {
225  cout << "Maximum number of packets: " << maxpkts << " exceeded." << endl;
226  // exit (1); // Liang Hong, 08/19/2020
227  }
228  // cout << apid << " " << npkts << endl;
229 
230  // Read next packet
231  read_packet(tfileStream, NULL, len, apid, endfile, isSPW);
232  if (endfile)
233  return 0;
234 
235  if (len > PKTSIZE) {
236  cout << "Packet size > " << PKTSIZE << " (" << len << ")" << endl;
237  uint8_t *dummy_buf = new uint8_t[len];
238  read_packet(tfileStream, dummy_buf, len, apid, endfile, isSPW);
239  delete[] dummy_buf;
240  // exit(1);
241  } else {
242  read_packet(tfileStream, packet, len, apid, endfile, isSPW);
243  apid = (packet[0] % 8) * 256 + packet[1];
244  }
245 
246  while (apid < apidmin || apid > apidmax) {
247  read_packet(tfileStream, NULL, len, apid, endfile, isSPW);
248  if (endfile)
249  return 0;
250 
251  if (len > PKTSIZE) {
252  cout << "Packet size > " << PKTSIZE << " (" << len << ")" << endl;
253  uint8_t *dummy_buf = new uint8_t[len];
254  read_packet(tfileStream, dummy_buf, len, apid, endfile, isSPW);
255  delete[] dummy_buf;
256  // exit(1);
257  } else {
258  read_packet(tfileStream, packet, len, apid, endfile, isSPW);
259  apid = (packet[0] % 8) * 256 + packet[1];
260  }
261  }
262 
263  // Get spin number
264  if (apid == apida) {
265  memcpy(&ui32, &packet[24], 4);
266  spn = SWAP_4(ui32);
267  } else if (apid == apidn[0] || apid == apidn[1]) { // ver 1.00.01
268  memcpy(&ui32, &packet[6], 4);
269  spn = SWAP_4(ui32);
270  // } else if (apid == apidt[0] || apid == apidt[1] || apid == apidt[2] || apid == apidt[3]) {
271  } else if (std::find(std::begin(apidt), std::end(apidt), apid) != std::end(apidt)) {
272  memcpy(&ui32, &packet[12], 4);
273  spnt = SWAP_4(ui32);
274  } // if (apid == apidn[0] || apid == apidn[1]) {
275 
276  } // while (spn <= spnum && spnt <= spnum)
277 
278  // Report spin out of order
279  // if (spn < spnum)
280  // cout << "Next packet spin number: " << spn
281  // << " out of order for spin: " << spnum << endl;
282  // cout << "npkts: " << npkts << endl;
283 
284  return 0;
285 }
286 
287 int anc_compare(uint8_t *apacket0, uint8_t *apacket) {
288  int anc_compare = 1;
289 
290  double stime;
291  int32_t iyear, iday;
292  uint32_t ui32;
293 
294  if (apacket[15] > 0)
295  return anc_compare; // ancillary packet without valid mode table, v0.20
296 
297  // Compare spatial data collection fields
298  int ioff = 36;
299  size_t ilen = 40;
300 
301  uint16_t aggi = 0;
302  for (size_t i = 0; i < ilen; i++)
303  if (apacket[ioff + i] != apacket0[ioff + i])
304  aggi++;
305 
306  if (aggi > 0) {
307  get_anc_packet_time(apacket, iyear, iday, stime);
308  memcpy(&ui32, &apacket[24], 4);
309  int32_t spn = SWAP_4(ui32);
310  uint16_t ih = (uint16_t)floor(stime / 3600);
311  uint16_t mn = (uint16_t)floor((stime - ih * 3600) / 60);
312  uint16_t sec = (uint16_t)floor(stime - ih * 3600 - mn * 60);
313  cout << "Spatial table change at: spin=" << spn << ", " << ih << ":" << mn << ":" << sec << endl;
314  anc_compare = 0;
315  }
316 
317  // Compare spectral data collection fields
318  // int joff = 76;
319  // size_t jlen = 16;
320  int joff = 80;
321  size_t jlen = 12;
322 
323  uint16_t aggj = 0;
324  for (size_t i = 0; i < jlen; i++)
325  if (apacket[joff + i] != apacket0[joff + i])
326  aggj++;
327 
328  if (aggj > 0) {
329  get_anc_packet_time(apacket, iyear, iday, stime);
330  memcpy(&ui32, &apacket[24], 4);
331  int32_t spn = SWAP_4(ui32);
332  uint16_t ih = (uint16_t)floor(stime / 3600);
333  uint16_t mn = (uint16_t)floor((stime - ih * 3600) / 60);
334  uint16_t sec = (uint16_t)floor(stime - ih * 3600 - mn * 60);
335  cout << "Spectral table change at: spin=" << spn << ", " << ih << ":" << mn << ":" << sec << endl;
336  anc_compare = 0;
337  }
338 
339  return anc_compare;
340 }
341 
342 int read_packet(L0Stream *tfileStream, uint8_t *packet, uint32_t &len, uint32_t &apid, int32_t &endfile,
343  bool isSPW) {
344  if (tfileStream->eof()) {
345  endfile = 1;
346  len = 0;
347  cout << "End of packet file" << endl;
348  return 0;
349  }
350 
351  // Read spacewire header if needed
352  if (isSPW) {
353  uint8_t spwhead[2];
354  tfileStream->read((char *)&spwhead, 2);
355  }
356 
357  // Read packet header
358  uint8_t phead[6];
359  if (packet == NULL) {
360  tfileStream->read((char *)&phead, 6);
361 
362  // Get length of packet body and APID
363  len = phead[4] * 256 + phead[5] + 1 + 6;
364  apid = (phead[0] % 8) * 256 + phead[1];
365 
366  if (tfileStream->tellg() == -1)
367  endfile = 1;
368  tfileStream->seekg(-6, ios_base::cur);
369  if (isSPW)
370  tfileStream->seekg(-2, ios_base::cur);
371  return 0;
372  }
373  tfileStream->read((char *)packet, len);
374  // cout << tfileStream->tellg() << endl;
375 
376  return 0;
377 }
378 
379 int get_swir_mode(uint8_t (*pbuffer)[PKTSIZE], uint32_t npkts, uint16_t &smode) {
380  // Look for a SWIR band packet
381 
382  // smode = 0;
383  for (size_t ipkt = 0; ipkt < npkts; ipkt++) {
384  uint32_t apid = (pbuffer[ipkt][0] % 8) * 256 + pbuffer[ipkt][1];
385  // if (npkts == 5) cout << "apid: " << apid << endl;
386  if (apid == 720) {
387  // Get mode from packet
388  uint8_t smeta = pbuffer[ipkt][12];
389  smode = (smeta % 64) / 16;
390  return 0;
391  }
392  }
393 
394  return 0;
395  // cout << "SWIR mode not found" << endl;
396  // exit(1);
397 }
L0Stream::tellg
int tellg()
Obtain the position of the read/write head from the beginning of the current filestream.
Definition: l0stream.cpp:146
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
read_packet
int read_packet(L0Stream *tfileStream, uint8_t *packet, uint32_t &len, uint32_t &apid, int32_t &endfile, bool isSPW)
Definition: common.cpp:342
NULL
#define NULL
Definition: decode_rs.h:63
itab::dtype
short dtype
Definition: common.h:11
L0Stream::seekg
void seekg(std::streamoff num, std::ios_base::seekdir dir)
Move the read/write head a number of bytes in a particular direction, switching streams as necessary.
Definition: l0stream.cpp:115
SWAP_4
#define SWAP_4(x)
Definition: common.h:18
L0Stream
A class to simulate a single stream object over multiple OCI L0 files.
Definition: l0stream.hpp:16
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
common.h
leapseconds_since_1993
int leapseconds_since_1993(double tai93time)
Definition: leapsecond.c:58
L0Stream::read
void read(char *buffer, size_t num)
Read from the current stream, switching streams as necessary.
Definition: l0stream.cpp:83
get_anc_packet_time
int get_anc_packet_time(uint8_t *apacket, int32_t &iyear, int32_t &iday, double &stime)
Definition: common.cpp:104
BAGGBASE
#define BAGGBASE
Definition: common.cpp:8
jd
Definition: jd.py:1
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
ih
int32_t ih
Definition: atrem_corl1.h:161
anc_compare
int anc_compare(uint8_t *apacket0, uint8_t *apacket)
Definition: common.cpp:287
ioff
int32 ioff
Definition: get_dem_height.c:17
dtype
dtype
Definition: DDataset.hpp:31
itab
Definition: common.h:10
L0Stream::eof
bool eof()
Indicate if fstream::eof is true.
Definition: l0stream.cpp:139
RAGGBASE
#define RAGGBASE
Definition: common.cpp:9
itab::iagg
short iagg
Definition: common.h:12
leap
logical function leap(YEAR)
Definition: leap.f:10
check_products.end
end
Definition: check_products.py:189
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
i
int i
Definition: decode_rs.h:71
PKTSIZE
#define PKTSIZE
Definition: common.h:8
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


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