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