OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
l3bin.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <math.h>
3 #include <time.h>
4 #include <libgen.h>
5 #include <sys/types.h>
6 
7 #include "netcdf.h"
8 
9 #include "l3bin_input.h"
10 #include <timeutils.h>
11 #include <genutils.h>
12 #include "hdf_bin.h"
13 #include "sensorInfo.h"
14 
15 #include <hdf.h>
16 #include <mfhdf.h>
17 
18 using namespace std;
19 
20 #define MAXNFILES 256
21 
22 #define BYTE unsigned char
23 
24 #define BINCHECK -1
25 
26 #define VERSION "5.13"
27 
28 static instr input;
29 
30 extern "C" int l3bin_input(int argc, char **argv, instr *input, const char* prog, const char* version);
31 
32 
33 int main(int argc, char **argv) {
34  intn i;
35  int retval = 0;
36 
37  int32_t irow;
38  int32_t kbin;
39  int32_t iprod;
40 
41  int32_t ifile, nfiles;
42 
43  int32_t nread;
44  int32_t offset;
45  int64_t offset_out;
46  int32_t offmin;
47  int32_t offmax;
48  int64_t bin_num;
49  int64_t bin_num_out;
50  int32_t row_write;
51  int32_t n_write_total = 0;
52  int32_t n_write = 0;
53  int32_t nprod[MAXNFILES];
54  int32_t ncols;
55  int32_t ncols_out;
56 
58  int16_t composite_outprod_index;
59 
60  int32_t nrows;
61 
62  int32_t reduce_fac;
63 
64  float wgt;
65  float *sort_array[MAXNVDATA];
66 
67  char buf[FILENAME_MAX];
68  char filename_L3[FILENAME_MAX];
69 
70  float *in_sum_buf[MAXNVDATA - 2];
71  float *out_sum_buf[MAXNVDATA - 2];
72  uint8 * in_qual_buf[MAXNFILES + 1], *uint8_buf;
73 
74  int status;
75 
76  char ptime[17];
77  char proc_con[2048];
78 
79  float f32;
80 
81  float minlon;
82  float maxlon;
83  float minlat;
84  float maxlat;
85 
86  float lat;
87 
88  time_t tnow;
89  struct tm *tmnow;
90 
91  FILE *fp, *fp_L3;
92 
93  void insertion_sort(float a[], int length);
94 
96 
97  printf("%s %s (%s %s)\n", "L3BIN", VERSION, __DATE__, __TIME__);
98 
99  if (l3bin_input(argc, argv, &input, "l3bin", VERSION) != 0) {
100  exit(1);
101  }
102 
103  get_time(ptime);
104 
105 
106  strcpy(proc_con, argv[0]);
107  for (i = 1; i < argc; i++) {
108  strcat(proc_con, " ");
109  strcat(proc_con, argv[i]);
110  }
111 
112  reduce_fac = input.reduce_fac;
113  if (reduce_fac != 1 &&
114  reduce_fac != 2 &&
115  reduce_fac != 4 &&
116  reduce_fac != 8 &&
117  reduce_fac != 16) {
118  printf("Reduction factor must be power of 2 less than 16\n");
119  exit(-1);
120  }
121 
122  if (input.loneast <= input.lonwest) {
123  printf("loneast: %f must be greater than lonwest: %f.\n",
124  input.loneast, input.lonwest);
125  exit(-1);
126  }
127 
128  if (input.latnorth <= input.latsouth) {
129  printf("latnorth: %f must be greater than latsouth: %f.\n",
130  input.latnorth, input.latsouth);
131  exit(-1);
132  }
133 
134  /* Get lon/lat limits */
135  minlon = input.lonwest;
136  maxlon = input.loneast;
137  minlat = input.latsouth;
138  maxlat = input.latnorth;
139 
140 
141  /* Determine number of input files */
142  /* ------------------------------- */
143  nfiles = 0;
144 
145  bool isHDF4 = false;
146  bool isHDF5 = false;
147  bool isCDF4 = false;
148 
149  Hdf::hdf_bin * input_binfile[MAXNFILES];
150 
151  /* Single HDF input */
152  /* ---------------- */
153  if (Hishdf(input.infile) == TRUE || H5Fis_hdf5(input.infile) == TRUE) {
154  printf("Single input file\n");
155  nfiles = 1;
156 
157  if (Hishdf(input.infile) == TRUE) {
158  isHDF4 = true;
159  input_binfile[0] = new Hdf::hdf4_bin;
160  }
161 
162  if (H5Fis_hdf5(input.infile) == TRUE) {
163  int ncid;
164  char nam_buf[256];
165  bzero(nam_buf, 256);
166  status = nc_open(input.infile, NC_NOWRITE, &ncid);
167  if (status != NC_NOERR) {
168  isHDF5 = true;
169  input_binfile[0] = new Hdf::hdf5_bin;
170  } else {
171  status = nc_get_att(ncid, NC_GLOBAL, "Mission", nam_buf);
172  if (status != NC_NOERR)
173  status = nc_get_att(ncid, NC_GLOBAL, "mission", nam_buf);
174 
175  if ((status == NC_NOERR) &&
176  ((strcmp(nam_buf, "SAC-D Aquarius") == 0) ||
177  (strcmp(nam_buf, "SMAP") == 0))) {
178  nc_close(ncid);
179  isHDF5 = true;
180  input_binfile[0] = new Hdf::hdf5_bin;
181  } else {
182  isCDF4 = true;
183  input_binfile[0] = new Hdf::cdf4_bin;
184  }
185  }
186  }
187 
188  input_binfile[0]->open(input.infile);
189  nprod[0] = input_binfile[0]->nprod();
190  } else {
191  // Work with a list of files
192  fp = fopen(input.infile, "r");
193  if (fp == NULL) {
194  printf("Input listing file: \"%s\" not found.\n", input.infile);
195  return -1;
196  }
197  while (fgets(buf, 256, fp) != NULL) {
198  // Verify that the file (in this list of files) does indeed exist
199  parse_file_name(buf, filename_L3);
200  fp_L3 = fopen(filename_L3, "r");
201  if (fp_L3 == NULL) {
202  printf("Input file: \"%s\" not found.\n", filename_L3);
203  return -1;
204  }
205  fclose(fp_L3);
206  // Increment count of files (in this list of files)
207  nfiles++;
208  }
209  fclose(fp);
210  printf("%d input files\n", nfiles);
211 
212 
213  /* Open L3 input files */
214  /* ------------------- */
215  fp = fopen(input.infile, "r");
216  for (ifile = 0; ifile < nfiles; ifile++) {
217  fgets(buf, 256, fp);
218  buf[strlen(buf) - 1] = 0;
219 
220  if (ifile == 0) {
221  if (Hishdf(buf) == TRUE) {
222  isHDF4 = true;
223  } else if (H5Fis_hdf5(buf) == TRUE) {
224  int ncid;
225  char nam_buf[256];
226  status = nc_open(buf, NC_NOWRITE, &ncid);
227  if (status != NC_NOERR) {
228  isHDF5 = true;
229  } else {
230  status = nc_get_att(ncid, NC_GLOBAL, "Mission", nam_buf);
231  if (status != NC_NOERR)
232  status = nc_get_att(ncid, NC_GLOBAL, "mission", nam_buf);
233  if (status == NC_NOERR) {
234  if ((strcmp(nam_buf, "SAC-D Aquarius") == 0) ||
235  (strcmp(nam_buf, "SMAP") == 0)) {
236  isHDF5 = true;
237  } else {
238  isCDF4 = true;
239  }
240  } else {
241  isCDF4 = true;
242  }
243  nc_close(ncid);
244  }
245  }
246  }
247  if (isHDF4)
248  input_binfile[ifile] = new Hdf::hdf4_bin;
249  else if (isHDF5)
250  input_binfile[ifile] = new Hdf::hdf5_bin;
251  else if (isCDF4)
252  input_binfile[ifile] = new Hdf::cdf4_bin;
253 
254  printf("%d %s\n", ifile, buf);
255  input_binfile[ifile]->open(buf);
256  nprod[ifile] = input_binfile[ifile]->nprod();
257 
258  //printf("open status: %d\n", status);
259 
260  } /* ifile loop */
261 
262  fclose(fp);
263  }
264 
265 
266  nrows = input_binfile[0]->nrows;
267 
268  /* Generate output product list from 1st input file if DEFAULT */
269  /* ----------------------------------------------------------- */
270  if (strcmp(input.out_parm, ":DEFAULT:") == 0) {
271 
272  // Read input file product list
273  input_binfile[0]->query(input.out_parm);
274 
275  // printf("out_parm: %s\n", &input.out_parm[1]);
276  } else {
277  strcpy(buf, &input.out_parm[1]);
278  buf[strlen(buf) - 1] = 0;
279  for (i = 0; i < (intn) strlen(buf); i++)
280  if (buf[i] == ':') buf[i] = ',';
281  strcpy(input.out_parm, buf);
282  }
283 
284  /* Determine active products */
285  /* ------------------------- */
286  for (ifile = 0; ifile < nfiles; ifile++) {
287  int status = input_binfile[ifile]->read(input.out_parm);
288  if (status == -1) {
289  printf("Not all output products found in input file %d (1-based)\n",
290  ifile);
291  exit(-1);
292  }
293  }
294 
295  /* Check whether Composite product exists in L3 inputfiles */
296  /* ------------------------------------------------------- */
297  if (input.composite_prod[0] != 0) {
298  for (ifile = 0; ifile < nfiles; ifile++) {
299  int tmpNumProducts = input_binfile[ifile]->n_active_prod;
300  int i;
301  for (i = 0; i < tmpNumProducts; i++)
302  if (strcmp(input.composite_prod,
303  input_binfile[0]->getActiveProdName(i)) == 0) break;
304 
306 
307  if (i == tmpNumProducts) {
308  printf("Composite product: \"%s\" not found in L3 dataset %d (1-based)\n",
309  input.composite_prod, ifile);
310  exit(-1);
311  }
312  }
313  }
314 
315 #if 0
316  /* Check whether NDVI product exists in L3 inputfiles */
317  /* -------------------------------------------------- */
318  if (input.composite_prod[0] != 0) {
319  for (ifile = 0; ifile < nfiles; ifile++) {
320  int tmpNumProducts = input_binfile[ifile]->n_active_prod;
321  int i;
322  for (i = 0; i < tmpNumProducts; i++)
323  if (strcmp("ndvi",
324  input_binfile[0]->getActiveProdName(i)) == 0) break;
325 
326  if (i == tmpNumProducts) {
327  printf("NDVI product not found in L3 dataset %d (1-based)\n", ifile);
328  exit(-1);
329  }
330  }
331  }
332 #endif
333 
334  /* Create output file */
335  /* ------------------ */
336  Hdf::hdf_bin *output_binfile;
337 
338  if (getFileFormatName(input.oformat) == NULL) {
339  if (isHDF4) strcpy(input.oformat, "HDF4");
340  if (isHDF5) strcpy(input.oformat, "HDF5");
341  if (isCDF4) strcpy(input.oformat, "netCDF4");
342  } else {
343  // normalize the oformat name
344  strcpy(input.oformat, getFileFormatName(input.oformat));
345  }
346 
347  strcpy(buf, input.ofile);
348 
349  if (strcmp(input.oformat, "HDF4") == 0) {
350  output_binfile = new Hdf::hdf4_bin;
351  if (input.noext == 0) {
352  output_binfile->hasNoext = false;
353  strcat(buf, ".main");
354  } else {
355  output_binfile->hasNoext = true;
356  }
357 
358  } else if (strcmp(input.oformat, "HDF5") == 0) {
359  output_binfile = new Hdf::hdf5_bin;
360  } else if (strcmp(input.oformat, "netCDF4") == 0) {
361  output_binfile = new Hdf::cdf4_bin;
362  output_binfile->deflate = input.deflate;
363  }
364 
365  int tmpNumProducts = input_binfile[0]->n_active_prod;
366  char* tmpProductNames[tmpNumProducts];
367  for (int i = 0; i < tmpNumProducts; i++) {
368  tmpProductNames[i] = (char*) input_binfile[0]->getActiveProdName(i);
369  }
370 
371 
372  output_binfile->setProductList(tmpNumProducts, tmpProductNames);
373  if (input_binfile[0]->has_qual()) {
374  output_binfile->hasQual = true;
375  }
376  output_binfile->create(buf, nrows / reduce_fac);
377 
378 
379  /* Allocate I/O buffers */
380  /* -------------------- */
381  ncols = 2 * nrows;
382  ncols_out = 2 * nrows / reduce_fac;
383 
384  for (iprod = 0; iprod < nprod[0]; iprod++) {
385  if (input_binfile[0]->active_data_prod[iprod] == true) {
386  in_sum_buf[iprod] = (float *) calloc(ncols, 2 * sizeof (float));
387  out_sum_buf[iprod] = (float *) calloc(ncols_out, 2 * sizeof (float));
388 
389  if (input.composite_prod[0] != 0) {
390  if (strcmp(input.composite_prod,
391  input_binfile[0]->getActiveProdName(iprod)) == 0)
392  composite_outprod_index = iprod;
393  }
394  }
395  } /* iprod loop */
396 
397 
398  /* Allocate quality buffer */
399  /* ----------------------- */
400  for (ifile = 0; ifile < nfiles; ifile++) {
401  in_qual_buf[ifile] = (uint8 *) calloc(ncols, sizeof (uint8));
402  }
403  in_qual_buf[nfiles] = (uint8 *) calloc(ncols, sizeof (uint8));
404 
405  uint8_buf = (uint8 *) calloc(ncols, sizeof (uint8));
406 
407 
408  /* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
409  /* For each scan ... (Main Loop) */
410  /* ----------------------------- */
411  for (irow = 0; irow < nrows; irow++) {
412 
413  if ((irow % 500) == 0) {
414  time(&tnow);
415  tmnow = localtime(&tnow);
416  printf("irow:%6d of %8d %s", irow, nrows, asctime(tmnow));
417  }
418 
419  int32_t max_out_kbin;
420  if ((irow % reduce_fac) == 0)
421  max_out_kbin = -1;
422 
423  // Get basebin and numbin for this input row
424  int64_t basebin = input_binfile[0]->get_basebin(irow);
425  int64_t numbin = input_binfile[0]->get_numbin(irow);
426 
427  double ratio = 1.0;
428  if (reduce_fac > 1)
429  ratio = 1.0 * input_binfile[0]->get_numbin(irow) /
430  output_binfile->get_numbin(irow / reduce_fac);
431 
432 
433  // If median allocate and initialize storage
434  if (input.median) {
435  for (iprod = 0; iprod < nprod[0]; iprod++) {
436  if (input_binfile[0]->active_data_prod[iprod] == true) {
437  sort_array[iprod] =
438  (float *) calloc(nfiles*numbin, sizeof (float));
439  for (int64_t i = 0; i < nfiles * numbin; i++) sort_array[iprod][i] = -999;
440  }
441  }
442  }
443 
444 
445  /* Clear output binlist, sum and quality buffers */
446  /* --------------------------------------------- */
447  if ((irow % reduce_fac) == 0) {
448  row_write = 0;
449 
450  output_binfile->clear_binlist();
451  for (iprod = 0; iprod < nprod[0]; iprod++) {
452  if (input_binfile[0]->active_data_prod[iprod] == true) {
453  memset(&out_sum_buf[iprod][0], 0, ncols_out * 2 * sizeof (float));
454  }
455  } /* iprod loop */
456 
457  memset(in_qual_buf[nfiles], 255, ncols);
458  for (ifile = 0; ifile < nfiles; ifile++) {
459  memset(in_qual_buf[ifile], 255, ncols);
460  }
461  }
462 
463  /* Get bin & qual info */
464  /* ------------------- */
465  for (ifile = 0; ifile < nfiles; ifile++) {
466  input_binfile[ifile]->readBinIndex(irow);
467 
468  int ext = input_binfile[ifile]->get_ext();
469 
470  /* Read BinList */
471  /* ------------ */
472  // Get current binlist pointer
473  int32_t list_ptr = input_binfile[ifile]->get_list_ptr();
474 
475  if (ext > 0) {
476  nread = input_binfile[ifile]->readBinList(ext);
477 
478  if (nread == -1) {
479  printf("Unable to read bin numbers...: %d\n", ext);
480  }
481  }
482 
483  /* Read quality if vdata exists */
484  /* Note: in_qual_buf is "uncompressed" along row */
485  if (input_binfile[ifile]->has_qual()) {
486  if ((irow % reduce_fac) == 0) {
487  int32_t ext_qual = ext;
488  for (int32_t i = 0; i < reduce_fac; i++) {
489  switch (i) {
490  case 0:
491  input_binfile[ifile]->readQual(uint8_buf, ext_qual, list_ptr);
492  break;
493  default:
494  input_binfile[ifile]->readBinIndex(irow + i);
495  ext_qual = input_binfile[ifile]->get_ext();
496  if (ext_qual != 0) {
497  nread = input_binfile[ifile]->readBinList(ext_qual);
498  input_binfile[ifile]->readQual(uint8_buf, ext_qual);
499  }
500  } // switch
501 
502  // Update in_qual_buf for (irow+i)
503  if (ext_qual != 0) {
504  int64_t basebin = input_binfile[0]->get_basebin(irow + i);
505 
506  double ratio = 1.0 * input_binfile[0]->get_numbin(irow + i) /
507  output_binfile->get_numbin(irow / reduce_fac);
508 
509  for (kbin = 0; kbin < ext_qual; kbin++) {
510  bin_num = input_binfile[ifile]->get_bin_num(kbin);
511  offset = bin_num - basebin;
512  if (offset < 0) {
513  cout << "bin_num - basebin is negative for ifile: " << ifile;
514  cout << " irow: " << irow << " kbin: " << kbin << endl;
515  cout << "bin_num: " << bin_num << " basebin: " << basebin
516  << endl;
517  exit(1);
518  }
519 
520  int32_t j = reduce_fac * (int32_t) ((offset / ratio) + 0.0);
521  if ((j < ncols) && (in_qual_buf[ifile][j] > uint8_buf[kbin])) {
522  in_qual_buf[ifile][j] = uint8_buf[kbin];
523  }
524  }
525  } // ext_qual != 0
526  } // i (reduce_fac) loop
527 
528  // Reset to current row if reduce_fac != 1
529  if (reduce_fac != 1) {
530  input_binfile[ifile]->readBinIndex(irow);
531  nread = input_binfile[ifile]->readBinList(ext, list_ptr);
532  }
533 
534  } // if ( (irow % reduce_fac) == 0)
535  } // if ( input_binfile[ifile]->has_qual())
536  } /* ifile loop */
537 
538 
539  /* Find best quality */
540  /* ----------------- */
541  for (kbin = 0; kbin < ncols; kbin++) {
542  for (ifile = 0; ifile < nfiles; ifile++) {
543  int32_t j = reduce_fac * (int32_t) ((kbin / ratio) + 0.0);
544  if (j < ncols) {
545  if (in_qual_buf[ifile][j] < in_qual_buf[nfiles][j])
546  in_qual_buf[nfiles][j] = in_qual_buf[ifile][j];
547  }
548  }
549  }
550 
551  /* For each file ... */
552  /* ----------------- */
553  for (ifile = 0; ifile < nfiles; ifile++) {
554 
555  int64_t beg = input_binfile[ifile]->get_beg();
556  int ext = input_binfile[ifile]->get_ext();
557 
558  /* ========== If row has data ... ========== */
559  /* ----------------------------------------- */
560  if (beg != 0) {
561  /* printf("row has data: %d\n", irow);*/
562 
563  /* Determine lon kbin limits */
564  /* ------------------------- */
565  offmin =
566  (int32_t) ((minlon + 180) * (numbin / 360.0) + 0.5);
567  offmax =
568  (int32_t) ((maxlon + 180) * (numbin / 360.0) + 0.5);
569 
570 
571  /* Get data values (sum, sum_sq) for each filled bin in row */
572  /* -------------------------------------------------------- */
573  int nbins_to_read = ext;
574  for (iprod = 0; iprod < nprod[ifile]; iprod++) {
575  if (input_binfile[0]->active_data_prod[iprod] == true) {
576 
577  input_binfile[ifile]->readSums(&in_sum_buf[iprod][0],
578  nbins_to_read, iprod);
579  }
580  } /* iprod loop */
581  if (isHDF5 || isCDF4)
582  input_binfile[ifile]->setDataPtr(nbins_to_read);
583  // row_write = 1;
584 
585 
586  /* Skip row if not between minlat & maxlat */
587  lat = ((irow + 0.5) / nrows) * 180.0 - 90.0;
588  if (lat < minlat || lat > maxlat) {
589  // row_write = 0;
590  continue;
591  }
592 
593  /* Fill output buffers with input bin data */
594  /* --------------------------------------- */
595  for (kbin = 0; kbin < ext; kbin++) {
596 
597  /* Store bin number */
598  /* ---------------- */
599  bin_num = input_binfile[ifile]->get_bin_num(kbin);
600  offset = bin_num - basebin;
601  if (offset < 0) {
602  cout << "bin_num - basebin is negative for ifile: " << ifile;
603  cout << " irow: " << irow << " kbin: " << kbin << endl;
604  cout << "bin_num: " << bin_num << " basebin: " << basebin << endl;
605  exit(1);
606  }
607 
608  /* If bin outside lon range then skip */
609  /* ---------------------------------- */
610  if (offset < offmin || offset > offmax)
611  continue;
612 
613  float weights = input_binfile[ifile]->get_weights(kbin);
614  float time_rec = input_binfile[ifile]->get_time_rec(kbin);
615 
616  /* Skip if not good enough */
617  /* ----------------------- */
618  int32_t j = reduce_fac * (int32_t) ((offset / ratio) + 0.0);
619  if (j < ncols) {
620  if (in_qual_buf[ifile][j] > in_qual_buf[nfiles][j])
621  continue;
622  }
623 
624  // Assign output offset & bin number
625  if (reduce_fac > 1) {
626  offset_out = (int64_t) ((offset / ratio) + 0.0);
627  // if ( offset_out == ncols_out)
628  // offset_out = ncols_out - 1;
629  bin_num_out =
630  offset_out + output_binfile->get_basebin(irow / reduce_fac);
631  } else {
632  offset_out = offset;
633  bin_num_out = bin_num;
634  }
635 
636  if (offset_out >= ncols_out) {
637  printf("Bad write to BINLIST: %d %d %d %ld\n",
638  ifile, irow, ncols_out, (long int) offset_out);
639  exit(1);
640  }
641 
642  if (offset_out > max_out_kbin)
643  max_out_kbin = offset_out;
644 
645  output_binfile->set_bin_num(offset_out, bin_num_out);
646  row_write = 1;
647 
648  /* Sum & store number of observations,nscenes */
649  /* ------------------------------------------ */
650  int nobs, nscenes;
651  if (input.composite_prod[0] == 0) {
652  nobs = input_binfile[ifile]->get_nobs(kbin);
653  output_binfile->inc_nobs(offset_out, nobs);
654  nscenes = input_binfile[ifile]->get_nscenes(kbin);
655  output_binfile->inc_nscenes(offset_out, nscenes);
656 
657  /* Sum & store weights */
658  /* ------------------- */
659  if (input.unit_wgt || input.median) {
660  output_binfile->set_weights(offset_out, 1);
661  } else {
662  output_binfile->inc_weights(offset_out, weights);
663  }
664 
665  output_binfile->inc_time_rec(offset_out, time_rec);
666  } else {
667  if (output_binfile->get_nobs(offset_out) == 0) {
668  output_binfile->inc_nobs(offset_out, 1);
669  nscenes = input_binfile[ifile]->get_nscenes(kbin);
670  output_binfile->inc_nscenes(offset_out, nscenes);
671  output_binfile->set_weights(offset_out, 1);
672  } else {
673  f32 = in_sum_buf[composite_prod_index[ifile]][2 * kbin];
674  if (strcmp(input.composite_scheme, "max") == 0) {
675  if (f32 < out_sum_buf[composite_outprod_index][2 * offset_out]) continue;
676  } else {
677  if (f32 > out_sum_buf[composite_outprod_index][2 * offset_out]) continue;
678  }
679  }
680  }
681 
682  /* Product loop */
683  /* ------------ */
684  for (iprod = 0; iprod < nprod[ifile]; iprod++) {
685  if (input_binfile[ifile]->active_data_prod[iprod] == true) {
686 
687  if (input.unit_wgt) {
688  wgt = weights;
689  f32 = in_sum_buf[iprod][2 * kbin];
690  out_sum_buf[iprod][2 * offset_out] += f32 / wgt;
691  out_sum_buf[iprod][2 * offset_out + 1] += (f32 / wgt)*(f32 / wgt);
692  } else if (input.median) {
693  wgt = weights;
694  f32 = in_sum_buf[iprod][2 * kbin];
695  sort_array[iprod][nfiles * offset + ifile] = f32 / wgt;
696  } else if (input.composite_prod[0] != 0) {
697  f32 = in_sum_buf[iprod][2 * kbin];
698  out_sum_buf[iprod][2 * offset_out] = f32;
699  f32 = in_sum_buf[iprod][2 * kbin + 1];
700  out_sum_buf[iprod][2 * offset_out + 1] = f32;
701  } else {
702  /* Add new sum to accumulated sum & sum2 */
703  /* ------------------------------------- */
704  f32 = in_sum_buf[iprod][2 * kbin];
705  out_sum_buf[iprod][2 * offset_out] += f32;
706  f32 = in_sum_buf[iprod][2 * kbin + 1];
707  out_sum_buf[iprod][2 * offset_out + 1] += f32;
708  } /* input.unit_wgt */
709 
710  } /* active_prod */
711 
712  } /* product loop */
713 
714  } /* kbin loop */
715 
716  } /* ========== if row has data ========== */
717 
718  } /* ifile loop */
719 
720 
721  // Compute median
722  if (row_write && input.median) {
723  for (kbin = 0; kbin < numbin; kbin++) {
724  for (iprod = 0; iprod < nprod[0]; iprod++) {
725  if (input_binfile[0]->active_data_prod[iprod] == true) {
726  float *sort_buf = (float *) calloc(nfiles, sizeof (float));
727  int nsort = 0;
728  for (ifile = 0; ifile < nfiles; ifile++) {
729  f32 = sort_array[iprod][nfiles * kbin + ifile];
730  if (f32 != -999)
731  sort_buf[nsort++] = sort_array[iprod][nfiles * kbin + ifile];
732  }
733  // Call insertion sort
734  if (nsort > 0) {
735  insertion_sort(sort_buf, nsort);
736  out_sum_buf[iprod][2 * (kbin / reduce_fac)] = sort_buf[nsort / 2];
737  out_sum_buf[iprod][2 * (kbin / reduce_fac) + 1] = 1.0;
738  free(sort_buf);
739  }
740  }
741  }
742  }
743  }
744 
745  /* Write output vdatas */
746  /* ------------------- */
747  if (row_write && (irow % reduce_fac) == reduce_fac - 1) {
748 
749  n_write = 0;
750  for (kbin = 0; kbin <= max_out_kbin; kbin++) {
751 
752  int64_t bin_num = output_binfile->get_bin_num(kbin);
753  if (bin_num != 0) {
754 
755  /* Loop over data products */
756  /* ----------------------- */
757  for (iprod = 0; iprod < nprod[0]; iprod++) {
758  if (input_binfile[0]->active_data_prod[iprod] == true) {
759 
760  /* Remove "blank" bin records */
761  /* -------------------------- */
762  if (n_write != kbin)
763  memcpy(&out_sum_buf[iprod][2 * n_write],
764  &out_sum_buf[iprod][2 * kbin], 8);
765  }
766  } /* iprod loop */
767 
768  /* Remove "blank" bin records */
769  /* -------------------------- */
770  if (n_write != kbin)
771  output_binfile->copy_binlist(kbin, n_write);
772 
773  n_write++;
774  n_write_total++;
775 
776  } /* bin_num != 0 */
777  } /* kbin loop */
778 
779  /* Write BinList & Data Products */
780  /* ----------------------------- */
781  if (n_write > 0) {
782  output_binfile->writeBinList(n_write);
783 
784  for (iprod = 0; iprod < nprod[0]; iprod++) {
785  if (input_binfile[0]->active_data_prod[iprod] == true) {
786 
787  input_binfile[0]->get_prodname(iprod, buf);
788  output_binfile->writeSums(&out_sum_buf[iprod][0], n_write, buf);
789  }
790  }
791  if (strcmp(input.oformat, "HDF5") == 0 ||
792  strcmp(input.oformat, "netCDF4") == 0)
793  output_binfile->incNumRec(n_write);
794  // if ( isHDF5 || isCDF4) output_binfile->incNumRec( n_write);
795  }
796 
797  /* Write to Quality Vdata */
798  /* ---------------------- */
799  i = 0;
800  int64_t nbin;
801  if (irow < nrows / 2)
802  nbin = input_binfile[0]->get_numbin(irow);
803  else
804  nbin = input_binfile[0]->get_numbin(irow - (reduce_fac / 2));
805 
806  offmin =
807  (int32_t) ((minlon + 180) * (nbin / 360.0) + 0.5);
808  offmax =
809  (int32_t) ((maxlon + 180) * (nbin / 360.0) + 0.5);
810 
811  for (kbin = 0; kbin < ncols; kbin++) {
812  if (kbin < offmin || kbin > offmax) continue;
813  if (in_qual_buf[nfiles][kbin] != 255) {
814  in_qual_buf[nfiles][i++] = in_qual_buf[nfiles][kbin];
815  }
816  }
817 
818  if (input_binfile[0]->has_qual()) {
819  if ((i - n_write) > 2) {
820  cout << "Problem with Quality write: irow: " << irow << " i: " <<
821  i << " n_write: " << n_write << " " << max_out_kbin << endl;
822  exit(1);
823  }
824 
825  output_binfile->writeQual(in_qual_buf[nfiles], n_write);
826  } // has_qual() == true
827  } /* row_write = 1 */
828 
829 
830  // if median free storage arrays
831  if (input.median) {
832  for (iprod = 0; iprod < nprod[0]; iprod++) {
833  if (input_binfile[0]->active_data_prod[iprod] == true) {
834  free(sort_array[iprod]);
835  }
836  }
837  }
838 
839  } /* irow loop (Main loop) */
840  /* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
841 
842  for (ifile = 0; ifile <= nfiles; ifile++) {
843  free(in_qual_buf[ifile]);
844  }
845  free(uint8_buf);
846 
847  for (iprod = 0; iprod < nprod[0]; iprod++) {
848  if (input_binfile[0]->active_data_prod[iprod] == true) {
849  free(in_sum_buf[iprod]);
850  free(out_sum_buf[iprod]);
851  }
852  } /* iprod loop */
853 
854 
855  if(output_binfile->n_data_records > 0) {
856 
857  // Copy metadata from input to output binfile
858  output_binfile->copymeta(nfiles, &input_binfile[0]);
859  int sensorID = sensorName2SensorId(output_binfile->meta_l3b.sensor_name);
860  output_binfile->meta_l3b.sensorID = sensorID;
861 
862  // put in a fix to the missing Mission
863  if (strlen(output_binfile->meta_l3b.mission) == 0) {
864  if (sensorID != -1) {
866  }
867  }
868 
869  strcpy(buf, input.ofile);
870  strcpy(output_binfile->meta_l3b.product_name, buf);
871  strcpy(output_binfile->meta_l3b.pversion, input.pversion);
872  strcpy(output_binfile->meta_l3b.soft_name, "L3BIN");
873  strcpy(output_binfile->meta_l3b.soft_ver, VERSION);
874  if (isCDF4 == 1) {
875  // 1994-11-05T13:15:30Z
876  strcpy(output_binfile->meta_l3b.ptime, unix2isodate(tnow, 'G'));
877  } else {
878  // yyyydddhhmmssmmm
879  strcpy(output_binfile->meta_l3b.ptime, ydhmsf(tnow, 'L'));
880  }
881  strcpy(output_binfile->meta_l3b.ptime, ptime);
882  strcpy(output_binfile->meta_l3b.proc_con, proc_con);
883  strcpy(output_binfile->meta_l3b.input_parms, input.parms);
884 
885  char buf2[LG_ATTRSZ];
886  if (Hishdf(input.infile) == TRUE || H5Fis_hdf5(input.infile) == TRUE) {
887  strcpy(output_binfile->meta_l3b.infiles, basename(input.infile));
888 
889  } else {
890 
891  fp = fopen(input.infile, "r");
892  buf2[0] = 0;
893  for (ifile = 0; ifile < nfiles; ifile++) {
894  fgets(buf, 256, fp);
895  buf[strlen(buf) - 1] = 0;
896 
897  strcat(buf2, basename(buf));
898  strcat(buf2, ",");
899  } /* ifile loop */
900  fclose(fp);
901  buf2[strlen(buf2) - 1] = 0;
902  strncpy(output_binfile->meta_l3b.infiles, buf2, LG_ATTRSZ - 1);
903  }
904  } else {
905  retval = 110;
906  printf("No valid data was binned\n");
907  }
908 
909 
910  for (ifile = 0; ifile < nfiles; ifile++) {
911  input_binfile[ifile]->close();
912  delete input_binfile[ifile];
913  }
914 
915  output_binfile->close();
916  delete output_binfile;
917 
918  if(retval == 110) {
919  string cmd = "rm -f ";
920  cmd += input.ofile;
921  system(cmd.c_str());
922  }
923 
924  printf("Done\n");
925  return retval;
926 }
927 
928 
929 /* Copyright (c) 2009 the authors listed at the following URL, and/or
930 the authors of referenced articles or incorporated external code:
931 http://en.literateprograms.org/Insertion_sort_(C)?action=history&offset=20081205204844
932 
933 Permission is hereby granted, free of charge, to any person obtaining
934 a copy of this software and associated documentation files (the
935 "Software"), to deal in the Software without restriction, including
936 without limitation the rights to use, copy, modify, merge, publish,
937 distribute, sublicense, and/or sell copies of the Software, and to
938 permit persons to whom the Software is furnished to do so, subject to
939 the following conditions:
940 
941 The above copyright notice and this permission notice shall be
942 included in all copies or substantial portions of the Software.
943 
944 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
945 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
946 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
947 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
948 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
949 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
950 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
951 
952 Retrieved from: http://en.literateprograms.org/Insertion_sort_(C)?oldid=15530
953  */
954 
955 /* Sort an array of integers */
956 void insertion_sort(float a[], int length) {
957  int i;
958  for (i = 0; i < length; i++) {
959  /* Insert a[i] into the sorted sublist */
960  int j;
961  float v = a[i];
962 
963  for (j = i - 1; j >= 0; j--) {
964  if (a[j] <= v) break;
965  a[j + 1] = a[j];
966  }
967  a[j + 1] = v;
968 
969  }
970 }
971 
char * ydhmsf(double dtime, char zone)
Definition: ydhmsf.c:12
virtual int64_t get_bin_num(int kbin)=0
virtual int get_nobs(int kbin)=0
virtual int32_t get_list_ptr()
Definition: hdf_bin.h:112
char product_name[SM_ATTRSZ]
Definition: meta_l3b.h:16
uint8_t isHDF5
Definition: get_l3m.c:8
virtual const char * getActiveProdName(int prodNum) const
Definition: bin_io.cpp:2858
int j
Definition: decode_rs.h:73
int status
Definition: l1_czcs_hdf.c:32
virtual int readSums(float *sums, int32_t nbins_to_read, int iprod)=0
virtual void setProductList(int numProducts, char *prodNames[])
Definition: bin_io.cpp:2807
char infiles[LG_ATTRSZ]
Definition: meta_l3b.h:38
#define VERSION
Definition: l3bin.cpp:26
virtual int inc_weights(int offset, float weights)=0
#define NULL
Definition: decode_rs.h:63
char mission[SM_ATTRSZ]
Definition: meta_l3b.h:21
virtual int create(const char *l3b_filename, int32_t nrows)=0
char input_parms[LG_ATTRSZ]
Definition: meta_l3b.h:36
int sensorName2SensorId(const char *name)
Definition: sensorInfo.c:268
const char * sensorId2PlatformName(int sensorId)
Definition: sensorInfo.c:226
#define TRUE
Definition: rice.h:165
float * lat
float tm[MODELMAX]
virtual int close()=0
virtual int64_t get_numbin(int irow)
Definition: hdf_bin.h:72
void insertion_sort(float a[], int length)
Definition: l3bin.cpp:956
virtual int64_t get_beg()=0
virtual int copy_binlist(int src, int dest)=0
void bzero()
virtual int get_nscenes(int kbin)=0
int32_t nobs
Definition: atrem_cor.h:93
instr * input
virtual int readQual(uint8_t *qual, int32_t nbins_to_read)=0
int setlinebuf(FILE *stream)
virtual int get_ext()=0
virtual int query()
Definition: bin_io.cpp:2764
bool hasNoext
Definition: hdf_bin.h:127
int16_t composite_prod_index[MAXNFILES]
Definition: l2bin.cpp:101
virtual int set_weights(int offset, float weights)=0
virtual int copymeta(int32_t nfiles, Hdf::hdf_bin *input_binfile[])
Definition: bin_io.cpp:2874
void get_time(char *pr_time)
Definition: get_time.c:28
const char * getFileFormatName(const char *str)
virtual int64_t get_basebin(int irow)
Definition: hdf_bin.h:76
int32_t n_active_prod
Definition: hdf_bin.h:123
int32_t nrows
Definition: hdf_bin.h:119
int l3bin_input(int argc, char **argv, instr *input, const char *prog, const char *version)
Definition: l3bin_input.c:336
virtual int32_t nprod()
Definition: hdf_bin.h:108
bool hasQual
Definition: hdf_bin.h:126
meta_l3bType meta_l3b
Definition: hdf_bin.h:131
int32 nrows
virtual int open(const char *l3b_filename)=0
virtual int incNumRec(int n_write)=0
virtual int writeSums(float *sums, int32_t nbins_to_write, const char *prodname)=0
virtual int setDataPtr(int nbins_to_read)=0
virtual float get_time_rec(int kbin)=0
#define MAXNFILES
Definition: l3bin.cpp:20
void parse_file_name(const char *inpath, char *outpath)
int sensorID
Definition: meta_l3b.h:18
#define basename(s)
Definition: l0chunk_modis.c:29
virtual int get_prodname(int iprod, char *prodname)
Definition: bin_io.cpp:2802
char proc_con[MD_ATTRSZ]
Definition: meta_l3b.h:35
int16 ncols[121]
virtual int inc_nobs(int offset, int nobs)=0
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned resolving bug report Corrected to filter out Aqua attitude records with missing status helping resolve bug MOD_PR03 will still correctly write scan and pixel data that does not depend upon the start time
Definition: HISTORY.txt:248
char soft_ver[SM_ATTRSZ]
Definition: meta_l3b.h:33
char * unix2isodate(double dtime, char zone)
Definition: unix2isodate.c:10
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second and prod_ix set to PARAM_TYPE_INT name_prefix is compared with the beginning of the product name If name_suffix is not empty the it must match the end of the product name The characters right after the prefix are read as an integer and prod_ix is set to that number strncpy(l2prod->name_prefix, "myprod", UNITLEN)
int16_t * nscenes
Definition: l2bin.cpp:86
#define MAXNVDATA
virtual int set_bin_num(int offset, int64_t bin_num)=0
char ptime[SM_ATTRSZ]
Definition: meta_l3b.h:34
virtual int clear_binlist()=0
virtual int read(char *product_list)
Definition: bin_io.cpp:358
int32_t n_data_records
Definition: hdf_bin.h:122
char pversion[SM_ATTRSZ]
Definition: meta_l3b.h:30
virtual int readBinIndex(int row_num_to_read)=0
l2prod offset
virtual int readBinList(int nbins_to_read)=0
virtual int inc_time_rec(int offset, float time_rec)=0
int main(int argc, char **argv)
Definition: l3bin.cpp:33
char sensor_name[SM_ATTRSZ]
Definition: meta_l3b.h:19
int i
Definition: decode_rs.h:71
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
virtual int writeQual(uint8_t *qual, int32_t nbins_to_write)=0
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to a
Definition: HISTORY.txt:424
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:97
virtual int inc_nscenes(int offset, int nscenes)=0
#define LG_ATTRSZ
Definition: meta_l3b.h:9
char soft_name[SM_ATTRSZ]
Definition: meta_l3b.h:32
uint32_t deflate
Definition: hdf_bin.h:129
float32 f32
Definition: l2bin.cpp:104
version
Definition: setup.py:15
virtual float get_weights(int kbin)=0
virtual int writeBinList(int32_t nbins_to_write)=0