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
l1bgen_oci.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <string.h>
5 #include <unistd.h>
6 
7 #include <sstream>
8 #include <fstream>
9 #include <iomanip>
10 #include <getopt.h>
11 #include <libgen.h>
12 #include <dirent.h>
13 #include <sys/stat.h>
14 
15 #include "nc4utils.h"
16 #include "global_attrs.h"
17 #include "l1bgen_oci.h"
18 
19 #include <clo.h>
20 #include <genutils.h>
21 
22 using namespace std;
23 using namespace netCDF;
24 using namespace netCDF::exceptions;
25 
26 #define VERSION "1.2"
27 
28 // NOTE: we are comparing to the int val of solz (scale=0.01)
29 #define MAX_SOLZ (88 * 100)
30 
31 // Modification history:
32 // Programmer Organization Date Ver Description of change
33 // ---------- ------------ ---- --- ---------------------
34 // Joel Gales SAIC 04/29/20 0.01 Original development
35 // Joel Gales SAIC 01/13/21 0.02 Add support for SWIR
36 // Joel Gales SAIC 08/11/21 0.03 Add support for handling
37 // fill values in science data
38 // Joel Gales SAIC 08/12/21 0.04 Add support for hysteresis
39 // correction
40 // Joel Gales SAIC 09/23/21 0.045 Initialize uninitialized
41 // variables
42 // Joel Gales SAIC 05/29/23 0.060 Implemented F.Patt code changes
43 // (05/08/23) #272
44 // Joel Gales SAIC 06/21/23 0.061 Fix hysteresis bugs
45 // Joel Gales SAIC 06/23/23 0.062 Move get_agg_mat to common.cpp
46 // Split off get_nib_nbb
47 // Gwyn Fireman SAIC 07/10/23 0.063 Read global metadata from json file
48 // Joel Gales SAIC 07/24/23 0.063 Add support for saturation
49 // Joel Gales SAIC 08/09/23 0.065 Add check_scan_times
50 // Joel Gales SAIC 08/25/23 0.700 Add support for reflectance
51 // output
52 // Joel Gales SAIC 08/29/23 0.710 Convert thetap,thetas todeg
53 // Joel Gales SAIC 09/29/23 0.800 Call geolocation as function
54 // rather than run beforehand
55 // Joel Gales SAIC 10/16/23 0.810 Add planarity correction
56 // Joel Gales SAIC 10/19/23 0.820 Fix metadata issues
57 // Joel Gales SAIC 12/06/23 0.823 Add rhot description attribute
58 // Joel Gales SAIC 12/10/23 0.840 Add tilt_home
59 // Joel Gales SAIC 01/02/24 0.860 Fix encoder interpolation bug
60 // Joel Gales SAIC 02/09/24 0.870 Add scan angle & polarization
61 // coefficients
62 
63 ofstream tempOut;
64 
65 int main (int argc, char* argv[])
66 {
67 
68  cout << "l1bgen_oci " << VERSION << " ("
69  << __DATE__ << " " << __TIME__ << ")" << endl;
70 
71  clo_optionList_t *optionList = clo_createList();
72  clo_addOption(optionList, "ifile", CLO_TYPE_IFILE, NULL, "Input L1A file");
73  clo_addOption(optionList, "ofile", CLO_TYPE_OFILE, NULL, "Output L1B file");
74  clo_addOption(optionList, "cal_lut", CLO_TYPE_IFILE, NULL, "CAL LUT file");
75  clo_addOption(optionList, "geo_lut", CLO_TYPE_IFILE, NULL, "GEO LUT file");
76  clo_addOption(optionList, "doi", CLO_TYPE_STRING, NULL,
77  "Digital Object Identifier (DOI) string");
78  clo_addOption(optionList, "pversion", CLO_TYPE_STRING, "Unspecified",
79  "processing version string");
80  clo_addOption(optionList, "demfile", CLO_TYPE_STRING,
81  "$OCDATAROOT/common/gebco_ocssw_v2020.nc",
82  "Digital elevation map file");
83  clo_addOption(optionList, "radiance", CLO_TYPE_BOOL, "false",
84  "Generate radiances");
85  clo_addOption(optionList, "ephfile", CLO_TYPE_STRING, nullptr,
86  "Definitive ephemeris file name");
87 
88  string l1a_filename;
89  string l1b_filename;
90  string cal_lut_filename;
91  string geo_lut_filename;
92  string dem_file;
93  string doi;
94  string pversion;
95  string ephFile;
96 
97  bool radiance;
98 
99  if (argc == 1) {
100  clo_printUsage(optionList);
101  exit(EXIT_FAILURE);
102  }
103  clo_readArgs(optionList, argc, argv);
104 
105  // Grab the OPER LUTs
106  string CALDIR = getenv("OCVARROOT");
107  CALDIR.append("/oci/cal/OPER/");
108  vector<string> lut_files;
109 
110  DIR *caldir;
111  struct dirent *caldirptr;
112  if ((caldir = opendir(CALDIR.c_str())) != NULL) {
113  while ((caldirptr = readdir(caldir)) != NULL) {
114  lut_files.push_back(string(caldirptr->d_name));
115  }
116  closedir(caldir);
117  }
118 
119  if (clo_isSet(optionList, "ifile")) {
120  l1a_filename = clo_getString(optionList, "ifile");
121  } else {
122  clo_printUsage(optionList);
123  exit(EXIT_FAILURE);
124  }
125  if (clo_isSet(optionList, "ofile")) {
126  l1b_filename = clo_getString(optionList, "ofile");
127  } else {
128  clo_printUsage(optionList);
129  exit(EXIT_FAILURE);
130  }
131  struct stat fstat_buffer;
132  if (clo_isSet(optionList, "cal_lut")) {
133  cal_lut_filename = clo_getString(optionList, "cal_lut");
134  } else {
135  for(const string& lut : lut_files) {
136  if (lut.find("PACE_OCI_L1B_LUT") != std::string::npos) {
137  cal_lut_filename = CALDIR;
138  cal_lut_filename.append(lut);
139  break;
140  }
141  }
142  }
143  if (cal_lut_filename.empty() || (stat (cal_lut_filename.c_str(), &fstat_buffer) != 0)) {
144  cout << "Error: input CAL LUT file: " << cal_lut_filename.c_str() << " does not exist\n";
145  exit(EXIT_FAILURE);
146  }
147  if (clo_isSet(optionList, "geo_lut")) {
148  geo_lut_filename = clo_getString(optionList, "geo_lut");
149  } else {
150  for(const string& lut : lut_files) {
151  if (lut.find("PACE_OCI_GEO_LUT") != std::string::npos) {
152  geo_lut_filename = CALDIR;
153  geo_lut_filename.append(lut);
154  break;
155  }
156  }
157  }
158  if (geo_lut_filename.empty() || (stat (geo_lut_filename.c_str(), &fstat_buffer) != 0)) {
159  cout << "Error: input GEO LUT file: " << geo_lut_filename.c_str() << " does not exist\n";
160  exit(EXIT_FAILURE);
161  }
162  radiance = clo_getBool(optionList, "radiance");
163  if (clo_isSet(optionList, "doi")) {
164  doi = clo_getString(optionList, "doi");
165  if (doi == "None"){
166  doi.clear();
167  }
168  }
169 
170  if (clo_isSet(optionList, "ephfile")) {
171  ephFile = clo_getString(optionList, "ephfile");
172  }
173 
174  if (clo_isSet(optionList, "pvesion")) {
175  pversion = clo_getString(optionList, "pversion");
176  }
177 
178  char tmp_filename[FILENAME_MAX];
179  parse_file_name(clo_getString(optionList, "demfile"), tmp_filename);
180  dem_file = tmp_filename;
181  if ((stat (dem_file.c_str(), &fstat_buffer) != 0)) {
182  cout << "Error: DEM file: " << dem_file.c_str() << " does not exist\n";
183  exit(EXIT_FAILURE);
184  }
185  // ********************************* //
186  // *** Read calibration LUT file *** //
187  // ********************************* //
188  nc_set_chunk_cache(CHUNK_CACHE_SIZE, CHUNK_CACHE_NELEMS, CHUNK_CACHE_PREEMPTION);
189  NcFile *calLUTfile = new NcFile( cal_lut_filename, NcFile::read);
190 
191  NcGroup gidCommon, gidBlue, gidRed, gidSWIR;
192  gidCommon = calLUTfile->getGroup( "common");
193  gidBlue = calLUTfile->getGroup( "blue");
194  gidRed = calLUTfile->getGroup( "red");
195  gidSWIR = calLUTfile->getGroup( "SWIR");
196 
197  float bwave[NBWAVE];
198  float rwave[NRWAVE];
199  float swave[NIWAVE];// = {939.716, 1038.317, 1250.375, 1248.550, 1378.169,
200  // 1619.626, 1618.036, 2130.593, 2258.432};
201  float bf0[NBWAVE];
202  float rf0[NRWAVE];
203  float sf0[NIWAVE];// = {81.959, 67.021, 44.350, 44.490, 35.551, 23.438, 23.476,
204  // 9.122, 7.427};
205  float spass[NIWAVE] = {45, 80, 30, 30, 15, 75, 75, 50, 75};
206 
207  NcDim timeDim = calLUTfile->getDim("number_of_times");
208  if(timeDim.isNull()) {
209  cout << "Error: could not read number_of_times from " << cal_lut_filename << "\n";
210  exit(EXIT_FAILURE);
211  }
212  size_t numTimes = timeDim.getSize();
213 
214  double *K2t = (double*)malloc(numTimes * sizeof(double));
215 
216  NcVar var;
217 
218  var = gidCommon.getVar( "blue_wavelength");
219  var.getVar( bwave);
220 
221  var = gidCommon.getVar( "blue_F0");
222  var.getVar( bf0);
223 
224  var = gidCommon.getVar( "red_wavelength");
225  var.getVar( rwave);
226 
227  var = gidCommon.getVar( "red_F0");
228  var.getVar( rf0);
229 
230  var = gidCommon.getVar( "SWIR_wavelength");
231  var.getVar( swave);
232 
233  var = gidCommon.getVar( "SWIR_F0");
234  var.getVar( sf0);
235 
236  var = gidCommon.getVar( "K2t");
237  var.getVar( K2t);
238 
239  string tag;
240 
241  cal_lut_struct blue_lut;
242  cal_lut_struct red_lut;
243  cal_lut_struct swir_lut;
244 
245  uint32_t bbanddim, rbanddim, sbanddim, nldim, poldim;
246 
247  tag.assign("blue");
248  read_oci_cal_lut( calLUTfile, tag, gidBlue, bbanddim, 1, nldim, poldim,
249  blue_lut);
250 
251  NcDim K3Tdim = calLUTfile->getDim("number_of_temperatures");
252 
253  float *K3T = new float[K3Tdim.getSize()];
254  var = gidCommon.getVar( "K3T");
255  var.getVar( K3T);
256 
257  tag.assign("red");
258  read_oci_cal_lut( calLUTfile, tag, gidRed, rbanddim, 1, nldim, poldim,
259  red_lut);
260 
261  tag.assign("SWIR");
262  read_oci_cal_lut( calLUTfile, tag, gidSWIR, sbanddim, 2, nldim, poldim,
263  swir_lut);
264 
265  // Read hysterisis parameters
266  float hysttime[9][4];
267  float hystamp[9][4];
268 
269  var = gidSWIR.getVar( "hyst_time_const");
270  var.getVar( &hysttime[0][0]);
271  var = gidSWIR.getVar( "hyst_amplitude");
272  var.getVar( &hystamp[0][0]);
273 
274  calLUTfile->close();
275 
276  NcGroupAtt att;
277 
278  geo_struct geoLUT;
279  geolocate_oci( l1a_filename, geo_lut_filename, geoLUT, l1b_filename, dem_file,
280  radiance, doi, ephFile, pversion);
281 
282  static l1bFile outfile;
283  outfile.l1bfile = new NcFile( l1b_filename, NcFile::write);
284 
285  outfile.ncGrps[0] = outfile.l1bfile->getGroup( "sensor_band_parameters");
286  outfile.ncGrps[1] = outfile.l1bfile->getGroup( "scan_line_attributes");
287  outfile.ncGrps[2] = outfile.l1bfile->getGroup( "geolocation_data");
288  outfile.ncGrps[3] = outfile.l1bfile->getGroup( "navigation_data");
289  outfile.ncGrps[4] = outfile.l1bfile->getGroup( "observation_data");
290 
291 
292  // Append call sequence to existing history
293  string history = get_history(outfile.l1bfile); // if updating existing file
294  history.append(call_sequence(argc, argv));
295 
296  double distcorr;
297  att = outfile.l1bfile->getAtt("earth_sun_distance_correction");
298  att.getValues(&distcorr);
299 
300  NcDim nscanl1b_dim = outfile.l1bfile->getDim("number_of_scans");
301  uint32_t nscanl1b = nscanl1b_dim.getSize();
302 
303  NcDim ccd_pixels_dim = outfile.l1bfile->getDim("ccd_pixels");
304  uint32_t ccd_pixels = ccd_pixels_dim.getSize();
305 
306  double *evtime = new double[nscanl1b];
307  var = outfile.ncGrps[1].getVar( "time");
308  var.getVar( evtime);
309 
310  short *solz;
311  float *csolz;
312  unsigned char *qualFlag;
313 
314  // var = outfile.ncGrps[4].getVar( "rhot_blue");
315  // if ( !var.isNull()) reflectance=true; else reflectance=false;
316 
317  if (!radiance) {
318  solz = new short[nscanl1b*ccd_pixels];
319  var = outfile.ncGrps[2].getVar( "solar_zenith");
320  var.getVar( solz);
321 
322  csolz = new float[nscanl1b*ccd_pixels];
323  for (size_t i=0; i<nscanl1b*ccd_pixels; i++)
324  csolz[i] = cos(solz[i]*M_PI/180/100);
325 
326  qualFlag = new unsigned char[nscanl1b*ccd_pixels];
327  var = outfile.ncGrps[2].getVar("quality_flag");
328  var.getVar(qualFlag);
329  }
330 
331  // Open and read data from L1Afile
332  NcFile *l1afile = new NcFile( l1a_filename, NcFile::read);
333 
334  NcGroup ncGrps[4];
335 
336  ncGrps[0] = l1afile->getGroup( "scan_line_attributes");
337  ncGrps[1] = l1afile->getGroup( "spatial_spectral_modes");
338  ncGrps[2] = l1afile->getGroup( "engineering_data");
339  ncGrps[3] = l1afile->getGroup( "science_data");
340 
341  // Get date (change this when year and day are added to time field)
342  string tstart, tend;
343  att = l1afile->getAtt("time_coverage_start");
344  att.getValues(tstart);
345  cout << "time_coverage_start: " << tstart << endl;
346 
347  att = l1afile->getAtt("time_coverage_end");
348  att.getValues(tend);
349  cout << "time_coverage_end: " << tend << endl;
350 
351  uint16_t iyr, imn, idom;
352  istringstream iss;
353 
354  iss.str(tstart.substr(0,4));
355  iss >> iyr; iss.clear();
356  iss.str(tstart.substr(5,2));
357  iss >> imn; iss.clear();
358  iss.str(tstart.substr(8,2));
359  iss >> idom;
360  int32_t jd = jday(iyr, imn, idom);
361 
362  int32_t iyr32, idy32;
363  jdate( jd, &iyr32, &idy32);
364 
365  // Get numbers of blue and red bands
366  NcDim blue_dim = l1afile->getDim("blue_bands");
367  uint32_t bbands = blue_dim.getSize();
368  NcDim red_dim = l1afile->getDim("red_bands");
369  uint32_t rbands = red_dim.getSize();
370 
371  // Scan time, spin ID and HAM side
372  NcDim nscan_dim = l1afile->getDim("number_of_scans");
373  uint32_t nscan = nscan_dim.getSize();
374 
375  NcDim mcescan_dim = l1afile->getDim("number_of_mce_scans");
376  uint32_t nmcescan = mcescan_dim.getSize();
377 
378  double *sstime = new double[nscan];
379  var = ncGrps[0].getVar( "scan_start_time");
380  var.getVar( sstime);
381 
382  int32_t *spin = new int32_t[nscan];
383  var = ncGrps[0].getVar( "spin_ID");
384  var.getVar( spin);
385 
386  uint8_t *hside = new uint8_t[nscan];
387  var = ncGrps[0].getVar( "HAM_side");
388  var.getVar( hside);
389 
390  uint32_t nscan_good=0;
391  for (size_t i=0; i<nscan; i++) {
392  if ( spin[i] > 0) {
393  sstime[nscan_good] = sstime[i];
394  hside[nscan_good] = hside[i];
395  nscan_good++;
396  }
397  }
398 
399  // Check for and fill in missing scan times
400  short *tfl = new short[nscan_good]();
401  check_scan_times( nscan_good, sstime, tfl);
402 
403  // ******************************************** //
404  // *** Get spatial and spectral aggregation *** //
405  // ******************************************** //
406  NcDim ntaps_dim = l1afile->getDim("number_of_taps");
407  uint32_t ntaps = ntaps_dim.getSize();
408  NcDim spatzn_dim = l1afile->getDim("spatial_zones");
409  uint32_t spatzn = spatzn_dim.getSize();
410 
411  int16_t *dtype = new int16_t[spatzn];
412  var = ncGrps[1].getVar( "spatial_zone_data_type");
413  var.getVar( dtype);
414 
415  int16_t *lines = new int16_t[spatzn];
416  var = ncGrps[1].getVar( "spatial_zone_lines");
417  var.getVar( lines);
418 
419  int16_t *iagg = new int16_t[spatzn];
420  var = ncGrps[1].getVar( "spatial_aggregation");
421  var.getVar( iagg);
422 
423  int16_t *bagg = new int16_t[ntaps];
424  var = ncGrps[1].getVar( "blue_spectral_mode");
425  var.getVar( bagg);
426 
427  int16_t *ragg = new int16_t[ntaps];
428  var = ncGrps[1].getVar( "red_spectral_mode");
429  var.getVar( ragg);
430 
431 
432  // ********************************************************************* //
433  // *** Get # of EV lines/offset from scan start time to EV mid-time *** //
434  // ********************************************************************* //
435  // This will be done by geolocation when integrated into L1B
436  int32_t ppr_off;
437  double revpsec, secpline;
438 
439  int32_t *mspin = new int32_t[nmcescan];
440  int32_t *ot_10us = new int32_t[nmcescan];
441  uint8_t *enc_count = new uint8_t[nmcescan];
442  int16_t board_id;
443 
444  NcDim nenc_dim = l1afile->getDim("encoder_samples");
445  uint32_t nenc = nenc_dim.getSize();
446 
447  float **hamenc = new float *[nmcescan];
448  hamenc[0] = new float[nenc*nmcescan];
449  for (size_t i=1; i<nmcescan; i++) hamenc[i] = hamenc[i-1] + nenc;
450 
451  float **rtaenc = new float *[nmcescan];
452  rtaenc[0] = new float[nenc*nmcescan];
453  for (size_t i=1; i<nmcescan; i++) rtaenc[i] = rtaenc[i-1] + nenc;
454 
455  int16_t iret;
456  read_mce_tlm( l1afile, geoLUT, ncGrps[2], nmcescan, nenc, ppr_off, revpsec,
457  secpline, board_id, mspin, ot_10us, enc_count,
458  &hamenc[0], rtaenc, iret);
459 
460  float clines[32400], slines[4050];
461  uint16_t pcdim, psdim;
462  // int16_t iret;
463  double ev_toff, deltc[32400], delts[4050];
464  bool dark = false;
465  get_ev( secpline, dtype, lines, iagg, pcdim, psdim, ev_toff, clines, slines,
466  deltc, delts, dark, iret);
467  if ( iret < 0) {
468  cout << "No science collect in file: " << l1a_filename.c_str() << endl;
469  l1afile->close();
470  return 1;
471  }
472 
473  size_t ka;
474  for (size_t i=0; i<spatzn; i++) {
475  if ( dtype[i] != 0 && dtype[i] != 2 && dtype[i] != 10) {
476  ka = i;
477  break;
478  }
479  }
480 
481 
482  // *********************************************************** //
483  // *** Generate matrices for spectral and gain aggregation *** //
484  // *********************************************************** //
485  size_t *ia = new size_t[ntaps];
486  uint32_t ntb[16];
487 
488  // Blue bands
489  uint32_t bib=1, bbb=1;
490 
491  if (get_nib_nbb( ntaps, ia, ntb, bagg, bib, bbb) == 2)
492  cout << "All blue taps disabled" << endl;
493 
494  // Note: bgmat & rgmat not necessarily contiguous
495  float **bamat = new float*[bib];
496  float **bgmat = new float*[512];
497 
498  if (bib != 1)
499  get_agg_mat( ia, bagg, ntb, bib, bbb, bamat, bgmat);
500 
501  if (bib != bbands) {
502  cout << "Number of blue bands in file: " << l1a_filename.c_str() <<
503  " not consistent with spectral aggregation" << endl;
504  l1afile->close();
505  return 1;
506  } else if ( bib < 4) cout << "No blue bands in file: " <<
507  l1a_filename.c_str() << endl;
508 
509 
510  // Red bands
511  uint32_t rib=1, rbb=1;
512 
513  if (get_nib_nbb( ntaps, ia, ntb, ragg, rib, rbb) == 2)
514  cout << "All red taps disabled" << endl;
515 
516  float **ramat = new float*[rib];
517  float **rgmat = new float*[512];
518 
519  if (rib != 1)
520  get_agg_mat( ia, ragg, ntb, rib, rbb, ramat, rgmat);
521 
522  if (rib != rbands) {
523  cout << "Number of red bands in file: " << l1a_filename.c_str() <<
524  " not consistent with spectral aggregation" << endl;
525  l1afile->close();
526  return 1;
527  } else if ( rib < 4) cout << "No red bands in file: " << l1a_filename.c_str()
528  << endl;
529 
530  uint16_t swb = 9;
531 
532 
533  // ********************************* //
534  // *** Get dark collect location *** //
535  // ********************************* //
536  int16_t kd=-1;
537  for (size_t i=0; i<spatzn; i++) {
538  if ( dtype[i] == 2) {
539  kd = (int16_t) i;
540  }
541  }
542  if ( kd == -1) {
543  cout << "No dark collect in file: " << l1a_filename.c_str() << endl;
544  l1afile->close();
545  return 1;
546  }
547 
548  int16_t ldc=0, lds=0;
549  for (size_t i=0; i<(size_t) kd; i++) {
550  if ( dtype[i] != 0 && dtype [1] != 10) {
551  ldc += lines[i] / iagg[i];
552  lds += lines[i] / 8;
553  }
554  }
555  int16_t ndc = lines[kd] / iagg[kd];
556  int16_t nds = lines[kd] / 8;
557 
558 
559  // *********************************************************************** //
560  // *** Generate band gain structs from LUTs, date/time & gain matrices *** //
561  // *********************************************************************** //
562  gains_struct bgains;
563  gains_struct rgains;
564  gains_struct sgains;
565 
566  // Note: bgmat & rgmat not necessarily contiguous
567  // That's why we pass the location of the array of pointers
568  if ( bib >= 4)
569  make_oci_gains( bib, bbanddim, iyr, jd, evtime[0], numTimes, K2t, -1,
570  iagg[ka], bagg, blue_lut, &bgmat[0], bgains);
571 
572  if ( rib >= 4)
573  make_oci_gains( rib, rbanddim, iyr, jd, evtime[0], numTimes, K2t, -1,
574  iagg[ka], ragg, red_lut, &rgmat[0], rgains);
575 
576  float **sgmat = new float*[swb];
577  for (size_t i=0; i<swb; i++) {
578  sgmat[i] = new float[swb];
579  for (size_t j=0; j<swb; j++) {
580  if (i == j) sgmat[i][j] = 1.0; else sgmat[i][j] = 0.0;
581  }
582  }
583  make_oci_gains( swb, swb, iyr, jd, evtime[0], numTimes, K2t, board_id, -1, NULL,
584  swir_lut, &sgmat[0], sgains);
585 
586  // Compute bibf0,b1bf0
587  float *bibf0 = new float[bib];
588  for (size_t i=0; i<bib; i++) {
589  bibf0[i] = 0.0;
590  for (size_t j=0; j<512; j++) {
591  bibf0[i] += bf0[j]*bgmat[j][i];
592  }
593  }
594 
595  float *b1bf0 = new float[bbb];
596  for (size_t i=0; i<bbb; i++) {
597  b1bf0[i] = 0.0;
598  for (size_t j=0; j<bib; j++) {
599  b1bf0[i] += bibf0[j]*bamat[j][i];
600  }
601  }
602 
603  // Compute ribf0,r1bf0
604  float *ribf0 = new float[rib];
605  for (size_t i=0; i<rib; i++) {
606  ribf0[i] = 0.0;
607  for (size_t j=0; j<512; j++) {
608  ribf0[i] += rf0[j]*rgmat[j][i];
609  }
610  }
611 
612  float *r1bf0 = new float[rbb];
613  for (size_t i=0; i<rbb; i++) {
614  r1bf0[i] = 0.0;
615  for (size_t j=0; j<rib; j++) {
616  r1bf0[i] += ribf0[j]*ramat[j][i];
617  }
618  }
619 
620  // Read selected temperature fields and interpolate to scan times
621  uint16_t ntemps = bgains.ldims[1] + sgains.ldims[1];
622 
623  float **caltemps = new float *[nscan_good];
624  caltemps[0] = new float[ntemps*nscan_good];
625  for (size_t i=1; i<nscan_good; i++) caltemps[i] = caltemps[i-1] + ntemps;
626 
627  for (size_t i=0; i<nscan_good; i++)
628  for (size_t j=1; j<ntemps; j++)
629  caltemps[i][j] = 0.0;
630  get_oci_cal_temps( l1afile, ncGrps[2], ntemps, nscan_good, evtime, caltemps);
631  uint16_t nctemps = bgains.ldims[1];
632 
633 
634  // Read dark collects from science data arrays
635  vector<size_t> start, count;
636  start.push_back(0);
637  start.push_back(0);
638  start.push_back(ldc);
639 
640  size_t dims[3];
641  dims[0] = nscan_good; dims[1] = bib; dims[2] = ndc;
642  uint16_t ***bdark = make3dT<uint16_t>(dims);
643  dims[0] = nscan_good; dims[1] = rib; dims[2] = ndc;
644  uint16_t ***rdark = make3dT<uint16_t>(dims);
645 
646  uint16_t bfill;
647  uint16_t rfill;
648 
649  if ( bib > 4) {
650  count.push_back(nscan_good);
651  count.push_back(bib);
652  count.push_back(ndc);
653 
654  var = ncGrps[3].getVar( "sci_blue");
655  var.getVar( start, count, &bdark[0][0][0]);
656 
657  var.getAtt("_FillValue").getValues(&bfill);
658  }
659 
660  if ( rib > 4) {
661  count.clear();
662  count.push_back(nscan_good);
663  count.push_back(rib);
664  count.push_back(ndc);
665 
666  var = ncGrps[3].getVar( "sci_red");
667  var.getVar( start, count, &rdark[0][0][0]);
668 
669  var.getAtt("_FillValue").getValues(&rfill);
670  }
671 
672  dims[0] = nscan_good; dims[1] = swb; dims[2] = nds;
673  uint32_t ***sdark = make3dT<uint32_t>(dims);
674 
675  uint32_t sfill;
676 
677  start.clear();
678  start.push_back(0);
679  start.push_back(0);
680  start.push_back(lds);
681 
682  count.clear();
683  count.push_back(nscan_good);
684  count.push_back(swb);
685  count.push_back(nds);
686 
687  var = ncGrps[3].getVar( "sci_SWIR");
688  var.getVar( start, count, &sdark[0][0][0]);
689 
690  var.getAtt("_FillValue").getValues(&sfill);
691 
692  uint32_t fill32;
693 
694  // number of scans of dark data to average; will make this an input parameter
695  uint16_t ndsc = 1;
696  //number of dark pixels to skip; will make this an input parameter
697  uint16_t nskp = 0;
698 
699  uint32_t bicount[3] = {1,bib,(uint32_t) pcdim};
700  uint32_t ricount[3] = {1,rib,(uint32_t) pcdim};
701  uint32_t sicount[3] = {1,swb,(uint32_t) psdim};
702  uint32_t bcount[3] = {bbb,1,(uint32_t) pcdim};
703  uint32_t rcount[3] = {rbb,1,(uint32_t) pcdim};
704  uint32_t scount[3] = {swb,1,(uint32_t) psdim};
705 
706  // Calibrated data variables
707  float **bdn = new float *[bib];
708  bdn[0] = new float[pcdim*bib];
709  for (size_t i=1; i<bib; i++) bdn[i] = bdn[i-1] + pcdim;
710 
711  float **rdn = new float *[rib];
712  rdn[0] = new float[pcdim*rib];
713  for (size_t i=1; i<rib; i++) rdn[i] = rdn[i-1] + pcdim;
714 
715  float **sdn = new float *[swb];
716  sdn[0] = new float[psdim*swb];
717  for (size_t i=1; i<swb; i++) sdn[i] = sdn[i-1] + psdim;
718 
719  float **bcal = new float *[bib];
720  bcal[0] = new float[pcdim*bib];
721  for (size_t i=1; i<bib; i++) bcal[i] = bcal[i-1] + pcdim;
722 
723  float **rcal = new float *[rib];
724  rcal[0] = new float[pcdim*rib];
725  for (size_t i=1; i<rib; i++) rcal[i] = rcal[i-1] + pcdim;
726 
727  float **scal = new float *[swb];
728  scal[0] = new float[psdim*swb];
729  for (size_t i=1; i<swb; i++) scal[i] = scal[i-1] + psdim;
730 
731  // Saturation arrays
732  uint8_t *bsat = new uint8_t[bib];
733  uint8_t **bqual = new uint8_t *[bbb];
734  bqual[0] = new uint8_t[pcdim*bbb];
735  for (size_t i=1; i<bbb; i++) bqual[i] = bqual[i-1] + pcdim;
736 
737  uint8_t *rsat = new uint8_t[rib];
738  uint8_t **rqual = new uint8_t *[rbb];
739  rqual[0] = new uint8_t[pcdim*rbb];
740  for (size_t i=1; i<rbb; i++) rqual[i] = rqual[i-1] + pcdim;
741 
742  uint8_t *ssat = new uint8_t[swb];
743  uint8_t **squal = new uint8_t *[swb];
744  squal[0] = new uint8_t[psdim*swb];
745  for (size_t i=1; i<swb; i++) squal[i] = squal[i-1] + psdim;
746 
747 
748  double *thetap = new double[pcdim];
749  double *thetas = new double[psdim];
750 
751  float **pview = new float *[pcdim];
752  pview[0] = new float[3*pcdim];
753  for (size_t i=1; i<pcdim; i++) pview[i] = pview[i-1] + 3;
754 
755  float **sview = new float *[psdim];
756  sview[0] = new float[3*psdim];
757  for (size_t i=1; i<psdim; i++) sview[i] = sview[i-1] + 3;
758 
759  uint16_t **bsci = new uint16_t *[bib];
760  bsci[0] = new uint16_t[pcdim*bib];
761  for (size_t i=1; i<bib; i++) bsci[i] = bsci[i-1] + pcdim;
762 
763  uint16_t **rsci = new uint16_t *[rib];
764  rsci[0] = new uint16_t[pcdim*rib];
765  for (size_t i=1; i<rib; i++) rsci[i] = rsci[i-1] + pcdim;
766 
767  uint32_t **ssci = new uint32_t *[swb];
768  ssci[0] = new uint32_t[psdim*swb];
769  for (size_t i=1; i<swb; i++) ssci[i] = ssci[i-1] + psdim;
770 
771  float **bcalb = new float *[bbb];
772  bcalb[0] = new float[pcdim*bbb];
773  for (size_t i=1; i<bbb; i++) bcalb[i] = bcalb[i-1] + pcdim;
774 
775  float **rcalb = new float *[rbb];
776  rcalb[0] = new float[pcdim*rbb];
777  for (size_t i=1; i<rbb; i++) rcalb[i] = rcalb[i-1] + pcdim;
778 
780  // Main loop //
782 
783  // Read, calibrate and write science data
784  for (size_t iscn=0; iscn<nscan_good; iscn++) {
785 
786  if ((iscn % 50) == 0) cout << "Calibrating scan: " << iscn << endl;
787 
788  // Check for valid mirror side
789  if ( hside[iscn] == 0 || hside[iscn] == 1) {
790 
791  // Get scan angle
792  get_oci_vecs( nscan, pcdim, geoLUT.as_planarity, geoLUT.at_planarity,
793  geoLUT.rta_nadir, geoLUT.ham_ct_angles,
794  ev_toff, spin[iscn], hside[iscn],
795  clines, deltc, revpsec, ppr_off, board_id, nmcescan, mspin,
796  enc_count, &hamenc[0], &rtaenc[0], pview, thetap, iret);
797 
798  for (size_t k=0; k<pcdim; k++) thetap[k] *= 180/M_PI;
799 
800  get_oci_vecs( nscan, psdim, geoLUT.as_planarity, geoLUT.at_planarity,
801  geoLUT.rta_nadir, geoLUT.ham_ct_angles,
802  ev_toff, spin[iscn], hside[iscn],
803  slines, delts, revpsec, ppr_off, board_id, nmcescan, mspin,
804  enc_count, &hamenc[0], &rtaenc[0], sview, thetas, iret);
805 
806  for (size_t k=0; k<psdim; k++) thetas[k] *= 180/M_PI;
807 
808  // Write scan angles
809  start.clear();
810  start.push_back(iscn);
811  start.push_back(0);
812 
813  count.clear();
814  count.push_back(1);
815  count.push_back(pcdim);
816 
817  var = outfile.ncGrps[3].getVar( "CCD_scan_angles");
818  var.putVar( start, count, thetap);
819 
820  count[1] = psdim;
821  var = outfile.ncGrps[3].getVar( "SWIR_scan_angles");
822  var.putVar( start, count, thetas);
823 
824  // Blue bands
825  if (bib >= 4) {
826 
827  start.clear();
828  start.push_back(iscn);
829  start.push_back(0);
830  start.push_back(0);
831 
832  count.clear();
833  count.push_back(bicount[0]);
834  count.push_back(bicount[1]);
835  count.push_back(bicount[2]);
836 
837 
838  var = ncGrps[3].getVar( "sci_blue");
839  var.getVar( start, count, bsci[0]);
840 
841  // Compute dark offset, correct data, and apply absolute and
842  // temporal gain and temperature correction
843  float *bdc = new float[bib];
844 
845  fill32 = bfill;
846  int16_t iret;
847  get_oci_dark<uint16_t>( iscn, nscan_good, hside, ndsc, nskp,
848  iagg[ka], iagg[kd], ntaps, bagg, fill32,
849  ndc, bdark, bib, bdc, iret);
850 
851  float *k3 = new float[bib];
852  get_oci_temp_corr( bib, bgains, K3T, caltemps[iscn], nscan_good, k3);
853  for (size_t j=0; j<bib; j++) {
854  for (size_t k=0; k<pcdim; k++) {
855 
856  // Handle fill value
857  if (bsci[j][k] == bfill) {
858  bdn[j][k] = -32767;
859  bcal[j][k] = -32767;
860  continue;
861  }
862 
863  // Need to save dn for linearity correction
864  bdn[j][k] = bsci[j][k] - bdc[j];
865  bcal[j][k] = k3[j] * bgains.K1K2[j][hside[iscn]] * bdn[j][k];
866  }
867  }
868 
869  delete [] k3;
870  delete [] bdc;
871 
872  // Compute and apply RVS and linearity
873  float **k4 = new float *[bib];
874  k4[0] = new float[pcdim*bib];
875  for (size_t i=1; i<bib; i++) k4[i] = k4[i-1] + pcdim;
876  get_oci_rvs_corr( bib, pcdim, hside[iscn], bgains, thetap, k4);
877 
878  float **k5 = new float *[bib];
879  k5[0] = new float[pcdim*bib];
880  for (size_t i=1; i<bib; i++) k5[i] = k5[i-1] + pcdim;
881  get_oci_lin_corr( bib, pcdim, nldim, bgains, bdn, k5);
882 
883  for (size_t j=0; j<bib; j++) {
884  for (size_t k=0; k<pcdim; k++) {
885  if(bcal[j][k] != -32767)
886  bcal[j][k] *= k4[j][k] * k5[j][k];
887  }
888  }
889 
890  delete [] k4[0];
891  delete [] k4;
892  delete [] k5[0];
893  delete [] k5;
894 
895  // Aggregate to L1B bands
896  //bcalb = transpose(bamat#transpose(bcal))
897  for (size_t j=0; j<bbb; j++) {
898  for (size_t k=0; k<pcdim; k++) {
899  float sum = 0.0;
900  for (size_t l=0; l<bib; l++)
901  if (bcal[l][k] != -32767) sum += bamat[l][j]*bcal[l][k];
902  bcalb[j][k] = sum;
903  if (!radiance) {
904  if(solz[iscn*pcdim+k] < MAX_SOLZ) // NOTE: solz is short with scale of 0.01
905  bcalb[j][k] *= M_PI*distcorr/(b1bf0[j]*csolz[iscn*pcdim+k]);
906  else
907  bcalb[j][k] = -32767;
908  }
909  }
910  }
911 
912  // Check for saturation
913  for (size_t k=0; k<pcdim; k++) {
914  for (size_t j=0; j<bib; j++) {
915  bsat[j] = 0;
916  if ( bdn[j][k] >= bgains.sat_thres[j]) bsat[j] = 1;
917  }
918  for (size_t j=0; j<bbb; j++) {
919  float sum = 0.0;
920  bqual[j][k] = 0;
921  for (size_t l=0; l<bib; l++) sum += bamat[l][j]*bsat[l];
922  if ( sum > 0) bqual[j][k] = 1;
923  }
924  }
925 
926  start.clear();
927  start.push_back(0);
928  start.push_back(iscn);
929  start.push_back(0);
930 
931  count.clear();
932  count.push_back(bcount[0]);
933  count.push_back(bcount[1]);
934  count.push_back(bcount[2]);
935 
936  // Output to L1B file
937  if (!radiance)
938  var = outfile.ncGrps[4].getVar( "rhot_blue");
939  else
940  var = outfile.ncGrps[4].getVar( "Lt_blue");
941  var.putVar( start, count, &bcalb[0][0]);
942 
943  var = outfile.ncGrps[4].getVar( "qual_blue");
944  var.putVar( start, count, &bqual[0][0]);
945  } // End blue
946 
947 
948  // Red bands
949  if (rib >= 4) {
950 
951  start.clear();
952  start.push_back(iscn);
953  start.push_back(0);
954  start.push_back(0);
955 
956  count.clear();
957  count.push_back(ricount[0]);
958  count.push_back(ricount[1]);
959  count.push_back(ricount[2]);
960 
961  var = ncGrps[3].getVar( "sci_red");
962  var.getVar( start, count, rsci[0]);
963 
964  // Compute dark offset, correct data, and apply absolute and
965  // temporal gain and temperature correction
966  float *rdc = new float[rib];
967  fill32 = rfill;
968  int16_t iret;
969  get_oci_dark<uint16_t>( iscn, nscan_good, hside, ndsc, nskp,
970  iagg[ka], iagg[kd], ntaps, ragg, fill32,
971  ndc, rdark, rib, rdc, iret);
972 
973  float *k3 = new float[rib];
974  get_oci_temp_corr( rib, rgains, K3T, caltemps[iscn], nscan_good, k3);
975 
976  for (size_t j=0; j<rib; j++) {
977  for (size_t k=0; k<pcdim; k++) {
978 
979  // Handle fill value
980  if (rsci[j][k] == rfill) {
981  rdn[j][k] = -32767;
982  rcal[j][k] = -32767;
983  continue;
984  }
985 
986  // Need to save dn for linearity correction
987  rdn[j][k] = rsci[j][k] - rdc[j];
988  rcal[j][k] = k3[j] * rgains.K1K2[j][hside[iscn]] * rdn[j][k];
989  }
990  }
991 
992  delete [] k3;
993  delete [] rdc;
994 
995  // Compute and apply RVS and linearity
996  float **k4 = new float *[rib];
997  k4[0] = new float[pcdim*rib];
998  for (size_t i=1; i<rib; i++) k4[i] = k4[i-1] + pcdim;
999  get_oci_rvs_corr( rib, pcdim, hside[iscn], rgains, thetap, k4);
1000 
1001  float **k5 = new float *[rib];
1002  k5[0] = new float[pcdim*rib];
1003  for (size_t i=1; i<rib; i++) k5[i] = k5[i-1] + pcdim;
1004 
1005  get_oci_lin_corr( rib, pcdim, nldim, rgains, rdn, k5);
1006 
1007  for (size_t j=0; j<rib; j++) {
1008  for (size_t k=0; k<pcdim; k++) {
1009  if(rcal[j][k] != -32767)
1010  rcal[j][k] *= k4[j][k] * k5[j][k];
1011  }
1012  }
1013 
1014  delete [] k4[0];
1015  delete [] k4;
1016  delete [] k5[0];
1017  delete [] k5;
1018 
1019  // Aggregate to L1B bands
1020  for (size_t j=0; j<rbb; j++) {
1021  for (size_t k=0; k<pcdim; k++) {
1022  float sum = 0.0;
1023  for (size_t l=0; l<rib; l++)
1024  if (rcal[l][k] != -32767) sum += ramat[l][j]*rcal[l][k];
1025  rcalb[j][k] = sum;
1026  if (!radiance) {
1027  if(solz[iscn*pcdim+k] < MAX_SOLZ) // NOTE: solz is short with scale of 0.01
1028  rcalb[j][k] *= M_PI*distcorr/(r1bf0[j]*csolz[iscn*pcdim+k]);
1029  else
1030  rcalb[j][k] = -32767;
1031  }
1032  }
1033  }
1034 
1035  // Check for saturation
1036  for (size_t k=0; k<pcdim; k++) {
1037  for (size_t j=0; j<rib; j++) {
1038  rsat[j] = 0;
1039  if ( rdn[j][k] >= rgains.sat_thres[j]) rsat[j] = 1;
1040  }
1041  for (size_t j=0; j<rbb; j++) {
1042  float sum = 0.0;
1043  rqual[j][k] = 0;
1044  for (size_t l=0; l<rib; l++) sum += ramat[l][j]*rsat[l];
1045  if ( sum > 0) rqual[j][k] = 1;
1046  }
1047  }
1048 
1049  start.clear();
1050  start.push_back(0);
1051  start.push_back(iscn);
1052  start.push_back(0);
1053 
1054  count.clear();
1055  count.push_back(rcount[0]);
1056  count.push_back(rcount[1]);
1057  count.push_back(rcount[2]);
1058 
1059  // Output to L1B file
1060  if (!radiance)
1061  var = outfile.ncGrps[4].getVar( "rhot_red");
1062  else
1063  var = outfile.ncGrps[4].getVar( "Lt_red");
1064  var.putVar( start, count, &rcalb[0][0]);
1065 
1066  var = outfile.ncGrps[4].getVar( "qual_red");
1067  var.putVar( start, count, &rqual[0][0]);
1068 
1069  } // End red
1070 
1071  // SWIR bands
1072  start.clear();
1073  start.push_back(iscn);
1074  start.push_back(0);
1075  start.push_back(0);
1076 
1077  count.clear();
1078  count.push_back(sicount[0]);
1079  count.push_back(sicount[1]);
1080  count.push_back(sicount[2]);
1081 
1082  var = ncGrps[3].getVar( "sci_SWIR");
1083  var.getVar( start, count, ssci[0]);
1084 
1085  // Compute dark offset, correct data, and apply absolute and
1086  // temporal gain and temperature correction
1087  float *sdc = new float[swb];
1088  int16_t sagg = 1;
1089  int16_t iret;
1090 
1091  float *k3 = new float[swb];
1092 
1093  float **k4 = new float *[swb];
1094  k4[0] = new float[psdim*swb];
1095  for (size_t i=1; i<swb; i++) k4[i] = k4[i-1] + psdim;
1096 
1097  float **k5 = new float *[swb];
1098  k5[0] = new float[psdim*swb];
1099  for (size_t i=1; i<swb; i++) k5[i] = k5[i-1] + psdim;
1100 
1101 
1102  // initialize k4, k5 and sdn
1103  // for (int band=0; band<swb; band++) {
1104  // for (int pix=0; pix<psdim; pix++) {
1105  // k4[band][pix] = 0.0;
1106  // k5[band][pix] = 0.0;
1107  // sdn[band][pix] = 0.0;
1108  // }
1109  // }
1110  for (int pix=0; pix<psdim*swb; pix++) {
1111  k4[0][pix] = 0.0;
1112  k5[0][pix] = 0.0;
1113  sdn[0][pix] = 0.0;
1114  }
1115 
1116  get_oci_dark<uint32_t>( iscn, nscan_good, hside, ndsc, nskp,
1117  1, 1, 1, &sagg, sfill, nds, sdark,
1118  swb, sdc, iret);
1119 
1120  if ( iret != -1) {
1121  get_oci_temp_corr( swb, sgains, &K3T[nctemps], &caltemps[iscn][nctemps],
1122  nscan_good, k3);
1123 
1124  // Compute and apply RVS and linearity
1125  get_oci_rvs_corr( swb, psdim, hside[iscn], sgains, thetas, k4);
1126  get_oci_lin_corr( swb, psdim, nldim, sgains, sdn, k5);
1127 
1128 
1129  // SWIR band loop
1130  for (size_t j=0; j<swb; j++) {
1131 
1132  uint32_t goodcnt=0;
1133  int32_t *indx = new int32_t[psdim];
1134  for (size_t k=0; k<psdim; k++) {
1135 
1136  // radiance == true, do not check qualFlag bc it does not retrieve it
1137  if (radiance && ssci[j][k] != sfill) {
1138  indx[goodcnt++] = k;
1139  }
1140  // radiance == false, check quality flag
1141  else if (!radiance && ssci[j][k] != sfill && (int)qualFlag[iscn*pcdim+k] == 0) {
1142  indx[goodcnt++] = k;
1143  }
1144  else {
1145  // Handle fill value
1146  sdn[j][k] = -32767;
1147  scal[j][k] = -32767;
1148  }
1149  }
1150 
1151  // Hysteresis correction
1152  float hc_prev[4]={0,0,0,0};
1153  float hc[4]={0,0,0,0};
1154  float *hyst = new float[psdim];
1155 
1156  for (size_t k=0; k<goodcnt; k++) {
1157  // Need to save dn for linearity correction
1158  sdn[j][indx[k]] = ssci[j][indx[k]] - sdc[j];
1159 
1160  if ( k > 0) {
1161  hyst[k] = 0.0;
1162  for (size_t l=0; l<3; l++) {
1163  // Compute exponential decay constants
1164  float e = exp(-1.0/hysttime[j][l]);
1165 
1166  hc[l] = hc_prev[l]*e + sdn[j][indx[k-1]]*hystamp[j][l];
1167  hyst[k] += hc[l];
1168 
1169  hc_prev[l] = hc[l];
1170  } // l-loop
1171  } else {
1172  hyst[k] = 0.0;
1173  }
1174  } // (reduced) k-loop
1175 
1176  for (size_t k=0; k<goodcnt; k++) {
1177  scal[j][indx[k]] =
1178  k3[j] * sgains.K1K2[j][hside[iscn]] * (sdn[j][indx[k]] - hyst[k]);
1179 
1180  scal[j][indx[k]] *= k5[j][indx[k]]*k4[j][indx[k]];
1181 
1182  if (!radiance) {
1183  if(solz[iscn*pcdim+indx[k]] < MAX_SOLZ) // NOTE: solz is short with scale of 0.01
1184  scal[j][indx[k]] *= M_PI*distcorr/(sf0[j]*csolz[iscn*psdim+indx[k]]);
1185  else
1186  scal[j][indx[k]] = -32767;
1187  }
1188  } // (reduced) k-loop
1189 
1190  delete [] indx;
1191  delete [] hyst;
1192  } // j-loop (band)
1193  } // iret != -1
1194 
1195 
1196  // Check for saturation
1197  for (size_t k=0; k<psdim; k++) {
1198  for (size_t j=0; j<swb; j++) {
1199  ssat[j] = 0;
1200  if ( ssci[j][k] >= sgains.sat_thres[j]) ssat[j] = 1;
1201  }
1202  for (size_t j=0; j<swb; j++) {
1203  squal[j][k] = ssat[j];
1204  }
1205  }
1206 
1207  delete [] k3;
1208  delete [] sdc;
1209 
1210  delete [] k4[0];
1211  delete [] k4;
1212  delete [] k5[0];
1213  delete [] k5;
1214 
1215  start.clear();
1216  start.push_back(0);
1217  start.push_back(iscn);
1218  start.push_back(0);
1219 
1220  count.clear();
1221  count.push_back(scount[0]);
1222  count.push_back(scount[1]);
1223  count.push_back(scount[2]);
1224 
1225  // Output to L1B file
1226  if (!radiance)
1227  var = outfile.ncGrps[4].getVar( "rhot_SWIR");
1228  else
1229  var = outfile.ncGrps[4].getVar( "Lt_SWIR");
1230  var.putVar( start, count, &scal[0][0]);
1231 
1232  var = outfile.ncGrps[4].getVar( "qual_SWIR");
1233  var.putVar( start, count, &squal[0][0]);
1234 
1235  } else {
1236  cout << "No mirror side index for scan: " << iscn << endl;
1237  } // Check for valid mirror side
1238  } // Scan loop
1239 
1240  // End Main loop
1241 
1242  // Write spectral band information
1243  // Calculate band centers for aggregated hyperspectral bands
1244 
1245  if (bib >= 4) {
1246  // bgmat#bwave
1247  float *b1bwave = new float[bbb];
1248  float *sum = new float[bib];
1249  for (size_t i=0; i<bib; i++) {
1250  sum[i] = 0.0;
1251  for (size_t j=0; j<512; j++) {
1252  sum[i] += bwave[j]*bgmat[j][i];
1253  }
1254  }
1255 
1256  // bamat#sum
1257  for (size_t i=0; i<bbb; i++) {
1258  b1bwave[i] = 0.0;
1259  for (size_t j=0; j<bib; j++) {
1260  b1bwave[i] += bamat[j][i]*sum[j];
1261  }
1262  }
1263 
1264  start.clear();
1265  start.push_back(0);
1266 
1267  count.clear();
1268  count.push_back(bbb);
1269  var = outfile.ncGrps[0].getVar( "blue_wavelength");
1270  var.putVar( start, count, b1bwave);
1271 
1272  var = outfile.ncGrps[0].getVar( "blue_solar_irradiance");
1273  var.putVar( start, count, b1bf0);
1274 
1275  // Add m12_coef
1276  // b1bm12[ip,im,*] = bamat#bgmat#transpose(blue_lut.m12_coef[ip,im,*])
1277  // blue_m12_coef(blue_bands, HAM_sides, polarization_coefficients) ;
1278 
1279  dims[0] = bbb; dims[1] = 2; dims[2] = 3;
1280  float ***b1bm12 = make3dT<float>(dims);
1281  float ***b1bm13 = make3dT<float>(dims);
1282 
1283  for (size_t l=0; l<3; l++) {
1284  for (size_t m=0; m<2; m++) {
1285 
1286  for (size_t i=0; i<bib; i++) {
1287  sum[i] = 0.0;
1288  for (size_t j=0; j<NBWAVE; j++) {
1289  sum[i] += blue_lut.m12_coef[j][m][l]*bgmat[j][i];
1290  }
1291  }
1292 
1293  // bamat#sum
1294  for (size_t i=0; i<bbb; i++) {
1295  b1bm12[i][m][l] = 0.0;
1296  for (size_t j=0; j<bib; j++) {
1297  b1bm12[i][m][l] += bamat[j][i]*sum[j];
1298  }
1299  }
1300 
1301  for (size_t i=0; i<bib; i++) {
1302  sum[i] = 0.0;
1303  for (size_t j=0; j<NBWAVE; j++) {
1304  sum[i] += blue_lut.m13_coef[j][m][l]*bgmat[j][i];
1305  }
1306  }
1307 
1308  // bamat#sum
1309  for (size_t i=0; i<bbb; i++) {
1310  b1bm13[i][m][l] = 0.0;
1311  for (size_t j=0; j<bib; j++) {
1312  b1bm13[i][m][l] += bamat[j][i]*sum[j];
1313  }
1314  }
1315 
1316  }
1317  }
1318 
1319  start.clear();
1320  start.push_back(0); start.push_back(0); start.push_back(0);
1321 
1322  count.clear();
1323  count.push_back(bbb); count.push_back(2); count.push_back(3);
1324 
1325  var = outfile.ncGrps[0].getVar( "blue_m12_coef");
1326  var.putVar( start, count, &b1bm12[0][0][0]);
1327  var = outfile.ncGrps[0].getVar( "blue_m13_coef");
1328  var.putVar( start, count, &b1bm13[0][0][0]);
1329 
1330  delete [] b1bwave;
1331  delete [] sum;
1332 
1333  delete [] b1bm12[0][0];
1334  delete [] b1bm13[0][0];
1335  }
1336 
1337  if (rib >= 4) {
1338  float *r1bwave = new float[rbb];
1339  float *sum = new float[rib];
1340  for (size_t i=0; i<rib; i++) {
1341  sum[i] = 0.0;
1342  for (size_t j=0; j<512; j++) {
1343  sum[i] += rwave[j]*rgmat[j][i];
1344  }
1345  }
1346 
1347  for (size_t i=0; i<rbb; i++) {
1348  r1bwave[i] = 0.0;
1349  for (size_t j=0; j<rib; j++) {
1350  r1bwave[i] += ramat[j][i]*sum[j];
1351  }
1352  }
1353 
1354  start.clear();
1355  start.push_back(0);
1356 
1357  count.clear();
1358  count.push_back(rbb);
1359  var = outfile.ncGrps[0].getVar( "red_wavelength");
1360  var.putVar( start, count, r1bwave);
1361 
1362  var = outfile.ncGrps[0].getVar( "red_solar_irradiance");
1363  var.putVar( start, count, r1bf0);
1364 
1365  dims[0] = rbb; dims[1] = 2; dims[2] = 3;
1366  float ***r1bm12 = make3dT<float>(dims);
1367  float ***r1bm13 = make3dT<float>(dims);
1368 
1369  for (size_t l=0; l<3; l++) {
1370  for (size_t m=0; m<2; m++) {
1371 
1372  for (size_t i=0; i<rib; i++) {
1373  sum[i] = 0.0;
1374  for (size_t j=0; j<NRWAVE; j++) {
1375  sum[i] += red_lut.m12_coef[j][m][l]*rgmat[j][i];
1376  }
1377  }
1378 
1379  // bamat#sum
1380  for (size_t i=0; i<rbb; i++) {
1381  r1bm12[i][m][l] = 0.0;
1382  for (size_t j=0; j<rib; j++) {
1383  r1bm12[i][m][l] += ramat[j][i]*sum[j];
1384  }
1385  }
1386 
1387  for (size_t i=0; i<rib; i++) {
1388  sum[i] = 0.0;
1389  for (size_t j=0; j<NRWAVE; j++) {
1390  sum[i] += red_lut.m13_coef[j][m][l]*rgmat[j][i];
1391  }
1392  }
1393 
1394  // bamat#sum
1395  for (size_t i=0; i<rbb; i++) {
1396  r1bm13[i][m][l] = 0.0;
1397  for (size_t j=0; j<rib; j++) {
1398  r1bm13[i][m][l] += ramat[j][i]*sum[j];
1399  }
1400  }
1401 
1402  }
1403  }
1404 
1405  start.clear();
1406  start.push_back(0); start.push_back(0); start.push_back(0);
1407 
1408  count.clear();
1409  count.push_back(rbb); count.push_back(2); count.push_back(3);
1410 
1411  var = outfile.ncGrps[0].getVar( "red_m12_coef");
1412  var.putVar( start, count, &r1bm12[0][0][0]);
1413  var = outfile.ncGrps[0].getVar( "red_m13_coef");
1414  var.putVar( start, count, &r1bm13[0][0][0]);
1415 
1416  delete [] r1bwave;
1417  delete [] sum;
1418 
1419  delete [] r1bm12[0][0];
1420  delete [] r1bm13[0][0];
1421  }
1422 
1423  // SWIR wavelengths/bandpass
1424  start.clear();
1425  start.push_back(0);
1426  count.clear();
1427  count.push_back(NIWAVE);
1428 
1429  var = outfile.ncGrps[0].getVar( "SWIR_wavelength");
1430  var.putVar( start, count, swave);
1431  var = outfile.ncGrps[0].getVar( "SWIR_bandpass");
1432  var.putVar( start, count, spass);
1433 
1434  var = outfile.ncGrps[0].getVar( "SWIR_solar_irradiance");
1435  var.putVar( start, count, sf0);
1436 
1437  start.clear();
1438  start.push_back(0); start.push_back(0); start.push_back(0);
1439 
1440  count.clear();
1441  count.push_back(9); count.push_back(2); count.push_back(3);
1442 
1443  var = outfile.ncGrps[0].getVar( "SWIR_m12_coef");
1444  var.putVar( start, count, &swir_lut.m12_coef[0][0][0]);
1445  var = outfile.ncGrps[0].getVar( "SWIR_m13_coef");
1446  var.putVar( start, count, &swir_lut.m13_coef[0][0][0]);
1447 
1448 
1449  // l1b_filename.assign(argv[4]);
1450  outfile.write_granule_metadata( tstart, tend, l1b_filename);
1451 
1452  // write global attributes, including history and date_created
1453  set_global_attrs(outfile.l1bfile, history, doi, pversion);
1454  outfile.close();
1455 
1456  delete [] sstime;
1457  delete [] spin;
1458  delete [] hside;
1459  delete [] tfl;
1460  delete [] dtype;
1461  delete [] lines;
1462  delete [] iagg;
1463  delete [] bagg;
1464  delete [] ragg;
1465  delete [] mspin;
1466  delete [] ot_10us;
1467  delete [] enc_count;
1468  delete [] sgmat;
1469  delete [] thetap;
1470  delete [] thetas;
1471  delete [] ia;
1472  delete [] evtime;
1473 
1474  delete [] hamenc[0];
1475  delete [] hamenc;
1476  delete [] rtaenc[0];
1477  delete [] rtaenc;
1478  delete [] caltemps[0];
1479  delete [] caltemps;
1480  delete [] bdn[0];
1481  delete [] bdn;
1482  delete [] rdn[0];
1483  delete [] rdn;
1484  delete [] sdn[0];
1485  delete [] sdn;
1486  delete [] bcal[0];
1487  delete [] bcal;
1488  delete [] rcal[0];
1489  delete [] rcal;
1490  delete [] scal[0];
1491  delete [] scal;
1492  delete [] pview[0];
1493  delete [] pview;
1494  delete [] sview[0];
1495  delete [] sview;
1496 
1497  delete [] bsci[0];
1498  delete [] bsci;
1499  delete [] rsci[0];
1500  delete [] rsci;
1501  delete [] ssci[0];
1502  delete [] ssci;
1503 
1504  delete [] bcalb[0];
1505  delete [] bcalb;
1506  delete [] rcalb[0];
1507  delete [] rcalb;
1508 
1509  delete [] blue_lut.K1[0];
1510  delete [] blue_lut.K1;
1511  delete [] blue_lut.K5_coef[0];
1512  delete [] blue_lut.K5_coef;
1513 
1514  delete [] K3T;
1515 
1516  if (bgains.K1K2 != NULL) delete [] bgains.K1K2[0];
1517  if (bgains.K1K2 != NULL) delete [] bgains.K1K2;
1518  if (bgains.K5_coef != NULL) delete [] bgains.K5_coef[0];
1519  if (bgains.K5_coef != NULL) delete [] bgains.K5_coef;
1520 
1521  // delete [] arrT[0][0]; delete [] arrT[0]; delete [] arrT;
1522 
1523  delete [] red_lut.K1[0];
1524  delete [] red_lut.K1;
1525  delete [] red_lut.K5_coef[0];
1526  delete [] red_lut.K5_coef;
1527 
1528  if (rgains.K1K2 != NULL) delete [] rgains.K1K2[0];
1529  if (rgains.K1K2 != NULL) delete [] rgains.K1K2;
1530  if (rgains.K5_coef != NULL) delete [] rgains.K5_coef[0];
1531  if (rgains.K5_coef != NULL) delete [] rgains.K5_coef;
1532 
1533  delete [] swir_lut.K1[0];
1534  delete [] swir_lut.K1;
1535  delete [] swir_lut.K5_coef[0];
1536  delete [] swir_lut.K5_coef;
1537 
1538  delete [] sgains.K1K2[0];
1539  delete [] sgains.K1K2;
1540  delete [] sgains.K5_coef[0];
1541  delete [] sgains.K5_coef;
1542 
1543  // Add delete for dark arrays
1544  delete [] bdark[0][0];
1545 
1546  clo_deleteList(optionList);
1547 
1548  return 0;
1549 }
1550 
1551 
1552 int read_oci_cal_lut( NcFile *calLUTfile, string tag, NcGroup gidLUT,
1553  uint32_t& banddim, uint32_t mcedim, uint32_t& nldim,
1554  uint32_t& poldim, cal_lut_struct& cal_lut) {
1555 
1556  size_t dims[3],dims4[4];
1557  NcDim ncDIM;
1558  uint32_t timedim, tempdim, tcdim, rvsdim, msdim=2;
1559  string bandname;
1560  bandname.assign( tag);
1561  bandname.append( "_bands");
1562 
1563  ncDIM = calLUTfile->getDim(bandname.c_str());
1564  banddim = ncDIM.getSize();
1565 
1566  ncDIM = calLUTfile->getDim("number_of_times");
1567  timedim = ncDIM.getSize();
1568 
1569  if (tag.compare("blue") == 0 || tag.compare("red") == 0)
1570  ncDIM = calLUTfile->getDim("number_of_CCD_temperatures");
1571  else
1572  ncDIM = calLUTfile->getDim("number_of_SWIR_temperatures");
1573  tempdim = ncDIM.getSize();
1574 
1575  ncDIM = calLUTfile->getDim("number_of_T_coefficients");
1576  tcdim = ncDIM.getSize();
1577  ncDIM = calLUTfile->getDim("number_of_RVS_coefficients");
1578  rvsdim = ncDIM.getSize();
1579  ncDIM = calLUTfile->getDim("number_of_nonlinearity_coefficients");
1580  nldim = ncDIM.getSize();
1581  ncDIM = calLUTfile->getDim("number_of_polarization_coefficients");
1582  poldim = ncDIM.getSize();
1583 
1584  float**K1 = new float *[banddim];
1585  K1[0] = new float[msdim*banddim];
1586  for (size_t i=1; i<banddim; i++) K1[i] = K1[i-1] + msdim;
1587  cal_lut.K1 = K1;
1588 
1589  dims[0] = banddim; dims[1] = msdim; dims[2] = timedim;
1590  float ***K2 = make3dT<float>(dims);
1591  cal_lut.K2 = K2;
1592  dims[0] = banddim; dims[1] = tempdim; dims[2] = tcdim;
1593  float ***K3_coef = make3dT<float>(dims);
1594  cal_lut.K3_coef = K3_coef;
1595  dims4[0] = banddim; dims4[1] = msdim; dims4[2] = mcedim; dims4[3] = rvsdim;
1596  float ****K4_coef = make4dT<float>(dims4);
1597  cal_lut.K4_coef = K4_coef;
1598 
1599  double **K5_coef = new double *[banddim];
1600  K5_coef[0] = new double[nldim*banddim];
1601  for (size_t i=1; i<banddim; i++) K5_coef[i] = K5_coef[i-1] + nldim;
1602  cal_lut.K5_coef = K5_coef;
1603 
1604  uint32_t *sat_thres = new uint32_t [banddim];
1605  cal_lut.sat_thres = sat_thres;
1606 
1607  dims[0] = banddim; dims[1] = msdim; dims[2] = poldim;
1608  float ***m12_coef = make3dT<float>(dims);
1609  cal_lut.m12_coef = m12_coef;
1610  float ***m13_coef = make3dT<float>(dims);
1611  cal_lut.m13_coef = m13_coef;
1612 
1613  // Assign cal lut dimensions array
1614  cal_lut.ldims[0] = timedim; cal_lut.ldims[1] = tempdim;
1615  cal_lut.ldims[2] = tcdim; cal_lut.ldims[3] = rvsdim;
1616  cal_lut.ldims[4] = nldim; cal_lut.ldims[5] = msdim;
1617  cal_lut.ldims[6] = poldim;
1618 
1619  NcVar var;
1620  var = gidLUT.getVar( "K1");
1621  var.getVar( &cal_lut.K1[0][0]);
1622  var = gidLUT.getVar( "K2");
1623  var.getVar( &cal_lut.K2[0][0][0]);
1624  var = gidLUT.getVar( "K3_coef");
1625  var.getVar( &cal_lut.K3_coef[0][0][0]);
1626  var = gidLUT.getVar( "K4_coef");
1627  var.getVar( &cal_lut.K4_coef[0][0][0][0]);
1628  var = gidLUT.getVar( "K5_coef");
1629  var.getVar( &cal_lut.K5_coef[0][0]);
1630  var = gidLUT.getVar( "sat_thres");
1631  var.getVar( &cal_lut.sat_thres[0]);
1632 
1633  var = gidLUT.getVar( "m12_coef");
1634  var.getVar( &cal_lut.m12_coef[0][0][0]);
1635  var = gidLUT.getVar( "m13_coef");
1636  var.getVar( &cal_lut.m13_coef[0][0][0]);
1637 
1638  return 0;
1639 }
1640 
1641 // float K1K2[nib][msdim]: absolute and temporal gain
1642 // float K3_coef[nib][tempdim][tcdim]: temperature correction
1643 // float K4_coef[nib][msdim][mcedim][rvsdim]: RVS correction
1644 // double K5_coef[nib][nldim]: linearity correction
1645 
1646 int make_oci_gains( uint32_t nib, uint32_t banddim, uint16_t iyr, uint32_t jd,
1647  double stime, size_t numTimes, double *K2t, int16_t board_id,
1648  int16_t iagg, int16_t *jagg, cal_lut_struct& cal_lut,
1649  float **gmat, gains_struct& gains) {
1650 
1651  for (size_t i=0; i<6; i++) gains.ldims[i] = cal_lut.ldims[i];
1652 
1653  uint16_t timedim = gains.ldims[0];
1654  uint16_t tempdim = gains.ldims[1];
1655  uint16_t tcdim = gains.ldims[2];
1656  uint16_t rvsdim = gains.ldims[3];
1657  uint16_t nldim = gains.ldims[4];
1658  uint16_t msdim = gains.ldims[5];
1659 
1660  bool hyper = false;
1661  uint16_t bd_id;
1662  int16_t *iaf = NULL;
1663 
1664  // Hyperspectral bands
1665  if ( board_id == -1) {
1666  hyper = true;
1667  iaf = new int16_t[nib];
1668  int ib = 0;
1669  for (size_t i=0; i<16; i++) {
1670  if ( jagg[i] > 0) {
1671  uint32_t nb = 32/jagg[i];
1672  for (size_t j=0; j<nb; j++) {
1673  if (iagg*jagg[i] < 4)
1674  iaf[ib+j] = 4/(iagg*jagg[i]);
1675  else
1676  iaf[ib+j] = 4/4;
1677  }
1678  ib += nb;
1679  }
1680  }
1681  } else bd_id = board_id % 2;
1682 
1683  size_t dims[3];
1684 
1685  float **K1K2 = new float *[nib];
1686  K1K2[0] = new float[nib*msdim];
1687  for (size_t i=1; i<nib; i++) K1K2[i] = K1K2[i-1] + msdim;
1688  gains.K1K2 = K1K2;
1689  dims[0] = nib; dims[1] = tempdim; dims[2] = tcdim;
1690  float ***K3_coef = make3dT<float>(dims);
1691  gains.K3_coef = K3_coef;
1692  dims[0] = nib; dims[1] = msdim; dims[2] = rvsdim;
1693  float ***K4_coef = make3dT<float>(dims);
1694  gains.K4_coef = K4_coef;
1695 
1696  double **K5_coef = new double *[nib];
1697  K5_coef[0] = new double[nib*nldim];
1698  for (size_t i=1; i<nib; i++) K5_coef[i] = K5_coef[i-1] + nldim;
1699  gains.K5_coef = K5_coef;
1700 
1701  uint32_t *sat_thres = new uint32_t[nib];
1702  gains.sat_thres = sat_thres;
1703 
1704  // Mirror-side dependent gains
1705  double *K2 = new double[banddim];
1706  for (size_t ms=0; ms<msdim; ms++) {
1707 
1708  // Get temporal gain and combine with absolute gain
1709  double d2 = jd - 2451545 + stime/86400.0;
1710 
1711  size_t kd=0;
1712  for (size_t j=numTimes-1; j>=0; j--) {
1713  if ( d2 > K2t[j]) {
1714  kd = j;
1715  break;
1716  }
1717  }
1718  if ( kd < (size_t) (timedim-1)) {
1719  double ff = (d2 - K2t[kd]) / (K2t[kd+1] - K2t[kd]);
1720  for (size_t j=0; j<banddim; j++)
1721  K2[j] = cal_lut.K2[j][ms][kd]*(1.0-ff) + cal_lut.K2[j][ms][kd+1]*ff;
1722  } else {
1723  for (size_t j=0; j<banddim; j++) K2[j] = cal_lut.K2[j][ms][kd];
1724  }
1725 
1726  int16_t iaf0 = 1;
1727  for (size_t i=0; i<nib; i++) {
1728  gains.K1K2[i][ms] = 0;
1729  if (hyper) iaf0 = iaf[i];
1730  for (size_t j=0; j<banddim; j++)
1731  gains.K1K2[i][ms] += gmat[j][i] * cal_lut.K1[j][ms]*K2[j]*iaf0;
1732  }
1733 
1734  // Generate RVS coefficents
1735  for (size_t i=0; i<nib; i++) {
1736  for (size_t k=0; k<rvsdim; k++) {
1737  gains.K4_coef[i][ms][k] = 0;
1738  for (size_t j=0; j<banddim; j++) {
1739  if (hyper)
1740  gains.K4_coef[i][ms][k] +=
1741  gmat[j][i] * cal_lut.K4_coef[j][ms][0][k];
1742  else
1743  gains.K4_coef[i][ms][k] +=
1744  gmat[j][i] * cal_lut.K4_coef[j][ms][bd_id][k];
1745  }
1746  }
1747  }
1748  }
1749  delete [] K2;
1750 
1751  // Generate temperature coefficients
1752  for (size_t i=0; i<nib; i++) {
1753  for (size_t k=0; k<tcdim; k++) {
1754  for (size_t l=0; l<tempdim; l++) {
1755  gains.K3_coef[i][l][k] = 0;
1756  for (size_t j=0; j<banddim; j++)
1757  gains.K3_coef[i][l][k] += gmat[j][i] * cal_lut.K3_coef[j][l][k];
1758  }
1759  }
1760  }
1761 
1762  // Generate linearity coefficients
1763  for (size_t i=0; i<nib; i++) {
1764  for (size_t k=0; k<nldim; k++) {
1765  gains.K5_coef[i][k] = 0;
1766  for (size_t j=0; j<banddim; j++) {
1767  if (hyper)
1768  gains.K5_coef[i][k] +=
1769  gmat[j][i] * cal_lut.K5_coef[j][k] * powf(iaf[i], (float) i);
1770  else
1771  gains.K5_coef[i][k] +=
1772  gmat[j][i] * cal_lut.K5_coef[j][k];
1773  }
1774  }
1775  }
1776 
1777  // Generate saturation thresholds
1778  for (size_t i=0; i<nib; i++) {
1779  gains.sat_thres[i] = 0;
1780  for (size_t j=0; j<banddim; j++) {
1781  if (hyper)
1782  gains.sat_thres[i] += gmat[j][i] * cal_lut.sat_thres[j] / iaf[i];
1783  else
1784  gains.sat_thres[i] += gmat[j][i] * cal_lut.sat_thres[j];
1785  }
1786  }
1787 
1788  if (hyper) delete[] iaf;
1789 
1790  return 0;
1791 }
1792 
1793 
1794 template <typename T>
1795 int get_oci_dark( size_t iscn, uint32_t nscan, uint8_t *hside, uint16_t ndsc,
1796  uint16_t nskp, int16_t iags, int16_t iagd, uint32_t ntaps,
1797  int16_t *jagg, uint32_t dfill, int16_t ndc, T ***dark,
1798  uint32_t nib, float *dc, int16_t& iret) {
1799 
1800  // Program to generate dark corrections for OCI data by averaging the
1801  // dark collect data and correcting for bit shift/truncation if necessary
1802 
1803  // Determine number of bands per tap for hyperspectral data
1804  int16_t *nbndt = new int16_t[ntaps];
1805 
1806  if (ntaps == 16) {
1807  // hyperspectral bands
1808  for (size_t i=0; i<ntaps; i++)
1809  if ( jagg[i] > 0) nbndt[i] = 32 / jagg[i]; else nbndt[i] = 0;
1810  } else {
1811  for (size_t i=0; i<ntaps; i++) nbndt[i] = 9;
1812  }
1813  int16_t nbnd = 0;
1814  for (size_t i=0; i<ntaps; i++) nbnd += nbndt[i];
1815 
1816  // Select data for HAM side and determine scan indices
1817  int32_t *kh = new int32_t[nscan];
1818  int32_t nkh=0;
1819  for (size_t i=0; i<nscan; i++) {
1820  if ( hside[i] == hside[iscn]) {
1821  kh[i] = (int32_t) i;
1822  nkh++;
1823  } else {
1824  kh[i] = -1;
1825  }
1826  }
1827 
1828  int32_t js=0;
1829  for (size_t i=0; i<nscan; i++) {
1830  if ( kh[i] == (int32_t) iscn) {
1831  js = (int32_t) i;
1832  break;
1833  }
1834  }
1835 
1836  // Check for valid dark collect data within specified range
1837  uint16_t ndscl = ndsc;
1838  bool valid_dark_found = false;
1839  int32_t is1=js, is2=js;
1840 
1841  while (!valid_dark_found && ndscl <= nkh) {
1842  if ( ndsc > 1) {
1843  is1 = js - ndsc/2;
1844  is2 = js + ndsc/2;
1845  // Check for start or end of granule
1846  if (is1 < 0) {
1847  is1 = 0;
1848  is2 = ndsc - 1;
1849  }
1850  if (is2 >= nkh) {
1851  is1 = nkh - ndsc;
1852  is2 = nkh - 1;
1853  }
1854  }
1855 
1856  // If no valid dark data, expand scan range
1857  for (size_t i=is1; i<=(size_t) is2; i++) {
1858  for (size_t j=nskp; j<(size_t) ndc; j++) {
1859  if ( dark[kh[i]][0][j] != dfill) {
1860  valid_dark_found = true;
1861  break;
1862  }
1863  }
1864  }
1865  if ( !valid_dark_found) ndscl += 2;
1866  }
1867 
1868  if ( !valid_dark_found) {
1869  iret =-1;
1870  return 0;
1871  }
1872 
1873  // Loop through taps and compute dark correction
1874  int16_t ibnd=0;
1875  for (size_t i=0; i<ntaps; i++) {
1876  if ( jagg[i] > 0) {
1877  float ddiv = 1.0;
1878  float doff = 0.0;
1879  if ( iags*jagg[i] > 4) {
1880  ddiv = iagd * jagg[i] / 4.0;
1881  doff = (ddiv-1) / (2*ddiv);
1882  }
1883 
1884  for (size_t j=0; j<(size_t) nbndt[i]; j++) {
1885  float sum = 0.0;
1886  int nv = 0;
1887  for (size_t k=kh[is1]; k<=(size_t) kh[is2]; k++) {
1888  for (size_t l=nskp; l<(size_t) ndc; l++) {
1889  if (dark[k][ibnd+j][l] != dfill) {
1890  sum += dark[k][ibnd+j][l];
1891  nv++;
1892  }
1893  }
1894  }
1895  dc[ibnd+j] = ((sum/(nv))/ddiv) - doff;
1896  } // j loop
1897  } // if (
1898  ibnd += nbndt[i];
1899  } // i loop
1900 
1901  delete [] kh;
1902  delete [] nbndt;
1903 
1904  iret = 0;
1905  if (ndscl > ndsc) iret = 1;
1906 
1907  return 0;
1908 }
1909 
1910 
1911 int get_oci_temp_corr( uint32_t nib, gains_struct gains, float *K3T,
1912  float *caltemps, uint32_t nscan, float *k3) {
1913 
1914  uint16_t tempdim = gains.ldims[1];
1915  uint16_t tcdim = gains.ldims[2];
1916 
1917  for (size_t i=0; i<nib; i++) k3[i] = 1.0;
1918 
1919  for (size_t i=0; i<tempdim; i++) {
1920  float td = caltemps[i] - K3T[i];
1921  for (size_t j=0; j<tcdim; j++) {
1922  for (size_t k=0; k<nib; k++) {
1923  k3[k] -= gains.K3_coef[k][i][j] * powf(td, j+1);
1924  }
1925  }
1926  }
1927 
1928  return 0;
1929 }
1930 
1931 int get_oci_rvs_corr( uint32_t nib, uint16_t pdim, uint8_t hside,
1932  gains_struct gains, double *theta, float **k4) {
1933 
1934  // Program to compute RVS correction from coefficients and scan angle
1935 
1936  uint16_t rvsdim = gains.ldims[3];
1937 
1938  for (size_t i=0; i<nib; i++)
1939  for (size_t j=0; j<pdim; j++)
1940  k4[i][j] = 1.0;
1941 
1942  for (size_t i=0; i<rvsdim; i++)
1943  for (size_t j=0; j<nib; j++)
1944  for (size_t k=0; k<pdim; k++)
1945  k4[j][k] += gains.K4_coef[j][hside][i] * powf(theta[k], i+1);
1946 
1947  return 0;
1948 }
1949 
1950 int get_oci_lin_corr( uint32_t nib, uint16_t pdim, uint32_t nldim,
1951  gains_struct gains, float **dn, float **k5) {
1952 
1953  for (size_t i=0; i<pdim; i++) {
1954  // Zeroth-order correction
1955  for (size_t j=0; j<nib; j++) {
1956  k5[j][i] = gains.K5_coef[j][0];
1957 
1958  for (size_t k=1; k<nldim; k++) {
1959  k5[j][i] += gains.K5_coef[j][k]*powf(dn[j][i], k);
1960  }
1961  }
1962  }
1963 
1964  return 0;
1965 }
1966 
1967 
1968 int get_oci_cal_temps( NcFile *l1afile, NcGroup egid,
1969  uint16_t ntemps, uint32_t nscan, double *evtime,
1970  float **caltemps) {
1971 
1972  // Program to read temperatures used in calibration from L1A file and
1973  // interpolate to scan times
1974  // The order of temperatures is:
1975  // Lens housings: blue CCD side, blue grating side,
1976  // red CCD side, red grating side
1977  // CCDs: red right, red left,
1978  // blue right, blue left
1979  // SDA detector 1 - 16
1980  // AOB 1, 2, 3, 4, 7, 8
1981  // MOSB near MLA
1982 
1983  NcDim ncDIM;
1984  uint32_t tlmdim, daudim, icdudim;
1985  ncDIM = l1afile->getDim("tlm_packets");
1986  tlmdim = ncDIM.getSize();
1987  ncDIM = l1afile->getDim("DAUC_temps");
1988  daudim = ncDIM.getSize();
1989  ncDIM = l1afile->getDim("ICDU_therm");
1990  icdudim = ncDIM.getSize();
1991 
1992  double *dauctime = new double[tlmdim];
1993 
1994  float **dauctemp = new float *[tlmdim];
1995  dauctemp[0] = new float[daudim*tlmdim];
1996  for (size_t i=1; i<tlmdim; i++) dauctemp[i] = dauctemp[i-1] + daudim;
1997 
1998  double *icdutime = new double[tlmdim];
1999 
2000  float **icdutherm = new float *[tlmdim];
2001  icdutherm[0] = new float[icdudim*tlmdim];
2002  for (size_t i=1; i<tlmdim; i++) icdutherm[i] = icdutherm[i-1] + icdudim;
2003 
2004  NcVar var;
2005  var = egid.getVar( "DAUC_temp_time");
2006  var.getVar( dauctime);
2007  var = egid.getVar( "DAUC_temperatures");
2008  var.getVar( &dauctemp[0][0]);
2009  var = egid.getVar( "TC_tlm_time");
2010  var.getVar( icdutime);
2011  var = egid.getVar( "ICDU_thermisters");
2012  var.getVar( &icdutherm[0][0]);
2013 
2014  /*
2015 DAUCTEMP FLOAT = Array[69, 301]
2016 DAUCTIME DOUBLE = Array[301]
2017 ICDUTHERM FLOAT = Array[74, 301]
2018 ICDUTIME DOUBLE = Array[301]
2019 
2020 ICDU_therm = 74 ;
2021  */
2022 
2023  // Indices of required temperatures
2024  // MLA and lens housings (blue CCD, blue grating, red CCD, red grating)
2025  uint32_t iicdu[5] = {23,24,25,26,11};
2026  // Red and blue CCDs, SDA detectors, AOB 1, 2, 3, 4, 7, 8
2027  uint32_t idauc[26] = {5,6,12,13,
2028  14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,
2029  30,31,32,33,36,37};
2030 
2031  // Loop through temperatures
2032  double *dstime = new double[nscan];
2033  for (size_t i=0; i<nscan; i++) dstime[i] = evtime[i] - evtime[0];
2034 
2035  // ICDU thermistors
2036  double *ditime = new double[tlmdim];
2037  double *dummy = new double[tlmdim];
2038  double c0, c1, cov00, cov01, cov11, sumsq;
2039  for (size_t i=0; i<=3; i++) {
2040  uint32_t k=0;
2041  for (size_t j=0; j<tlmdim; j++)
2042  if ( icdutime[j] > 0) {
2043  ditime[k] = icdutime[j] - evtime[0];
2044  dummy[k++] = icdutherm[j][iicdu[i]];
2045  }
2046 
2047  gsl_fit_linear(ditime, 1, dummy, 1, k,
2048  &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
2049 
2050  for (size_t j=0; j<nscan; j++) caltemps[j][i] = c0 + c1*dstime[j];
2051  }
2052 
2053  size_t k=0;
2054  for (size_t j=0; j<tlmdim; j++)
2055  if ( icdutime[j] > 0)
2056  dummy[k++] = icdutherm[j][iicdu[4]];
2057 
2058  gsl_fit_linear(ditime, 1, dummy, 1, k,
2059  &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
2060  for (size_t j=0; j<nscan; j++) caltemps[j][30] = c0 + c1*dstime[j];
2061 
2062  // DAUC temperatures
2063  double *ddtime = new double[tlmdim];
2064  for (size_t i=0; i<=25; i++) {
2065  size_t k=0;
2066  for (size_t j=0; j<tlmdim; j++)
2067  if ( dauctime[j] > 0) {
2068  ddtime[k] = dauctime[j] - evtime[0];
2069  dummy[k++] = dauctemp[j][idauc[i]];
2070  }
2071 
2072  gsl_fit_linear(ddtime, 1, dummy, 1, k,
2073  &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
2074 
2075  for (size_t j=0; j<nscan; j++) caltemps[j][i+4] = c0 + c1*dstime[j];
2076  }
2077 
2078  delete [] dauctime;
2079  delete [] dauctemp[0];
2080  delete [] dauctemp;
2081 
2082  delete [] icdutime;
2083  delete [] icdutherm[0];
2084  delete [] icdutherm;
2085 
2086  delete [] ditime;
2087  delete [] ddtime;
2088 
2089  return 0;
2090 }
2091 
2092 
2093 int l1bFile::write_granule_metadata( std::string tstart, std::string tend,
2094  std::string l1b_name) {
2095 
2096  l1bfile->putAtt("time_coverage_start", tstart.c_str());
2097  l1bfile->putAtt("time_coverage_end", tend.c_str());
2098 
2099  // Write product file name
2100  l1bfile->putAtt("product_name", l1b_name.c_str());
2101 
2102  return 0;
2103 }
2104 
2105 
2106 int l1bFile::close() {
2107 
2108  try {
2109  l1bfile->close();
2110  }
2111  catch ( NcException& e) {
2112  cout << e.what() << endl;
2113  cerr << "Failure closing: " + fileName << endl;
2114  exit(1);
2115  }
2116 
2117  return 0;
2118 }
2119 
2120 
2121 template <typename T>
2122 T*** make3dT( size_t dims[3]) {
2123 
2124  T ***arr3d = new T **[dims[0]];
2125 
2126  arr3d[0] = new T*[dims[0]*dims[1]];
2127  arr3d[0][0] = new T[dims[0]*dims[1]*dims[2]];
2128 
2129  for (size_t i=1; i<dims[0]; i++) arr3d[i] = arr3d[i-1] + dims[1];
2130 
2131  for (size_t i=0; i<dims[0]; i++) {
2132  if ( i > 0) arr3d[i][0] = arr3d[i-1][0] + dims[1]*dims[2];
2133  for (size_t j=1; j<dims[1]; j++)
2134  arr3d[i][j] = arr3d[i][j-1] + dims[2];
2135  }
2136 
2137  return arr3d;
2138 }
2139 
2140 template <typename T>
2141 T**** make4dT( size_t dims[4]) {
2142 
2143  T ****arr4d = new T ***[dims[0]];
2144 
2145  arr4d[0] = new T**[dims[0]*dims[1]];
2146  arr4d[0][0] = new T*[dims[0]*dims[1]*dims[2]];
2147  arr4d[0][0][0] = new T[dims[0]*dims[1]*dims[2]*dims[3]];
2148 
2149  for (size_t i=0; i<dims[0]; i++) {
2150  if ( i > 0) {
2151  arr4d[i] = arr4d[i-1] + dims[1];
2152  arr4d[i][0] = arr4d[i-1][0] + dims[1]*dims[2];
2153  arr4d[i][0][0] = arr4d[i-1][0][0] + dims[1]*dims[2]*dims[3];
2154  }
2155  for (size_t j=0; j<dims[1]; j++) {
2156  if ( j > 0) {
2157  arr4d[i][j] = arr4d[i][j-1] + dims[2];
2158  arr4d[i][j][0] = arr4d[i][j-1][0] + dims[2]*dims[3];
2159  }
2160  for (size_t k=1; k<dims[2]; k++) {
2161  arr4d[i][j][k] = arr4d[i][j][k-1] + dims[3];
2162  }
2163  }
2164  }
2165 
2166  return arr4d;
2167 }
CHUNK_CACHE_NELEMS
#define CHUNK_CACHE_NELEMS
Definition: l1agen_oci.h:16
read_oci_cal_lut
int read_oci_cal_lut(NcFile *calLUTfile, string tag, NcGroup gidLUT, uint32_t &banddim, uint32_t mcedim, uint32_t &nldim, uint32_t &poldim, cal_lut_struct &cal_lut)
Definition: l1bgen_oci.cpp:1552
imn
int32_t imn
Definition: atrem_corl1.h:161
clo_addOption
clo_option_t * clo_addOption(clo_optionList_t *list, const char *key, enum clo_dataType_t dataType, const char *defaultVal, const char *desc)
Definition: clo.c:684
get_oci_lin_corr
int get_oci_lin_corr(uint32_t nib, uint16_t pdim, uint32_t nldim, gains_struct gains, float **dn, float **k5)
Definition: l1bgen_oci.cpp:1950
l1bgen_oci.h
j
int j
Definition: decode_rs.h:73
geo_struct
Definition: common.h:25
main
int main(int argc, char *argv[])
Definition: l1bgen_oci.cpp:65
cal_lut_struct::K5_coef
double ** K5_coef
Definition: l1bgen_oci.h:26
clo_getString
char * clo_getString(clo_optionList_t *list, const char *key)
Definition: clo.c:1357
clo_readArgs
void clo_readArgs(clo_optionList_t *list, int argc, char *argv[])
Definition: clo.c:2103
make_oci_gains
int make_oci_gains(uint32_t nib, uint32_t banddim, uint16_t iyr, uint32_t jd, double stime, size_t numTimes, double *K2t, int16_t board_id, int16_t iagg, int16_t *jagg, cal_lut_struct &cal_lut, float **gmat, gains_struct &gains)
Definition: l1bgen_oci.cpp:1646
NULL
#define NULL
Definition: decode_rs.h:63
spin
void spin(double st, double *pos1, double *pos2)
check_scan_times
int check_scan_times(uint32_t nscan, double *sstime, short *sfl)
Definition: common.cpp:863
sumsq
unsigned long long sumsq(signed short *in, int cnt)
clo_isSet
int clo_isSet(clo_optionList_t *list, const char *key)
Definition: clo.c:2270
make4dT
T **** make4dT(size_t dims[4])
Definition: l1bgen_oci.cpp:2141
geo_struct::as_planarity
double as_planarity[5]
Definition: common.h:44
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
nscan
int32 nscan
Definition: l1_czcs_hdf.c:19
string
@ string
Definition: GEO_read_param_file.c:64
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
CLO_TYPE_BOOL
@ CLO_TYPE_BOOL
Definition: clo.h:81
nbnd
int nbnd
Definition: get_cmp.c:30
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
names::if
character(len=1000) if
Definition: names.f90:13
K1
#define K1
Definition: get_cal_swf.c:89
gains_struct::K3_coef
float *** K3_coef
Definition: l1bgen_oci.h:35
CHUNK_CACHE_SIZE
#define CHUNK_CACHE_SIZE
Definition: l1agen_oci.h:15
ff
float ff(float u, float v)
Definition: par_utils.c:1461
clo_optionList_t
Definition: clo.h:126
cal_lut_struct::K1
float ** K1
Definition: l1bgen_oci.h:22
cal_lut_struct
Definition: l1bgen_oci.h:20
gains_struct
Definition: l1bgen_oci.h:32
clo_createList
clo_optionList_t * clo_createList()
Definition: clo.c:532
geolocate_oci
int geolocate_oci(std::string l1a_filename, std::string geo_lut_filename, geo_struct &geoLUT, std::string l1b_filename, std::string dem_file, bool radiance, std::string doi, const std::string ephFile, std::string pversion)
call_sequence
string call_sequence(int argc, char *argv[])
Definition: global_attrs.cpp:15
get_oci_vecs
int get_oci_vecs(uint32_t nscan, uint16_t pdim, double as_planarity[5], double at_planarity[5], int32_t *rta_nadir, double ham_ct_angles[2], double ev_toff, int32_t spin, uint8_t hside, float *clines, double *delt, double revpsec, int32_t ppr_off, int16_t board_id, uint32_t nmcescan, int32_t *mspin, uint8_t *enc_count, float **hamenc, float **rtaenc, float **pview, double *theta, int16_t &iret)
Definition: common.cpp:195
get_history
string get_history(int argc, char *argv[])
Definition: DDOptions.cpp:496
get_oci_temp_corr
int get_oci_temp_corr(uint32_t nib, gains_struct gains, float *K3T, float *caltemps, uint32_t nscan, float *k3)
Definition: l1bgen_oci.cpp:1911
set_global_attrs
void set_global_attrs(NcFile *outfile, string history, string doi, string pversion)
Definition: global_attrs.cpp:40
gains_struct::sat_thres
uint32_t * sat_thres
Definition: l1bgen_oci.h:38
hktgen_pace.att
int att
Definition: hktgen_pace.py:28
get_oci_rvs_corr
int get_oci_rvs_corr(uint32_t nib, uint16_t pdim, uint8_t hside, gains_struct gains, double *theta, float **k4)
Definition: l1bgen_oci.cpp:1931
cal_lut_struct::m13_coef
float *** m13_coef
Definition: l1bgen_oci.h:29
clo_printUsage
void clo_printUsage(clo_optionList_t *list)
Definition: clo.c:1988
oci_hdr_strip.argc
argc
Definition: oci_hdr_strip.py:7
jd
Definition: jd.py:1
geo_struct::ham_ct_angles
double ham_ct_angles[2]
Definition: common.h:38
cal_lut_struct::ldims
uint16_t ldims[7]
Definition: l1bgen_oci.h:21
CLO_TYPE_IFILE
@ CLO_TYPE_IFILE
Definition: clo.h:87
tempOut
ofstream tempOut
Definition: l1bgen_oci.cpp:63
CLO_TYPE_OFILE
@ CLO_TYPE_OFILE
Definition: clo.h:88
cal_lut_struct::K2
float *** K2
Definition: l1bgen_oci.h:23
oci_hdr_strip.tmp_filename
string tmp_filename
Definition: oci_hdr_strip.py:9
get_nib_nbb
int get_nib_nbb(uint32_t ntaps, size_t *ia, uint32_t ntb[16], int16_t jagg[16], uint32_t &nib, uint32_t &nbb)
Definition: common.cpp:742
l1bFile::write_granule_metadata
int write_granule_metadata(std::string tstart, std::string tend, std::string l1b_name)
Definition: l1bgen_oci.cpp:2093
run_local.outfile
outfile
Definition: run_local.py:76
ncattredit.history
string history
Definition: ncattredit.py:30
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
NIWAVE
#define NIWAVE
Definition: l1bgen_oci.h:16
cal_lut_struct::K4_coef
float **** K4_coef
Definition: l1bgen_oci.h:25
get_oci_cal_temps
int get_oci_cal_temps(NcFile *l1afile, NcGroup egid, uint16_t ntemps, uint32_t nscan, double *evtime, float **caltemps)
Definition: l1bgen_oci.cpp:1968
optionList
clo_optionList_t * optionList
Definition: main_l2extract.c:68
parse_file_name
void parse_file_name(const char *inpath, char *outpath)
Definition: parse_file_name.c:23
gains_struct::K1K2
float ** K1K2
Definition: l1bgen_oci.h:34
fix_mac_rpath.lines
lines
Definition: fix_mac_rpath.py:35
dtype
dtype
Definition: DDataset.hpp:31
get_oci_dark
int get_oci_dark(size_t iscn, uint32_t nscan, uint8_t *hside, uint16_t ndsc, uint16_t nskp, int16_t iags, int16_t iagd, uint32_t ntaps, int16_t *jagg, uint32_t dfill, int16_t ndc, T ***dark, uint32_t nib, float *dc, int16_t &iret)
Definition: l1bgen_oci.cpp:1795
global_attrs.h
clo_deleteList
void clo_deleteList(clo_optionList_t *list)
Definition: clo.c:875
cal_lut_struct::K3_coef
float *** K3_coef
Definition: l1bgen_oci.h:24
gains_struct::K4_coef
float *** K4_coef
Definition: l1bgen_oci.h:36
geo_struct::at_planarity
double at_planarity[5]
Definition: common.h:45
read_mce_tlm
int read_mce_tlm(NcFile *l1afile, geo_struct &geo_lut, NcGroup egid, uint32_t nmcescan, uint32_t nenc, int32_t &ppr_off, double &revpsec, double &secpline, int16_t &board_id, int32_t *mspin, int32_t *ot_10us, uint8_t *enc_count, float **hamenc, float **rtaenc, int16_t &iret)
Definition: common.cpp:8
MAX_SOLZ
#define MAX_SOLZ
Definition: l1bgen_oci.cpp:29
l1bFile::close
int close()
Definition: l1bgen_oci.cpp:2106
NBWAVE
#define NBWAVE
Definition: l1bgen_oci.h:14
gains_struct::ldims
uint16_t ldims[6]
Definition: l1bgen_oci.h:33
CHUNK_CACHE_PREEMPTION
#define CHUNK_CACHE_PREEMPTION
Definition: l1agen_oci.h:17
clo_getBool
int clo_getBool(clo_optionList_t *list, const char *key)
Definition: clo.c:1375
cal_lut_struct::m12_coef
float *** m12_coef
Definition: l1bgen_oci.h:28
VERSION
#define VERSION
Definition: l1bgen_oci.cpp:26
NRWAVE
#define NRWAVE
Definition: l1bgen_oci.h:15
i
int i
Definition: decode_rs.h:71
get_ev
int get_ev(double secpline, int16_t *dtype, int16_t *lines, int16_t *iagg, uint16_t &pcdim, uint16_t &psdim, double &ev_toff, float *clines, float *slines, double *deltc, double *delts, bool dark, int16_t &iret)
Definition: common.cpp:142
CLO_TYPE_STRING
@ CLO_TYPE_STRING
Definition: clo.h:86
geo_struct::rta_nadir
int32_t rta_nadir[2]
Definition: common.h:42
iyr
int32_t iyr
Definition: atrem_corl1.h:161
k
int k
Definition: decode_rs.h:73
make3dT
T *** make3dT(size_t dims[3])
Definition: l1bgen_oci.cpp:2122
gains_struct::K5_coef
double ** K5_coef
Definition: l1bgen_oci.h:37
nb
int32_t nb
Definition: atrem_corl1.h:132
get_agg_mat
int get_agg_mat(size_t *ia, int16_t jagg[16], uint32_t ntb[16], uint32_t nib, uint32_t nbb, float **amat, float **gmat)
Definition: common.cpp:777
cal_lut_struct::sat_thres
uint32_t * sat_thres
Definition: l1bgen_oci.h:27
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