NASA Logo
Ocean Color Science Software

ocssw V2022
l2bin.cpp
Go to the documentation of this file.
1 #include <cstdlib>
2 #include <cstdint>
3 #include <libgen.h>
4 #include <cfloat>
5 #include <cmath>
6 #include <cstdio>
7 #include <cstring>
8 #include <ctime>
9 #include <iostream>
10 #include <vector>
11 #include <algorithm>
12 #include <string>
13 #include <netcdf>
14 #include <boost/algorithm/string.hpp>
15 #include <boost/geometry.hpp>
16 #include <boost/geometry/geometries/polygon.hpp>
17 #include <boost/geometry/geometries/box.hpp>
18 #include <boost/foreach.hpp>
19 #include <boost/algorithm/string/trim.hpp>
20 #include <get_dataday.hpp>
21 #include <meta_l3b.h>
22 #include <seaproto.h>
23 #include <readL2scan.h>
24 #include <timeutils.h>
25 #include <genutils.h>
26 #include <setupflags.h>
27 #include <sensorInfo.h>
28 #include <nc4utils.h>
29 #include <ncdfbin_utils.h>
30 #include <chrono>
31 #include <hdf.h>
32 #include "l2bin_input.h"
33 #include "L3Shape.h"
34 #include "L3ShapeIsine.h"
35 #include <get_geospatial.hpp>
36 #include "expand3D.hpp"
37 #include "find_variable.hpp"
38 
39 // using namespace std;
40 // using namespace netCDF;
41 // using namespace netCDF::exceptions;
42 
43 namespace bg = boost::geometry;
44 
45 // #define DEBUG_PRINT
46 // #define DEBUG_PRINT2
47 // extern void cdata_();
48 #ifdef DEBUG_PRINT
49 bool enableDebugPrint = true;
50 #endif
51 
52 #define MTILT_DIMS_2 20
53 #define LTILT_DIMS_2 2
54 #define MAXALLOCPERBIN 20
55 
56 
57 static constexpr int max_l3b_products = MAXNUMBERPRODUCTS;
58 
59 /* Global variables */
60 static instr input;
61 static l2_prod l2_str[MAXNFILES];
62 // bin file shape
63 static int32_t nrows = -1;
64 static l3::L3Shape *shape;
65 static int64_t *basebin;
66 
67 // global variables for the group
68 int32_t n_allocperbin; // number of obs to allocate to make another chunk
69 int16_t *allocated_space; // number of obs allocated for each bin
70 static int16_t *nobs; // number of observations in each bin
71 
72 static float32 **data_values;
73 static double **data_areas; // in fractions of a bin
74 static int16_t **file_index;
75 static uint8_t **data_quality;
76 static float64 **time_value;
77 static double time_tai93; // TAI 93 time for the current L2 line
78 
79 int32_t *bin_flag;
80 int16_t *tilt, *qual, *nscenes, *lastfile;
81 
82 int32_t n_bins_in_group; // number of bins in group
83 int32_t n_rows_in_group = -1; // number of rows in group
84 int32_t krow; // first bin row number of group
85 
86 // other globals
87 int32_t tiltstate = 0;
88 int32_t l3b_nprod; // number of products in bin file
89 static char buf[LG_ATTRSZ];
90 
92 
97 
98 float32 f32;
99 std::vector<int32_t> thirdDimId;
100 std::vector<float32> min_value, max_value;
101 
102 typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree>> Point_t;
103 typedef bg::model::polygon<Point_t> Polygon_t;
104 typedef bg::model::box<Point_t> Box_t;
105 
106 #define VERSION "7.1.0"
107 #define PROGRAM "l2bin"
108 
109 int32_t get_l2prod_index(const l2_prod &l2, const char *prodname)
110 {
111  int32_t index;
112  for (index = 0; index < l2.nprod; index++)
113  if (strcmp(prodname, l2.prodname[index]) == 0)
114  break;
115  if (index == l2.nprod)
116  index = -1;
117  return index;
118 }
119 
120 void addPixelToBin(int32_t ifile, int32_t ipixl, uint64_t bin, bool is_l2_flags_defined = true, double areaFrac = 1.0)
121 {
122 
123  int32_t ibin = bin - basebin[krow];
124 
125  /* if bin not within bin row group then continue */
126  /* --------------------------------------------- */
127  if (ibin < 0 ||
128  ibin >= n_bins_in_group ||
129  (int64_t)bin >= basebin[krow + n_rows_in_group] ||
130  (int64_t)bin < basebin[krow])
131  return;
132 
133  /* GOOD OBSERVATION FOUND */
134  /* ---------------------- */
135 #ifdef MALLINFO
136  if (input.dcinfo)
137  {
138  if ((l2_str[ifile].longitude[ipixl] <= -160) ||
139  (l2_str[ifile].longitude[ipixl] >= +160))
140  {
141  printf("DC: %10ld %12d %8.2f %8.2f\n",
142  (long int)bin,
143  (int)dbldate,
144  l2_str[ifile].longitude[ipixl],
145  l2_str[ifile].latitude[ipixl]);
146  }
147  }
148 #endif
149 
150  /* "OR" flags in swath pixel & set tilt & increment nscenes */
151  /* -------------------------------------------------------- */
152  if (is_l2_flags_defined)
153  bin_flag[ibin] = bin_flag[ibin] | l2_str[ifile].l2_flags[ipixl];
154  tilt[ibin] = tiltstate;
155  if (ifile != lastfile[ibin])
156  {
157  nscenes[ibin]++;
158  lastfile[ibin] = ifile;
159  }
160 
161  /* Allocate space for file index & bin data values */
162  /* ----------------------------------------------- */
163  if (file_index[ibin] == NULL)
164  {
165  file_index[ibin] = (int16_t *)
166  calloc(n_allocperbin, sizeof(int16_t));
167 
168  data_values[ibin] = (float32 *)
169  calloc(n_allocperbin * l3b_nprod, sizeof(float32));
170  if (data_values[ibin] == 0x0)
171  {
172  perror(buf);
173  printf("Allocation failed for data_values[ibin]: %d %s\n",
174  ibin, buf);
175  exit(EXIT_FAILURE);
176  }
177 
178  data_areas[ibin] = (double *)calloc(n_allocperbin, sizeof(double));
179  if (data_areas[ibin] == 0x0)
180  {
181  perror(buf);
182  printf("Allocation failed for data_areas[ibin]: %d %s\n",
183  ibin, buf);
184  exit(EXIT_FAILURE);
185  }
186 
187  data_quality[ibin] = (uint8_t *)
188  calloc(n_allocperbin, sizeof(uint8_t));
189  if (data_quality[ibin] == 0x0)
190  {
191  perror(buf);
192  printf("Allocation failed for data_quality[ibin]: %d %s\n",
193  ibin, buf);
194  exit(EXIT_FAILURE);
195  }
196 
197  time_value[ibin] = (float64 *)
198  calloc(n_allocperbin, sizeof(float64));
199  if (time_value[ibin] == 0x0)
200  {
201  perror(buf);
202  printf("Allocation failed for time_value[ibin]: %d %s\n",
203  ibin, buf);
204  exit(EXIT_FAILURE);
205  }
206 
208  }
209 
210  /* Set file_index for each observation */
211  /* ----------------------------------- */
212  file_index[ibin][nobs[ibin]] = ifile;
213 
214  /* Get data quality */
215  /* ---------------- */
216  if (input.qual_prod[0] != 0)
217  {
218  data_quality[ibin][nobs[ibin]] = l2_str[ifile].l2_data[qual_prod_index[ifile]][ipixl];
219  // a temporary fix to mask the last 4 pixels of a MODIS Terra scan
220  // for SST product
221  if (strncmp(input.suite, "SST", 3) == 0 && sensorID[ifile] == MODIST && ipixl > 1349)
222  {
223  data_quality[ibin][nobs[ibin]] = 3;
224  }
225  }
226 
227  /* Store time_value (TAI93) */
228  /* ---------------------- */
229  time_value[ibin][nobs[ibin]] = time_tai93;
230 
231  /* Get composite data */
232  /* ------------------- */
233  if (input.composite_prod[0] != 0)
234  {
235  int idx = composite_prod_index[ifile];
236  f32 = l2_str[ifile].l2_data[idx][ipixl];
237  if (f32 == -32767)
238  return;
239 
240  if (nobs[ibin] != 0)
241  {
242  if (strcmp(input.composite_scheme, "max") == 0)
243  {
244  if (f32 < data_values[ibin][composite_l3prod_index])
245  return;
246  }
247  else
248  {
249  if (f32 > data_values[ibin][composite_l3prod_index])
250  return;
251  }
252  nobs[ibin] = 0;
253  }
254  }
255 
256  /* Get data area for pixel */
257  /* ----------------------------------- */
258  data_areas[ibin][nobs[ibin]] = areaFrac;
259 
260  /* Get data values for all L3 products */
261  /* ----------------------------------- */
262  for (int32_t l3_iprod = 0; l3_iprod < l3b_nprod; l3_iprod++)
263  {
264 
265  int32_t l2_iprod = numer[ifile][l3_iprod];
266  f32 = l2_str[ifile].l2_data[l2_iprod][ipixl * l2_str[ifile].thirdDim[l2_iprod] + thirdDimId[l3_iprod]]; // probably the best way is to pass wavelength
267 
268  /* Set -32767 value to "bad" quality */
269  if (f32 == -32767)
270  if (input.qual_prod[0] != 0)
271  data_quality[ibin][nobs[ibin]] = 4;
272 
273  if (input.composite_prod[0] != 0)
274  {
275  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = f32;
276  }
277  else
278  {
279  if (denom[ifile][l3_iprod] == -1 && f32 >= min_value[l3_iprod] && f32 <= max_value[l3_iprod])
280  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = f32;
281 
282  if (denom[ifile][l3_iprod] == -1 && f32 < min_value[l3_iprod])
283  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = min_value[l3_iprod];
284 
285  if (denom[ifile][l3_iprod] == -1 && f32 > max_value[l3_iprod])
286  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = max_value[l3_iprod];
287 
288  if (denom[ifile][l3_iprod] == -2 && is_l2_flags_defined)
289  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] =
290  (l2_str[ifile].l2_flags[ipixl] >>
291  numer[ifile][l3_iprod]) &
292  1;
293  }
294 
295  /* ratio product */
296  /* ------------- */
297  if (denom[ifile][l3_iprod] >= 0)
298  {
299 
300  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] = f32;
301 
302  l2_iprod = denom[ifile][l3_iprod];
303  f32 = l2_str[ifile].l2_data[l2_iprod][ipixl * l2_str[ifile].thirdDim[l2_iprod] + thirdDimId[l3_iprod]];
304 
305  if (f32 >= min_value[l3_iprod] && f32 <= max_value[l3_iprod])
306  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] /= f32;
307  else if (f32 < min_value[l3_iprod])
308  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] /= min_value[l3_iprod];
309  else
310  data_values[ibin][l3b_nprod * nobs[ibin] + l3_iprod] /= max_value[l3_iprod];
311  }
312 
313  } /* iprod loop */
314 
315  /* Increment number of observations in bin */
316  /* --------------------------------------- */
317  nobs[ibin]++;
318  if (input.composite_prod[0] != 0)
319  nobs[ibin] = 1;
320 
321  /* Reallocate if necessary */
322  /* ----------------------- */
323  if (nobs[ibin] == allocated_space[ibin])
324  {
325 
326  file_index[ibin] = (int16_t *)
327  realloc(file_index[ibin],
328  (nobs[ibin] + n_allocperbin) * sizeof(int16_t));
329 
330  data_values[ibin] = (float32 *)
331  realloc(data_values[ibin],
332  (nobs[ibin] + n_allocperbin) * l3b_nprod * sizeof(float32));
333  if (data_values[ibin] == 0x0)
334  {
335  perror(buf);
336  printf("Reallocation failed for data_values[ibin]: %d %s\n",
337  ibin, buf);
338  exit(EXIT_FAILURE);
339  }
340 
341  data_areas[ibin] = (double *)realloc(data_areas[ibin],
342  (nobs[ibin] + n_allocperbin) * sizeof(double));
343  if (data_areas[ibin] == 0x0)
344  {
345  perror(buf);
346  printf("Reallocation failed for data_areas[ibin]: %d %s\n",
347  ibin, buf);
348  exit(EXIT_FAILURE);
349  }
350 
351  data_quality[ibin] = (uint8_t *)
352  realloc(data_quality[ibin],
353  (nobs[ibin] + n_allocperbin) * sizeof(uint8_t));
354  if (data_quality[ibin] == 0x0)
355  {
356  perror(buf);
357  printf("Reallocation failed for data_quality[ibin]: %d %s\n",
358  ibin, buf);
359  exit(EXIT_FAILURE);
360  }
361 
362  time_value[ibin] = (float64 *)
363  realloc(time_value[ibin],
364  (nobs[ibin] + n_allocperbin) * sizeof(float64));
365  if (time_value[ibin] == 0x0)
366  {
367  perror(buf);
368  printf("Reallocation failed for time_value[ibin]: %d %s\n",
369  ibin, buf);
370  exit(EXIT_FAILURE);
371  }
372 
374  } /* end reallocate */
375 }
376 
377 bool binIntersectsPixel(int32_t row, int32_t col, Box_t &pixelBox, double &areaFrac)
378 {
379  bool result = false;
380  double n, s, e, w;
381  shape->rowcol2bounds(row, col, n, s, e, w);
382 
383 #ifdef DEBUG_PRINT2
384  // debug plot the bin
385  if (enableDebugPrint)
386  {
387  printf("lat = [%f, %f, %f, %f, %f]\n", s, n, n, s, s);
388  printf("lon = [%f, %f, %f, %f, %f]\n", w, w, e, e, w);
389  printf("plt.plot(lon, lat)\n\n");
390  }
391 #endif
392 
393  Box_t box(Point_t(w, s), Point_t(e, n));
394  areaFrac = 0;
395 
396  if (!bg::disjoint(box, pixelBox))
397  {
398  Box_t output;
399  if (bg::intersection(box, pixelBox, output))
400  {
401  double intersectArea = bg::area(output);
402  if (intersectArea > 0)
403  {
404  result = true;
405  double binArea = (n - s) * (e - w);
406  areaFrac = intersectArea / binArea;
407  }
408  }
409  }
410  return result;
411 }
412 
413 bool binIntersectsPixel(int32_t row, int32_t col, Polygon_t &pixelPoly, double &areaFrac)
414 {
415  bool result = false;
416  double n, s, e, w;
417  shape->rowcol2bounds(row, col, n, s, e, w);
418 
419 #ifdef DEBUG_PRINT2
420  // debug plot the bin
421  if (enableDebugPrint)
422  {
423  printf("lat = [%f, %f, %f, %f, %f]\n", s, n, n, s, s);
424  printf("lon = [%f, %f, %f, %f, %f]\n", w, w, e, e, w);
425  printf("plt.plot(lon, lat)\n\n");
426  }
427 #endif
428 
429  Box_t box(Point_t(w, s), Point_t(e, n));
430  areaFrac = 0;
431 
432  if (!bg::disjoint(box, pixelPoly))
433  {
434  std::deque<Polygon_t> output;
435  if (bg::intersection(box, pixelPoly, output))
436  {
437  double binArea = (n - s) * (e - w);
438  BOOST_FOREACH (Polygon_t const &p, output)
439  {
440  double intersectArea = bg::area(p);
441  if (intersectArea > 0.0)
442  {
443  result = true;
444  areaFrac += intersectArea / binArea;
445  }
446  }
447  }
448  }
449  return result;
450 }
451 
452 template <class T>
453 bool getBinsFromRow(double lat, double lon, T &pixelPoly, std::map<uint64_t, double> &areas)
454 {
455  int32_t row0, col0;
456  int32_t col;
457  bool result = false;
458  double areaFrac;
459  uint64_t bin;
460 
461  shape->latlon2rowcol(lat, lon, row0, col0);
462 
463  // look right
464  col = col0;
465  while (binIntersectsPixel(row0, col, pixelPoly, areaFrac))
466  {
467  result = true;
468  bin = shape->rowcol2bin(row0, col);
469  areas.emplace(bin, areaFrac);
470  col++;
471  shape->constrainRowCol(row0, col);
472  if (col == col0)
473  break;
474  }
475 
476  // look left
477  col = col0 - 1;
478  while (binIntersectsPixel(row0, col, pixelPoly, areaFrac))
479  {
480  result = true;
481  bin = shape->rowcol2bin(row0, col);
482  areas.emplace(bin, areaFrac);
483  col--;
484  shape->constrainRowCol(row0, col);
485  if (col == col0 - 1)
486  break;
487  }
488 
489  return result;
490 }
491 
492 void getBins(int32_t ifile, int32_t ipixl, std::map<uint64_t, double> &areas)
493 {
494 
495  areas.clear();
496 
497  double lat0 = l2_str[ifile].latitude[ipixl];
498  double lon0 = l2_str[ifile].longitude[ipixl];
499  int32_t row0;
500  double lat;
501 
502 #ifdef DEBUG_PRINT
503  // debug plotting
504  // plot center point
505  if (enableDebugPrint)
506  {
507  printf("lat=%f\n", lat0);
508  printf("lon=%f\n", lon0);
509  printf("plt.plot(lon, lat, \"ro\")\n\n");
510  }
511 #endif
512 
513  // use the center of the bin lat for search
514  row0 = shape->lat2row(lat0);
515  lat0 = shape->row2lat(row0);
516 
517  double deltaLat = 180.0 / shape->getNumRows();
518 
519  if (input.area_weighting == 1)
520  {
521  Box_t pixelBox;
522 
523  // use the pixel box
524  // lat1/lon1 contain pixel delta
525  pixelBox.min_corner().set<0>(l2_str[ifile].longitude[ipixl] - l2_str[ifile].lon1[ipixl]);
526  pixelBox.min_corner().set<1>(l2_str[ifile].latitude[ipixl] - l2_str[ifile].lat1[ipixl]);
527  pixelBox.max_corner().set<0>(l2_str[ifile].longitude[ipixl] + l2_str[ifile].lon1[ipixl]);
528  pixelBox.max_corner().set<1>(l2_str[ifile].latitude[ipixl] + l2_str[ifile].lat1[ipixl]);
529 
530 #ifdef DEBUG_PRINT
531  // debug
532  // plot pixel box
533  if (enableDebugPrint)
534  {
535  printf("lat = [%f, %f, %f, %f, %f]\n",
536  pixelBox.min_corner().y(), pixelBox.max_corner().y(), pixelBox.max_corner().y(),
537  pixelBox.min_corner().y(), pixelBox.min_corner().y());
538  printf("lon = [%f, %f, %f, %f, %f]\n",
539  pixelBox.min_corner().x(), pixelBox.min_corner().x(), pixelBox.max_corner().x(),
540  pixelBox.max_corner().x(), pixelBox.min_corner().x());
541  printf("plt.plot(lon, lat)\n\n");
542  }
543 #endif
544 
545  // look up
546  lat = lat0;
547  while (getBinsFromRow(lat, lon0, pixelBox, areas))
548  {
549  lat += deltaLat;
550  if (lat > 90.0)
551  break;
552  }
553 
554  // look down
555  lat = lat0 - deltaLat;
556  while (getBinsFromRow(lat, lon0, pixelBox, areas))
557  {
558  lat -= deltaLat;
559  if (lat < -90.0)
560  break;
561  }
562  }
563  else if (input.area_weighting == 2)
564  {
565  Box_t pixelBox;
566 
567  // use the pixel bounding box
568  float latMax = std::max(l2_str[ifile].lat1[ipixl], l2_str[ifile].lat1[ipixl + 1]);
569  latMax = std::max(latMax, l2_str[ifile].lat2[ipixl]);
570  latMax = std::max(latMax, l2_str[ifile].lat2[ipixl + 1]);
571 
572  float latMin = std::min(l2_str[ifile].lat1[ipixl], l2_str[ifile].lat1[ipixl + 1]);
573  latMin = std::min(latMin, l2_str[ifile].lat2[ipixl]);
574  latMin = std::min(latMin, l2_str[ifile].lat2[ipixl + 1]);
575 
576  float lonMax = std::max(l2_str[ifile].lon1[ipixl], l2_str[ifile].lon1[ipixl + 1]);
577  lonMax = std::max(lonMax, l2_str[ifile].lon2[ipixl]);
578  lonMax = std::max(lonMax, l2_str[ifile].lon2[ipixl + 1]);
579 
580  float lonMin = std::min(l2_str[ifile].lon1[ipixl], l2_str[ifile].lon1[ipixl + 1]);
581  lonMin = std::min(lonMin, l2_str[ifile].lon2[ipixl]);
582  lonMin = std::min(lonMin, l2_str[ifile].lon2[ipixl + 1]);
583 
584  if ((lonMax - lonMin) > 180)
585  {
586  float tmpF = lonMin;
587  lonMin = lonMax;
588  lonMax = tmpF + 360.0;
589  }
590 
591  pixelBox.min_corner().set<0>(lonMin);
592  pixelBox.min_corner().set<1>(latMin);
593  pixelBox.max_corner().set<0>(lonMax);
594  pixelBox.max_corner().set<1>(latMax);
595 
596 #ifdef DEBUG_PRINT
597  // debug
598  // plot pixel box
599  if (enableDebugPrint)
600  {
601  printf("lat = [%f, %f, %f, %f, %f]\n",
602  pixelBox.min_corner().y(), pixelBox.max_corner().y(), pixelBox.max_corner().y(),
603  pixelBox.min_corner().y(), pixelBox.min_corner().y());
604  printf("lon = [%f, %f, %f, %f, %f]\n",
605  pixelBox.min_corner().x(), pixelBox.min_corner().x(), pixelBox.max_corner().x(),
606  pixelBox.max_corner().x(), pixelBox.min_corner().x());
607  printf("plt.plot(lon, lat)\n\n");
608  }
609 #endif
610 
611  // look up
612  lat = lat0;
613  while (getBinsFromRow(lat, lon0, pixelBox, areas))
614  {
615  lat += deltaLat;
616  if (lat > 90.0)
617  break;
618  }
619 
620  // look down
621  lat = lat0 - deltaLat;
622  while (getBinsFromRow(lat, lon0, pixelBox, areas))
623  {
624  lat -= deltaLat;
625  if (lat < -90.0)
626  break;
627  }
628  }
629  else
630  {
631  // use the exact pixel polygon
632  Polygon_t pixelPoly;
633 
634 #ifdef DEBUG_PRINT
635  // debug
636  // plot pixel polygon
637  if (enableDebugPrint)
638  {
639  printf("lat = [%f, %f, %f, %f, %f]\n",
640  l2_str[ifile].lat1[ipixl],
641  l2_str[ifile].lat2[ipixl],
642  l2_str[ifile].lat2[ipixl + 1],
643  l2_str[ifile].lat1[ipixl + 1],
644  l2_str[ifile].lat1[ipixl]);
645  printf("lon = [%f, %f, %f, %f, %f]\n",
646  l2_str[ifile].lon1[ipixl],
647  l2_str[ifile].lon2[ipixl],
648  l2_str[ifile].lon2[ipixl + 1],
649  l2_str[ifile].lon1[ipixl + 1],
650  l2_str[ifile].lon1[ipixl]);
651  printf("plt.plot(lon, lat)\n\n");
652  }
653 #endif
654 
655  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon1[ipixl], l2_str[ifile].lat1[ipixl]));
656  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon2[ipixl], l2_str[ifile].lat2[ipixl]));
657  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon2[ipixl + 1], l2_str[ifile].lat2[ipixl + 1]));
658  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon1[ipixl + 1], l2_str[ifile].lat1[ipixl + 1]));
659  bg::append(pixelPoly.outer(), Point_t(l2_str[ifile].lon1[ipixl], l2_str[ifile].lat1[ipixl]));
660 
661  // make sure the polygon is defined properly
662  bg::correct(pixelPoly);
663 
664  // look up
665  lat = lat0;
666  while (getBinsFromRow(lat, lon0, pixelPoly, areas))
667  {
668  lat += deltaLat;
669  if (lat > 90.0)
670  break;
671  }
672 
673  // look down
674  lat = lat0 - deltaLat;
675  while (getBinsFromRow(lat, lon0, pixelPoly, areas))
676  {
677  lat -= deltaLat;
678  if (lat < -90.0)
679  break;
680  }
681  }
682 }
683 
692 bool skip_DL(float lon, int side, int night_flag, time_t end_day, time_t beg_day) {
693  if (night_flag == 1) {
694  if ((side == -1) && (beg_day == -1) && (lon < 0))
695  return true;
696  if ((side == +1) && (end_day == 0) && (lon > 0))
697  return true;
698  } else {
699  if ((side == -1) && (beg_day <= 0) && (lon < 0))
700  return true;
701  if ((side == +1) && (end_day >= 0) && (lon > 0))
702  return true;
703  }
704  return false;
705 }
706 
707 int main(int argc, char **argv)
708 {
709  int i, j, k;
710  int status;
711  intn ret_status = 0;
712 
713  int32_t index;
714  int32_t ifile, jsrow, ipixl;
715  uint64_t bin;
716  int32_t ibin;
717  int32_t nfiles;
718  int *fileused;
719  int32_t n_active_files;
720 
721  int32_t within_flag;
722 
723  int32_t n_filled_bins;
724  int64_t total_filled_bins = 0;
725  int32_t date;
726 
727  int16_t brk_scan[MAXNFILES];
728 
729  float32 *slat = NULL;
730  float32 *elat = NULL;
731  float32 dlat;
732  double latbin = 0.0;
733  double lonbin = 0.0;
734 
735  int32_t igroup, ngroup;
736 
737  static int32_t *bscan_row[MAXNFILES];
738  static int32_t *escan_row[MAXNFILES];
739  static unsigned char *scan_in_rowgroup[MAXNFILES];
740 
741  int32_t fileid_w;
742  int32_t out_grp;
743 
744  int64_t *beg;
745  int32_t *ext;
746  int64_t *binnum_data;
747  int32_t i32;
748  int32_t iprod;
749 
750  time_t diffday_beg, diffday_end;
751  int32_t syear, sday, eyear, eday;
752 
753  int32_t ntilts;
754  int16_t tilt_flags[MTILT_DIMS_2];
755  int16_t tilt_ranges[LTILT_DIMS_2][MTILT_DIMS_2];
756 
757  uint32_t flagusemask;
758  uint32_t required;
759  uint32_t flagcheck;
760 
761  int32_t proc_day_beg, proc_day_end;
762 
763  uint8_t qual_max_allowed;
764  std::vector<uint8_t> best_qual;
765  double *sum_bin;
766  double *sum2_bin;
767  double wgt;
768  float32 northmost = -90.0, southmost = 90.0, eastmost = -180.0, westmost = 180.0;
769 
770 
771  double time_rec = 0;
772  time_t tnow;
773  int32_t dataday0, dataday1, startdate, enddate;
774  double dbldate;
775  struct tm *tmnow;
776  static meta_l2Type meta_l2;
777  static meta_l3bType meta_l3b;
778 
779  static float lonRanges[4];
780  lonRanges[0] = 0; // min, max in Western hemisphere
781  lonRanges[1] = -180;
782  lonRanges[2] = 180; // min, max in Eastern hemisphere
783  lonRanges[3] = 0;
784 
785  char units[1024];
786  std::vector<std::string> l2_prodname;
787  std::vector<std::string> l3_prodname;
788 
789  /* Function Prototypes */
790  int64_t getbinnum(int32_t, int32_t);
791 
792  FILE *fp2 = NULL;
793 
794  binListStruct_nc *binList32nc = NULL;
795  binListStruct64_nc *binList64nc = NULL;
796 
797  binIndexStruct_nc *binIndex32nc = NULL;
798  binIndexStruct64_nc *binIndex64nc = NULL;
799 
800 #ifdef MALLINFO
801  struct mallinfo minfo;
802 #endif
803 
805 
806  /* From Fred Patt
807 
808  sum(data) sum(data)*sqrt(n)
809  s = --------- = ----------------- = avg(data)*sqrt(n)
810  sqrt(n) n
811 
812  */
813 
815  printf("%s %s (%s %s)\n", PROGRAM, VERSION, __DATE__, __TIME__);
816 
817  l2bin_input(argc, argv, &input, PROGRAM, VERSION);
818 
819  // set dictionary for file names
820  std::map<std::string,std::string> output_l3_filenames;
821  if(input.output_product_names[0]){
822  std::vector<std::string> output_pairs_l2_l3_names;
823  boost::split(output_pairs_l2_l3_names,std::string(input.output_product_names),boost::is_any_of(","));
824  for(const auto & pair_l2_l3 : output_pairs_l2_l3_names){
825  std::vector<std::string> l2_l3_name;
826  boost::split(l2_l3_name,pair_l2_l3,boost::is_any_of(":"));
827  boost::trim(l2_l3_name[0]);
828  boost::trim(l2_l3_name[1]);
829  output_l3_filenames[l2_l3_name[0]] = l2_l3_name[1];
830  }
831  }
832 
833  // check to make sure area weighting AND composite_prod are not both set
834  if(input.composite_prod[0] != 0 && input.area_weighting != 0) {
835  printf("-E- Arguments \"composite_prod\" and \"area_weighting\" are incompatable. Do not set both.\n");
836  exit(EXIT_FAILURE);
837  }
838 
839  switch (input.area_weighting)
840  {
841  case 0:
842  enableL2PixelArea(L2PixelOff); // no corner calc
843  break;
844  case 1:
845  enableL2PixelArea(L2PixelDelta); // pixel deltas
846  break;
847  default:
848  enableL2PixelArea(L2PixelCorner); // pixel corners
849  break;
850  }
851 
852  // Get times for each file
853  nfiles = input.files.size();
854  printf("%d input files\n", nfiles);
855 
856  bool flags_l2_use = set_l2_flags_use(input.flaguse);
857  std::string products_requested_l3_temp,products_requested_l3;
858  std::string products_requested;
859  for (ifile = 0; ifile < nfiles; ifile++)
860  {
861  std::string product_list_temp;
862  std::string product_list;
863  products_requested = input.l3bprod;
864  if (products_requested == "ALL")
865  {
866  products_requested = "";
867  find_variables_geo_physical(input.files[ifile],products_requested);
868  }
869 
870  // check for quality products and composite products
871  if (input.composite_prod[0])
872  {
873  if (products_requested.find(input.composite_prod) == std::string::npos)
874  {
875  products_requested.append(",");
876  products_requested.append(input.composite_prod);
877  }
878  }
879  if (input.qual_prod[0])
880  {
881  if (products_requested.find(input.qual_prod) == std::string::npos)
882  {
883  products_requested.append(",");
884  products_requested.append(input.qual_prod);
885  }
886  }
887 
888  // set the list of requested products
889  std::vector<std::string> products_requested_separated;
890  set_mapper(input.files[ifile], products_requested, products_requested_separated,input.output_wavelengths);
891 
892  // sanitize the input product list: remove duplicates
893  std::unordered_set<std::string> products_l2_unique;
894  for (size_t i = 0; i < products_requested_separated.size(); i++)
895  {
896  const auto &prod = products_requested_separated.at(i);
897  const auto name_to_pass = prod; // return_l2_name(prod);
898  if (products_l2_unique.count(name_to_pass) > 0)
899  continue;
900  if (!product_list_temp.empty())
901  product_list_temp += "," + name_to_pass;
902  else
903  product_list_temp += name_to_pass; //
904  if(name_to_pass!=input.qual_prod)
905  {
906  if (!product_list.empty())
907  product_list += "," + name_to_pass;
908  else
909  product_list += name_to_pass;
910  }//
911 
912  products_l2_unique.insert(name_to_pass);
913  }
914  products_requested_l3_temp = product_list_temp;
915  products_requested_l3 = product_list;
916  status = openL2(input.files[ifile].c_str(), products_requested_l3_temp.c_str(), &l2_str[ifile]);
917  // l2_stimes[ifile] = (time_t)yds2unix(l2_str[ifile].syear,
918  // l2_str[ifile].sday,
919  // l2_str[ifile].smsec / 1000); // 86400;
920  closeL2(&l2_str[ifile], ifile);
921  }
922 
923  if ((fileused = (int *)calloc(nfiles, sizeof(int))) == NULL)
924  {
925  printf("-E- Problem allocating memory for fileused element of metadata structure\n");
926  exit(EXIT_FAILURE);
927  }
928 
929  proc_day_beg = input.sday;
930  proc_day_end = input.eday;
931  syear = (int32_t)input.sday / 1000.;
932  sday = input.sday - 1000 * syear;
933  startdate = (int32_t)(yds2unix(syear, sday, 0) / 86400);
934  eyear = (int32_t)input.eday / 1000.;
935  eday = input.eday - 1000 * eyear;
936  enddate = (int32_t)(yds2unix(eyear, eday, 0) / 86400);
937 
938  resolve2binRows(input.resolve, nrows, meta_l3b.resolution);
939 
940  if (nrows == -1)
941  {
942  printf("Grid resolution not defined.\n");
943  exit(EXIT_FAILURE);
944  }
945  dlat = 180. / nrows;
946  bool is64bit = (nrows > 2160 * 16);
947 
948  qual_max_allowed = input.qual_max;
949 
950  printf("Resolution: %s\n", input.resolve);
951  printf("Max Qual Allowed: %d\n", input.qual_max);
952 
953  /* Setup flag mask */
954  /* --------------- */
955  if (flags_l2_use)
956  {
957  strcpy(buf, l2_str[0].flagnames);
958  setupflags(buf, input.flaguse, &flagusemask, &required, &status,l2_str[0].l2_bits);
959  }
960 
961  // Parse L3 Product list
962  std::string delim1 = ", "; // product delimiters
963  std::string l3bprod = products_requested_l3;
964  boost::trim_if(l3bprod, boost::is_any_of(delim1));
965  std::vector<std::string> prodparam;
966  boost::algorithm::split(prodparam, l3bprod,
967  boost::is_any_of(delim1));
968 
969  // set min, max and thirdim
970  set_prodname_3d_to_l2(prodparam, l2_str[0], l2_prodname, l3_prodname, thirdDimId, min_value, max_value);
971  l3b_nprod = l3_prodname.size();
972  if (l3b_nprod > max_l3b_products)
973  {
974  printf("-E- Number of output products is greater than %d\n", max_l3b_products);
975  exit(EXIT_FAILURE);
976  }
977 
978  /* Initialize bscan_row, escan_row, numer, denom */
979  /* --------------------------------------------- */
980  for (i = 0; i < MAXNFILES; i++)
981  {
982  bscan_row[i] = NULL;
983  escan_row[i] = NULL;
984  numer[i] = NULL;
985  denom[i] = NULL;
986  scan_in_rowgroup[i] = NULL;
987  }
988 
989  /* Check each L2 file for required products */
990  /* ---------------------------------------- */
991  for (ifile = 0; ifile < nfiles; ifile++)
992  {
993 
994  numer[ifile] = (int16_t *)calloc(l3b_nprod, sizeof(int16_t));
995  denom[ifile] = (int16_t *)calloc(l3b_nprod, sizeof(int16_t));
996 
997  // Check whether L3 products exist in L2
998  for (iprod = 0; iprod < l3b_nprod; iprod++)
999  {
1000 
1001  std::vector<std::string> operands;
1002  boost::algorithm::split(operands, l2_prodname[iprod],
1003  boost::is_any_of("/"));
1004 
1005  // find numerator, or solo product
1006  index = get_l2prod_index(l2_str[ifile], operands[0].c_str());
1007  if (index < 0)
1008  {
1009  printf("L3 product: \"%s\" not found in L2 file \"%s\".\n",
1010  operands[0].c_str(), l2_str[ifile].filename);
1011  exit(EXIT_FAILURE);
1012  }
1013  else
1014  numer[ifile][iprod] = index;
1015 
1016  // find denominator product (as needed)
1017  denom[ifile][iprod] = -1;
1018  if (operands.size() > 1)
1019  {
1020  boost::replace_all(l3_prodname[iprod], "/", "_");
1021  index = get_l2prod_index(l2_str[ifile], operands[1].c_str());
1022  if (index < 0)
1023  {
1024  printf("L3 product: \"%s\" not found in L2 file \"%s\".\n",
1025  operands[1].c_str(), l2_str[ifile].filename);
1026  exit(EXIT_FAILURE);
1027  }
1028  else
1029  denom[ifile][iprod] = index;
1030  }
1031 
1032  } // iprod loop
1033 
1034  // Check whether Quality product exists in L2
1035  if (input.qual_prod[0] != 0)
1036  {
1037  index = get_l2prod_index(l2_str[ifile], input.qual_prod);
1038  if (index < 0)
1039  {
1040  printf("Quality product: \"%s\" not found in L2 file \"%s\".\n.See line %d in %s\n",
1041  input.qual_prod, l2_str[ifile].filename, __LINE__, __FILE__);
1042  exit(EXIT_FAILURE);
1043  }
1044  else
1046  }
1047 
1048  // Check whether Composite product exists in L2
1049  if (input.composite_prod[0] != 0)
1050  {
1051  index = get_l2prod_index(l2_str[ifile], input.composite_prod);
1052  if (index < 0)
1053  {
1054  printf("Composite product: \"%s\" not found in L2 file \"%s\".\n",
1055  input.composite_prod, l2_str[ifile].filename);
1056  exit(EXIT_FAILURE);
1057  }
1058  else
1060  }
1061 
1062  } // ifile loop
1063 
1064  /* Check whether composite product exists in L3 product list */
1065  /* --------------------------------------------------------- */
1066  if (input.composite_prod[0] != 0)
1067  {
1068  for (iprod = 0; iprod < l3b_nprod; iprod++)
1069  {
1070  if (l3_prodname[iprod] == input.composite_prod)
1071  {
1072  composite_l3prod_index = iprod;
1073  break;
1074  }
1075  }
1076  if (composite_l3prod_index == -1)
1077  {
1078  printf("Composite product: \"%s\" not found in L3 product list.\n",
1079  input.composite_prod);
1080  exit(EXIT_FAILURE);
1081  }
1082  }
1083 
1084  /* Find begin and end scan latitudes for each swath row */
1085  /* ---------------------------------------------------- */
1086  for (ifile = 0; ifile < nfiles; ifile++)
1087  {
1088  size_t nrec = l2_str[ifile].nrec;
1089  slat = (float32 *)calloc(l2_str[ifile].nrec, sizeof(float32));
1090  elat = (float32 *)calloc(l2_str[ifile].nrec, sizeof(float32));
1091  bscan_row[ifile] = (int32_t *)calloc(l2_str[ifile].nrec, sizeof(int32_t));
1092  escan_row[ifile] = (int32_t *)calloc(l2_str[ifile].nrec, sizeof(int32_t));
1094  try {
1095  netCDF::NcFile nc_input(input.files[ifile], netCDF::NcFile::read);
1096  Geospatialbounds geo_bounds(nc_input);
1097  instrument = geo_bounds.get_instrument();
1098  platform = geo_bounds.get_platform();
1099  std::copy(geo_bounds.get_slat(),geo_bounds.get_slat() + nrec,slat);
1100  std::copy(geo_bounds.get_elat(),geo_bounds.get_elat() + nrec,elat);
1101  nc_input.close();
1102  }
1103  catch (netCDF::exceptions::NcException const &e)
1104  {
1105  std::cout << e.what() << std::endl;
1106  exit(EXIT_FAILURE);
1107  }
1108  catch (std::exception const &e)
1109  {
1110  std::cerr << "-X- " << __FILE__ " line " << __LINE__ << ": "
1111  << e.what() << std::endl;
1112  exit(EXIT_FAILURE);
1113  }
1114 
1115  /* Get minimum/maximu value from product.xml */
1116  /* ---------------------------------- */
1118  platform.c_str());
1119  /* Note: bscan > escan */
1120 
1121  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++)
1122  {
1123  escan_row[ifile][jsrow] = (int32_t)((90 + elat[jsrow]) / dlat);
1124  bscan_row[ifile][jsrow] = (int32_t)((90 + slat[jsrow]) / dlat);
1125 
1126  if (escan_row[ifile][jsrow] > bscan_row[ifile][jsrow])
1127  {
1128  k = escan_row[ifile][jsrow];
1129  escan_row[ifile][jsrow] = bscan_row[ifile][jsrow];
1130  bscan_row[ifile][jsrow] = k;
1131  }
1132  escan_row[ifile][jsrow] -= 10;
1133  bscan_row[ifile][jsrow] += 10;
1134  }
1135 
1136  free(slat);
1137  free(elat);
1138 
1139  } /* ifile loop */
1140 
1141  /* Find begin & end scans for each input file */
1142  /* ------------------------------------------ */
1143  n_active_files = nfiles;
1144  for (ifile = 0; ifile < nfiles; ifile++)
1145  {
1146 
1147  /* Determine brk_scan value */
1148  /* ------------------------ */
1149  brk_scan[ifile] = 0;
1150 
1151  /* Regional Product */
1152  /* ---------------- */
1153  if (strcasecmp(input.prodtype, "regional") == 0)
1154  {
1155  // printf("%s brk:%5d %5d %3d %6d\n",
1156  // buf, brk_scan[ifile],
1157  // l2_str[ifile].nrec, l2_str[ifile].sday,
1158  // l2_str[ifile].smsec / 1000);
1159  continue;
1160  }
1161  // get data days
1162 
1163  elat = (float32 *)calloc(l2_str[ifile].nrec, sizeof(float32));
1164  slat = (float32 *)calloc(l2_str[ifile].nrec, sizeof(float32));
1165  size_t nrec = l2_str[ifile].nrec;
1166  try {
1167  netCDF::NcFile nc_input(input.files[ifile], netCDF::NcFile::read);
1168  Geospatialbounds geo_bounds(nc_input);
1169  std::copy(geo_bounds.get_slat(),geo_bounds.get_slat() + nrec,slat);
1170  std::copy(geo_bounds.get_elat(),geo_bounds.get_elat() + nrec,elat);
1171  std::pair<int32_t,int32_t> dates_0_1 = get_datadays(nc_input,input.deltaeqcross,input.night);
1172  dataday0 = dates_0_1.first;
1173  dataday1 = dates_0_1.second;
1174  nc_input.close();
1175  }
1176  catch (netCDF::exceptions::NcException const &e)
1177  {
1178  std::cout << e.what();
1179  exit(EXIT_FAILURE);
1180  }
1181  catch (std::exception const &e)
1182  {
1183  std::cerr << "-X- " << __FILE__ " line " << __LINE__ << ": "
1184  << e.what() << std::endl;
1185  exit(EXIT_FAILURE);
1186  }
1187 
1188  date = (time_t)yds2unix(l2_str[ifile].syear,
1189  l2_str[ifile].sday,
1190  l2_str[ifile].smsec / 1000) /
1191  86400;
1192  diffday_beg = date - startdate;
1193  diffday_end = date - enddate;
1194  if (dataday1 == dataday0)
1195  {
1196  if (dataday0 < startdate || dataday1 > enddate)
1197  brk_scan[ifile] = -9999;
1198  else
1199  brk_scan[ifile] = 0; // -9999;
1200  }
1201  else
1202  {
1203  if (dataday1 < startdate) // startdate is dataday conversion of input sday
1204  brk_scan[ifile] = -9999;
1205  else if (dataday0 > enddate)
1206  { // enddate is dataday conversion of input eday
1207  brk_scan[ifile] = -9999;
1208  }
1209  else
1210  {
1211  if (dataday1 > enddate)
1212  brk_scan[ifile] = 1;
1213  else
1214  brk_scan[ifile] = -1;
1215  }
1216  }
1217  if (brk_scan[ifile] == -9999)
1218  n_active_files--;
1219  free(elat);
1220  free(slat);
1221  } /* ifile loop */
1222 
1223 
1224  shape = new l3::L3ShapeIsine(nrows);
1225 
1226  /* Compute basebin array (Starting bin of each row [1-based]) */
1227  /* ---------------------------------------------------------- */
1228  basebin = (int64_t *)calloc(nrows + 1, sizeof(int64_t));
1229  for (i = 0; i <= nrows; i++)
1230  {
1231  basebin[i] = shape->getBaseBin(i);
1232  }
1233  basebin[nrows] = shape->getNumBins() + 1;
1234 
1235  printf("total number of bins: %ld\n", (long int)basebin[nrows] - 1);
1236 
1237  /*
1238  * Create output file
1239  */
1240  strcpy(buf, input.ofile);
1241  status = nc_create(buf, NC_NETCDF4, &fileid_w);
1242  check_err(status, __LINE__, __FILE__);
1243 
1244  status = nc_def_grp(fileid_w, "level-3_binned_data", &out_grp);
1245  check_err(status, __LINE__, __FILE__);
1246 
1247  idDS ds_id;
1248  ds_id.fid = fileid_w;
1249  ds_id.sid = -1;
1250  ds_id.fftype = DS_NCDF; // FMT_L2NCDF
1251 
1252  if (is64bit)
1253  status = defineBinList64_nc(input.deflate, out_grp);
1254  else
1255  status = defineBinList_nc(input.deflate, out_grp);
1256  if (status)
1257  {
1258  printf("-E- %s:%d Could not define binList variable in output file\n",
1259  __FILE__, __LINE__);
1260  exit(EXIT_FAILURE);
1261  }
1262 
1263  char **prodnames = (char **)malloc(l3_prodname.size() * sizeof(char *));
1264  for (size_t size = 0; size < l3_prodname.size(); size++)
1265  {
1266  // look for user requested product names
1267  if(output_l3_filenames.find(l3_prodname[size]) != output_l3_filenames.end())
1268  prodnames[size] = strdup(output_l3_filenames[l3_prodname[size]].c_str());
1269  else
1270  prodnames[size] = strdup(l3_prodname[size].c_str());
1271  }
1272  status = defineBinData_nc(input.deflate, out_grp, l3b_nprod, prodnames);
1273  // set wave float attribute here
1274  for (size_t size = 0; size < l3_prodname.size(); size++) {
1275  float wave = BAD_FLT;
1276  wave = l3_attr(l3_prodname[size]);
1277  if(wave !=BAD_FLT) {
1278  int varid;
1279  nc_inq_varid(out_grp, prodnames[size] , &varid);
1280  nc_put_att_float(out_grp,varid,"wavelength",NC_FLOAT,1,&wave);
1281  }
1282  }
1283 
1284  if (status)
1285  {
1286  printf("-E- %s:%d Could not define binData variable in output file\n",
1287  __FILE__, __LINE__);
1288  exit(EXIT_FAILURE);
1289  }
1290  for (size_t i = 0; i < l3_prodname.size(); i++)
1291  {
1292  free(prodnames[i]);
1293  }
1294  free(prodnames);
1295 
1296  if (is64bit)
1297  status = defineBinIndex64_nc(input.deflate, out_grp);
1298  else
1299  status = defineBinIndex_nc(input.deflate, out_grp);
1300  if (status)
1301  {
1302  printf("-E- %s:%d Could not define binIndex variable in output file\n",
1303  __FILE__, __LINE__);
1304  exit(EXIT_FAILURE);
1305  }
1306 
1307  if (input.qual_prod[0] != 0)
1308  {
1309  status = defineQuality_nc(input.deflate, out_grp);
1310  if (status)
1311  {
1312  printf("-E- %s:%d Could not define quality variable in output file\n",
1313  __FILE__, __LINE__);
1314  exit(EXIT_FAILURE);
1315  }
1316  }
1317 
1318  /* Allocate Arrays for Bin Index */
1319  /* ----------------------------- */
1320  beg = (int64_t *)calloc(nrows, sizeof(int64_t));
1321  ext = (int32_t *)calloc(nrows, sizeof(int32_t));
1322 
1323  if (is64bit)
1324  binIndex64nc = (binIndexStruct64_nc *)calloc(nrows, sizeof(binIndexStruct64_nc));
1325  else
1326  binIndex32nc = (binIndexStruct_nc *)calloc(nrows, sizeof(binIndexStruct_nc));
1327 
1328  /* Initialize bin_indx array */
1329  /* ------------------------- */
1330  for (i = 0; i < nrows; i++)
1331  {
1332 
1333  i32 = i;
1334 
1335  if (i32 < 0 || i32 >= nrows)
1336  {
1337  printf("%d %d\n", i, nrows);
1338  exit(-1);
1339  }
1340 
1341  if (is64bit)
1342  {
1343  binIndex64nc[i].begin = beg[i32];
1344  binIndex64nc[i].extent = ext[i32];
1345  binIndex64nc[i].max = shape->getNumCols(i32);
1346  }
1347  else
1348  {
1349  binIndex32nc[i].begin = beg[i32];
1350  binIndex32nc[i].extent = ext[i32];
1351  binIndex32nc[i].max = shape->getNumCols(i32);
1352  }
1353  }
1354 
1355  // Row Group
1356  n_rows_in_group = input.rowgroup;
1357  if (n_rows_in_group <= 0)
1358  {
1359  printf("row_group not defined, using 270.\n");
1360  n_rows_in_group = 270;
1361  }
1362  if (input.verbose)
1363  printf("%d %d %d\n", proc_day_beg, proc_day_end, n_rows_in_group);
1364 
1365  /* Find row_group that divides nrows */
1366  /* --------------------------------- */
1367  for (i = nrows; i > 0; i--)
1368  {
1369  if ((nrows % i) == 0)
1370  {
1371  if (i <= n_rows_in_group)
1372  {
1373  n_rows_in_group = i;
1374  break;
1375  }
1376  }
1377  }
1378  if (input.rowgroup != n_rows_in_group)
1379  {
1380  printf("Input row_group: %d Actual row_group: %d\n",
1381  input.rowgroup, n_rows_in_group);
1382  }
1383  ngroup = nrows / n_rows_in_group;
1384 
1385  /* Process each group of bin rows (Main Loop) */
1386  /* ========================================== */
1387  int64_t total_count = 0;
1388  for (krow = 0, igroup = 0; igroup < ngroup; igroup++)
1389  {
1390 
1391  if (shape->row2lat(krow + n_rows_in_group) < input.latsouth)
1392  {
1393  krow += n_rows_in_group;
1394  continue;
1395  }
1396  if (shape->row2lat(krow) > input.latnorth)
1397  {
1398  krow += n_rows_in_group;
1399  continue;
1400  }
1401 
1402  /* Print info on rowgroup */
1403  /* ---------------------- */
1404  time(&tnow);
1405  tmnow = localtime(&tnow);
1406  printf("krow:%6d out of %6d (%6.2f to %6.2f) ",
1407  krow, nrows,
1408  ((float32)(krow) / nrows) * 180 - 90,
1409  ((float32)(krow + n_rows_in_group) / nrows) * 180 - 90);
1410  printf("%s", asctime(tmnow));
1411 
1412  n_bins_in_group = basebin[krow + n_rows_in_group] - basebin[krow];
1413  within_flag = 0;
1414 
1415  /* Determine relevant swath rows for this bin row group for each file */
1416  /* ------------------------------------------------------------------ */
1417  for (ifile = 0; ifile < nfiles; ifile++)
1418  {
1419 
1420  /* add an extra 0 to the end of scan_in_rowgroup so the caching
1421  * code never reads past the end of the file */
1422  scan_in_rowgroup[ifile] = (unsigned char *)
1423  calloc(l2_str[ifile].nrec + 1, sizeof(unsigned char));
1424 
1425  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++)
1426  {
1427  scan_in_rowgroup[ifile][jsrow] = 1;
1428  if (bscan_row[ifile][jsrow] < krow ||
1429  escan_row[ifile][jsrow] >= (krow + n_rows_in_group - 1))
1430  {
1431  scan_in_rowgroup[ifile][jsrow] = 255;
1432  }
1433  } /* jsrow loop */
1434 
1435  /* Determine if within bin row group */
1436  /* --------------------------------- */
1437  if (!within_flag)
1438  {
1439  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++)
1440  {
1441  if (scan_in_rowgroup[ifile][jsrow] == 1)
1442  {
1443  within_flag = 1;
1444  break;
1445  }
1446  } /* scan row loop */
1447  }
1448 
1449  } /* ifile loop */
1450 
1451  /* If no swath rows within group then continue to next group */
1452  /* --------------------------------------------------------- */
1453  if (within_flag == 0)
1454  {
1455  for (ifile = 0; ifile < nfiles; ifile++)
1456  {
1457  if (scan_in_rowgroup[ifile] != NULL)
1458  {
1459  free(scan_in_rowgroup[ifile]);
1460  scan_in_rowgroup[ifile] = NULL;
1461  }
1462  }
1463  krow += n_rows_in_group;
1464  continue;
1465  }
1466 
1467  /* Allocate # pixels in bin, bin_flag, tilt, qual, & nscenes arrays */
1468  /* ---------------------------------------------------------------- */
1469  n_filled_bins = 0;
1470  bin_flag = (int32_t *)calloc(n_bins_in_group, sizeof(int32_t));
1471  tilt = (int16_t *)calloc(n_bins_in_group, sizeof(int16_t));
1472  qual = (int16_t *)calloc(n_bins_in_group, sizeof(int16_t));
1473  nscenes = (int16_t *)calloc(n_bins_in_group, sizeof(int16_t));
1474  lastfile = (int16_t *)calloc(n_bins_in_group, sizeof(int16_t));
1475 
1476  for (i = 0; i < n_bins_in_group; i++)
1477  {
1478  tilt[i] = -1;
1479  qual[i] = 3;
1480  lastfile[i] = -1;
1481  }
1482 
1483  /* Allocate bin accumulator & data value arrays */
1484  /* -------------------------------------------- */
1485  nobs = (int16_t *)calloc(n_bins_in_group, sizeof(int16_t));
1486  allocated_space = (int16_t *)calloc(n_bins_in_group, sizeof(int16_t));
1487  data_values = (float32 **)calloc(n_bins_in_group, sizeof(float32 *));
1488  data_areas = (double **)calloc(n_bins_in_group, sizeof(double *));
1489  file_index = (int16_t **)calloc(n_bins_in_group, sizeof(int16_t *));
1490  data_quality = (uint8_t **)calloc(n_bins_in_group, sizeof(uint8_t *));
1491  time_value = (float64 **)calloc(n_bins_in_group, sizeof(float64 *));
1492 
1493  /* Initialize bin counters */
1494  /* ----------------------- */
1495  n_allocperbin =
1496  n_active_files * l2_str[0].nrec * l2_str[0].nsamp / 50000000;
1497 
1498  if (n_allocperbin < 2)
1499  n_allocperbin = 2;
1502 
1503  if (input.verbose)
1504  printf("%-20s:%8d\n\n", "# allocated per bin", n_allocperbin);
1505 
1506  for (i = 0; i < n_bins_in_group; i++)
1507  {
1508  nobs[i] = 0;
1509  allocated_space[i] = 0;
1510  lastfile[i] = -1;
1511  }
1512 
1513  /* Read L2 files and fill data_values (L3b) array */
1514  /* ++++++++++++++++++++++++++++++++++++++++++++++ */
1515  for (ifile = 0; ifile < nfiles; ifile++)
1516  {
1517 
1519 
1520  /* if "early" or "late" input file then skip */
1521  /* ----------------------------------------- */
1522  if (brk_scan[ifile] == -9999)
1523  continue;
1524 
1525  /* if no scans in rowgroup for this file then skip */
1526  /* ----------------------------------------------- */
1527  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++)
1528  {
1529  if (scan_in_rowgroup[ifile][jsrow] == 1)
1530  {
1531  break;
1532  }
1533  }
1534  if (jsrow == l2_str[ifile].nrec)
1535  {
1536  continue;
1537  }
1538 
1539  reopenL2(ifile, &l2_str[ifile]);
1540 
1541  /* Get tilt flags & ranges */
1542  /* ----------------------- */
1543  ntilts = l2_str[ifile].ntilts;
1544  for (i = 0; i < ntilts; i++)
1545  {
1546  tilt_flags[i] = l2_str[ifile].tilt_flags[i];
1547  tilt_ranges[0][i] = l2_str[ifile].tilt_ranges[0][i];
1548  tilt_ranges[1][i] = l2_str[ifile].tilt_ranges[1][i];
1549  }
1550 
1551  /* Get date stuff */
1552  /* -------------- */
1553  date = (time_t)yds2unix(l2_str[ifile].syear, l2_str[ifile].sday,
1554  l2_str[ifile].smsec / 1000) /
1555  86400;
1556  diffday_beg = date - startdate;
1557  diffday_end = date - enddate;
1558 
1559  /* Loop over swath rows */
1560  /* ^^^^^^^^^^^^^^^^^^^^ */
1561  for (jsrow = 0; jsrow < l2_str[ifile].nrec; jsrow++)
1562  {
1563 
1564  /* if swath row not within group then continue */
1565  /* ------------------------------------------- */
1566  if (scan_in_rowgroup[ifile][jsrow] != 1)
1567  continue;
1568 
1569  /* Read swath record from L2 */
1570  /* ------------------------- */
1572  status = readL2(&l2_str[ifile], ifile, jsrow, -1,
1573  scan_in_rowgroup[ifile]);
1575  auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
1576  auto count_ = duration.count();
1577  total_count += count_;
1578  if (status == 5)
1579  continue;
1580  int scan_crosses_dateline = 0;
1581  float scan_lon = BAD_FLT;
1582  if (brk_scan[ifile] != 0) {
1583  int npixls = l2_str[ifile].nsamp - 1;
1584  float *ptr_lon = l2_str[ifile].longitude;
1585  float slon = BAD_FLT;
1586  float elon = BAD_FLT;
1587 
1588  for (int i_p = 0; i_p <= npixls; i_p++) {
1589  if (std::abs(slon) <= 180.0)
1590  break;
1591  slon = ptr_lon[i_p];
1592  }
1593  for (int i_p = npixls; i_p >= 0; i_p--) {
1594  if (std::abs(elon) <= 180.0)
1595  break;
1596  elon = ptr_lon[i_p];
1597  }
1598  if (std::abs(elon) > 180.0 || std::abs(slon) > 180.0)
1599  continue;
1600  scan_lon = elon;
1601  if (slon * elon < 0)
1602  scan_crosses_dateline = 1;
1603  }
1604  if (scan_crosses_dateline == 0 && (strcasecmp(input.prodtype, "regional") != 0)) {
1605  if (skip_DL(scan_lon, brk_scan[ifile], input.night, diffday_end, diffday_beg))
1606  continue;
1607  }
1608 
1609  /* Check tilt state */
1610  /* ---------------- */
1611  for (i = 0; i < ntilts; i++)
1612  {
1613  if ((jsrow + 1) <= tilt_ranges[1][i])
1614  {
1615  tiltstate = (tilt_flags[i] & 0xFF);
1616  break;
1617  }
1618  }
1619 
1620  if (l2_str[ifile].nsamp == 0)
1621  continue;
1622 
1623  if ((jsrow % 100) == 0 && input.verbose)
1624  {
1625  printf("ifile:%4d jsrow:%6d nsamp:%8d\n", ifile, jsrow, l2_str[ifile].nsamp);
1626  }
1627 
1628  // save the tai93 date for this line, so is can be saved in the bins
1629  dbldate = yds2unix(l2_str[ifile].year, l2_str[ifile].day, l2_str[ifile].msec / 1000.0);
1630  time_tai93 = unix_to_tai93(dbldate);
1631 
1632  /* ##### Loop over L2 pixels ##### */
1633  /* ------------------------------- */
1634  for (ipixl = 0; ipixl < l2_str[ifile].nsamp; ipixl++)
1635  {
1636 
1637 #ifdef DEBUG_PRINT
1638  if (ipixl >= 510 && ipixl < 512)
1639  enableDebugPrint = true;
1640  else
1641  enableDebugPrint = false;
1642 #endif
1643 
1644  /* if bin flagged then continue */
1645  /* ---------------------------- */
1646  if (flags_l2_use)
1647  {
1648  flagcheck = (l2_str[ifile].l2_flags[ipixl]);
1649  if ((flagcheck & flagusemask) != 0)
1650  continue;
1651  if ((flagcheck & required) != required)
1652  continue;
1653  }
1654 
1655  /* Check for dateline crossing */
1656  /* --------------------------- */
1657  if (scan_crosses_dateline == 1) {
1658  if (skip_DL(l2_str[ifile].longitude[ipixl], brk_scan[ifile], input.night, diffday_end,
1659  diffday_beg))
1660  continue;
1661  }
1662  /* Check for bad value in any of the products */
1663  /* ------------------------------------------ */
1664  int bad_value = 0;
1665  for (iprod = 0; iprod < l3b_nprod; iprod++) {
1666  int32_t l2_iprod = numer[ifile][iprod];
1667  for (int i = 0; i < l2_str[ifile].thirdDim[l2_iprod]; i++) {
1668  f32 =
1669  l2_str[ifile].l2_data[l2_iprod][ipixl * l2_str[ifile].thirdDim[l2_iprod] + i];
1670  if (f32 == BAD_FLT) {
1671  bad_value = 1;
1672  break;
1673  }
1674  }
1675  if (bad_value)
1676  break;
1677  }
1678  if (bad_value == 1)
1679  continue;
1680 
1681  /* Check if within longitude boundaries */
1682  /* ------------------------------------ */
1683  if (input.lonwest != 0.0 || input.loneast != 0.0)
1684  {
1685  if (l2_str[ifile].longitude[ipixl] < input.lonwest)
1686  continue;
1687  if (l2_str[ifile].longitude[ipixl] > input.loneast)
1688  continue;
1689  }
1690  /* Check if within latitude boundaries */
1691  /* ------------------------------------ */
1692  if (l2_str[ifile].latitude[ipixl] < input.latsouth)
1693  continue;
1694  if (l2_str[ifile].latitude[ipixl] > input.latnorth)
1695  continue;
1696  if(std::abs(l2_str[ifile].latitude[ipixl]) > 90)
1697  continue;
1698  if(std::abs(l2_str[ifile].longitude[ipixl]) > 180)
1699  continue;
1700  if (input.area_weighting) {
1701  if (std::abs(l2_str[ifile].lat1[ipixl]) > 90)
1702  continue;
1703  if (std::abs(l2_str[ifile].lon1[ipixl]) > 180)
1704  continue;
1705  if (std::abs(l2_str[ifile].lat1[ipixl + 1]) > 90)
1706  continue;
1707  if (std::abs(l2_str[ifile].lon1[ipixl + 1]) > 180)
1708  continue;
1709  }
1710  if (input.area_weighting >= 2) {
1711  if (std::abs(l2_str[ifile].lat2[ipixl]) > 90)
1712  continue;
1713  if (std::abs(l2_str[ifile].lon2[ipixl]) > 180)
1714  continue;
1715  if (std::abs(l2_str[ifile].lat2[ipixl + 1]) > 90)
1716  continue;
1717  if (std::abs(l2_str[ifile].lon2[ipixl + 1]) > 180)
1718  continue;
1719  }
1720  if (input.area_weighting)
1721  {
1722  /* Get map of bin numbers and intersection area that intersect pixel */
1723  std::map<uint64_t, double> areas;
1724  getBins(ifile, ipixl, areas);
1725  for (auto const &val : areas)
1726  {
1727  addPixelToBin(ifile, ipixl, val.first, flags_l2_use, val.second);
1728  }
1729  }
1730  else
1731  {
1732  /* Get Bin Number for Pixel */
1733  /* ------------------------ */
1734  bin = getbinnum(ifile, ipixl); // bin is 1-based
1735  addPixelToBin(ifile, ipixl, bin, flags_l2_use);
1736  }
1737 
1738  } /* ##### i (ipixl) loop ##### */
1739 
1740  } /* ^^^^^^^^^^ jsrow loop ^^^^^^^^^^ */
1741 
1742  closeL2(&l2_str[ifile], ifile);
1743 
1744 #ifdef MALLINFO
1745  if (input.meminfo)
1746  {
1747  int32_t total_alloc_space;
1748  /* malloc_stats();*/
1749  minfo = mallinfo();
1750  total_alloc_space = 0;
1751  for (i = 0; i < n_bins_in_group; i++)
1752  {
1753  total_alloc_space += allocated_space[i];
1754  }
1755  printf("Used space: %10d\n", minfo.uordblks);
1756  printf("Allo space: %10d\n", total_alloc_space * (2 + l3b_nprod * 4));
1757  }
1758 #endif
1759 
1760  } /* ++++++++++ ifile loop ++++++++++ */
1761 
1762  if (input.verbose)
1763  {
1764  time(&tnow);
1765  tmnow = localtime(&tnow);
1766  printf("krow:%5d After data_value fill: %s\n", krow, asctime(tmnow));
1767  }
1768 
1769  /* Compute Total # of filled bins */
1770  /* ------------------------------ */
1771  for (ibin = 0; ibin < n_bins_in_group; ibin++)
1772  {
1773  if (nobs[ibin] > 0 && nobs[ibin] < input.minobs)
1774  nobs[ibin] = 0;
1775 
1776  if (nobs[ibin] != 0)
1777  n_filled_bins++;
1778  } /* ibin loop */
1779 
1780  best_qual = std::vector<uint8_t>(n_bins_in_group, 255);
1781  // memset(best_qual, 255, n_bins_in_group * sizeof(uint8_t));
1782 
1783  /* ********** If filled bins ********** */
1784  /* ------------------------------------ */
1785  if (n_filled_bins > 0)
1786  {
1787 
1788  /* Fill "Bin List" vdata array */
1789  /* --------------------------- */
1790  if (is64bit)
1791  binList64nc = (binListStruct64_nc *)
1792  calloc(n_filled_bins, sizeof(binListStruct64_nc));
1793  else
1794  binList32nc = (binListStruct_nc *)
1795  calloc(n_filled_bins, sizeof(binListStruct_nc));
1796 
1797  i = 0;
1798  for (ibin = 0; ibin < n_bins_in_group; ibin++)
1799  {
1800 
1801  /* Adjust for bins with "bad" quality values */
1802  /* ----------------------------------------- */
1803  if (input.qual_prod[0] != 0 && nobs[ibin] > 0)
1804  {
1805  best_qual[ibin] = 255;
1806  for (j = 0; j < nobs[ibin]; j++)
1807  if (data_quality[ibin][j] < best_qual[ibin])
1808  best_qual[ibin] = data_quality[ibin][j];
1809 
1810  k = 0;
1811  for (j = 0; j < nobs[ibin]; j++)
1812  {
1813  if ((data_quality[ibin][j] <= best_qual[ibin]) &&
1814  (data_quality[ibin][j] <= qual_max_allowed))
1815  {
1816  if (k < j)
1817  {
1818  data_areas[ibin][k] = data_areas[ibin][j];
1819  for (iprod = 0; iprod < l3b_nprod; iprod++)
1820  {
1821  data_values[ibin][k * l3b_nprod + iprod] =
1822  data_values[ibin][j * l3b_nprod + iprod];
1823  }
1824  }
1825  k++;
1826  }
1827  }
1828  nobs[ibin] = k;
1829 
1830  if (nobs[ibin] == 0)
1831  n_filled_bins--;
1832  }
1833 
1834  if (nobs[ibin] != 0)
1835  {
1836  bin = (uint64_t)ibin + basebin[krow];
1837 
1838  if (is64bit)
1839  {
1840  binList64nc[i].binnum = bin;
1841  binList64nc[i].nobs = nobs[ibin];
1842  binList64nc[i].nscenes = nscenes[ibin];
1843  }
1844  else
1845  {
1846  binList32nc[i].binnum = bin;
1847  binList32nc[i].nobs = nobs[ibin];
1848  binList32nc[i].nscenes = nscenes[ibin];
1849  }
1850 
1851  /* weights {=sqrt(# of L2 files in given bin)} */
1852  /* ------------------------------------------- */
1853  wgt = 0.0;
1854  for (ifile = 0; ifile <= nfiles; ifile++)
1855  {
1856  double area = 0.0;
1857  for (j = 0; j < nobs[ibin]; j++)
1858  {
1859  if (file_index[ibin][j] == ifile)
1860  area += data_areas[ibin][j];
1861  }
1862  wgt += sqrt(area);
1863  }
1864 
1865  time_rec = 0.0;
1866  for (j = 0; j < nobs[ibin]; j++)
1867  time_rec += time_value[ibin][j];
1868 
1869  if (is64bit)
1870  {
1871  binList64nc[i].weights = wgt;
1872  binList64nc[i].time_rec = time_rec;
1873  }
1874  else
1875  {
1876  binList32nc[i].weights = wgt;
1877  binList32nc[i].time_rec = time_rec;
1878  }
1879 
1880  i++;
1881 
1882  /* Update Max/Min Lon/Lat */
1883  /* ---------------------- */
1884  shape->bin2latlon(bin, latbin, lonbin);
1885 
1886  if (latbin > northmost)
1887  northmost = latbin;
1888  if (latbin < southmost)
1889  southmost = latbin;
1890 
1891  float minNeg = lonRanges[0]; // min, max in Western hemisphere
1892  float maxNeg = lonRanges[1];
1893  float minPos = lonRanges[2]; // min, max in Eastern hemisphere
1894  float maxPos = lonRanges[3];
1895  if (lonbin < 0)
1896  { // Western hemisphere
1897  minNeg = fmin(minNeg, lonbin);
1898  maxNeg = fmax(maxNeg, lonbin);
1899  }
1900  else
1901  { // Eastern hemisphere
1902  minPos = fmin(minPos, lonbin);
1903  maxPos = fmax(maxPos, lonbin);
1904  }
1905  // Adjust east and west granule bounding coordinates
1906  float max_angle = 90.0;
1907  int8_t lon000 = ((minPos - maxNeg) < max_angle);
1908  int8_t lon180 = ((maxPos - minNeg) > (360.0 - max_angle));
1909  int8_t polar = (lon000 && lon180) ||
1910  (90.0 - fmax(northmost, -1 * southmost) <= FLT_EPSILON);
1911 
1912  if (minNeg >= maxNeg)
1913  { // Entirely in Eastern hemisphere
1914  eastmost = maxPos;
1915  westmost = minPos;
1916  }
1917  else if (minPos >= maxPos)
1918  { // Entirely in Western hemisphere
1919  eastmost = maxNeg;
1920  westmost = minNeg;
1921  }
1922  else if (polar)
1923  { // Polar granule
1924  eastmost = 180;
1925  westmost = -180;
1926  if (northmost > 0)
1927  northmost = 90; // North pole
1928  else
1929  southmost = -90; // South pole
1930  }
1931  else if (lon000)
1932  { // Granule crosses Longitude 0
1933  eastmost = maxPos;
1934  westmost = minNeg;
1935  }
1936  else if (lon180)
1937  { // Granule crosses Longitude 180
1938  eastmost = maxNeg;
1939  westmost = minPos;
1940  }
1941  lonRanges[0] = minNeg;
1942  lonRanges[1] = maxNeg;
1943  lonRanges[2] = minPos;
1944  lonRanges[3] = maxPos;
1945  } /* nobs[ibin] != 0 */
1946  } /* ibin loop */
1947 
1948  /* if no good obs left than bail */
1949  /* ----------------------------- */
1950  if (n_filled_bins == 0)
1951  goto freemem;
1952 
1953  /* Print info on filled row group */
1954  /* ------------------------------ */
1955  printf("%-20s:%8d\n", "# bins in row group", n_bins_in_group);
1956  printf("%-20s:%8d\n", "# filled bins", n_filled_bins);
1957 
1958  /* Write "Bin List" vdata */
1959  /* ---------------------- */
1960  if (is64bit)
1961  {
1962  writeBinList_nc(out_grp, n_filled_bins, (VOIDP)binList64nc);
1963  free(binList64nc);
1964  }
1965  else
1966  {
1967  writeBinList_nc(out_grp, n_filled_bins, (VOIDP)binList32nc);
1968  free(binList32nc);
1969  }
1970 
1971  /* Allocate sum & sum-squared arrays */
1972  /* --------------------------------- */
1973  sum_bin = (double *)calloc(n_bins_in_group, sizeof(double));
1974  sum2_bin = (double *)calloc(n_bins_in_group, sizeof(double));
1975 
1976  /* Loop over all L3 products to fill sum arrays */
1977  /* -------------------------------------------- */
1978  for (iprod = 0; iprod < l3b_nprod; iprod++)
1979  {
1980  memset(sum_bin, 0, n_bins_in_group * sizeof(double));
1981  memset(sum2_bin, 0, n_bins_in_group * sizeof(double));
1982 
1983  /* Process bins */
1984  /* ------------ */
1985  for (ibin = 0; ibin < n_bins_in_group; ibin++)
1986  {
1987  if (nobs[ibin] == 0)
1988  continue;
1989 
1990  /* Process data file by file */
1991  /* ------------------------- */
1992  int32_t npix_file; // for weighting spatial average
1993  double sum_file; // accumulators for each file
1994  double sum2_file;
1995  double area_file;
1996  int16_t prev_file; // index of previous file considered
1997 
1998  // initialize sum for this bin with first file
1999  j = 0;
2000  float32 pixval = data_values[ibin][j * l3b_nprod + iprod];
2001  double pixarea = data_areas[ibin][j];
2002  npix_file = 1;
2003  area_file = pixarea;
2004  sum_file = pixval * pixarea;
2005  sum2_file = pixval * pixarea * pixval * pixarea;
2006  prev_file = file_index[ibin][j];
2007 
2008  // add weighted data for rest of observations (files)
2009  for (j = 1; j < nobs[ibin]; j++)
2010  {
2011  pixval = data_values[ibin][j * l3b_nprod + iprod];
2012  pixarea = data_areas[ibin][j];
2013 
2014  if (file_index[ibin][j] == prev_file)
2015  { // same file
2016  npix_file += 1;
2017  area_file += pixarea;
2018  sum_file += pixval * pixarea;
2019  sum2_file += pixval * pixarea * pixval * pixarea;
2020  }
2021  else
2022  { // new file
2023 
2024  // finalize weighted sum for previous file
2025  sum_bin[ibin] += (sum_file / sqrt(area_file));
2026  sum2_bin[ibin] += (sum2_file / sqrt(area_file));
2027 
2028  // initialize sum for current file
2029  npix_file = 1;
2030  area_file = pixarea;
2031  sum_file = pixval * pixarea;
2032  sum2_file = pixval * pixarea * pixval * pixarea;
2033  prev_file = file_index[ibin][j];
2034  }
2035  } /* observation loop */
2036 
2037  // add data from final file
2038  sum_bin[ibin] += (sum_file / sqrt(area_file));
2039  sum2_bin[ibin] += (sum2_file / sqrt(area_file));
2040 
2041  } /* ibin loop */
2042 
2043  /* Write Product Vdatas */
2044  /* -------------------- */
2045  float *productData = (float *)calloc(2 * n_filled_bins, sizeof(float));
2046 
2047  /* Fill bin data array */
2048  /* ------------------- */
2049  i = 0;
2050  for (ibin = 0; ibin < n_bins_in_group; ibin++)
2051  {
2052  if (nobs[ibin] != 0)
2053  {
2054  productData[i * 2] = sum_bin[ibin];
2055  productData[i * 2 + 1] = sum2_bin[ibin];
2056  // memcpy(&productData[i * 8], &sum_bin[ibin], 4);
2057  // memcpy(&productData[i * 8 + 4], &sum2_bin[ibin], 4);
2058  i++;
2059  }
2060  }
2061 
2062  // char_ptr1 = strchr(prodname[iprod], '/');
2063  // if (char_ptr1 != NULL) *char_ptr1 = '_';
2064  writeBinData_nc(out_grp, i, iprod, productData);
2065  // if (char_ptr1 != NULL) *char_ptr1 = '/';
2066 
2067  free(productData);
2068 
2069  } /* iprod loop */
2070 
2071  /* Write Quality vdata */
2072  /* ------------------- */
2073  if (input.qual_prod[0] != 0)
2074  {
2075  uint8_t *qualityData = (uint8_t *)calloc(n_filled_bins, 1);
2076 
2077  i = 0;
2078  for (ibin = 0; ibin < n_bins_in_group; ibin++)
2079  {
2080  if (nobs[ibin] != 0)
2081  {
2082  memcpy(&qualityData[i], &best_qual[ibin], 1);
2083  i++;
2084  }
2085  }
2086 
2087  writeQuality_nc(out_grp, n_filled_bins, (VOIDP)qualityData);
2088  free(qualityData);
2089  }
2090 
2091  /* Free dynamic memory */
2092  /* ------------------- */
2093  if (sum_bin != NULL)
2094  free(sum_bin);
2095  if (sum2_bin != NULL)
2096  free(sum2_bin);
2097 
2098  /* Compute "begin" & "extent" vdata entries */
2099  /* ---------------------------------------- */
2100  binnum_data = (int64_t *)calloc(n_filled_bins, sizeof(int64_t));
2101 
2102  i = 0;
2103  for (ibin = 0; ibin < n_bins_in_group; ibin++)
2104  {
2105  if (nobs[ibin] != 0)
2106  {
2107  binnum_data[i] = (int64_t)ibin + basebin[krow];
2108 
2109  if (i < 0 || i >= n_filled_bins)
2110  {
2111  printf("Error: %d %d %d %d\n",
2112  i, ibin, n_filled_bins, n_bins_in_group);
2113  }
2114  i++;
2115  }
2116  }
2117 
2118  get_beg_ext(n_filled_bins, binnum_data, basebin, nrows, beg, ext);
2119 
2120  free(binnum_data);
2121 
2122  total_filled_bins += n_filled_bins;
2123 
2124  } /* ********** n_filled_bin > 0 ********** */
2125 
2126  // free(best_qual);
2127 
2128  if (input.verbose)
2129  {
2130  time(&tnow);
2131  tmnow = localtime(&tnow);
2132  printf("krow:%5d After bin processing: %s", krow, asctime(tmnow));
2133  }
2134 
2135  /* Fill BinIndex Vdata */
2136  /* ------------------- */
2137  for (i = 0; i < n_rows_in_group; i++)
2138  {
2139 
2140  i32 = i + krow;
2141  if (i32 < 0 || i32 >= nrows)
2142  {
2143  printf("Error: %d %d\n", i, krow);
2144  exit(-1);
2145  }
2146 
2147  if (is64bit)
2148  {
2149  binIndex64nc[i + krow].start_num = basebin[i32];
2150  binIndex64nc[i + krow].begin = beg[i32];
2151  binIndex64nc[i + krow].extent = ext[i32];
2152  binIndex64nc[i + krow].max = shape->getNumCols(i32);
2153  }
2154  else
2155  {
2156  binIndex32nc[i + krow].start_num = basebin[i32];
2157  binIndex32nc[i + krow].begin = beg[i32];
2158  binIndex32nc[i + krow].extent = ext[i32];
2159  binIndex32nc[i + krow].max = shape->getNumCols(i32);
2160  }
2161 
2162  } /* row_group loop */
2163 
2164  /* End Bin Index Fill */
2165 
2166  freemem:
2167  /* Free Dynamic Memory */
2168  /* ------------------- */
2169  if (bin_flag != NULL)
2170  free(bin_flag);
2171  if (tilt != NULL)
2172  free(tilt);
2173  if (qual != NULL)
2174  free(qual);
2175  if (nscenes != NULL)
2176  free(nscenes);
2177  if (lastfile != NULL)
2178  free(lastfile);
2179 
2180  for (ifile = 0; ifile < nfiles; ifile++)
2181  {
2182  if (scan_in_rowgroup[ifile] != NULL)
2183  free(scan_in_rowgroup[ifile]);
2184  }
2185 
2186  if (input.verbose)
2187  {
2188  time(&tnow);
2189  tmnow = localtime(&tnow);
2190  printf("krow:%5d Befre free dynic mem: %s", krow, asctime(tmnow));
2191  }
2192 
2193  for (i = 0; i < n_bins_in_group; i++)
2194  {
2195  if (file_index[i] != NULL)
2196  free(file_index[i]);
2197  if (data_values[i] != NULL)
2198  free(data_values[i]);
2199  if (data_areas[i] != NULL)
2200  free(data_areas[i]);
2201  if (data_quality[i] != NULL)
2202  free(data_quality[i]);
2203  if (time_value[i] != NULL)
2204  free(time_value[i]);
2205  }
2206 
2207  if (input.verbose)
2208  {
2209  time(&tnow);
2210  tmnow = localtime(&tnow);
2211  printf("krow:%5d After free dynic mem: %s", krow, asctime(tmnow));
2212  }
2213 
2214  free(data_values);
2215  free(data_areas);
2216  free(data_quality);
2217  free(time_value);
2218  free(nobs);
2219  free(allocated_space);
2220  free(file_index);
2221 
2222  krow += n_rows_in_group;
2223 
2224  } /* ========== End krow (Main) loop ========== */
2225 
2226  if (input.verbose)
2227  {
2228  time(&tnow);
2229  tmnow = localtime(&tnow);
2230  printf("krow:%5d %s", krow, asctime(tmnow));
2231  }
2232  printf("total_filled_bins: %ld\n", (long int)total_filled_bins);
2233 
2234  // {
2235  // std::cout << "I/O (readL2) takes " << total_count/1000.0 << " seconds" << std::endl;
2236  // }
2237 
2238  if (total_filled_bins == 0)
2239  {
2240  strcpy(buf, "rm -f ");
2241  strcat(buf, input.ofile);
2242  strcat(buf, "*");
2243  printf("%s\n", buf);
2244  system(buf);
2245  ret_status = 110;
2246  exit(ret_status);
2247  }
2248 
2249  /* Write and Close BinIndex Vdata */
2250  /* ------------------------------ */
2251  if (is64bit)
2252  writeBinIndex_nc(out_grp, nrows, binIndex64nc);
2253  else
2254  writeBinIndex_nc(out_grp, nrows, binIndex32nc);
2255 
2256  /*
2257  * determine list of files actually used in the bin output
2258  */
2259  if (input.fileuse[0] != 0)
2260  {
2261  fp2 = fopen(input.fileuse, "w");
2262  for (ifile = 0; ifile < nfiles; ifile++)
2263  {
2264  if (brk_scan[ifile] != -9999)
2265  {
2266  fileused[ifile] = 1;
2267  fprintf(fp2, "%s\n", l2_str[ifile].filename);
2268  }
2269  }
2270  fclose(fp2);
2271  }
2272  else
2273  {
2274  for (ifile = 0; ifile < nfiles; ifile++)
2275  {
2276  if (brk_scan[ifile] != -9999)
2277  {
2278  fileused[ifile] = 1;
2279  }
2280  }
2281  }
2282 
2283  /* Read and write global attributes */
2284  /* -------------------------------- */
2285  if (input.verbose)
2286  printf("Writing Global Attributes\n");
2287 
2288  status = reopenL2(0, &l2_str[0]);
2289  status = readL2meta(&meta_l2, 0);
2290 
2291  strcpy(meta_l3b.product_name, input.ofile);
2292  if (meta_l2.sensor_name != NULL)
2293  strcpy(meta_l3b.sensor_name, meta_l2.sensor_name);
2294  strcpy(meta_l3b.sensor, meta_l2.sensor);
2295  meta_l3b.sensorID = sensorID[0];
2296  strcpy(meta_l3b.title, meta_l3b.sensor_name);
2297  strcat(meta_l3b.title, " Level-3 Binned Data");
2298 
2299  if (input.suite) {
2300  strcat(meta_l3b.title, " ");
2301  strncat(meta_l3b.title, input.suite, sizeof(meta_l3b.title) - strlen(meta_l3b.title) - strlen(input.suite));
2302  }
2303 
2304  if (meta_l2.mission != NULL)
2305  strcpy(meta_l3b.mission, meta_l2.mission);
2306  strcpy(meta_l3b.prod_type, input.prodtype);
2307 
2308  strcat(meta_l3b.pversion, input.pversion);
2309  strcat(meta_l3b.soft_name, "l2bin");
2310  strcat(meta_l3b.soft_ver, VERSION);
2311 
2312  /*
2313  * loop through the input files
2314  *
2315  * set some ridiculous start and end times to get the ball rolling...
2316  * if this code is still in use in the 22nd century, damn we're good!
2317  * ...but this line will need tweaking.
2318  *
2319  */
2320  meta_l3b.end_orb = 0;
2321  meta_l3b.start_orb = 1000000;
2322  double startTime = yds2unix(2100, 1, 0);
2323  double endTime = yds2unix(1900, 1, 0);
2324  double tmpTime;
2325  strcpy(meta_l3b.infiles, basename(l2_str[0].filename));
2326  for (i = 0; i < nfiles; i++)
2327  {
2328  if (fileused[i] == 1)
2329  {
2330  // update orbit numbers
2331  if (l2_str[i].orbit > meta_l3b.end_orb)
2332  meta_l3b.end_orb = l2_str[i].orbit;
2333  if (l2_str[i].orbit < meta_l3b.start_orb)
2334  meta_l3b.start_orb = l2_str[i].orbit;
2335 
2336  // update start/end times
2337  tmpTime = yds2unix(l2_str[i].syear,
2338  l2_str[i].sday,
2339  (float64)(l2_str[i].smsec / 1000.0));
2340  if (tmpTime < startTime)
2341  startTime = tmpTime;
2342  tmpTime = yds2unix(l2_str[i].eyear,
2343  l2_str[i].eday,
2344  (float64)(l2_str[i].emsec / 1000.0));
2345  if (tmpTime > endTime)
2346  endTime = tmpTime;
2347  }
2348  // append to the file list - all files input, not just those used in the binning.
2349  if (i > 0)
2350  {
2351  strcat(meta_l3b.infiles, ",");
2352  strcat(meta_l3b.infiles, basename(l2_str[i].filename));
2353  }
2354  }
2355 
2356  meta_l3b.startTime = startTime;
2357  meta_l3b.endTime = endTime;
2358 
2359  strcpy(meta_l3b.binning_scheme, "Integerized Sinusoidal Grid");
2360 
2361  meta_l3b.north = northmost;
2362  meta_l3b.south = southmost;
2363  meta_l3b.east = eastmost;
2364  meta_l3b.west = westmost;
2365 
2366  strcpy(buf, argv[0]);
2367  for (i = 1; i < argc; i++)
2368  {
2369  strcat(buf, " ");
2370  strcat(buf, argv[i]);
2371  }
2372  strcpy(meta_l3b.proc_con, buf);
2373  strcpy(meta_l3b.input_parms, input.parms);
2374  strcpy(meta_l3b.flag_names, input.flaguse);
2375 
2376  buf[0] = 0;
2377  for (iprod = 0; iprod < l3b_nprod; iprod++)
2378  {
2379  char *prodName = strdup(l3_prodname[iprod].c_str());
2380  getL3units(&l2_str[0], 0, prodName, units);
2381  strcat(buf, prodName);
2382  free(prodName);
2383  strcat(buf, ":");
2384  strcat(buf, units);
2385  strcat(buf, ",");
2386  if (strlen(buf) >= MD_ATTRSZ)
2387  {
2388  printf("-E- units metadata is too long. Remove some products from product list.\n");
2389  exit(EXIT_FAILURE);
2390  }
2391  }
2392  buf[strlen(buf) - 1] = 0;
2393  strcpy(meta_l3b.units, buf);
2394 
2395  meta_l3b.data_bins = total_filled_bins;
2396  int64_t totalNumBins;
2397  totalNumBins = basebin[nrows] - 1;
2398  meta_l3b.pct_databins = (float)(total_filled_bins) / totalNumBins * 100.0;
2399 
2400  strcpy(meta_l3b.doi, input.doi);
2401 
2402  const char *keywordStr = getGCMDKeywords(input.suite);
2403  if (keywordStr)
2404  strcpy(meta_l3b.keywords, keywordStr);
2405 
2406  time(&tnow);
2407  strcpy(meta_l3b.ptime, unix2isodate(tnow, 'G'));
2409  write_l3b_meta_netcdf4(ds_id, &meta_l3b, is64bit);
2410 
2411  freeL2meta(&meta_l2);
2412 
2413  /* Free Dynamic Memory */
2414  /* ------------------- */
2415  if (input.verbose)
2416  printf("Freeing Dynamic Memory\n");
2417  free(fileused);
2418 
2419  if (is64bit)
2420  free(binIndex64nc);
2421  else
2422  free(binIndex32nc);
2423 
2424  free(ext);
2425  free(beg);
2426 
2427  free(basebin);
2428 
2429  for (ifile = 0; ifile < nfiles; ifile++)
2430  {
2431  if (bscan_row[ifile] != NULL)
2432  free(bscan_row[ifile]);
2433  if (escan_row[ifile] != NULL)
2434  free(escan_row[ifile]);
2435 
2436  if (numer[ifile] != NULL)
2437  free(numer[ifile]);
2438  if (denom[ifile] != NULL)
2439  free(denom[ifile]);
2440  }
2441 
2442  /* Close L2 input files and free allocated L2 memory */
2443  /* ------------------------------------------------- */
2444  if (input.verbose)
2445  printf("Freeing L2 arrays\n");
2446  closeL2(&l2_str[0], 0);
2447  for (ifile = 0; ifile < nfiles; ifile++)
2448  {
2449  freeL2(&l2_str[ifile]);
2450  }
2451  freeL2(NULL);
2452 
2453  /* Close L3B output file */
2454  /* --------------------- */
2455  status = nc_close(fileid_w);
2456 
2457  return ret_status;
2458 }
2459 
2460 int64_t getbinnum(int32_t ifile, int32_t ipixl)
2461 {
2462 
2463  return shape->latlon2bin(l2_str[ifile].latitude[ipixl], l2_str[ifile].longitude[ipixl]);
2464 }
virtual int64_t getNumBins() const
Definition: L3Shape.cpp:63
void setupflags(char *flagdef, char *flaguse, uint32_t *flagusemask, uint32_t *required, int *status, int32_t *BITS)
Definition: setupflags.c:7
@ L2PixelOff
Definition: readL2scan.h:163
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
void addPixelToBin(int32_t ifile, int32_t ipixl, uint64_t bin, bool is_l2_flags_defined=true, double areaFrac=1.0)
Definition: l2bin.cpp:120
int16_t qual_prod_index[MAXNFILES]
Definition: l2bin.cpp:94
#define MAXALLOCPERBIN
Definition: l2bin.cpp:54
virtual void constrainRowCol(int32_t &row, int32_t &col) const =0
virtual int64_t latlon2bin(double lat, double lon) const =0
int16 eday
Definition: l1_czcs_hdf.c:17
int32_t reopenL2(int32_t fileindex, l2_prod *l2_str)
Definition: readL2scan.c:1050
char units[MD_ATTRSZ]
Definition: meta_l3b.h:28
char binning_scheme[SM_ATTRSZ]
Definition: meta_l3b.h:48
int j
Definition: decode_rs.h:73
const char * getGCMDKeywords(const char *suite)
bg::model::box< Point_t > Box_t
Definition: l2bin.cpp:104
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:405
const std::string & get_platform() const
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:707
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:93
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:163
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:453
int32_t * bin_flag
Definition: l2bin.cpp:79
int defineBinData_nc(int32_t deflate, int32_t grpid, int32_t nprods, char **prodnames)
int16_t * qual
Definition: l2bin.cpp:80
#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:2151
bg::model::point< double, 2, bg::cs::geographic< bg::degree > > Point_t
Definition: l2bin.cpp:102
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
#define LTILT_DIMS_2
Definition: l2bin.cpp:53
char keywords[SM_ATTRSZ]
Definition: meta_l3b.h:55
bg::model::polygon< Point_t > Polygon_t
Definition: l2bin.cpp:103
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:301
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:107
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:99
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:2325
int32_t tiltstate
Definition: l2bin.cpp:87
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:1336
ds_format_t fftype
Definition: dfutils.h:31
int32 smsec
Definition: l1_czcs_hdf.c:16
bool skip_DL(float lon, int side, int night_flag, time_t end_day, time_t beg_day)
checks if lon should be skipped
Definition: l2bin.cpp:692
@ string
int32_t closeL2(l2_prod *l2_str, int32_t ifile)
Definition: readL2scan.c:1995
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:163
void init_rowgroup_cache()
Definition: readL2scan.c:287
int setlinebuf(FILE *stream)
int16_t composite_l3prod_index
Definition: l2bin.cpp:96
int defineQuality_nc(int32_t deflate, int32_t grpid)
const float * get_elat()
int16_t composite_prod_index[MAXNFILES]
Definition: l2bin.cpp:95
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
std::multimap< std::string, netCDF::NcVar > find_variables_geo_physical(const std::string &nc_path, std::string &prod_list)
int32_t krow
Definition: l2bin.cpp:84
void enableL2PixelArea(enum L2PixelMode_t val)
Definition: readL2scan.c:213
double endTime
Definition: meta_l3b.h:40
bool binIntersectsPixel(int32_t row, int32_t col, Box_t &pixelBox, double &areaFrac)
Definition: l2bin.cpp:377
int32_t openL2(const char *fname, const char *plist, l2_prod *l2_str)
Definition: readL2scan.c:316
int32_t n_bins_in_group
Definition: l2bin.cpp:82
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:492
int32_t n_rows_in_group
Definition: l2bin.cpp:83
bg::model::polygon< Point_t > Polygon_t
Definition: get_dataday.cpp:25
void free_rowgroup_cache()
Definition: readL2scan.c:259
virtual int64_t getBaseBin(int32_t row) const =0
#define MAXNUMBERPRODUCTS
Definition: readL2scan.h:10
double resolution
Definition: meta_l3b.h:53
int defineBinIndex_nc(int32_t deflate, int32_t grpid)
int32_t resolve2binRows(const char *resolve)
convert l2bin style resolve string to number line in L3bin file
bg::model::point< double, 2, bg::cs::spherical_equatorial< bg::degree > > Point_t
Definition: get_dataday.cpp:23
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:2460
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
void set_prodname_3d_to_l2(const std::vector< std::string > &prodparam, l2_prod &l2_str, std::vector< std::string > &l2_prodname, std::vector< std::string > &l3_prodname, std::vector< int32_t > &thirdDimId, std::vector< float > &min_value, std::vector< float > &max_value)
Definition: expand3D.cpp:383
const std::string & get_instrument() const
integer, parameter double
float north
Definition: meta_l3b.h:49
@ DS_NCDF
Definition: dfutils.h:20
std::vector< float32 > min_value
Definition: l2bin.cpp:100
string & trim(string &s, const string &delimiters)
Definition: EnvsatUtil.cpp:29
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
bool set_l2_flags_use(const std::string &flagsuse)
Definition: expand3D.cpp:407
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:68
#define basename(s)
Definition: l0chunk_modis.c:29
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)
std::vector< int32_t > thirdDimId
Definition: l2bin.cpp:99
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:88
virtual int32_t getNumRows() const
Definition: L3Shape.cpp:55
#define BAD_FLT
Definition: jplaeriallib.h:19
float south
Definition: meta_l3b.h:50
subroutine geometry
int16_t * numer[MAXNFILES]
Definition: l2bin.cpp:93
char proc_con[MD_ATTRSZ]
Definition: meta_l3b.h:35
int16_t * side
Definition: l1a_seawifs.c:88
#define MD_ATTRSZ
Definition: meta_l3b.h:8
const float * get_slat()
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
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 the sector rotation may actually occur during one of the scans earlier than the one where it is first reported As a the b1 values are about the LOCALGRANULEID metadata should have an extension NRT It is requested to identify the NRT to fill pixels affected by dead subframes with a special value Output the metadata of noisy and dead subframe Dead Subframe EV and Detector Quality Flag2 Removed the function call of Fill_Dead_Detector_SI to stop interpolating SI values for dead but also for all downstream products for science test only Changes from v5 which will affect scientific to conform to MODIS requirements Removed the Mixed option from the ScanType in the code because the L1A Scan Type is never Mixed Changed for ANSI C compliance and comments to better document the fact that when the HDF_EOS metadata is stricly the and products are off by and in the track respectively Corrected some misspelling of RCS swir_oob_sending_detector to the Reflective LUTs to enable the SWIR OOB correction detector so that if any of the sending detectors becomes noisy or non near by good detectors from the same sending band can be specified as the substitute in the new look up table Code change for adding an additional dimension of mirror side to the Band_21_b1 LUT to separate the coefficient of the two mirror sides for just like other thermal emissive so that the L1B code can calibrate Band scan to scan with mirror side dependency which leads better calibration result Changes which do not affect scientific when the EV data are not provided in this Crosstalk Correction will not be performed to the Band calibration data Changes which do not affect scientific and BB_500m in L1A Logic was added to turn off the or to spatial aggregation processes and the EV_250m_Aggr1km_RefSB and EV_500m_Aggr1km_RefSB fields were set to fill values when SDSs EV_250m and EV_500m are absent in L1A file Logic was added to skip the processing and turn off the output of the L1B QKM and HKM EV data when EV_250m and EV_500m are absent from L1A In this the new process avoids accessing and reading the and L1A EV skips and writing to the L1B and EV omits reading and subsampling SDSs from geolocation file and writing them to the L1B and omits writing metadata to L1B and EV and skips closing the L1A and L1B EV and SDSs Logic was added to turn off the L1B OBC output when the high resolution OBC SDSs are absent from L1A This is accomplished by skipping the openning the writing of metadata and the closing of the L1B OBC hdf which is Bit in the scan by scan bit QA has been changed Until now
Definition: HISTORY.txt:361
char * unix2isodate(double dtime, char zone)
Definition: unix2isodate.c:10
void copy(double **aout, double **ain, int n)
float pct_databins
Definition: meta_l3b.h:47
int16_t * nscenes
Definition: l2bin.cpp:80
virtual double row2lat(int32_t row) const =0
Definition: dfutils.h:28
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:69
data_t s[NROOTS]
Definition: decode_rs.h:75
char pversion[SM_ATTRSZ]
Definition: meta_l3b.h:30
float l3_attr(const std::string &inp3)
Definition: expand3D.cpp:421
int write_l3b_meta_netcdf4(idDS ds_id, meta_l3bType *meta_l3b, int write64bit)
int32_t getL3units(l2_prod *l2_str, int32_t ifile, char *l3b_prodname, char *units)
Definition: readL2scan.c:2346
float east
Definition: meta_l3b.h:51
int64_t data_bins
Definition: meta_l3b.h:46
int16_t * tilt
Definition: l2bin.cpp:80
char sensor_name[SM_ATTRSZ]
Definition: meta_l3b.h:19
#define abs(a)
Definition: misc.h:90
int i
Definition: decode_rs.h:71
std::vector< float32 > max_value
Definition: l2bin.cpp:100
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:91
int32_t freeL2(l2_prod *l2_str)
Definition: readL2scan.c:2083
#define LG_ATTRSZ
Definition: meta_l3b.h:9
int32_t end_orb
Definition: meta_l3b.h:43
#define VERSION
Definition: l2bin.cpp:106
char soft_name[SM_ATTRSZ]
Definition: meta_l3b.h:32
int k
Definition: decode_rs.h:73
#define MTILT_DIMS_2
Definition: l2bin.cpp:52
float p[MODELMAX]
Definition: atrem_corl1.h:131
float32 f32
Definition: l2bin.cpp:98
void set_mapper(const std::string &input_file_name, std::string &products_requested, std::vector< std::string > &products_requested_separated, const std::string &requested_wavelengths)
readL2, openL2 uses the maps from this module. If the module is not initialized, the function read an...
Definition: expand3D.cpp:235
bg::model::box< Point_t > Box_t
Definition: l1c.cpp:69
l2prod max
int16_t * lastfile
Definition: l2bin.cpp:80
int32_t get_l2prod_index(const l2_prod &l2, const char *prodname)
Definition: l2bin.cpp:109