OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
l2bin.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdint.h>
3 #include <libgen.h>
4 #include <float.h>
5 #include <math.h>
6 #include <stdio.h>
7 #include <string.h>
8 #include <sys/types.h>
9 #include <time.h>
10 #include <iostream>
11 #include <vector>
12 #include <algorithm>
13 #include <string>
14 #include <netcdf>
15 #include <boost/algorithm/string.hpp>
16 #include <boost/geometry.hpp>
17 #include <boost/geometry/geometries/point_xy.hpp>
18 #include <boost/geometry/geometries/polygon.hpp>
19 #include <boost/geometry/geometries/box.hpp>
20 #include <boost/foreach.hpp>
21 
22 #include <meta_l3b.h>
23 #include <seaproto.h>
24 #include <readL2scan.h>
25 #include <timeutils.h>
26 #include <genutils.h>
27 #include <setupflags.h>
28 #include <sensorInfo.h>
29 #include <sensorDefs.h>
30 #include <nc4utils.h>
31 #include <ncdfbin_utils.h>
32 
33 #include <hdf.h>
34 #include <mfhdf.h>
35 
36 #include "dataday.h"
37 #include "l2bin_input.h"
38 #include "L3Shape.h"
39 #include "L3ShapeIsine.h"
40 
41 using namespace std;
42 using namespace netCDF;
43 using namespace netCDF::exceptions;
44 
45 namespace bg = boost::geometry;
46 
47 
48 //#define DEBUG_PRINT
49 //#define DEBUG_PRINT2
50 
51 #ifdef DEBUG_PRINT
52 bool enableDebugPrint=true;
53 #endif
54 
55 
56 #define PI 3.141592653589793
57 #define MTILT_DIMS_2 20
58 #define LTILT_DIMS_2 2
59 #define MAXALLOCPERBIN 20
60 #define EARTH_RADIUS 6378.14
61 
62 static const int max_l3b_products = 400;
63 
64 /* Global variables */
65 static instr input;
66 static l2_prod l2_str[MAXNFILES];
67 
68 // bin file shape
69 static int32_t nrows = -1;
70 static l3::L3Shape *shape;
71 static int64_t *basebin;
72 
73 // global variables for the group
74 int32_t n_allocperbin; // number of obs to allocate to make another chunk
75 int16_t *allocated_space; // number of obs allocated for each bin
76 static int16_t *nobs; // number of observations in each bin
77 
78 static float32 **data_values;
79 static double **data_areas; // in fractions of a bin
80 static int16_t **file_index;
81 static uint8_t **data_quality;
82 static float64 **time_value;
83 static double time_tai93; // TAI 93 time for the current L2 line
84 
85 int32_t *bin_flag;
86 int16_t *tilt, *qual, *nscenes, *lastfile;
87 
88 int32_t n_bins_in_group; // number of bins in group
89 int32_t n_rows_in_group = -1; // number of rows in group
90 int32_t krow; // first bin row number of group
91 
92 // other globals
93 int32_t tiltstate = 0;
94 int32_t l3b_nprod; // number of products in bin file
95 static char buf[LG_ATTRSZ];
96 
98 
103 
104 float32 f32;
105 vector<int32_t> thirdDimId;
106 vector<float32> min_value;
107 
108 typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree>> Point_t;
109 typedef bg::model::polygon<Point_t> Polygon_t;
110 typedef bg::model::box<Point_t> Box_t;
111 
112 #define VERSION "7.0.3"
113 #define PROGRAM "l2bin"
114 
115 /*
116  *
117  Revision 7.0.1 2021-04-22
118  Read bandsPerPixel from openL2
119  J. Gales
120 
121  Revision 7.0.0 2021-01-29
122  Add support the 3D L2 data products
123  J. Gales
124 
125  Revision 4.4.1 2018-09-12
126  Remove support for HEALPIX
127  S. Bailey, G. Fireman
128 
129  Revision 4.4.0 11/29/2016
130  Added clo library to replace command line processing functions
131  R. Healy
132 
133  Revision 4.3.0 11/28/2016
134  Added in xml code to read product min/max value from product.xml file. Changed version to match l2bin64
135  R. Healy
136 
137  Revision 4.1.2 08/04/2016
138  Add composite_scheme parameter for land composite binning.
139  J. Gales
140 
141  Revision 4.1.1 07/02/2016
142  Add support for land composite binning.
143  J. Gales
144 
145  Revision 4.1.0 06/16/2016
146  Allocate and clear dolat/dolon arrays in dataday code to avoid time bombs.
147  J. Gales
148 
149  Revision 4.0.9 01/06/2015
150  Fixed issue introduced with 4.0.5 with default sday/eday which resulted
151  in a change in behavior
152  * NOTE: if this program is still in use in 2038...may need a tweak :)
153  Modified test for "regional" prodtype to be case insensitive
154  S. Bailey
155 
156  Revision 4.0.8 11/30/2015
157  Fixed opening group id for multiple files.
158  D. Shea
159 
160  Revision 4.0.7 10/26/2015
161  Fix reading of year,day,time values for time_rec
162  D. Shea
163 
164  Revision 4.0.6 10/26/2015
165  Accumulate time_rec sums at double precision, store at float
166  nature of crossing the pole.
167  J. Gales
168 
169  Revision 4.0.5 10/01/2015
170  Added logic to handle orbit based files that cross the dateline by
171  nature of crossing the pole.
172  S. Bailey
173 
174  Revision 4.0.4 09/01/15
175  Add support for time_rec field in the binlist records
176  J. Gales
177 
178  *
179  * Revision 4.0
180  * Replaced dateline crossing (dataday) with functions
181  * based on Norman Kuring's brstobins7
182  * Replaced diffday calculations using unix time
183  * Added sensorID dependent equator crossing times
184  * R. Healy
185 
186  Revision 3.1.2 03/10/14
187  Remove paths from input filenames
188  J. Gales
189 
190  Revision 3.1.1 03/06/14
191  Add support for binlist and bindata compression
192  J. Gales
193 
194  Revision 3.1.0 03/05/14
195  Add support for CF metadata
196  J. Gales
197 
198  Revision 2.5.1 05/13/13
199  Add support for HEALPIX midaverage (despeckling)
200  J. Gales
201 
202  Revision 2.5.0 04/08/13
203  Add support for NETCDF4
204  J. Gales
205 
206  Revision 2.4.9 10/14/12
207  Put MERIS p1hr dataday parameter back to 19
208  J. Gales
209 
210  Revision 2.4.8 10/14/12
211  Adjust MERIS dataday parameters
212  J. Gales
213 
214  Revision 2.4.7 10/10/12
215  Switch MERIS from TERRA-like to SEAWIFS-like
216  J. Gales
217 
218  Revision 2.4.6 07/30/12
219  Incorporate cache allocation fix in readL2scan
220  J. Gales
221 
222  Revision 2.4.5 09/12/11
223  Put back "platformInformation" & "instrumentInformation" metadata
224  J. Gales
225 
226  Revision 2.4.4 06/27/11
227  Add MERIS to "TERRA-like" dataday
228  J. Gales
229 
230  Revision 2.4.3 08/05/10
231  Exit on non-existent parm/suite file
232  J. Gales
233 
234  Revision 2.4.2 02/04/10
235  Change parameter for cde=5 (dataday) to 0.92
236  J. Gales
237 
238  Revision 2.4.0 02/04/10
239  Add caching for better i/o performance.
240  D. Shea
241 
242  Revision 2.3.2 08/31/09
243  Determine bad value by testing for nan set in readL2 function.
244  J. Gales
245 
246  Revision 2.3.2 08/25/09
247  Add suite parameter
248  Modify l2bin_input code to handle suite defaults
249  J. Gales
250 
251  Revision 2.3.1 08/11/09
252  Add pversion parameter
253  Remove "Replacement Flag" metadata
254  J. Gales
255 
256  Revision 2.3.0 11/15/07
257  Add ability to "flag" on bad value defined in SDS metadata
258  J. Gales
259 
260  Revision 2.2.3 05/02/07
261  Check that no more than one delimiter is specified.
262  J. Gales
263 
264  Revision 2.2.2 04/10/07
265  Tweak ssec limit for TERRA night
266  previous day/no scancross granules to 10.1*60*60
267  (T2000138100000.L2_LAC_SST)
268  J. Gales
269 
270  Revision 2.2.1 04/06/07
271  Set bins with -32767 l2 pixel values to data_quality = 4
272  J. Gales
273 
274  Revision 2.2.0 04/03/07
275  Only throw out bad geolocated scans, not entire granules
276  J. Gales
277 
278  Revision 2.1.9 04/02/07
279  Fix memory leak when no good bins are left after qual check.
280  J. Gales
281 
282  Revision 2.1.8 03/27/07
283  Add check for non-navigatable file
284  J. Gales
285 
286  Revision 2.1.7 03/02/07
287  Remove previous revision (2.1.6)
288  Use spline check in libl2 to catch bad lon/lat
289  Skip pixels with bad lon/lat
290  J. Gales
291 
292  Revision 2.1.6 03/01/07
293  Skip swath scan lines with bad lon/lat
294  J. Gales
295 
296  Revision 2.1.5 02/09/07
297  Disable INTERP parameter
298  Fix missing and incorrect entries in Input Parameters Attribute
299  J. Gales
300 
301  Revision 2.1.4 12/12/06
302  Add ',' and ' ' as product delimiters.
303  Add '=' as min value delimiters.
304  Trap bad minimum values.
305  J. Gales
306 
307  Revision 2.1.3 12/01/06
308  Add VERBOSE input parameter
309  J. Gales
310 
311  Revision 2.1.2 11/10/06
312  Fix interp_distance value definition
313  J. Gales
314 
315  Revision 2.1.1 09/25/06
316  Set MAXNFILES to 544 (Fix made in readL2scan.h)
317  J. Gales
318 
319  Revision 2.1.0 09/05/06
320  Fix problem with prodtype=regional and fileuse
321  J. Gales
322 
323  Revision 2.0.9 05/12/06
324  Fix problem when escan_row > bscan_row.
325  J. Gales
326 
327  Revision 2.0.8 04/27/06
328  Fix dataday problems for day granules
329  J. Gales
330 
331  Revision 2.0.7 04/10/06
332  Support for TERRA SST
333  J. Gales
334 
335  Revision 2.0.6 03/31/06
336  Add longitude boundary check
337  (input.lonwest & input.loneast input parameters)
338  J. Gales
339 
340  Revision 2.0.5 03/16/06
341  Added support for HMODIST
342  J. Gales
343 */
344 
345 int32_t get_l2prod_index(const l2_prod& l2, const char* prodname) {
346  int32_t index;
347  for (index = 0; index < l2.nprod; index++)
348  if (strcmp(prodname, l2.prodname[index]) == 0)
349  break;
350  if (index == l2.nprod)
351  index = -1;
352  return index;
353 }
354 
355 void addPixelToBin(int32_t ifile, int32_t ipixl, uint64_t bin, double areaFrac = 1.0) {
356 
357  int32_t ibin = bin - basebin[krow];
358 
359  /* if bin not within bin row group then continue */
360  /* --------------------------------------------- */
361  if (ibin < 0 ||
362  ibin >= n_bins_in_group ||
363  (int64_t) bin >= basebin[krow + n_rows_in_group] ||
364  (int64_t) bin < basebin[krow])
365  return;
366 
367  /* GOOD OBSERVATION FOUND */
368  /* ---------------------- */
369 #ifdef MALLINFO
370  if (input.dcinfo) {
371  if ((l2_str[ifile].longitude[ipixl] <= -160) ||
372  (l2_str[ifile].longitude[ipixl] >= +160)) {
373  printf("DC: %10ld %12d %8.2f %8.2f\n",
374  (long int) bin,
375  (int) dbldate,
376  l2_str[ifile].longitude[ipixl],
377  l2_str[ifile].latitude[ipixl]);
378  }
379  }
380 #endif
381 
382  /* "OR" flags in swath pixel & set tilt & increment nscenes */
383  /* -------------------------------------------------------- */
384  bin_flag[ibin] = bin_flag[ibin] | l2_str[ifile].l2_flags[ipixl];
385  tilt[ibin] = tiltstate;
386  if (ifile != lastfile[ibin]) {
387  nscenes[ibin]++;
388  lastfile[ibin] = ifile;
389  }
390 
391  /* Allocate space for file index & bin data values */
392  /* ----------------------------------------------- */
393  if (file_index[ibin] == NULL) {
394  file_index[ibin] = (int16_t *)
395  calloc(n_allocperbin, sizeof (int16_t));
396 
397  data_values[ibin] = (float32 *)
398  calloc(n_allocperbin * l3b_nprod, sizeof (float32));
399  if (data_values[ibin] == 0x0) {
400  perror(buf);
401  printf("Allocation failed for data_values[ibin]: %d %s\n",
402  ibin, buf);
403  exit(EXIT_FAILURE);
404  }
405 
406  data_areas[ibin] = (double *) calloc(n_allocperbin, sizeof (double));
407  if (data_areas[ibin] == 0x0) {
408  perror(buf);
409  printf("Allocation failed for data_areas[ibin]: %d %s\n",
410  ibin, buf);
411  exit(EXIT_FAILURE);
412  }
413 
414  data_quality[ibin] = (uint8_t *)
415  calloc(n_allocperbin, sizeof (uint8_t));
416  if (data_quality[ibin] == 0x0) {
417  perror(buf);
418  printf("Allocation failed for data_quality[ibin]: %d %s\n",
419  ibin, buf);
420  exit(EXIT_FAILURE);
421  }
422 
423  time_value[ibin] = (float64 *)
424  calloc(n_allocperbin, sizeof (float64));
425  if (time_value[ibin] == 0x0) {
426  perror(buf);
427  printf("Allocation failed for time_value[ibin]: %d %s\n",
428  ibin, buf);
429  exit(EXIT_FAILURE);
430  }
431 
433  }
434 
435  /* Set file_index for each observation */
436  /* ----------------------------------- */
437  file_index[ibin][nobs[ibin]] = ifile;
438 
439  /* Get data quality */
440  /* ---------------- */
441  if (input.qual_prod[0] != 0) {
442  data_quality[ibin][nobs[ibin]] = l2_str[ifile].l2_data[qual_prod_index[ifile]][ipixl];
443  // a temporary fix to mask the last 4 pixels of a MODIS Terra scan
444  // for SST product
445  if (strncmp(input.suite, "SST", 3) == 0
446  && sensorID[ifile] == MODIST && ipixl > 1349) {
447  data_quality[ibin][nobs[ibin]] = 3;
448  }
449  }
450 
451  /* Store time_value (TAI93) */
452  /* ---------------------- */
453  time_value[ibin][nobs[ibin]] = time_tai93;
454 
455  /* Get composite data */
456  /* ------------------- */
457  if (input.composite_prod[0] != 0) {
458  int idx = composite_prod_index[ifile];
459  f32 = l2_str[ifile].l2_data[idx][ipixl];
460  if (f32 == -32767)
461  return;
462 
463  if (nobs[ibin] != 0) {
464  if (strcmp(input.composite_scheme, "max") == 0) {
465  if (f32 < data_values[ibin][composite_l3prod_index])
466  return;
467  } else {
468  if (f32 > data_values[ibin][composite_l3prod_index])
469  return;
470  }
471  nobs[ibin] = 0;
472  }
473  }
474 
475  /* Get data area for pixel */
476  /* ----------------------------------- */
477  data_areas[ibin][nobs[ibin]] = areaFrac;
478 
479  /* Get data values for all L3 products */
480  /* ----------------------------------- */
481  for (int32_t l3_iprod = 0; l3_iprod < l3b_nprod; l3_iprod++) {
482 
483  int32_t l2_iprod = numer[ifile][l3_iprod];
484  f32 = l2_str[ifile].l2_data[l2_iprod][ipixl*l2_str[ifile].thirdDim[l2_iprod] + thirdDimId[l3_iprod]];
485 
486  /* Set -32767 value to "bad" quality */
487  if (f32 == -32767)
488  if (input.qual_prod[0] != 0)
489  data_quality[ibin][nobs[ibin]] = 4;
490 
491  if (input.composite_prod[0] != 0) {
492  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = f32;
493  } else {
494  if (denom[ifile][l3_iprod] == -1 && f32 >= min_value[l3_iprod])
495  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = f32;
496 
497  if (denom[ifile][l3_iprod] == -1 && f32 < min_value[l3_iprod])
498  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = min_value[l3_iprod];
499 
500  if (denom[ifile][l3_iprod] == -2)
501  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] =
502  (l2_str[ifile].l2_flags[ipixl] >>
503  numer[ifile][l3_iprod]) & 1;
504  }
505 
506  /* ratio product */
507  /* ------------- */
508  if (denom[ifile][l3_iprod] >= 0) {
509 
510  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = f32;
511 
512  l2_iprod = denom[ifile][l3_iprod];
513  f32 = l2_str[ifile].l2_data[l2_iprod][ipixl*l2_str[ifile].thirdDim[l2_iprod] + thirdDimId[l3_iprod]];
514 
515  if (f32 >= min_value[l3_iprod])
516  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] /= f32;
517  else
518  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] /= min_value[l3_iprod];
519  }
520 
521  } /* iprod loop */
522 
523  /* Increment number of observations in bin */
524  /* --------------------------------------- */
525  nobs[ibin]++;
526  if (input.composite_prod[0] != 0) nobs[ibin] = 1;
527 
528  /* Reallocate if necessary */
529  /* ----------------------- */
530  if (nobs[ibin] == allocated_space[ibin]) {
531 
532  file_index[ibin] = (int16_t *)
533  realloc(file_index[ibin],
534  (nobs[ibin] + n_allocperbin)
535  * sizeof (int16_t));
536 
537  data_values[ibin] = (float32 *)
538  realloc(data_values[ibin],
539  (nobs[ibin] + n_allocperbin)
540  * l3b_nprod * sizeof (float32));
541  if (data_values[ibin] == 0x0) {
542  perror(buf);
543  printf("Reallocation failed for data_values[ibin]: %d %s\n",
544  ibin, buf);
545  exit(EXIT_FAILURE);
546  }
547 
548  data_areas[ibin] = (double *) realloc(data_areas[ibin],
549  (nobs[ibin] + n_allocperbin) * sizeof (double));
550  if (data_areas[ibin] == 0x0) {
551  perror(buf);
552  printf("Reallocation failed for data_areas[ibin]: %d %s\n",
553  ibin, buf);
554  exit(EXIT_FAILURE);
555  }
556 
557  data_quality[ibin] = (uint8_t *)
558  realloc(data_quality[ibin],
559  (nobs[ibin] + n_allocperbin)
560  * sizeof (uint8_t));
561  if (data_quality[ibin] == 0x0) {
562  perror(buf);
563  printf("Reallocation failed for data_quality[ibin]: %d %s\n",
564  ibin, buf);
565  exit(EXIT_FAILURE);
566  }
567 
568  time_value[ibin] = (float64 *)
569  realloc(time_value[ibin],
570  (nobs[ibin] + n_allocperbin)
571  * sizeof (float64));
572  if (time_value[ibin] == 0x0) {
573  perror(buf);
574  printf("Reallocation failed for time_value[ibin]: %d %s\n",
575  ibin, buf);
576  exit(EXIT_FAILURE);
577  }
578 
580  } /* end reallocate */
581 
582 }
583 
584 bool binIntersectsPixel(int32_t row, int32_t col, Box_t &pixelBox, double &areaFrac) {
585  bool result = false;
586  double n,s,e,w;
587  shape->rowcol2bounds(row, col, n, s, e, w);
588 
589 
590 #ifdef DEBUG_PRINT2
591  // debug plot the bin
592  if(enableDebugPrint) {
593  printf("lat = [%f, %f, %f, %f, %f]\n", s, n, n, s, s);
594  printf("lon = [%f, %f, %f, %f, %f]\n", w, w, e, e, w);
595  printf("plt.plot(lon, lat)\n\n");
596  }
597 #endif
598 
599 
600  Box_t box(Point_t(w,s), Point_t(e,n));
601  areaFrac = 0;
602 
603  if(!bg::disjoint(box, pixelBox)) {
604  Box_t output;
605  if(bg::intersection(box, pixelBox, output)) {
606  double intersectArea = bg::area(output);
607  if(intersectArea > 0) {
608  result = true;
609  double binArea = (n-s) * (e-w);
610  areaFrac = intersectArea / binArea;
611  }
612  }
613  }
614  return result;
615 }
616 
617 bool binIntersectsPixel(int32_t row, int32_t col, Polygon_t &pixelPoly, double &areaFrac) {
618  bool result = false;
619  double n,s,e,w;
620  shape->rowcol2bounds(row, col, n, s, e, w);
621 
622 
623 #ifdef DEBUG_PRINT2
624  // debug plot the bin
625  if(enableDebugPrint) {
626  printf("lat = [%f, %f, %f, %f, %f]\n", s, n, n, s, s);
627  printf("lon = [%f, %f, %f, %f, %f]\n", w, w, e, e, w);
628  printf("plt.plot(lon, lat)\n\n");
629  }
630 #endif
631 
632 
633  Box_t box(Point_t(w,s), Point_t(e,n));
634  areaFrac = 0;
635 
636  if(!bg::disjoint(box, pixelPoly)) {
637  std::deque<Polygon_t> output;
638  if(bg::intersection(box, pixelPoly, output)) {
639  double binArea = (n-s) * (e-w);
640  BOOST_FOREACH(Polygon_t const& p, output)
641  {
642  double intersectArea = bg::area(p);
643  if(intersectArea > 0.0) {
644  result = true;
645  areaFrac += intersectArea / binArea;
646  }
647  }
648  }
649  }
650  return result;
651 }
652 
653 template<class T>
654 bool getBinsFromRow(double lat, double lon, T &pixelPoly, std::map<uint64_t, double> &areas) {
655  int32_t row0, col0;
656  int32_t col;
657  bool result = false;
658  double areaFrac;
659  uint64_t bin;
660 
661  shape->latlon2rowcol(lat, lon, row0, col0);
662 
663  // look right
664  col = col0;
665  while(binIntersectsPixel(row0, col, pixelPoly, areaFrac)) {
666  result = true;
667  bin = shape->rowcol2bin(row0, col);
668  areas.emplace(bin, areaFrac);
669  col++;
670  shape->constrainRowCol(row0, col);
671  if(col == col0)
672  break;
673  }
674 
675  // look left
676  col = col0 - 1;
677  while(binIntersectsPixel(row0, col, pixelPoly, areaFrac)) {
678  result = true;
679  bin = shape->rowcol2bin(row0, col);
680  areas.emplace(bin, areaFrac);
681  col--;
682  shape->constrainRowCol(row0, col);
683  if(col == col0 - 1)
684  break;
685  }
686 
687  return result;
688 }
689 
690 void getBins(int32_t ifile, int32_t ipixl, std::map<uint64_t, double> &areas) {
691 
692  areas.clear();
693 
694  double lat0 = l2_str[ifile].latitude[ipixl];
695  double lon0 = l2_str[ifile].longitude[ipixl];
696  int32_t row0;
697  double lat;
698 
699 
700 #ifdef DEBUG_PRINT
701  // debug plotting
702  // plot center point
703  if(enableDebugPrint) {
704  printf("lat=%f\n", lat0);
705  printf("lon=%f\n", lon0);
706  printf("plt.plot(lon, lat, \"ro\")\n\n");
707  }
708 #endif
709 
710 
711  // use the center of the bin lat for search
712  row0 = shape->lat2row(lat0);
713  lat0 = shape->row2lat(row0);
714 
715  double deltaLat = 180.0 / shape->getNumRows();
716 
717 
718  if(input.area_weighting == 1) {
719  Box_t pixelBox;
720 
721  // use the pixel box
722  // lat1/lon1 contain pixel delta
723  pixelBox.min_corner().set<0>(l2_str[ifile].longitude[ipixl] - l2_str[ifile].lon1[ipixl]);
724  pixelBox.min_corner().set<1>(l2_str[ifile].latitude[ipixl] - l2_str[ifile].lat1[ipixl]);
725  pixelBox.max_corner().set<0>(l2_str[ifile].longitude[ipixl] + l2_str[ifile].lon1[ipixl]);
726  pixelBox.max_corner().set<1>(l2_str[ifile].latitude[ipixl] + l2_str[ifile].lat1[ipixl]);
727 
728 #ifdef DEBUG_PRINT
729  // debug
730  // plot pixel box
731  if(enableDebugPrint) {
732  printf("lat = [%f, %f, %f, %f, %f]\n",
733  pixelBox.min_corner().y(), pixelBox.max_corner().y(), pixelBox.max_corner().y(),
734  pixelBox.min_corner().y(), pixelBox.min_corner().y());
735  printf("lon = [%f, %f, %f, %f, %f]\n",
736  pixelBox.min_corner().x(), pixelBox.min_corner().x(), pixelBox.max_corner().x(),
737  pixelBox.max_corner().x(), pixelBox.min_corner().x());
738  printf("plt.plot(lon, lat)\n\n");
739  }
740 #endif
741 
742  // look up
743  lat = lat0;
744  while(getBinsFromRow(lat, lon0, pixelBox, areas)) {
745  lat += deltaLat;
746  if(lat > 90.0)
747  break;
748  }
749 
750  // look down
751  lat = lat0 - deltaLat;
752  while(getBinsFromRow(lat, lon0, pixelBox, areas)) {
753  lat -= deltaLat;
754  if(lat < -90.0)
755  break;
756  }
757 
758 
759  } else if(input.area_weighting == 2) {
760  Box_t pixelBox;
761 
762  // use the pixel bounding box
763  float latMax = std::max(l2_str[ifile].lat1[ipixl], l2_str[ifile].lat1[ipixl+1]);
764  latMax = std::max(latMax, l2_str[ifile].lat2[ipixl]);
765  latMax = std::max(latMax, l2_str[ifile].lat2[ipixl+1]);
766 
767  float latMin = std::min(l2_str[ifile].lat1[ipixl], l2_str[ifile].lat1[ipixl+1]);
768  latMin = std::min(latMin, l2_str[ifile].lat2[ipixl]);
769  latMin = std::min(latMin, l2_str[ifile].lat2[ipixl+1]);
770 
771  float lonMax = std::max(l2_str[ifile].lon1[ipixl], l2_str[ifile].lon1[ipixl+1]);
772  lonMax = std::max(lonMax, l2_str[ifile].lon2[ipixl]);
773  lonMax = std::max(lonMax, l2_str[ifile].lon2[ipixl+1]);
774 
775  float lonMin = std::min(l2_str[ifile].lon1[ipixl], l2_str[ifile].lon1[ipixl+1]);
776  lonMin = std::min(lonMin, l2_str[ifile].lon2[ipixl]);
777  lonMin = std::min(lonMin, l2_str[ifile].lon2[ipixl+1]);
778 
779  if((lonMax - lonMin) > 180) {
780  float tmpF = lonMin;
781  lonMin = lonMax;
782  lonMax = tmpF + 360.0;
783  }
784 
785  pixelBox.min_corner().set<0>(lonMin);
786  pixelBox.min_corner().set<1>(latMin);
787  pixelBox.max_corner().set<0>(lonMax);
788  pixelBox.max_corner().set<1>(latMax);
789 
790 
791 #ifdef DEBUG_PRINT
792  // debug
793  // plot pixel box
794  if(enableDebugPrint) {
795  printf("lat = [%f, %f, %f, %f, %f]\n",
796  pixelBox.min_corner().y(), pixelBox.max_corner().y(), pixelBox.max_corner().y(),
797  pixelBox.min_corner().y(), pixelBox.min_corner().y());
798  printf("lon = [%f, %f, %f, %f, %f]\n",
799  pixelBox.min_corner().x(), pixelBox.min_corner().x(), pixelBox.max_corner().x(),
800  pixelBox.max_corner().x(), pixelBox.min_corner().x());
801  printf("plt.plot(lon, lat)\n\n");
802  }
803 #endif
804 
805  // look up
806  lat = lat0;
807  while(getBinsFromRow(lat, lon0, pixelBox, areas)) {
808  lat += deltaLat;
809  if(lat > 90.0)
810  break;
811  }
812 
813  // look down
814  lat = lat0 - deltaLat;
815  while(getBinsFromRow(lat, lon0, pixelBox, areas)) {
816  lat -= deltaLat;
817  if(lat < -90.0)
818  break;
819  }
820 
821  } else {
822  // use the exact pixel polygon
823  Polygon_t pixelPoly;
824 
825 
826 #ifdef DEBUG_PRINT
827  // debug
828  // plot pixel polygon
829  if(enableDebugPrint) {
830  printf("lat = [%f, %f, %f, %f, %f]\n",
831  l2_str[ifile].lat1[ipixl],
832  l2_str[ifile].lat2[ipixl],
833  l2_str[ifile].lat2[ipixl+1],
834  l2_str[ifile].lat1[ipixl+1],
835  l2_str[ifile].lat1[ipixl]);
836  printf("lon = [%f, %f, %f, %f, %f]\n",
837  l2_str[ifile].lon1[ipixl],
838  l2_str[ifile].lon2[ipixl],
839  l2_str[ifile].lon2[ipixl+1],
840  l2_str[ifile].lon1[ipixl+1],
841  l2_str[ifile].lon1[ipixl]);
842  printf("plt.plot(lon, lat)\n\n");
843  }
844 #endif
845 
846 
847  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon1[ipixl], l2_str[ifile].lat1[ipixl]));
848  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon2[ipixl], l2_str[ifile].lat2[ipixl]));
849  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon2[ipixl+1], l2_str[ifile].lat2[ipixl+1]));
850  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon1[ipixl+1], l2_str[ifile].lat1[ipixl+1]));
851  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon1[ipixl], l2_str[ifile].lat1[ipixl]));
852 
853  // make sure the polygon is defined properly
854  bg::correct(pixelPoly);
855 
856  // look up
857  lat = lat0;
858  while(getBinsFromRow(lat, lon0, pixelPoly, areas)) {
859  lat += deltaLat;
860  if(lat > 90.0)
861  break;
862  }
863 
864  // look down
865  lat = lat0 - deltaLat;
866  while(getBinsFromRow(lat, lon0, pixelPoly, areas)) {
867  lat -= deltaLat;
868  if(lat < -90.0)
869  break;
870  }
871 
872  }
873 
874 }
875 
876 int main(int argc, char **argv) {
877  int i, j, k;
878  int status;
879  intn ret_status = 0;
880  int plusday = 0;
881 
882  int32_t index;
883  int32_t ifile, jsrow, ipixl;
884  uint64_t bin;
885  int32_t ibin;
886  int32_t nfiles;
887  int *fileused;
888  int32_t n_active_files;
889 
890  int32_t within_flag;
891 
892  int32_t n_filled_bins;
893  int64_t total_filled_bins = 0;
894  int32_t date;
895 
896  int16_t brk_scan[MAXNFILES];
897 
898  float32 *slat = NULL;
899  float32 *elat = NULL;
900  float32 *clat = NULL;
901  float32 *slon = NULL;
902  float32 *elon = NULL;
903  float32 *clon = NULL;
904 
905  float32 dlat;
906  double latbin = 0.0;
907  double lonbin = 0.0;
908 
909  int32_t igroup, ngroup;
910 
911  static int32_t * bscan_row[MAXNFILES];
912  static int32_t * escan_row[MAXNFILES];
913  static unsigned char *scan_in_rowgroup[MAXNFILES];
914 
915  int32_t fileid_w;
916  int32_t out_grp;
917 
918  static ptrdiff_t stride[] = { 20, 5, 1 };
919  int64_t *beg;
920  int32_t *ext;
921  int64_t *binnum_data;
922  int32_t i32;
923  int32_t iprod;
924 
925  time_t diffday_beg, diffday_end;
926  int32_t syear, sday, eyear, eday;
927 
928  int32_t ntilts;
929  int16_t tilt_flags[MTILT_DIMS_2];
930  int16_t tilt_ranges[LTILT_DIMS_2][MTILT_DIMS_2];
931 
932  uint32_t flagusemask;
933  uint32_t required;
934  uint32_t flagcheck;
935 
936  int32_t proc_day_beg, proc_day_end;
937 
938  uint8_t *best_qual, qual_max_allowed;
939 
940  double *sum_bin;
941  double *sum2_bin;
942  double wgt;
943  float32 northmost = -90.0, southmost = 90.0, eastmost = -180.0, westmost = 180.0;
944 
945  int32_t bad_value;
946 
947  double time_rec = 0;
948  time_t tnow;
949  time_t l2_stimes[MAXNFILES];
950 
951  int8_t **dorn;
952  float onorth, osouth, owest, oeast, eqcross;
953  int dateline;
954  int32_t dataday0, dataday1, startdate, enddate;
955  double ddstime, ddetime, dbldate;
956 
957  float **dolat, **dolon, *lonArray, *latArray;
958  float *dnlat[2], *dnlon[2];
959  int32_t n[2];
960  int32_t *years, *days, *msecs;
961  struct tm *tmnow;
962 
963  vector<NcDim> dims_nc;
964  int32_t dims[3];
965  int32_t small_dims[3];
966 
967  static meta_l2Type meta_l2;
968  static meta_l3bType meta_l3b;
969 
970  static float lonRanges[4];
971  lonRanges[0] = 0; // min, max in Western hemisphere
972  lonRanges[1] = -180;
973  lonRanges[2] = 180; // min, max in Eastern hemisphere
974  lonRanges[3] = 0;
975 
976  char units[1024];
977  vector<string> l2_prodname;
978  vector<string> l3_prodname;
979 
980  /* Function Prototypes */
981  int64_t getbinnum(int32_t, int32_t);
982 
983  FILE *fp2 = NULL;
984 
985  int dnStatus;
986 
987  binListStruct_nc* binList32nc = NULL;
988  binListStruct64_nc* binList64nc = NULL;
989 
990  binIndexStruct_nc* binIndex32nc = NULL;
991  binIndexStruct64_nc* binIndex64nc = NULL;
992 
993 #ifdef MALLINFO
994  struct mallinfo minfo;
995 #endif
996 
998 
999  /* From Fred Patt
1000 
1001  sum(data) sum(data)*sqrt(n)
1002  s = --------- = ----------------- = avg(data)*sqrt(n)
1003  sqrt(n) n
1004 
1005  */
1006 
1007  setlinebuf(stdout);
1008  printf("%s %s (%s %s)\n", PROGRAM, VERSION, __DATE__, __TIME__);
1009 
1010  l2bin_input(argc, argv, &input, PROGRAM, VERSION);
1011 
1012  cdata_(); // initialize global FORTRAN common block data for l_sun call
1013 
1014  switch(input.area_weighting) {
1015  case 0:
1016  enableL2PixelArea(L2PixelOff); // no corner calc
1017  break;
1018  case 1:
1019  enableL2PixelArea(L2PixelDelta); // pixel deltas
1020  break;
1021  default:
1022  enableL2PixelArea(L2PixelCorner); // pixel corners
1023  break;
1024  }
1025 
1026  // Get times for each file
1027  nfiles = input.files.size();
1028  printf("%d input files\n", nfiles);
1029  for (ifile = 0; ifile < nfiles; ifile++) {
1030  status = openL2(input.files[ifile].c_str(), 0x0, &l2_str[ifile]);
1031  l2_stimes[ifile] = (time_t) yds2unix(l2_str[ifile].syear,
1032  l2_str[ifile].sday,
1033  l2_str[ifile].smsec / 1000); // 86400;
1034  closeL2(&l2_str[ifile], ifile);
1035  }
1036 
1037  if ((fileused = (int *) calloc(nfiles, sizeof(int))) == NULL) {
1038  printf("-E- Problem allocating memory for fileused element of metadata structure\n");
1039  exit(EXIT_FAILURE);
1040  }
1041 
1042  proc_day_beg = input.sday;
1043  proc_day_end = input.eday;
1044  syear = (int32_t) input.sday / 1000.;
1045  sday = input.sday - 1000 * syear;
1046  startdate = (int32_t) (yds2unix(syear, sday, 0) / 86400);
1047  eyear = (int32_t) input.eday / 1000.;
1048  eday = input.eday - 1000 * eyear;
1049  enddate = (int32_t) (yds2unix(eyear, eday, 0) / 86400);
1050 
1051  if (strcmp(input.resolve, "1D") == 0) {
1052  nrows = 180;
1053  meta_l3b.resolution = str2resolution("1deg");
1054  } else if (strcmp(input.resolve, "HD") == 0) {
1055  nrows = 360;
1056  meta_l3b.resolution = str2resolution("0.5deg");
1057  } else if (strcmp(input.resolve, "QD") == 0) {
1058  nrows = 720;
1059  meta_l3b.resolution = str2resolution("0.25deg");
1060  } else if (strcmp(input.resolve, "36") == 0) {
1061  nrows = 2160 / 4;
1062  meta_l3b.resolution = str2resolution("36km");
1063  } else if (strcmp(input.resolve, "18") == 0) {
1064  nrows = 2160 / 2;
1065  meta_l3b.resolution = str2resolution("18km");
1066  } else if (strcmp(input.resolve, "9") == 0) {
1067  nrows = 2160;
1068  meta_l3b.resolution = str2resolution("9km");
1069  } else if (strcmp(input.resolve, "4") == 0) {
1070  nrows = 2160 * 2;
1071  meta_l3b.resolution = str2resolution("4km");
1072  } else if (strcmp(input.resolve, "2") == 0) {
1073  nrows = 2160 * 4;
1074  meta_l3b.resolution = str2resolution("2km");
1075  } else if (strcmp(input.resolve, "1") == 0) {
1076  nrows = 2160 * 8;
1077  meta_l3b.resolution = str2resolution("1km");
1078  } else if (strcmp(input.resolve, "H") == 0) {
1079  nrows = 2160 * 16;
1080  meta_l3b.resolution = str2resolution("hkm");
1081  } else if (strcmp(input.resolve, "Q") == 0) {
1082  nrows = 2160 * 32;
1083  meta_l3b.resolution = str2resolution("qkm");
1084  } else if (strcmp(input.resolve, "HQ") == 0) {
1085  nrows = 2160 * 80;
1086  meta_l3b.resolution = str2resolution("hqkm");
1087  } else if (strcmp(input.resolve, "HH") == 0) {
1088  nrows = 2160 * 160;
1089  meta_l3b.resolution = str2resolution("hhkm");
1090  }
1091  if (nrows == -1) {
1092  printf("Grid resolution not defined.\n");
1093  exit(EXIT_FAILURE);
1094  }
1095  dlat = 180. / nrows;
1096  bool is64bit = (nrows > 2160 * 16);
1097 
1098  qual_max_allowed = input.qual_max;
1099 
1100  printf("Resolution: %s\n", input.resolve);
1101  printf("Max Qual Allowed: %d\n", input.qual_max);
1102 
1103  /* Setup flag mask */
1104  /* --------------- */
1105  strcpy(buf, l2_str[0].flagnames);
1106  setupflags(buf, input.flaguse, &flagusemask, &required, &status);
1107  printf("flagusemask: %d\n", flagusemask);
1108  printf("required: %d\n", required);
1109 
1110  // If l3bprod="ALL", replace with product list from first file
1111  if (boost::iequals(input.l3bprod, "ALL")) {
1112  vector<string> l2prods(l2_str[0].prodname, l2_str[0].prodname+l2_str[0].nprod);
1113  string joined = boost::algorithm::join(l2prods, ":");
1114  strcpy(input.l3bprod, joined.c_str());
1115 
1116  // fix l3bprod entry in Input Parameters attribute
1117  string parms = input.parms;
1118  boost::ireplace_all(parms, "l3bprod = ALL", "l3bprod = " + joined);
1119  strcpy(input.parms, parms.c_str());
1120  }
1121 
1122  // Parse L3 Product list
1123  string delim1 = ":, "; // product delimiters
1124  string l3bprod = input.l3bprod;
1125  boost::trim_if(l3bprod, boost::is_any_of(delim1));
1126  vector<string> prodparam;
1127  boost::algorithm::split(prodparam, l3bprod,
1128  boost::is_any_of(delim1));
1129 
1130  // expand 3D products
1131  string delim2 = ";="; // minval delimiters
1132  string wave_string=std::string{","} + std::string{input.output_wavelengths} + ",";
1133 
1134  for (size_t iprod = 0; iprod < prodparam.size(); iprod++) {
1135 
1136  // get prodname and min value
1137  vector<string> thisprod;
1138  boost::algorithm::split(thisprod, prodparam[iprod],
1139  boost::is_any_of(delim2));
1140 
1141  int32_t l2_iprod = get_l2prod_index(l2_str[0], thisprod[0].c_str());
1142  if(l2_iprod == -1) {
1143  printf("-E- Product: %s was not found in the L2 file: %s\n", thisprod[0].c_str(), l2_str[0].filename);
1144  exit(EXIT_FAILURE);
1145  }
1146  int32_t tmpThirdDim = l2_str[0].thirdDim[l2_iprod];
1147  if (tmpThirdDim == 1) {
1148  l2_prodname.push_back(thisprod[0]);
1149  l3_prodname.push_back(thisprod[0]);
1150  thirdDimId.push_back(0);
1151  if (thisprod.size() > 1) {
1152  min_value.push_back((float32) strtod(thisprod[1].c_str(), NULL));
1153  } else {
1154  min_value.push_back(BAD_FLT);
1155  }
1156  } else {
1157  for (int32_t ibpp=0; ibpp < tmpThirdDim; ibpp++) {
1158 
1159  string subprodname = thisprod[0];
1160  string tag = to_string(l2_str[0].wavelength[ibpp]);
1161  subprodname.append("_");
1162  subprodname.append(tag);
1163 
1164  if (wave_string != ",ALL,") {
1165  size_t found = wave_string.find(","+tag+",");
1166  if (found == string::npos) {
1167  continue;
1168  }
1169  }
1170 
1171  l2_prodname.push_back(thisprod[0]);
1172  l3_prodname.push_back(subprodname);
1173  thirdDimId.push_back(ibpp);
1174  if (thisprod.size() > 1) {
1175  min_value.push_back((float32) strtod(thisprod[1].c_str(), NULL));
1176  } else {
1177  min_value.push_back(BAD_FLT);
1178  }
1179  } /* ibpp loop */
1180  }
1181  } /* iprod loop */
1182 
1183  l3b_nprod = l3_prodname.size();
1184  if(l3b_nprod > max_l3b_products) {
1185  printf("-E- Number of output products is greater than %d\n", max_l3b_products);
1186  exit(EXIT_FAILURE);
1187  }
1188 
1189  /* Initialize bscan_row, escan_row, numer, denom */
1190  /* --------------------------------------------- */
1191  for (i = 0; i < MAXNFILES; i++) {
1192  bscan_row[i] = NULL;
1193  escan_row[i] = NULL;
1194  numer[i] = NULL;
1195  denom[i] = NULL;
1196  scan_in_rowgroup[i] = NULL;
1197  }
1198 
1199  /* Check each L2 file for required products */
1200  /* ---------------------------------------- */
1201  for (ifile = 0; ifile < nfiles; ifile++) {
1202 
1203  numer[ifile] = (int16_t *) calloc(l3b_nprod, sizeof(int16_t));
1204  denom[ifile] = (int16_t *) calloc(l3b_nprod, sizeof(int16_t));
1205 
1206  // Check whether L3 products exist in L2
1207  for (iprod = 0; iprod < l3b_nprod; iprod++) {
1208 
1209  vector<string> operands;
1210  boost::algorithm::split(operands, l2_prodname[iprod],
1211  boost::is_any_of("/"));
1212 
1213  // find numerator, or solo product
1214  index = get_l2prod_index(l2_str[ifile], operands[0].c_str());
1215  if (index < 0) {
1216  printf("L3 product: \"%s\" not found in L2 file \"%s\".\n",
1217  operands[0].c_str(), l2_str[ifile].filename);
1218  exit(EXIT_FAILURE);
1219  } else
1220  numer[ifile][iprod] = index;
1221 
1222 
1223  // find denominator product (as needed)
1224  denom[ifile][iprod] = -1;
1225  if (operands.size() > 1) {
1226  boost::replace_all(l3_prodname[iprod], "/", "_");
1227  index = get_l2prod_index(l2_str[ifile], operands[1].c_str());
1228  if (index < 0) {
1229  printf("L3 product: \"%s\" not found in L2 file \"%s\".\n",
1230  operands[1].c_str(), l2_str[ifile].filename);
1231  exit(EXIT_FAILURE);
1232  } else
1233  denom[ifile][iprod] = index;
1234  }
1235 
1236  } // iprod loop
1237 
1238  // Check whether Quality product exists in L2
1239  if (input.qual_prod[0] != 0) {
1240  index = get_l2prod_index(l2_str[ifile], input.qual_prod);
1241  if (index < 0) {
1242  printf("Quality product: \"%s\" not found in L2 file \"%s\".\n",
1243  input.qual_prod, l2_str[ifile].filename);
1244  exit(EXIT_FAILURE);
1245  } else
1247  }
1248 
1249  // Check whether Composite product exists in L2
1250  if (input.composite_prod[0] != 0) {
1251  index = get_l2prod_index(l2_str[ifile], input.composite_prod);
1252  if (index < 0) {
1253  printf("Composite product: \"%s\" not found in L2 file \"%s\".\n",
1254  input.composite_prod, l2_str[ifile].filename);
1255  exit(EXIT_FAILURE);
1256  } else
1258  }
1259 
1260  } // ifile loop
1261 
1262 
1263  /* Check whether composite product exists in L3 product list */
1264  /* --------------------------------------------------------- */
1265  if (input.composite_prod[0] != 0) {
1266  for (iprod = 0; iprod < l3b_nprod; iprod++) {
1267  if(l3_prodname[iprod] == input.composite_prod) {
1268  composite_l3prod_index = iprod;
1269  break;
1270  }
1271  }
1272  if (composite_l3prod_index == -1) {
1273  printf("Composite product: \"%s\" not found in L3 product list.\n",
1274  input.composite_prod);
1275  exit(EXIT_FAILURE);
1276  }
1277  }
1278 
1279  /* Find begin and end scan latitudes for each swath row */
1280  /* ---------------------------------------------------- */
1281  for (ifile = 0; ifile < nfiles; ifile++) {
1282 
1283  slat = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1284  elat = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1285 
1286  bscan_row[ifile] = (int32_t *) calloc(l2_str[ifile].nrec, sizeof(int32_t));
1287  escan_row[ifile] = (int32_t *) calloc(l2_str[ifile].nrec, sizeof(int32_t));
1288 
1289  string instrument, platform;
1290  try {
1291  NcFile nc_input(input.files[ifile], NcFile::read);
1292  nc_input.getAtt("instrument").getValues(instrument);
1293  nc_input.getAtt("platform").getValues(platform);
1294  NcGroup scan_grp = nc_input.getGroup("scan_line_attributes");
1295  scan_grp.getVar("slat").getVar(slat);
1296  scan_grp.getVar("elat").getVar(elat);
1297  nc_input.close();
1298 
1299  } catch (NcException const & e) {
1300  e.what();
1301  exit(EXIT_FAILURE);
1302  } catch (exception const & e) {
1303  cerr << "-X- " << __FILE__ " line " << __LINE__ << ": "
1304  << e.what() << endl;
1305  exit(EXIT_FAILURE);
1306  }
1307 
1308  /* Get minimum value from product.xml */
1309  /* ---------------------------------- */
1311  platform.c_str());
1312  productInfo_t* info = allocateProductInfo();
1313  for (iprod = 0; iprod < l3b_nprod; iprod++) {
1314  if (min_value[iprod] == BAD_FLT) {
1315  if (findProductInfo(l2_prodname[iprod].c_str(), sensorID[ifile], info)) {
1316  min_value[iprod] = info->validMin;
1317  } else {
1318  min_value[iprod] = 0;
1319  }
1320  }
1321  }
1322  freeProductInfo(info);
1323 
1324  /* Note: bscan > escan */
1325 
1326  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++) {
1327  escan_row[ifile][jsrow] = (int32_t) ((90 + elat[jsrow]) / dlat);
1328  bscan_row[ifile][jsrow] = (int32_t) ((90 + slat[jsrow]) / dlat);
1329 
1330  if (escan_row[ifile][jsrow] > bscan_row[ifile][jsrow]) {
1331  k = escan_row[ifile][jsrow];
1332  escan_row[ifile][jsrow] = bscan_row[ifile][jsrow];
1333  bscan_row[ifile][jsrow] = k;
1334  }
1335  escan_row[ifile][jsrow] -= 10;
1336  bscan_row[ifile][jsrow] += 10;
1337  }
1338 
1339  free(slat);
1340  free(elat);
1341 
1342  } /* ifile loop */
1343 
1344  /* Find begin & end scans for each input file */
1345  /* ------------------------------------------ */
1346  n_active_files = nfiles;
1347  for (ifile = 0; ifile < nfiles; ifile++) {
1348  if (input.deltaeqcross == 0.0) {
1349  switch (sensorID[ifile]) {
1350  case HICO:
1351  eqcross = 12.0; // Made this up - RJH
1352  break;
1353  case MODISA:
1354  case VIIRSN:
1355  case VIIRSJ1:
1356  eqcross = 13.5;
1357  break;
1358  case MODIST:
1359  case MERIS:
1360  case OCTS:
1361  case HAWKEYE:
1362  eqcross = 10.5;
1363  break;
1364  case OLCIS3A:
1365  case OLCIS3B:
1366  eqcross = 10.0;
1367  break;
1368  case CZCS:
1369  eqcross = 12.0;
1370  break;
1371  case SEAWIFS:
1372  eqcross = 12.0;
1373  if (l2_str[ifile].syear > 2002) {
1374  // The constant 10957 makes d the number of days since 1 Jan. 2000
1375  int d = l2_stimes[ifile] / 86400 - 10957;
1376  /*
1377  * On 10 July 2010 (d=3843) OrbView-2/SeaWiFS was nearing the
1378  * end of its orbit-raising maneuvers. Before these maneuvers
1379  * the node-crossing time was drifting further into the afternoon
1380  * (first equation below). After the orbit raising, the node-crossing
1381  * time was drifting back towards noon (second equation).
1382 
1383  * Correction equations provided by Fred Patt.
1384  */
1385  double deg;
1386  if (d < 3843) {
1387  deg = 7.7517951e-10 * d * d * d
1388  - 2.1692192e-06 * d * d
1389  + 0.00706692410 * d
1390  - 4.1300585;
1391  } else {
1392  /*
1393  * This equation may need to be replaced with a polynomial
1394  * once we get more orbit data under our belts. (16-Jul-2010 N.Kuring)
1395  * ...note...the above comment was written before the demise of SeaWiFS
1396  * on 11-Dec-2010 ...no further modifications necessary...
1397  */
1398  deg = -0.024285181 * d + 128.86093;
1399  }
1400 
1401  // The above polynomials yield degrees; convert to hours.
1402  float hours = (float) (deg / 15.);
1403 
1404  eqcross = 12 + hours;
1405  }
1406  break;
1407 
1408  default:
1409  eqcross = 12.0;
1410 
1411  fprintf(stderr,
1412  "-Warning- Unknown equator crossing time for sensorID=%d file=%s\n ",
1413  sensorID[ifile], input.files[ifile].c_str());
1414  fprintf(stderr, "using crossing time of %f\n",eqcross);
1415  break;
1416 
1417  }
1418  } else {
1419  eqcross = 12.0 + (input.deltaeqcross / 60.);
1420  }
1421  // Shift the equatorial crossing time by 12 hours if doing night time binning
1422  // This is done so that if crossing the dateline the reference hour is being used
1423  // to move into a dataday in the same direction in function get_datadays
1424  if (input.night) {
1425  eqcross = eqcross - 12;
1426  if (eqcross < 0) {
1427  eqcross = eqcross + 24;
1428  plusday = 1; // if we do this, we're effectively going back
1429  // a day, and we so we need to add a day to the
1430  // output of get_datadays...
1431  }
1432  }
1433 
1434  slon = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1435  elon = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1436  clon = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1437  elat = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1438  slat = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1439  clat = (float32 *) calloc(l2_str[ifile].nrec, sizeof(float32));
1440  years = (int32_t *) calloc(l2_str[ifile].nrec, sizeof(int32_t));
1441  days = (int32_t *) calloc(l2_str[ifile].nrec, sizeof(int32_t));
1442  msecs = (int32_t *) calloc(l2_str[ifile].nrec, sizeof(int32_t));
1443  lonArray = (float *) malloc(sizeof(float)*l2_str[ifile].nrec * l2_str[ifile].nsamp);
1444  latArray = (float *) malloc(sizeof(float)*l2_str[ifile].nrec * l2_str[ifile].nsamp);
1445 
1446  try {
1447  NcFile nc_input(input.files[ifile], NcFile::read);
1448 
1449  NcGroup scan_grp = nc_input.getGroup("scan_line_attributes");
1450  scan_grp.getVar("slon").getVar(slon);
1451  scan_grp.getVar("elon").getVar(elon);
1452  scan_grp.getVar("clon").getVar(clon);
1453  scan_grp.getVar("slat").getVar(slat);
1454  scan_grp.getVar("elat").getVar(elat);
1455  scan_grp.getVar("clat").getVar(clat);
1456  scan_grp.getVar("year").getVar(years);
1457  scan_grp.getVar("day").getVar(days);
1458  scan_grp.getVar("msec").getVar(msecs);
1459 
1460  NcGroup nav_grp = nc_input.getGroup("navigation_data");
1461  nav_grp.getVar("longitude").getVar(lonArray);
1462  nav_grp.getVar("latitude").getVar(latArray);
1463 
1464  dims_nc = nav_grp.getVar("latitude").getDims();
1465  for (i = 0; i < 2; i++) {
1466  dims[i] = (int32_t) dims_nc[i].getSize();
1467  if (dims[i] / stride[i] < 5)
1468  stride[i] = 1;
1469  small_dims[i] = dims[i] / stride[i];
1470  }
1471 
1472  nc_input.close();
1473 
1474  } catch (NcException const & e) {
1475  e.what();
1476  exit(EXIT_FAILURE);
1477  } catch (exception const & e) {
1478  cerr << "-X- " << __FILE__ " line " << __LINE__ << ": "
1479  << e.what() << endl;
1480  exit(EXIT_FAILURE);
1481  }
1482 
1483  dolat = (float **) malloc(sizeof(float *)*small_dims[0]);
1484  dolon = (float **) malloc(sizeof(float *)*small_dims[0]);
1485 
1486  // Allocate and clear dolat/dolon arrays
1487  for (i = 0; i < small_dims[0]; i++) {
1488  dolat[i] = (float *) calloc(small_dims[1], sizeof(float));
1489  dolon[i] = (float *) calloc(small_dims[1], sizeof(float));
1490  }
1491 
1492  k = 0;
1493 
1494  for (j = 0; j < small_dims[0]; j++) {
1495  k = dims[1] * (j * stride[0]);
1496  for (i = 0; i < small_dims[1]; i++) {
1497  dolat[j][i] = latArray[k];
1498  dolon[j][i] = lonArray[k];
1499  k += stride[1];
1500  }
1501  }
1502 
1503  // reset stride
1504  stride[0] = 20;
1505  stride[1] = 5;
1506 
1507  /* Determine brk_scan value */
1508  /* ------------------------ */
1509  brk_scan[ifile] = 0;
1510 
1511  /* Regional Product */
1512  /* ---------------- */
1513  if (strcasecmp(input.prodtype, "regional") == 0) {
1514  printf("%s brk:%5d %5d %3d %6d\n",
1515  buf, brk_scan[ifile],
1516  l2_str[ifile].nrec, l2_str[ifile].sday,
1517  l2_str[ifile].smsec / 1000);
1518 
1519  free(slon);
1520  free(elon);
1521  free(clon);
1522  free(elat);
1523  free(slat);
1524  free(clat);
1525  for (i = 0; i < small_dims[0]; i++) {
1526  free(dolat[i]);
1527  free(dolon[i]);
1528  }
1529  free(dolat);
1530  free(dolon);
1531 
1532  free(years);
1533  free(days);
1534  free(msecs);
1535  free(lonArray);
1536  free(latArray);
1537 
1538  continue;
1539  }
1540 
1541  /* Allocate space for day or night array. */
1542  if ((dorn = (int8_t **) malloc(small_dims[0] * sizeof(int8_t *))) == NULL) {
1543  printf("%s -Error: Cannot allocate memory to dorn array\n",
1544  __FILE__);
1545  exit(EXIT_FAILURE);
1546  }
1547  for (i = 0; i < small_dims[0]; i++) {
1548  if ((dorn[i] = (int8_t *) calloc(small_dims[1], sizeof(int8_t))) == NULL) {
1549  printf("%s -Error: Cannot allocate memory to dorn array\n",
1550  __FILE__);
1551  exit(EXIT_FAILURE);
1552  }
1553  }
1554 
1555  date = (time_t) yds2unix(l2_str[ifile].syear,
1556  l2_str[ifile].sday,
1557  l2_str[ifile].smsec / 1000) / 86400;
1558  diffday_beg = date - startdate;
1559  diffday_end = date - enddate;
1560 
1561  // get the vertices of the outline of the swath into one array for lat/lon
1562  // to pass into get_coord_extrema, which calculates whether the dateline is crossed
1563  dnlat[0] = NULL;
1564  dnlat[1] = NULL;
1565  dnlon[0] = NULL;
1566  dnlon[1] = NULL;
1567  dnStatus = daynight_outlines(years, days, msecs,
1568  small_dims[1], small_dims[0],
1569  dolat,
1570  dolon,
1571  n, dnlat, dnlon, dorn);
1572 
1573  if (dnStatus != 0) {
1574  fprintf(stderr,
1575  "Consider QC fail for file: %s\n...look at the file "
1576  "though, as I might be lying...\n",
1577  l2_str[ifile].filename);
1578  if (ret_status == 0) {
1579  ret_status = 120;
1580 
1581  // if fileuse is set, open a corresponding qcfail file
1582  if (input.fileuse[0] != 0) {
1583  char *qcFailFile = (char*) malloc(strlen(input.fileuse) + 10);
1584  strcpy(qcFailFile, input.fileuse);
1585  strcat(qcFailFile, ".qcfail");
1586  fp2 = fopen(qcFailFile, "w");
1587  free(qcFailFile);
1588  }
1589 
1590  }
1591  fprintf(fp2, "%s\n", l2_str[ifile].filename);
1592 
1593  }
1594 
1595  // First argument, isnight, set to 0. We want to use the dateline all the time
1596  if (dnStatus == 0 && n[input.night] > 0) {
1597  get_coord_extrema(0, n[input.night], dnlat[input.night], dnlon[input.night],
1598  &onorth, &osouth, &owest, &oeast, &dateline);
1599  // Third argument, isnight, set to 0. We want to use the dateline all the time
1600  // Determine dataday0 and dataday1 based on dateline determined from get_coord_extrema
1601  get_datadays(l2_stimes[ifile], eqcross, input.night, dateline, owest, oeast,
1602  &dataday0, &dataday1);
1603  dataday0 += plusday;
1604  dataday1 += plusday;
1605  if (dataday1 == dataday0) {
1606  if (dataday0 < startdate || dataday1 > enddate)
1607  brk_scan[ifile] = -9999;
1608  else
1609  brk_scan[ifile] = 0; // -9999;
1610  } else {
1611  if (dataday1 < startdate) // startdate is dataday conversion of input sday
1612  brk_scan[ifile] = -9999;
1613  else if (dataday0 > enddate) { // enddate is dataday conversion of input eday
1614  brk_scan[ifile] = -9999;
1615  } else {
1616  if (dataday1 > enddate)
1617  brk_scan[ifile] = 1;
1618  else
1619  brk_scan[ifile] = -1;
1620 
1621  // file is within sday and eday, brk_scan tells whether
1622  // it went from the data day to the next day over the dateline (1),
1623  // or data day to a previous day over the dateline (-1)
1624  }
1625  }
1626 
1627  } else {
1628  brk_scan[ifile] = -9999;
1629  }
1630 
1631  for (i = 0; i < small_dims[0]; i++) {
1632  free(dorn[i]);
1633  }
1634  free(dorn);
1635  if (dnlat[0])
1636  free(dnlat[0]);
1637  if (dnlat[1])
1638  free(dnlat[1]);
1639  if (dnlon[0])
1640  free(dnlon[0]);
1641  if (dnlon[1])
1642  free(dnlon[1]);
1643 
1644  if (brk_scan[ifile] == -9999) n_active_files--;
1645 
1646  free(slon);
1647  free(elon);
1648  free(clon);
1649  free(elat);
1650  free(slat);
1651  free(clat);
1652 
1653  free(lonArray);
1654  free(latArray);
1655 
1656  for (i = 0; i < small_dims[0]; i++) {
1657  free(dolat[i]);
1658  free(dolon[i]);
1659  }
1660  free(dolat);
1661  free(dolon);
1662  free(years);
1663  free(days);
1664  free(msecs);
1665 
1666  } /* ifile loop */
1667  if (ret_status == 120)
1668  fclose(fp2);
1669 
1670  shape = new l3::L3ShapeIsine(nrows);
1671 
1672  /* Compute basebin array (Starting bin of each row [1-based]) */
1673  /* ---------------------------------------------------------- */
1674  basebin = (int64_t *) calloc(nrows + 1, sizeof(int64_t));
1675  for (i = 0; i <= nrows; i++) {
1676  basebin[i] = shape->getBaseBin(i);
1677  }
1678  basebin[nrows] = shape->getNumBins() + 1;
1679 
1680  printf("total number of bins: %ld\n", (long int) basebin[nrows] - 1);
1681 
1682  /*
1683  * Create output file
1684  */
1685  strcpy(buf, input.ofile);
1686  status = nc_create(buf, NC_NETCDF4, &fileid_w);
1687  check_err(status, __LINE__, __FILE__);
1688 
1689  status = nc_def_grp(fileid_w, "level-3_binned_data", &out_grp);
1690  check_err(status, __LINE__, __FILE__);
1691 
1692  idDS ds_id;
1693  ds_id.fid = fileid_w;
1694  ds_id.sid = -1;
1695  ds_id.fftype = DS_NCDF; // FMT_L2NCDF
1696 
1697  if (is64bit)
1698  status = defineBinList64_nc(input.deflate, out_grp);
1699  else
1700  status = defineBinList_nc(input.deflate, out_grp);
1701  if (status) {
1702  printf("-E- %s:%d Could not define binList variable in output file\n",
1703  __FILE__, __LINE__);
1704  exit(EXIT_FAILURE);
1705  }
1706 
1707  char** prodnames = (char**) malloc(l3_prodname.size() * sizeof(char*));
1708  for(size_t i=0; i<l3_prodname.size(); i++) {
1709  prodnames[i] = strdup(l3_prodname[i].c_str());
1710  }
1711  status = defineBinData_nc(input.deflate, out_grp, l3b_nprod, prodnames);
1712  if (status) {
1713  printf("-E- %s:%d Could not define binData variable in output file\n",
1714  __FILE__, __LINE__);
1715  exit(EXIT_FAILURE);
1716  }
1717  for(size_t i=0; i<l3_prodname.size(); i++) {
1718  free(prodnames[i]);
1719  }
1720  free(prodnames);
1721 
1722  if (is64bit)
1723  status = defineBinIndex64_nc(input.deflate, out_grp);
1724  else
1725  status = defineBinIndex_nc(input.deflate, out_grp);
1726  if (status) {
1727  printf("-E- %s:%d Could not define binIndex variable in output file\n",
1728  __FILE__, __LINE__);
1729  exit(EXIT_FAILURE);
1730  }
1731 
1732  if (input.qual_prod[0] != 0) {
1733  status = defineQuality_nc(input.deflate, out_grp);
1734  if (status) {
1735  printf("-E- %s:%d Could not define quality variable in output file\n",
1736  __FILE__, __LINE__);
1737  exit(EXIT_FAILURE);
1738  }
1739  }
1740 
1741  /* Allocate Arrays for Bin Index */
1742  /* ----------------------------- */
1743  beg = (int64_t *) calloc(nrows, sizeof(int64_t));
1744  ext = (int32_t *) calloc(nrows, sizeof(int32_t));
1745 
1746  if (is64bit)
1747  binIndex64nc = (binIndexStruct64_nc *) calloc(nrows, sizeof(binIndexStruct64_nc));
1748  else
1749  binIndex32nc = (binIndexStruct_nc *) calloc(nrows, sizeof(binIndexStruct_nc));
1750 
1751  /* Initialize bin_indx array */
1752  /* ------------------------- */
1753  for (i = 0; i < nrows; i++) {
1754 
1755  i32 = i;
1756 
1757  if (i32 < 0 || i32 >= nrows) {
1758  printf("%d %d\n", i, nrows);
1759  exit(-1);
1760  }
1761 
1762  if (is64bit) {
1763  binIndex64nc[i].begin = beg[i32];
1764  binIndex64nc[i].extent = ext[i32];
1765  binIndex64nc[i].max = shape->getNumCols(i32);
1766  } else {
1767  binIndex32nc[i].begin = beg[i32];
1768  binIndex32nc[i].extent = ext[i32];
1769  binIndex32nc[i].max = shape->getNumCols(i32);
1770  }
1771  }
1772 
1773  // Row Group
1774  n_rows_in_group = input.rowgroup;
1775  if (n_rows_in_group <= 0) {
1776  printf("row_group not defined, using 270.\n");
1777  n_rows_in_group = 270;
1778  }
1779  if (input.verbose)
1780  printf("%d %d %d\n", proc_day_beg, proc_day_end, n_rows_in_group);
1781 
1782  /* Find row_group that divides nrows */
1783  /* --------------------------------- */
1784  for (i = nrows; i > 0; i--) {
1785  if ((nrows % i) == 0) {
1786  if (i <= n_rows_in_group) {
1787  n_rows_in_group = i;
1788  break;
1789  }
1790  }
1791  }
1792  if (input.rowgroup != n_rows_in_group) {
1793  printf("Input row_group: %d Actual row_group: %d\n",
1794  input.rowgroup, n_rows_in_group);
1795  }
1796  ngroup = nrows / n_rows_in_group;
1797 
1798  /* Process each group of bin rows (Main Loop) */
1799  /* ========================================== */
1800  for (krow = 0, igroup = 0; igroup < ngroup; igroup++) {
1801 
1802  if (shape->row2lat(krow + n_rows_in_group) < input.latsouth) {
1803  krow += n_rows_in_group;
1804  continue;
1805  }
1806  if (shape->row2lat(krow) > input.latnorth) {
1807  krow += n_rows_in_group;
1808  continue;
1809  }
1810 
1811  /* Print info on rowgroup */
1812  /* ---------------------- */
1813  time(&tnow);
1814  tmnow = localtime(&tnow);
1815  printf("krow:%6d out of %6d (%6.2f to %6.2f) ",
1816  krow, nrows,
1817  ((float32) (krow) / nrows) * 180 - 90,
1818  ((float32) (krow + n_rows_in_group) / nrows) * 180 - 90);
1819  printf("%s", asctime(tmnow));
1820 
1821  n_bins_in_group = basebin[krow + n_rows_in_group] - basebin[krow];
1822  within_flag = 0;
1823 
1824  /* Determine relevant swath rows for this bin row group for each file */
1825  /* ------------------------------------------------------------------ */
1826  for (ifile = 0; ifile < nfiles; ifile++) {
1827 
1828  /* add an extra 0 to the end of scan_in_rowgroup so the caching
1829  * code never reads past the end of the file */
1830  scan_in_rowgroup[ifile] = (unsigned char *)
1831  calloc(l2_str[ifile].nrec + 1, sizeof (unsigned char));
1832 
1833  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++) {
1834  scan_in_rowgroup[ifile][jsrow] = 1;
1835  if (bscan_row[ifile][jsrow] < krow ||
1836  escan_row[ifile][jsrow] >= (krow + n_rows_in_group - 1)) {
1837  scan_in_rowgroup[ifile][jsrow] = 255;
1838  }
1839  } /* jsrow loop */
1840 
1841  /* Determine if within bin row group */
1842  /* --------------------------------- */
1843  if(!within_flag) {
1844  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++) {
1845  if (scan_in_rowgroup[ifile][jsrow] == 1) {
1846  within_flag = 1;
1847  break;
1848  }
1849  } /* scan row loop */
1850  }
1851 
1852  } /* ifile loop */
1853 
1854  /* If no swath rows within group then continue to next group */
1855  /* --------------------------------------------------------- */
1856  if (within_flag == 0) {
1857  for (ifile = 0; ifile < nfiles; ifile++) {
1858  if (scan_in_rowgroup[ifile] != NULL) {
1859  free(scan_in_rowgroup[ifile]);
1860  scan_in_rowgroup[ifile] = NULL;
1861  }
1862  }
1863  krow += n_rows_in_group;
1864  continue;
1865  }
1866 
1867  /* Allocate # pixels in bin, bin_flag, tilt, qual, & nscenes arrays */
1868  /* ---------------------------------------------------------------- */
1869  n_filled_bins = 0;
1870  bin_flag = (int32_t *) calloc(n_bins_in_group, sizeof(int32_t));
1871  tilt = (int16_t *) calloc(n_bins_in_group, sizeof(int16_t));
1872  qual = (int16_t *) calloc(n_bins_in_group, sizeof(int16_t));
1873  nscenes = (int16_t *) calloc(n_bins_in_group, sizeof(int16_t));
1874  lastfile = (int16_t *) calloc(n_bins_in_group, sizeof(int16_t));
1875 
1876  for (i = 0; i < n_bins_in_group; i++) {
1877  tilt[i] = -1;
1878  qual[i] = 3;
1879  lastfile[i] = -1;
1880  }
1881 
1882  /* Allocate bin accumulator & data value arrays */
1883  /* -------------------------------------------- */
1884  nobs = (int16_t *) calloc(n_bins_in_group, sizeof(int16_t));
1885  allocated_space = (int16_t *) calloc(n_bins_in_group, sizeof(int16_t));
1886  data_values = (float32 **) calloc(n_bins_in_group, sizeof(float32 *));
1887  data_areas = (double **) calloc(n_bins_in_group, sizeof(double *));
1888  file_index = (int16_t **) calloc(n_bins_in_group, sizeof(int16_t *));
1889  data_quality = (uint8_t **) calloc(n_bins_in_group, sizeof(uint8_t *));
1890  time_value = (float64 **) calloc(n_bins_in_group, sizeof(float64 *));
1891 
1892  /* Initialize bin counters */
1893  /* ----------------------- */
1894  n_allocperbin =
1895  n_active_files * l2_str[0].nrec * l2_str[0].nsamp / 50000000;
1896 
1897  if (n_allocperbin < 2)
1898  n_allocperbin = 2;
1901 
1902  if (input.verbose)
1903  printf("%-20s:%8d\n\n", "# allocated per bin", n_allocperbin);
1904 
1905  for (i = 0; i < n_bins_in_group; i++) {
1906  nobs[i] = 0;
1907  allocated_space[i] = 0;
1908  lastfile[i] = -1;
1909  }
1910 
1911  /* Read L2 files and fill data_values (L3b) array */
1912  /* ++++++++++++++++++++++++++++++++++++++++++++++ */
1913  for (ifile = 0; ifile < nfiles; ifile++) {
1914 
1916 
1917  /* if "early" or "late" input file then skip */
1918  /* ----------------------------------------- */
1919  if (brk_scan[ifile] == -9999) continue;
1920 
1921  /* if no scans in rowgroup for this file then skip */
1922  /* ----------------------------------------------- */
1923  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++) {
1924  if (scan_in_rowgroup[ifile][jsrow] == 1) {
1925  break;
1926  }
1927  }
1928  if (jsrow == l2_str[ifile].nrec) {
1929  continue;
1930  }
1931 
1932  reopenL2(ifile, &l2_str[ifile]);
1933 
1934  /* Get tilt flags & ranges */
1935  /* ----------------------- */
1936  ntilts = l2_str[ifile].ntilts;
1937  for (i = 0; i < ntilts; i++) {
1938  tilt_flags[i] = l2_str[ifile].tilt_flags[i];
1939  tilt_ranges[0][i] = l2_str[ifile].tilt_ranges[0][i];
1940  tilt_ranges[1][i] = l2_str[ifile].tilt_ranges[1][i];
1941  }
1942 
1943  /* Get date stuff */
1944  /* -------------- */
1945  date = (time_t) yds2unix(l2_str[ifile].syear, l2_str[ifile].sday,
1946  l2_str[ifile].smsec / 1000) / 86400;
1947  diffday_beg = date - startdate;
1948  diffday_end = date - enddate;
1949 
1950  /* Loop over swath rows */
1951  /* ^^^^^^^^^^^^^^^^^^^^ */
1952  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++) {
1953 
1954  /* if swath row not within group then continue */
1955  /* ------------------------------------------- */
1956  if (scan_in_rowgroup[ifile][jsrow] != 1) continue;
1957 
1958  /* Read swath record from L2 */
1959  /* ------------------------- */
1960  status = readL2(&l2_str[ifile], ifile, jsrow, -1,
1961  scan_in_rowgroup[ifile]);
1962 
1963  if (status == 5) continue;
1964  /*
1965  * The following bits were added to address the issue of orbit files
1966  * (e.g. SeaWiFS) that cross the dateline because they cross the pole
1967  * The intent is to prevent the dateline crossing code from being as
1968  * draconian with the data it excludes for the portion of the orbit
1969  * that doesn't cross the dateline.
1970  */
1971  dbldate = yds2unix(l2_str[ifile].year, l2_str[ifile].day,
1972  l2_str[ifile].msec / 1000.0);
1973  ddstime = yds2unix(syear, sday - plusday,
1974  ((double) eqcross - 12) * 3600.);
1975  ddetime = yds2unix(eyear, eday - plusday,
1976  ((double) eqcross + 12) * 3600.);
1977 
1978  int scan_crosses_dateline = 0;
1979  if (brk_scan[ifile] != 0) {
1980  int npixls = l2_str[ifile].nsamp - 1;
1981  float slat = l2_str[ifile].longitude[0];
1982  float elat = l2_str[ifile].longitude[npixls];
1983  if (abs(slat - elat) > 180)
1984  scan_crosses_dateline = 1;
1985  }
1986 
1987  if (scan_crosses_dateline == 0
1988  && (strcasecmp(input.prodtype, "regional") != 0)) {
1989  // the bounds of the dataday(s), move along...
1990  if (dbldate < ddstime) continue;
1991  if (dbldate > ddetime) continue;
1992  }
1993 
1994  /* Check tilt state */
1995  /* ---------------- */
1996  for (i = 0; i < ntilts; i++) {
1997  if ((jsrow + 1) <= tilt_ranges[1][i]) {
1998  tiltstate = (tilt_flags[i] & 0xFF);
1999  break;
2000  }
2001  }
2002 
2003  if (l2_str[ifile].nsamp == 0) continue;
2004 
2005  if ((jsrow % 100) == 0 && input.verbose) {
2006  printf("ifile:%4d jsrow:%6d nsamp:%8d\n", ifile, jsrow, l2_str[ifile].nsamp);
2007  }
2008 
2009  // save the tai93 date for this line, so is can be saved in the bins
2010  time_tai93 = unix_to_tai93(dbldate);
2011 
2012  /* ##### Loop over L2 pixels ##### */
2013  /* ------------------------------- */
2014  for (ipixl = 0; ipixl < l2_str[ifile].nsamp; ipixl++) {
2015 
2016 #ifdef DEBUG_PRINT
2017  if(ipixl>=510 && ipixl<512)
2018  enableDebugPrint = true;
2019  else
2020  enableDebugPrint = false;
2021 #endif
2022 
2023  /* if bin flagged then continue */
2024  /* ---------------------------- */
2025  flagcheck = (l2_str[ifile].l2_flags[ipixl]);
2026  if ((flagcheck & flagusemask) != 0)
2027  continue;
2028  if ((flagcheck & required) != required)
2029  continue;
2030 
2031  /* Check for dateline crossing */
2032  /* --------------------------- */
2033  if (scan_crosses_dateline == 1) {
2034  // if the scan does cross the dateline, decide which
2035  // pixels to bin based on longitude and which side of
2036  // the dateline we should keep
2037  if (input.night == 1) {
2038 
2039  if ((brk_scan[ifile] == -1) &&
2040  (diffday_beg == -1) &&
2041  (l2_str[ifile].longitude[ipixl] < 0)) continue;
2042 
2043  if ((brk_scan[ifile] == +1) &&
2044  (diffday_end == 0) &&
2045  (l2_str[ifile].longitude[ipixl] > 0)) continue;
2046 
2047  } else {
2048 
2049  if ((brk_scan[ifile] == -1) &&
2050  (diffday_beg <= 0) &&
2051  (l2_str[ifile].longitude[ipixl] < 0)) continue;
2052 
2053  if ((brk_scan[ifile] == +1) &&
2054  (diffday_end >= 0) &&
2055  (l2_str[ifile].longitude[ipixl] > 0)) continue;
2056  }
2057  }
2058  /* Check for bad value in any of the products */
2059  /* ------------------------------------------ */
2060  bad_value = 0;
2061  for (iprod = 0; iprod < l3b_nprod; iprod++) {
2062  int32_t l2_iprod = numer[ifile][iprod];
2063  for(int i=0; i<l2_str[ifile].thirdDim[l2_iprod]; i++) {
2064  f32 = l2_str[ifile].l2_data[l2_iprod][ipixl*l2_str[ifile].thirdDim[l2_iprod]+i];
2065  if (isnan(f32)) {
2066  bad_value = 1;
2067  break;
2068  }
2069  }
2070  if(bad_value)
2071  break;
2072  }
2073  if (bad_value == 1) continue;
2074 
2075  /* Check if within longitude boundaries */
2076  /* ------------------------------------ */
2077  if (input.lonwest != 0.0 || input.loneast != 0.0) {
2078  if (l2_str[ifile].longitude[ipixl] < input.lonwest) continue;
2079  if (l2_str[ifile].longitude[ipixl] > input.loneast) continue;
2080  }
2081 
2082  if(input.area_weighting) {
2083  /* Get map of bin numbers and intersection area that intersect pixel */
2084  std::map<uint64_t, double> areas;
2085  getBins(ifile, ipixl, areas);
2086  for(auto const& val : areas) {
2087  addPixelToBin(ifile, ipixl, val.first, val.second);
2088  }
2089  } else {
2090  /* Get Bin Number for Pixel */
2091  /* ------------------------ */
2092  bin = getbinnum(ifile, ipixl); // bin is 1-based
2093  addPixelToBin(ifile, ipixl, bin);
2094  }
2095 
2096  } /* ##### i (ipixl) loop ##### */
2097 
2098  } /* ^^^^^^^^^^ jsrow loop ^^^^^^^^^^ */
2099 
2100  closeL2(&l2_str[ifile], ifile);
2101 
2102 #ifdef MALLINFO
2103  if (input.meminfo) {
2104  int32_t total_alloc_space;
2105  /* malloc_stats();*/
2106  minfo = mallinfo();
2107  total_alloc_space = 0;
2108  for (i = 0; i < n_bins_in_group; i++) {
2109  total_alloc_space += allocated_space[i];
2110  }
2111  printf("Used space: %10d\n", minfo.uordblks);
2112  printf("Allo space: %10d\n", total_alloc_space * (2 + l3b_nprod * 4));
2113  }
2114 #endif
2115 
2116  } /* ++++++++++ ifile loop ++++++++++ */
2117 
2118  if (input.verbose) {
2119  time(&tnow);
2120  tmnow = localtime(&tnow);
2121  printf("krow:%5d After data_value fill: %s\n", krow, asctime(tmnow));
2122  }
2123 
2124  /* Compute Total # of filled bins */
2125  /* ------------------------------ */
2126  for (ibin = 0; ibin < n_bins_in_group; ibin++) {
2127  if (nobs[ibin] > 0 && nobs[ibin] < input.minobs)
2128  nobs[ibin] = 0;
2129 
2130  if (nobs[ibin] != 0) n_filled_bins++;
2131  } /* ibin loop */
2132 
2133  best_qual = (uint8_t *) calloc(n_bins_in_group, sizeof(uint8_t));
2134  memset(best_qual, 255, n_bins_in_group * sizeof(uint8_t));
2135 
2136  /* ********** If filled bins ********** */
2137  /* ------------------------------------ */
2138  if (n_filled_bins > 0) {
2139 
2140  /* Fill "Bin List" vdata array */
2141  /* --------------------------- */
2142  if (is64bit)
2143  binList64nc = (binListStruct64_nc*)
2144  calloc(n_filled_bins, sizeof(binListStruct64_nc));
2145  else
2146  binList32nc = (binListStruct_nc*)
2147  calloc(n_filled_bins, sizeof(binListStruct_nc));
2148 
2149  i = 0;
2150  for (ibin = 0; ibin < n_bins_in_group; ibin++) {
2151 
2152  /* Adjust for bins with "bad" quality values */
2153  /* ----------------------------------------- */
2154  if (input.qual_prod[0] != 0 && nobs[ibin] > 0) {
2155  best_qual[ibin] = 255;
2156  for (j = 0; j < nobs[ibin]; j++)
2157  if (data_quality[ibin][j] < best_qual[ibin])
2158  best_qual[ibin] = data_quality[ibin][j];
2159 
2160  k = 0;
2161  for (j = 0; j < nobs[ibin]; j++) {
2162  if ((data_quality[ibin][j] <= best_qual[ibin]) &&
2163  (data_quality[ibin][j] <= qual_max_allowed)) {
2164  if (k < j) {
2165  data_areas[ibin][k] = data_areas[ibin][j];
2166  for (iprod = 0; iprod < l3b_nprod; iprod++) {
2167  data_values[ibin][k * l3b_nprod + iprod] =
2168  data_values[ibin][j * l3b_nprod + iprod];
2169  }
2170  }
2171  k++;
2172  }
2173  }
2174  nobs[ibin] = k;
2175 
2176  if (nobs[ibin] == 0) n_filled_bins--;
2177  }
2178 
2179  if (nobs[ibin] != 0) {
2180  bin = (uint64_t) ibin + basebin[krow];
2181 
2182  if (is64bit) {
2183  binList64nc[i].binnum = bin;
2184  binList64nc[i].nobs = nobs[ibin];
2185  binList64nc[i].nscenes = nscenes[ibin];
2186  } else {
2187  binList32nc[i].binnum = bin;
2188  binList32nc[i].nobs = nobs[ibin];
2189  binList32nc[i].nscenes = nscenes[ibin];
2190  }
2191 
2192  /* weights {=sqrt(# of L2 files in given bin)} */
2193  /* ------------------------------------------- */
2194  wgt = 0.0;
2195  for (ifile = 0; ifile <= nfiles; ifile++) {
2196  double area = 0.0;
2197  for (j = 0; j < nobs[ibin]; j++) {
2198  if (file_index[ibin][j] == ifile)
2199  area += data_areas[ibin][j];
2200  }
2201  wgt += sqrt(area);
2202  }
2203 
2204  time_rec = 0.0;
2205  for (j = 0; j < nobs[ibin]; j++)
2206  time_rec += time_value[ibin][j];
2207 
2208  if (is64bit) {
2209  binList64nc[i].weights = wgt;
2210  binList64nc[i].time_rec = time_rec;
2211  } else {
2212  binList32nc[i].weights = wgt;
2213  binList32nc[i].time_rec = time_rec;
2214  }
2215 
2216  i++;
2217 
2218  /* Update Max/Min Lon/Lat */
2219  /* ---------------------- */
2220  shape->bin2latlon(bin, latbin, lonbin);
2221 
2222  if (latbin > northmost) northmost = latbin;
2223  if (latbin < southmost) southmost = latbin;
2224 
2225  float minNeg = lonRanges[0]; // min, max in Western hemisphere
2226  float maxNeg = lonRanges[1];
2227  float minPos = lonRanges[2]; // min, max in Eastern hemisphere
2228  float maxPos = lonRanges[3];
2229  if (lonbin < 0) { // Western hemisphere
2230  minNeg = fmin(minNeg, lonbin);
2231  maxNeg = fmax(maxNeg, lonbin);
2232 
2233  } else { // Eastern hemisphere
2234  minPos = fmin(minPos, lonbin);
2235  maxPos = fmax(maxPos, lonbin);
2236  }
2237  // Adjust east and west granule bounding coordinates
2238  float max_angle = 90.0;
2239  int8_t lon000 = ((minPos - maxNeg) < max_angle);
2240  int8_t lon180 = ((maxPos - minNeg) > (360.0 - max_angle));
2241  int8_t polar = (lon000 && lon180) ||
2242  (90.0 - fmax(northmost, -1 * southmost) <= FLT_EPSILON);
2243 
2244  if (minNeg >= maxNeg) { // Entirely in Eastern hemisphere
2245  eastmost = maxPos;
2246  westmost = minPos;
2247  } else if (minPos >= maxPos) { // Entirely in Western hemisphere
2248  eastmost = maxNeg;
2249  westmost = minNeg;
2250  } else if (polar) { // Polar granule
2251  eastmost = 180;
2252  westmost = -180;
2253  if (northmost > 0)
2254  northmost = 90; // North pole
2255  else
2256  southmost = -90; // South pole
2257  } else if (lon000) { // Granule crosses Longitude 0
2258  eastmost = maxPos;
2259  westmost = minNeg;
2260  } else if (lon180) { // Granule crosses Longitude 180
2261  eastmost = maxNeg;
2262  westmost = minPos;
2263  }
2264  lonRanges[0] = minNeg;
2265  lonRanges[1] = maxNeg;
2266  lonRanges[2] = minPos;
2267  lonRanges[3] = maxPos;
2268  } /* nobs[ibin] != 0 */
2269  } /* ibin loop */
2270 
2271  /* if no good obs left than bail */
2272  /* ----------------------------- */
2273  if (n_filled_bins == 0) goto freemem;
2274 
2275  /* Print info on filled row group */
2276  /* ------------------------------ */
2277  printf("%-20s:%8d\n", "# bins in row group", n_bins_in_group);
2278  printf("%-20s:%8d\n", "# filled bins", n_filled_bins);
2279 
2280  /* Write "Bin List" vdata */
2281  /* ---------------------- */
2282  if (is64bit) {
2283  writeBinList_nc(out_grp, n_filled_bins, (VOIDP) binList64nc);
2284  free(binList64nc);
2285  } else {
2286  writeBinList_nc(out_grp, n_filled_bins, (VOIDP) binList32nc);
2287  free(binList32nc);
2288  }
2289 
2290  /* Allocate sum & sum-squared arrays */
2291  /* --------------------------------- */
2292  sum_bin = (double *) calloc(n_bins_in_group, sizeof(double));
2293  sum2_bin = (double *) calloc(n_bins_in_group, sizeof(double));
2294 
2295  /* Loop over all L3 products to fill sum arrays */
2296  /* -------------------------------------------- */
2297  for (iprod = 0; iprod < l3b_nprod; iprod++) {
2298  memset(sum_bin, 0, n_bins_in_group * sizeof(double));
2299  memset(sum2_bin, 0, n_bins_in_group * sizeof(double));
2300 
2301  /* Process bins */
2302  /* ------------ */
2303  for (ibin = 0; ibin < n_bins_in_group; ibin++) {
2304  if (nobs[ibin] == 0) continue;
2305 
2306  /* Process data file by file */
2307  /* ------------------------- */
2308  int32_t npix_file; // for weighting spatial average
2309  double sum_file; // accumulators for each file
2310  double sum2_file;
2311  double area_file;
2312  int16_t prev_file; // index of previous file considered
2313 
2314  // initialize sum for this bin with first file
2315  j = 0;
2316  float32 pixval = data_values[ibin][j * l3b_nprod + iprod];
2317  double pixarea = data_areas[ibin][j];
2318  npix_file = 1;
2319  area_file = pixarea;
2320  sum_file = pixval * pixarea;
2321  sum2_file = pixval * pixarea * pixval * pixarea;
2322  prev_file = file_index[ibin][j];
2323 
2324  // add weighted data for rest of observations (files)
2325  for (j = 1; j < nobs[ibin]; j++) {
2326  pixval = data_values[ibin][j * l3b_nprod + iprod];
2327  pixarea = data_areas[ibin][j];
2328 
2329  if (file_index[ibin][j] == prev_file) { // same file
2330  npix_file += 1;
2331  area_file += pixarea;
2332  sum_file += pixval * pixarea;
2333  sum2_file += pixval * pixarea * pixval * pixarea;
2334 
2335  } else { // new file
2336 
2337  // finalize weighted sum for previous file
2338  sum_bin[ibin] += (sum_file / sqrt(area_file));
2339  sum2_bin[ibin] += (sum2_file / sqrt(area_file));
2340 
2341  // initialize sum for current file
2342  npix_file = 1;
2343  area_file = pixarea;
2344  sum_file = pixval * pixarea;
2345  sum2_file = pixval * pixarea * pixval * pixarea;
2346  prev_file = file_index[ibin][j];
2347  }
2348  } /* observation loop */
2349 
2350  // add data from final file
2351  sum_bin[ibin] += (sum_file / sqrt(area_file));
2352  sum2_bin[ibin] += (sum2_file / sqrt(area_file));
2353 
2354  } /* ibin loop */
2355 
2356  /* Write Product Vdatas */
2357  /* -------------------- */
2358  float *productData = (float *) calloc(2 * n_filled_bins, sizeof(float));
2359 
2360  /* Fill bin data array */
2361  /* ------------------- */
2362  i = 0;
2363  for (ibin = 0; ibin < n_bins_in_group; ibin++) {
2364  if (nobs[ibin] != 0) {
2365  productData[i * 2] = sum_bin[ibin];
2366  productData[i * 2 + 1] = sum2_bin[ibin];
2367  //memcpy(&productData[i * 8], &sum_bin[ibin], 4);
2368  //memcpy(&productData[i * 8 + 4], &sum2_bin[ibin], 4);
2369  i++;
2370  }
2371  }
2372 
2373  //char_ptr1 = strchr(prodname[iprod], '/');
2374  //if (char_ptr1 != NULL) *char_ptr1 = '_';
2375  writeBinData_nc(out_grp, i, iprod, productData);
2376  //if (char_ptr1 != NULL) *char_ptr1 = '/';
2377 
2378  free(productData);
2379 
2380  } /* iprod loop */
2381 
2382  /* Write Quality vdata */
2383  /* ------------------- */
2384  if (input.qual_prod[0] != 0) {
2385  uint8_t *qualityData = (uint8_t *) calloc(n_filled_bins, 1);
2386 
2387  i = 0;
2388  for (ibin = 0; ibin < n_bins_in_group; ibin++) {
2389  if (nobs[ibin] != 0) {
2390  memcpy(&qualityData[i], &best_qual[ibin], 1);
2391  i++;
2392  }
2393  }
2394 
2395  writeQuality_nc(out_grp, n_filled_bins, (VOIDP) qualityData);
2396  free(qualityData);
2397  }
2398 
2399  /* Free dynamic memory */
2400  /* ------------------- */
2401  if (sum_bin != NULL) free(sum_bin);
2402  if (sum2_bin != NULL) free(sum2_bin);
2403 
2404  /* Compute "begin" & "extent" vdata entries */
2405  /* ---------------------------------------- */
2406  binnum_data = (int64_t *) calloc(n_filled_bins, sizeof(int64_t));
2407 
2408  i = 0;
2409  for (ibin = 0; ibin < n_bins_in_group; ibin++) {
2410  if (nobs[ibin] != 0) {
2411  binnum_data[i] = (int64_t) ibin + basebin[krow];
2412 
2413  if (i < 0 || i >= n_filled_bins) {
2414  printf("Error: %d %d %d %d\n",
2415  i, ibin, n_filled_bins, n_bins_in_group);
2416  }
2417  i++;
2418  }
2419  }
2420 
2421  get_beg_ext(n_filled_bins, binnum_data, basebin, nrows, beg, ext);
2422 
2423  free(binnum_data);
2424 
2425  total_filled_bins += n_filled_bins;
2426 
2427  } /* ********** n_filled_bin > 0 ********** */
2428 
2429  free(best_qual);
2430 
2431  if (input.verbose) {
2432  time(&tnow);
2433  tmnow = localtime(&tnow);
2434  printf("krow:%5d After bin processing: %s", krow, asctime(tmnow));
2435  }
2436 
2437  /* Fill BinIndex Vdata */
2438  /* ------------------- */
2439  for (i = 0; i < n_rows_in_group; i++) {
2440 
2441  i32 = i + krow;
2442  if (i32 < 0 || i32 >= nrows) {
2443  printf("Error: %d %d\n", i, krow);
2444  exit(-1);
2445  }
2446 
2447  if (is64bit) {
2448  binIndex64nc[i + krow].start_num = basebin[i32];
2449  binIndex64nc[i + krow].begin = beg[i32];
2450  binIndex64nc[i + krow].extent = ext[i32];
2451  binIndex64nc[i + krow].max = shape->getNumCols(i32);
2452  } else {
2453  binIndex32nc[i + krow].start_num = basebin[i32];
2454  binIndex32nc[i + krow].begin = beg[i32];
2455  binIndex32nc[i + krow].extent = ext[i32];
2456  binIndex32nc[i + krow].max = shape->getNumCols(i32);
2457  }
2458 
2459  } /* row_group loop */
2460 
2461  /* End Bin Index Fill */
2462 
2463  freemem:
2464  /* Free Dynamic Memory */
2465  /* ------------------- */
2466  if (bin_flag != NULL) free(bin_flag);
2467  if (tilt != NULL) free(tilt);
2468  if (qual != NULL) free(qual);
2469  if (nscenes != NULL) free(nscenes);
2470  if (lastfile != NULL) free(lastfile);
2471 
2472  for (ifile = 0; ifile < nfiles; ifile++) {
2473  if (scan_in_rowgroup[ifile] != NULL) free(scan_in_rowgroup[ifile]);
2474  }
2475 
2476  if (input.verbose) {
2477  time(&tnow);
2478  tmnow = localtime(&tnow);
2479  printf("krow:%5d Befre free dynic mem: %s", krow, asctime(tmnow));
2480  }
2481 
2482  for (i = 0; i < n_bins_in_group; i++) {
2483  if (file_index[i] != NULL) free(file_index[i]);
2484  if (data_values[i] != NULL) free(data_values[i]);
2485  if (data_areas[i] != NULL) free(data_areas[i]);
2486  if (data_quality[i] != NULL) free(data_quality[i]);
2487  if (time_value[i] != NULL) free(time_value[i]);
2488  }
2489 
2490  if (input.verbose) {
2491  time(&tnow);
2492  tmnow = localtime(&tnow);
2493  printf("krow:%5d After free dynic mem: %s", krow, asctime(tmnow));
2494  }
2495 
2496  free(data_values);
2497  free(data_areas);
2498  free(data_quality);
2499  free(time_value);
2500  free(nobs);
2501  free(allocated_space);
2502  free(file_index);
2503 
2504  krow += n_rows_in_group;
2505 
2506  } /* ========== End krow (Main) loop ========== */
2507 
2508  if (input.verbose) {
2509  time(&tnow);
2510  tmnow = localtime(&tnow);
2511  printf("krow:%5d %s", krow, asctime(tmnow));
2512  }
2513  printf("total_filled_bins: %ld\n", (long int) total_filled_bins);
2514 
2515  if (total_filled_bins == 0) {
2516  strcpy(buf, "rm -f ");
2517  strcat(buf, input.ofile);
2518  strcat(buf, "*");
2519  printf("%s\n", buf);
2520  system(buf);
2521  ret_status = 110;
2522  exit(ret_status);
2523  }
2524 
2525  /* Write and Close BinIndex Vdata */
2526  /* ------------------------------ */
2527  if (is64bit)
2528  writeBinIndex_nc(out_grp, nrows, binIndex64nc);
2529  else
2530  writeBinIndex_nc(out_grp, nrows, binIndex32nc);
2531 
2532  /*
2533  * determine list of files actually used in the bin output
2534  */
2535  if (input.fileuse[0] != 0) {
2536  fp2 = fopen(input.fileuse, "w");
2537  for (ifile = 0; ifile < nfiles; ifile++) {
2538  if (brk_scan[ifile] != -9999) {
2539  fileused[ifile] = 1;
2540  fprintf(fp2, "%s\n", l2_str[ifile].filename);
2541  }
2542  }
2543  fclose(fp2);
2544  } else {
2545  for (ifile = 0; ifile < nfiles; ifile++) {
2546  if (brk_scan[ifile] != -9999) {
2547  fileused[ifile] = 1;
2548  }
2549  }
2550  }
2551 
2552  /* Read and write global attributes */
2553  /* -------------------------------- */
2554  if (input.verbose)
2555  printf("Writing Global Attributes\n");
2556 
2557  status = reopenL2(0, &l2_str[0]);
2558  status = readL2meta(&meta_l2, 0);
2559 
2560  strcpy(meta_l3b.product_name, input.ofile);
2561  if (meta_l2.sensor_name != NULL)
2562  strcpy(meta_l3b.sensor_name, meta_l2.sensor_name);
2563  strcpy(meta_l3b.sensor, meta_l2.sensor);
2564  meta_l3b.sensorID = sensorID[0];
2565  strcpy(meta_l3b.title, meta_l3b.sensor_name);
2566  strcat(meta_l3b.title, " Level-3 Binned Data");
2567  if (meta_l2.mission != NULL)
2568  strcpy(meta_l3b.mission, meta_l2.mission);
2569  strcpy(meta_l3b.prod_type, input.prodtype);
2570 
2571  strcat(meta_l3b.pversion, input.pversion);
2572  strcat(meta_l3b.soft_name, "l2bin");
2573  strcat(meta_l3b.soft_ver, VERSION);
2574 
2575  /*
2576  * loop through the input files
2577  *
2578  * set some ridiculous start and end times to get the ball rolling...
2579  * if this code is still in use in the 22nd century, damn we're good!
2580  * ...but this line will need tweaking.
2581  *
2582  */
2583  meta_l3b.end_orb = 0;
2584  meta_l3b.start_orb = 1000000;
2585  double startTime = yds2unix(2100, 1, 0);
2586  double endTime = yds2unix(1900, 1, 0);
2587  double tmpTime;
2588  strcpy(meta_l3b.infiles, basename(l2_str[0].filename));
2589  for (i = 0; i < nfiles; i++) {
2590  if (fileused[i] == 1) {
2591  // update orbit numbers
2592  if (l2_str[i].orbit > meta_l3b.end_orb)
2593  meta_l3b.end_orb = l2_str[i].orbit;
2594  if (l2_str[i].orbit < meta_l3b.start_orb)
2595  meta_l3b.start_orb = l2_str[i].orbit;
2596 
2597  // update start/end times
2598  tmpTime = yds2unix(l2_str[i].syear,
2599  l2_str[i].sday,
2600  (float64) (l2_str[i].smsec / 1000.0));
2601  if (tmpTime < startTime)
2602  startTime = tmpTime;
2603  tmpTime = yds2unix(l2_str[i].eyear,
2604  l2_str[i].eday,
2605  (float64) (l2_str[i].emsec / 1000.0));
2606  if (tmpTime > endTime)
2607  endTime = tmpTime;
2608  }
2609  // append to the file list - all files input, not just those used in the binning.
2610  if (i > 0) {
2611  strcat(meta_l3b.infiles, ",");
2612  strcat(meta_l3b.infiles, basename(l2_str[i].filename));
2613  }
2614  }
2615 
2616  meta_l3b.startTime = startTime;
2617  meta_l3b.endTime = endTime;
2618 
2619  strcpy(meta_l3b.binning_scheme, "Integerized Sinusoidal Grid");
2620 
2621  meta_l3b.north = northmost;
2622  meta_l3b.south = southmost;
2623  meta_l3b.east = eastmost;
2624  meta_l3b.west = westmost;
2625 
2626  strcpy(buf, argv[0]);
2627  for (i = 1; i < argc; i++) {
2628  strcat(buf, " ");
2629  strcat(buf, argv[i]);
2630  }
2631  strcpy(meta_l3b.proc_con, buf);
2632  strcpy(meta_l3b.input_parms, input.parms);
2633  strcpy(meta_l3b.flag_names, input.flaguse);
2634 
2635  buf[0] = 0;
2636  for (iprod = 0; iprod < l3b_nprod; iprod++) {
2637  char* prodName = strdup(l3_prodname[iprod].c_str());
2638  getL3units(&l2_str[0], 0, prodName, units);
2639  strcat(buf, prodName);
2640  free(prodName);
2641  strcat(buf, ":");
2642  strcat(buf, units);
2643  strcat(buf, ",");
2644  if (strlen(buf) >= MD_ATTRSZ) {
2645  printf("-E- units metadata is too long. Remove some products from product list.\n");
2646  exit(EXIT_FAILURE);
2647  }
2648  }
2649  buf[strlen(buf) - 1] = 0;
2650  strcpy(meta_l3b.units, buf);
2651 
2652  meta_l3b.data_bins = total_filled_bins;
2653  int64_t totalNumBins;
2654  totalNumBins = basebin[nrows] - 1;
2655  meta_l3b.pct_databins = (float) (total_filled_bins) / totalNumBins * 100.0;
2656 
2657  const char* doiStr = getDOI(meta_l3b.mission, meta_l3b.sensor,
2658  "L3B", input.suite, input.pversion);
2659  if (doiStr)
2660  strcpy(meta_l3b.doi, doiStr);
2661 
2662  const char* keywordStr = getGCMDKeywords(input.suite);
2663  if (keywordStr)
2664  strcpy(meta_l3b.keywords, keywordStr);
2665 
2666  time(&tnow);
2667  strcpy(meta_l3b.ptime, unix2isodate(tnow, 'G'));
2669  write_l3b_meta_netcdf4(ds_id, &meta_l3b, is64bit);
2670 
2671  freeL2meta(&meta_l2);
2672 
2673  /* Free Dynamic Memory */
2674  /* ------------------- */
2675  if (input.verbose)
2676  printf("Freeing Dynamic Memory\n");
2677  free(fileused);
2678 
2679  if (is64bit)
2680  free(binIndex64nc);
2681  else
2682  free(binIndex32nc);
2683 
2684  free(ext);
2685  free(beg);
2686 
2687  free(basebin);
2688 
2689  for (ifile = 0; ifile < nfiles; ifile++) {
2690  if (bscan_row[ifile] != NULL) free(bscan_row[ifile]);
2691  if (escan_row[ifile] != NULL) free(escan_row[ifile]);
2692 
2693  if (numer[ifile] != NULL) free(numer[ifile]);
2694  if (denom[ifile] != NULL) free(denom[ifile]);
2695  }
2696 
2697  /* Close L2 input files and free allocated L2 memory */
2698  /* ------------------------------------------------- */
2699  if (input.verbose)
2700  printf("Freeing L2 arrays\n");
2701  closeL2(&l2_str[0], 0);
2702  for (ifile = 0; ifile < nfiles; ifile++) {
2703  freeL2(&l2_str[ifile]);
2704  }
2705  freeL2(NULL);
2706 
2707  /* Close L3B output file */
2708  /* --------------------- */
2709  status = nc_close(fileid_w);
2710 
2711  return ret_status;
2712 }
2713 
2714 int64_t getbinnum(int32_t ifile, int32_t ipixl) {
2715 
2716  return shape->latlon2bin(l2_str[ifile].latitude[ipixl], l2_str[ifile].longitude[ipixl]);
2717 
2718 }
virtual int64_t getNumBins() const
Definition: L3Shape.cpp:63
@ L2PixelOff
Definition: readL2scan.h:154
int32_t openL2(const char *fname, char *plist, l2_prod *l2_str)
Definition: readL2scan.c:296
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
char product_name[SM_ATTRSZ]
Definition: meta_l3b.h:16
int16_t qual_prod_index[MAXNFILES]
Definition: l2bin.cpp:100
#define MAXALLOCPERBIN
Definition: l2bin.cpp:59
virtual void constrainRowCol(int32_t &row, int32_t &col) const =0
virtual int64_t latlon2bin(double lat, double lon) const =0
#define OLCIS3A
Definition: sensorDefs.h:32
int16 eday
Definition: l1_czcs_hdf.c:17
int32_t reopenL2(int32_t fileindex, l2_prod *l2_str)
Definition: readL2scan.c:975
char units[MD_ATTRSZ]
Definition: meta_l3b.h:28
char binning_scheme[SM_ATTRSZ]
Definition: meta_l3b.h:48
void freeProductInfo(productInfo_t *info)
int j
Definition: decode_rs.h:73
const char * getGCMDKeywords(const char *suite)
bg::model::box< Point_t > Box_t
Definition: l2bin.cpp:110
int32_t day
int status
Definition: l1_czcs_hdf.c:32
double yds2unix(int16_t year, int16_t day, double secs)
Definition: yds2unix.c:7
int instrumentPlatform2SensorId(const char *instrument, const char *platform)
Definition: sensorInfo.c:302
virtual int32_t getNumCols(int32_t row) const =0
void check_err(const int stat, const int line, const char *file)
Definition: nc4utils.c:35
int main(int argc, char **argv)
Definition: l2bin.cpp:876
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:99
These are used to scale the SD before writing it to the HDF4 file The default is and which means the product is not scaled at all Since the product is usually stored as a float inside of this is a way to write the float out as a integer l2prod min
int defineBinList64_nc(int32_t deflate, int32_t grpid)
Definition: ncdfbin_utils.c:99
*********************************************HISTORY FILE for MOD_PR02AQUA **Version ****Point of no algorithm change is made in this so the and when the blackbody temperature is above a threshold Since the LWIR FPA temperature oscillates on orbit
Definition: HISTORY.txt:106
@ L2PixelCorner
Definition: readL2scan.h:154
char infiles[LG_ATTRSZ]
Definition: meta_l3b.h:38
real(single), dimension(:,:), allocatable longitude
bool getBinsFromRow(double lat, double lon, T &pixelPoly, std::map< uint64_t, double > &areas)
Definition: l2bin.cpp:654
int32_t * bin_flag
Definition: l2bin.cpp:85
int defineBinData_nc(int32_t deflate, int32_t grpid, int32_t nprods, char **prodnames)
int16_t * qual
Definition: l2bin.cpp:86
#define NULL
Definition: decode_rs.h:63
char mission[SM_ATTRSZ]
Definition: meta_l3b.h:21
int32_t readL2meta(meta_l2Type *meta_l2, int32_t ifile)
Definition: readL2scan.c:2081
bg::model::point< double, 2, bg::cs::geographic< bg::degree > > Point_t
Definition: l2bin.cpp:108
#define LTILT_DIMS_2
Definition: l2bin.cpp:58
#define VIIRSN
Definition: sensorDefs.h:23
char keywords[SM_ATTRSZ]
Definition: meta_l3b.h:55
bg::model::polygon< Point_t > Polygon_t
Definition: l2bin.cpp:109
int l2bin_input(int argc, char **argv, instr *input, const char *prog, const char *version)
real(single), dimension(:,:), allocatable latitude
int32_t start_orb
Definition: meta_l3b.h:42
char input_parms[LG_ATTRSZ]
Definition: meta_l3b.h:36
const char * sensorId2PlatformName(int sensorId)
Definition: sensorInfo.c:226
#define MERIS
Definition: sensorDefs.h:22
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 wavelength
float * lat
int32 * msec
Definition: l1_czcs_hdf.c:31
#define MODIST
Definition: sensorDefs.h:18
float tm[MODELMAX]
int writeBinIndex_nc(int32_t grpid, int32_t n_write, const void *data)
#define PROGRAM
Definition: l2bin.cpp:113
int16 eyear
Definition: l1_czcs_hdf.c:17
virtual int32_t lat2row(double lat) const =0
int syear
Definition: l1_czcs_hdf.c:15
virtual void rowcol2bounds(int32_t row, int32_t col, double &north, double &south, double &east, double &west) const =0
string prodName
Definition: l3mapgen.cpp:55
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
int32_t freeL2meta(meta_l2Type *meta_l2)
Definition: readL2scan.c:2255
int32_t tiltstate
Definition: l2bin.cpp:93
char flag_names[SM_ATTRSZ]
Definition: meta_l3b.h:37
int32_t readL2(l2_prod *l2_str, int32_t ifile, int32_t recnum, int32_t iprod, unsigned char *scan_in_rowgroup)
Definition: readL2scan.c:1250
ds_format_t fftype
Definition: dfutils.h:31
int32 smsec
Definition: l1_czcs_hdf.c:16
@ string
int32_t closeL2(l2_prod *l2_str, int32_t ifile)
Definition: readL2scan.c:1936
int64_t get_beg_ext(int32_t n_bins_write, int64_t *binnum_data, int64_t *basebin, int32_t nrows, int64_t *beg, int32_t *ext)
Definition: get_beg_ext.c:83
int32_t nobs
Definition: atrem_cor.h:93
int writeBinList_nc(int32_t grpid, int32_t nbins_to_write, const void *data)
int sday
Definition: l1_czcs_hdf.c:15
@ L2PixelDelta
Definition: readL2scan.h:154
void init_rowgroup_cache()
Definition: readL2scan.c:265
int setlinebuf(FILE *stream)
int16_t composite_l3prod_index
Definition: l2bin.cpp:102
int defineQuality_nc(int32_t deflate, int32_t grpid)
int16_t composite_prod_index[MAXNFILES]
Definition: l2bin.cpp:101
char sensor[SM_ATTRSZ]
Definition: meta_l3b.h:23
virtual void latlon2rowcol(double lat, double lon, int32_t &row, int32_t &col) const =0
double unix_to_tai93(double unixtime)
Definition: yds2tai.c:11
int32_t krow
Definition: l2bin.cpp:90
void enableL2PixelArea(enum L2PixelMode_t val)
Definition: readL2scan.c:193
double endTime
Definition: meta_l3b.h:40
bool binIntersectsPixel(int32_t row, int32_t col, Box_t &pixelBox, double &areaFrac)
Definition: l2bin.cpp:584
void setupflags(char *flagdef, char *flaguse, uint32_t *flagusemask, uint32_t *required, int *status)
Definition: setupflags.c:5
productInfo_t * allocateProductInfo()
int32_t n_bins_in_group
Definition: l2bin.cpp:88
char * strdup(const char *)
char prod_type[SM_ATTRSZ]
Definition: meta_l3b.h:29
void getBins(int32_t ifile, int32_t ipixl, std::map< uint64_t, double > &areas)
Definition: l2bin.cpp:690
int32_t n_rows_in_group
Definition: l2bin.cpp:89
bg::model::point< double, 2, bg::cs::spherical_equatorial< bg::degree > > Point_t
Definition: get_dataday.cpp:23
#define HAWKEYE
Definition: sensorDefs.h:39
void free_rowgroup_cache()
Definition: readL2scan.c:236
void cdata_()
virtual int64_t getBaseBin(int32_t row) const =0
double resolution
Definition: meta_l3b.h:53
int defineBinIndex_nc(int32_t deflate, int32_t grpid)
virtual void bin2latlon(int64_t bin, double &lat, double &lon)
Definition: L3Shape.cpp:99
int64_t getbinnum(int32_t ifile, int32_t ipixl)
Definition: l2bin.cpp:2714
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
integer, parameter double
float north
Definition: meta_l3b.h:49
@ DS_NCDF
Definition: dfutils.h:20
char title[SM_ATTRSZ]
Definition: meta_l3b.h:17
virtual int64_t rowcol2bin(int32_t row, int32_t col) const =0
#define MAXNFILES
Definition: l3bin.cpp:20
vector< int32_t > thirdDimId
Definition: l2bin.cpp:105
int writeQuality_nc(int32_t grpid, int32_t nbins_to_write, const void *data)
float west
Definition: meta_l3b.h:52
int32_t sid
Definition: dfutils.h:30
int sensorID
Definition: meta_l3b.h:18
int32_t n_allocperbin
Definition: l2bin.cpp:74
#define basename(s)
Definition: l0chunk_modis.c:29
int findProductInfo(const char *productName, int sensorId, productInfo_t *info)
char doi[SM_ATTRSZ]
Definition: meta_l3b.h:54
int writeBinData_nc(int32_t grpid, int32_t nbins_to_write, int32_t iprod, const void *data)
const char * getDOI(const char *platform, const char *sensor, const char *level, const char *suite, const char *version)
Definition: getDOI.cpp:16
int32 emsec
Definition: l1_czcs_hdf.c:18
int defineBinList_nc(int32_t deflate, int32_t grpid)
Definition: ncdfbin_utils.c:42
int32_t l3b_nprod
Definition: l2bin.cpp:94
virtual int32_t getNumRows() const
Definition: L3Shape.cpp:55
#define BAD_FLT
Definition: jplaeriallib.h:19
#define OCTS
Definition: sensorDefs.h:14
float south
Definition: meta_l3b.h:50
subroutine geometry
void get_datadays(time_t starttime, float equatorialCrossingTime, DL dateLineCrossed, float west, float east, int32_t *dataday0, int32_t *dataday1)
Definition: get_dataday.cpp:65
int16_t * numer[MAXNFILES]
Definition: l2bin.cpp:99
char proc_con[MD_ATTRSZ]
Definition: meta_l3b.h:35
#define MD_ATTRSZ
Definition: meta_l3b.h:8
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
double startTime
Definition: meta_l3b.h:39
char soft_ver[SM_ATTRSZ]
Definition: meta_l3b.h:33
int32_t fid
Definition: dfutils.h:29
char * unix2isodate(double dtime, char zone)
Definition: unix2isodate.c:10
#define CZCS
Definition: sensorDefs.h:17
float pct_databins
Definition: meta_l3b.h:47
int16_t * nscenes
Definition: l2bin.cpp:86
virtual double row2lat(int32_t row) const =0
Definition: dfutils.h:28
vector< float32 > min_value
Definition: l2bin.cpp:106
float * lon
char ptime[SM_ATTRSZ]
Definition: meta_l3b.h:34
int defineBinIndex64_nc(int32_t deflate, int32_t grpid)
int16_t * allocated_space
Definition: l2bin.cpp:75
data_t s[NROOTS]
Definition: decode_rs.h:75
char pversion[SM_ATTRSZ]
Definition: meta_l3b.h:30
int write_l3b_meta_netcdf4(idDS ds_id, meta_l3bType *meta_l3b, int write64bit)
#define OLCIS3B
Definition: sensorDefs.h:41
int32_t getL3units(l2_prod *l2_str, int32_t ifile, char *l3b_prodname, char *units)
Definition: readL2scan.c:2274
float east
Definition: meta_l3b.h:51
#define HICO
Definition: sensorDefs.h:25
int64_t data_bins
Definition: meta_l3b.h:46
bg::model::polygon< Point_t > Polygon_t
Definition: get_dataday.cpp:25
void get_coord_extrema(int isnight, int n, float *olat, float *olon, float *north, float *south, float *west, float *east, int *dateline)
Definition: dataday.c:65
#define SEAWIFS
Definition: sensorDefs.h:12
int16_t * tilt
Definition: l2bin.cpp:86
void addPixelToBin(int32_t ifile, int32_t ipixl, uint64_t bin, double areaFrac=1.0)
Definition: l2bin.cpp:355
char sensor_name[SM_ATTRSZ]
Definition: meta_l3b.h:19
#define abs(a)
Definition: misc.h:90
int i
Definition: decode_rs.h:71
msiBandIdx val
Definition: l1c_msi.cpp:34
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:97
int32_t freeL2(l2_prod *l2_str)
Definition: readL2scan.c:2026
#define MODISA
Definition: sensorDefs.h:19
#define LG_ATTRSZ
Definition: meta_l3b.h:9
int32_t end_orb
Definition: meta_l3b.h:43
#define VIIRSJ1
Definition: sensorDefs.h:37
#define VERSION
Definition: l2bin.cpp:112
char soft_name[SM_ATTRSZ]
Definition: meta_l3b.h:32
int k
Definition: decode_rs.h:73
double str2resolution(char const *resolutionStr)
Pixel resolution string to meters.
#define MTILT_DIMS_2
Definition: l2bin.cpp:57
float p[MODELMAX]
Definition: atrem_corl1.h:131
float32 f32
Definition: l2bin.cpp:104
int daynight_outlines(int32_t *year, int32_t *dayOfYear, int32_t *msecondOfDay, int32_t wid, int32_t hgt, float **lat, float **lon, int32_t n[2], float *olat[2], float *olon[2], int8_t **dorn)
Definition: dataday.c:244
bg::model::box< Point_t > Box_t
Definition: l1c.cpp:73
l2prod max
int16_t * lastfile
Definition: l2bin.cpp:86
int32_t get_l2prod_index(const l2_prod &l2, const char *prodname)
Definition: l2bin.cpp:345