NASA Logo
Ocean Color Science Software

ocssw V2022
sst_cloud_mask_utils.hpp
Go to the documentation of this file.
1 #ifndef SST_CLOUD_MASK_H
2 #define SST_CLOUD_MASK_H
3 
4 #include "l2_struc.h"
5 #include "l1q_struc.h"
6 #include <vector>
7 #include <iostream>
8 #include "l2_struc.h"
9 #include "l1q_struc.h"
10 #include <boost/variant.hpp>
11 #include <string>
12 #include <utility>
13 #include "flags_sst.h"
14 #include <assert.h>
15 #include <algorithm>
16 #include <fstream>
17 #include <cstdlib>
18 #include <rapidjson/document.h>
19 #include <rapidjson/istreamwrapper.h>
20 #include <rapidjson/writer.h>
21 #include <rapidjson/stringbuffer.h>
22 #include <rapidjson/ostreamwrapper.h>
23 #include <unordered_set>
24 #include <chrono>
25 #include <unordered_map>
26 #include <memory>
27 #include <netcdf>
28 #include "sst_dsdi.h"
29 #include <boost/algorithm/string.hpp>
30 #include <boost/unordered_map.hpp>
31 // parameters
32 namespace envset {
37 
41  const std::unordered_map<std::string, std::string (*)(const instr *)> get_sst_coeffs_path = {{std::string("SST"),
42  [](const instr *input) {
43  return std::string(
44  input->sstcoeffile);
45  }},
46  {std::string("SST4"),
47  [](const instr *input) {
48  return std::string(
49  input->sst4coeffile);
50  }},
51  {std::string("SST3"),
52  [](const instr *input) {
53  return std::string(
54  input->sst3coeffile);
55  }}};
56 
60  const std::unordered_map<std::string, std::string (*)(const instr *)> get_sses_coeffs_path = {{std::string("SST"),
61  [](const instr *input) {
62  return std::string(
63  input->sstssesfile);
64  }},
65  {std::string("SST4"),
66  [](const instr *input) {
67  return std::string(
68  input->sst4ssesfile);
69  }},
70  {std::string("SST3"),
71  [](const instr *input) {
72  return std::string(
73  input->sst3ssesfile);
74  }}};
75 
76 }
81 namespace cldmsk {
82  const float cldthresh = -1.0;
83  const std::unordered_map<std::string, float> cldthresh_list = {{"modis", 0.01},
84  {"viirs", 0.04}}; // threshold for cold SST test
85  const std::unordered_map<std::string, const std::unordered_map<std::string, int>> bands_set = // bands WV for VIIRS/MODIS needed to produce SST/Cloud mask
86  {{"modis",
87  {{"ibred", 678},
88  {"ib07", 748},
89  {"ib16", 1640},
90  {"ib37", 3750},
91  {"ib39", 3959},
92  {"ib40", 4050},
93  {"ib67", 6715},
94  {"ib73", 7325},
95  {"ib85", 8550},
96  {"ib11", 11000},
97  {"ib12", 12000}}},
98  {"viirs",
99  {{"ibred", 672},
100  {"ib07", 748},
101  {"ib16", 1601},
102  {"ib37", 3750},
103  {"ib40", 4050},
104  {"ib85", 8550},
105  {"ib11", 11000},
106  {"ib12", 12000}}}};
107  const std::unordered_map<std::string, size_t> bt_box_sizes = {{"modis", 3},
108  {"viirs", 5}}; // bt box size
109  const float hisenz = 55.0; // senz thresholds
110  const float vhisenz = 75.0;
111  const float vhisenzv2 = 65.0;
112  const float solznight = 85.0;
113  const float Btmin = -4.0; // BT ranges for BT range cloud mask tests
114  const float Btmax = 37.0;
115  const float Btmax40 = 35.0;
116  const float SSTmin = -1.8; // SST ranges for SST cloud mask tests
117  const float SSTmax = 40.0;
118  const float SSTmaxn = 37.0;
119  const float glintmax = 0.005;
120  const float dBtmin = 0.0; // dBT ranges for the cloud mask tests
121  const float dBtmax = 3.6;
122  const float dBt4min = 0.0;
123  const float dBt4max = 8.0;
124  const float SSTdiff = 3.0;
125  const float SSTvdiff = 5.0; // dSST ranges (thresholds) for the cloud mask tests
126  const float SST4diff1 = -0.8;
127  const float SST4diff2 = -1.0;
128  const float SST3diff1 = 0.8;
129  const float SST3diff2 = 1.0;
130  const float Bt11unif1 = 0.7; // BT thresholds for BT uniformity tests
131  const float Bt12unif1 = 0.7;
132  const float Bt11unif2 = 1.9;
133  const float Bt12unif2 = 1.9;
134  const float Bt37unif1 = 0.7;
135  const float Bt37unif2 = 1.9;
136  const float Bt39unif1 = 0.7;
137  const float Bt40unif1 = 0.7;
138  const float Bt39unif2 = 1.9;
139  const float Bt40unif2 = 1.9;
140  const float dBtrefmin = -1.1;
141  const float dBtrefmax = 10.0;
142  const float equatorialNorth = 30.0; // the ranges where SST diff tests are applied
143  const float equatorialSouth = -10.0;
144  const float equatorialWest = -105.0;
145  const float equatorialEast = 105.0;
146  const float invalid_val = (BAD_FLT + 0.1); // bad vals ranges
147  const float max_bt = BT_HI - 0.1;
148  const float min_bt = BT_LO + 0.1;
149  const float min_lt = BAD_FLT + 1.0;
150  const float max_lt = BT_HI - 1.0;
151  const std::unordered_set<int> modis_sensors = {MODIST,
152  MODISA}; // sensors/platforms INT values and corresponding string values
153  const std::unordered_set<int> viirs_sensors = {VIIRSJ1, VIIRSJ2, VIIRSN};
154  const std::unordered_map<std::string, std::unordered_set<int>> all_sensors = {{"modis", modis_sensors},
155  {"viirs", viirs_sensors}};
156  const std::unordered_map<int, std::string> platforms = {{MODIST, "terra"},
157  {MODISA, "aqua"},
158  {VIIRSN, "npp"},
159  {VIIRSJ1, "j1"},
160  {VIIRSJ2, "j2"}};
161  const double scan_time_modis_t_day1 = 972777600; // 29 Oct 2000, MODIS TERRA correction dates
162  const double scan_time_modis_t_day2 = 993945600; // 1 Jul 2001
163  const float el_corr_modis_t_1 = 0.4452; // MODIS TERRA correction values
164  const float el_corr_modis_t_2 = 0.2593;
165 
170  std::vector<float> & month_data();
171 
178 
187  template<class T>
188  std::ostream &operator<<(std::ostream &stream, const std::vector<T> &data) {
189  for (const auto &el: data) {
190  stream << el << " ";
191  }
192  stream << "\n";
193  return stream;
194  }
195 
206  float btrefdiffv6(int32_t ip, float BT39, float BT40, const l1str *l1rec, int fullscanpix);
207 
212  };
213 
214  // Functions ptr below get the values of BT/RHOs. Compute a value of BT/RHO or its mask at a given pixel
215  typedef bool (*get_valid)(const l1str &, int, int, int); //
216  typedef float (*get_value)(const l1str &, int, int, int);
217 
218  const get_valid cirrus_mask = [](const l1str &l1str_, int pixel, int nbands, int ib) {
219  return l1str_.rho_cirrus[pixel] > invalid_val;
220  };
221  const get_value cirrus_value = [](const l1str &l1str_, int pixel, int nbands,
222  int ib) { return l1str_.rho_cirrus[pixel]; };
223  // bts
224  const get_valid bt_mask = [](const l1str &l1str_, int pixel, int nbands, int ib) {
225  const float bt = l1str_.Bt[pixel * nbands + ib];
226  return bt > min_bt && bt < max_bt;
227  };
228  const get_value bt_value = [](const l1str &l1str_, int pixel, int nbands, int ib) {
229  return l1str_.Bt[pixel * nbands + ib];
230  };
231  // rho (RSMAS)
232  const get_valid rho_mask = [](const l1str &l1str_, int pixel, int nbands, int ib) {
233  const float lt = l1str_.Lt[pixel * nbands + ib];
234  return lt > min_lt && lt < max_lt;
235  };
236  const get_value rho_value = [](const l1str &l1str_, int pixel, int nbands, int ib) {
237  const float lt = l1str_.Lt[pixel * nbands + ib];
238  const float fo = l1str_.Fo[ib];
239  const float csolz = l1str_.csolz[pixel];
240  const float rho = M_PI * lt / fo / csolz;
241  return rho;
242  };
243 
244  const get_value cldrh_value = [](const l1str &l1str_, int pixel, int nbands, int ib) {
245  const float lt = l1str_.Lt[pixel * nbands + ib];
246  const float fo = l1str_.Fo[ib];
247  // const float csolz = l1str_.csolz[pixel];
248  const float tg_sol = l1str_.tg_sol[pixel * nbands + ib];
249  const float tg_sen = l1str_.tg_sen[pixel * nbands + ib];
250  const float t_sen = l1str_.t_sen[pixel * nbands + ib];
251  const float t_sol = l1str_.t_sol[pixel * nbands + ib];
252  const float cldrh = M_PI * lt / fo / tg_sol / tg_sen / t_sol / t_sen; // /csolz
253  return cldrh;
254  };
255  const get_valid cldrh_mask = [](const l1str &l1str_, int pixel, int nbands, int ib) {
256  const float lt = l1str_.Lt[pixel * nbands + ib];
257  return lt > 0.0 && lt < max_lt;
258  };
259  // other vars (not bandd)
260  const std::unordered_map<std::string, float *(*)(const l1str &)> get_non_BT_vars =
261  {
262  {"solz", [](const l1str &l1rec) { return l1rec.solz; }},
263  {"senz", [](const l1str &l1rec) { return l1rec.senz; }},
264  {"wv", [](const l1str &l1rec) { return l1rec.wv; }},
265  {"glintcoef", [](const l1str &l1rec) { return l1rec.glint_coef; }},
266  {"lat", [](const l1str &l1rec) { return l1rec.lat; }},
267  {"lon", [](const l1str &l1rec) { return l1rec.lon; }},
268  {"month", [](const l1str &l1rec) { return month_data().data(); }}};
269 
280  void get_var_mask(int *mask, const l1str &l1str, size_t npix, int nbands, int ib, get_valid get_mask);
281 
285  void get_var_vals(float *BT, const l1str &l1str, size_t npix, int nbands, int ib, get_value get_val);
286 
296  void get_window_1D_max(float *maxs_1d, const float *inp1d, const int *mask1d, size_t npix, size_t radius);
297 
308  void
309  get_total_2d_max(float *max_global, const float *inp2d, size_t npix, size_t nscan, size_t center, size_t radius);
310 
314  void get_window_1D_min(float *mins_1d, const float *inp1d, const int *mask1d, size_t npix, size_t radius);
315 
319  void
320  get_total_2d_min(float *min_global, const float *inp2d, size_t npix, size_t nscan, size_t center, size_t radius);
321 
335  void get_std_box(const float *inp2d, const int *mask2d, size_t npix, size_t nscan, size_t radius_x, size_t radius_y,
336  float *out, size_t center, l1qstr *l1qrec);
337 
343  struct SSESData {
344  size_t nsst;
345  size_t nday;
346  size_t nquar;
347  size_t nsenz;
348  size_t ndiff;
349  size_t nlat;
350  size_t nqual;
351  size_t total_size;
352  size_t npix;
353  std::vector<float> sst;
354  std::vector<float> senz;
355  std::vector<float> diff;
356  std::vector<float> lat;
357  std::vector<float> qual;
358  // LUTS
359  std::vector<float> bias;
360  std::vector<float> stdv;
361  std::vector<float> bias_mean;
362  std::vector<int16_t> counts;
363 
364  std::vector<float> bias_out;
365  std::vector<float> stdv_out;
366  std::vector<float> bias_mean_out;
367  std::vector<int16_t> counts_out;
368 
375  SSESData(const std::string &nc_path_LUTs, size_t npix) : npix(npix) {
376  try {
377  netCDF::NcFile sst_sses_lut_bias(nc_path_LUTs, netCDF::NcFile::read);
378  {
379  std::cout << "Reading SSES LUT netCDF file : " << nc_path_LUTs << std::endl;
380  }
381  const auto vars = sst_sses_lut_bias.getVars();
382  for (const auto &var: vars) {
383  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar(var.first);
384  auto dims = var_sses.getDims();
385  for (const auto &dim: dims) {
386  if (dim.getName() == "nlat") {
387  nlat = dim.getSize();
388  }
389  if (dim.getName() == "nsatz") {
390  nsenz = dim.getSize();
391  }
392  if (dim.getName() == "ndiff") {
393  ndiff = dim.getSize();
394  }
395  if (dim.getName() == "nquar") {
396  nquar = dim.getSize();
397  }
398  if (dim.getName() == "nqual") {
399  nqual = dim.getSize();
400  }
401  if (dim.getName() == "nsst") {
402  nsst = dim.getSize();
403  }
404  if (dim.getName() == "nday") {
405  nday = dim.getSize();
406  }
407  }
408  }
409  {
410  total_size = nday * nsst * nquar * nqual * nlat * ndiff * nsenz;
411  }
412  bias.resize(total_size);
413  stdv.resize(total_size);
414 
415  sst.resize(nsst);
416  senz.resize(nsenz);
417  diff.resize(ndiff);
418  lat.resize(nlat);
419  qual.resize(nqual);
420  {
421  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("bias");
422  var_sses.getVar(bias.data());
423  }
424  {
425  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("bias_mean");
426  if (!var_sses.isNull()) {
427  bias_mean.resize(total_size);
428  var_sses.getVar(bias_mean.data());
429  }
430  }
431  {
432  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("stdv");
433  var_sses.getVar(stdv.data());
434  }
435  {
436  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("counts");
437  if (!var_sses.isNull()) {
438  counts.resize(total_size);
439  var_sses.getVar(counts.data());
440  }
441  }
442  {
443  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("sst");
444  var_sses.getVar(sst.data());
445  }
446  {
447  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("senz");
448  var_sses.getVar(senz.data());
449  }
450  {
451  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("lat");
452  var_sses.getVar(lat.data());
453  }
454  {
455  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("qual");
456  var_sses.getVar(qual.data());
457  }
458  {
459  netCDF::NcVar var_sses = sst_sses_lut_bias.getVar("BTdiff");
460  var_sses.getVar(diff.data());
461  }
462  {
463  bias_out = std::vector<float>(npix, BAD_FLT);
464  stdv_out = std::vector<float>(npix, BAD_FLT);
465 
466  if (!counts.empty())
467  counts_out = std::vector<int16_t>(npix, 0);
468  if (!bias_mean.empty())
469  bias_mean_out = std::vector<float>(npix, BAD_FLT);
470  }
471  }
472  catch (...) {
473  std::cerr << "Error opening and reading SSES LUT file : " << nc_path_LUTs << std::endl;
474  std::cerr << "Exiting ... " << nc_path_LUTs << std::endl;
475  exit(EXIT_FAILURE);
476  }
477  };
478 
479  size_t
480  operator()(size_t iqual, size_t ilat, size_t idiff, size_t isenz, size_t iquar, size_t iday, size_t isst) {
481  size_t index = isst + iday * nsst + iquar * nday * nsst + isenz * nquar * nsst * nday +
482  idiff * nday * nsst * nquar * nsenz + ilat * nday * nsst * nquar * nsenz * ndiff +
483  iqual * nday * nsst * nquar * nsenz * nlat * ndiff;
484  return index;
485  }
486 
500  void
501  calculate_sses_bias_stats(const float *diff, const float *sst, const float *solz, const float *senz, size_t day,
502  const float *lat, const int8_t *qual, const int16_t *flags, const float *glint) {
503  if (!bias_mean_out.empty())
504  std::fill(bias_mean_out.begin(), bias_mean_out.end(), BAD_FLT);
505  std::fill(stdv_out.begin(), stdv_out.end(), BAD_FLT);
506  std::fill(bias_out.begin(), bias_out.end(), BAD_FLT);
507  if (!counts_out.empty())
508  std::fill(counts_out.begin(), counts_out.end(), 0);
509  size_t iquar, iday, ilat, isenz, isst, idiff;
510  size_t i;
511  for (size_t i_p = 0; i_p < npix; i_p++) {
512  size_t iqual = *(qual + i_p);
513  if (iqual >= nqual)
514  continue;
515  if (iqual == 1 && (glint[i_p] > glintmax) && ((flags[i_p] & SSTF_BTNONUNIF) == 0) &&
516  ((flags[i_p] & SSTF_HISENZ) == 0)) {
517  iqual = 0;
518  }
519  iquar = day / 91;
520  iday = 0;
521  if (solz[i_p] < solznight) {
522  iday = 1;
523  }
524  if (iday > nday - 1)
525  continue;
526  for (i = 1; i < nlat; i++)
527  if (*(lat + i_p) < this->lat.at(i))
528  break;
529  ilat = i - 1;
530  for (i = 1; i < nsenz; i++)
531  if (*(senz + i_p) < this->senz.at(i))
532  break;
533  isenz = i - 1;
534  for (i = 1; i < nsst; i++)
535  if (*(sst + i_p) < this->sst.at(i))
536  break;
537  isst = i - 1;
538  for (i = 1; i < ndiff; i++)
539  if (*(diff + i_p) <= this->diff.at(i))
540  break;
541  idiff = i - 1;
542  const size_t index = this->operator()(iqual, ilat, idiff, isenz, iquar, iday, isst);
543  bias_out.at(i_p) = bias.at(index);
544  if (!bias_mean_out.empty())
545  bias_mean_out.at(i_p) = bias_mean.at(index);
546  stdv_out.at(i_p) = stdv.at(index);
547  if (!counts_out.empty())
548  counts_out.at(i_p) = counts.at(index);
549  }
550  }
551  };
552 
553  void read_sst_bands(const std::string &sat, std::unordered_map<std::string, int> &bands);
554 }
555 
560 namespace bstats {
565  struct PairOfVars {
566  std::vector<float> difference;
567  std::vector<float> ratio;
568  size_t npix;
569  int current_scan = 0;
570 
572  }
573 
574  PairOfVars(size_t npix) : npix(npix) {
575  }
576 
577  float *operator()() {
578  return difference.data();
579  }
580 
581  float *get_Tdeflong() {
582  return ratio.data();
583  }
584  };
585 
589  struct PairOfVarsBT : public PairOfVars {
590  int ib1, ib2;
591  float offset = 0.00001; // in case if BT_2 == 0
593 
594  };
595 
596  PairOfVarsBT(size_t npix, int ib1, int ib2, l1str &l1str, bool has_difference = false, bool has_ratio = false)
597  : PairOfVars(npix), ib1(ib1), ib2(ib2) {
598  if (has_difference) {
599  difference = std::vector<float>(npix, BAD_FLT);
600  get_difference(l1str);
601  }
602  if (has_ratio) {
603  ratio = std::vector<float>(npix, BAD_FLT);
604  get_ratio(l1str);
605  }
606  }
607 
608  void get_difference(l1str &l1str) {
609  for (size_t i_p = 0; i_p < npix; i_p++) {
610  difference.at(i_p) =
611  cldmsk::bt_value(l1str, i_p, NBANDSIR, ib1) - cldmsk::bt_value(l1str, i_p, NBANDSIR, ib2);
612  }
613  if (!ratio.empty())
614  get_ratio(l1str);
615  }
616 
617  void get_ratio(l1str &l1str) {
618  if (ratio.empty())
619  ratio = std::vector<float>(npix, BAD_FLT);
620  for (size_t i_p = 0; i_p < npix; i_p++) {
621  if (cldmsk::bt_value(l1str, i_p, NBANDSIR, ib1) != 0)
622  ratio.at(i_p) = 1.0 - cldmsk::bt_value(l1str, i_p, NBANDSIR, ib2) /
623  cldmsk::bt_value(l1str, i_p, NBANDSIR, ib1);
624  else if (cldmsk::bt_value(l1str, i_p, NBANDSIR, ib2) == 0) {
625  ratio.at(i_p) = 0.0;
626  } else {
627  ratio.at(i_p) = 1.0 - cldmsk::bt_value(l1str, i_p, NBANDSIR, ib2) / offset;
628  }
629  }
630  }
631  };
632 
636  struct PairOfSST : public PairOfVars {
637  float *sst1;
638  float *sst2;
639 
641  }
642 
643  PairOfSST(size_t npix, float *sst1, float *sst2, bool has_difference) : PairOfVars(npix), sst1(sst1),
644  sst2(sst2) {
645  if (has_difference) {
646  difference = std::vector<float>(npix, BAD_FLT);
647  get_difference();
648  }
649  }
650 
651  void get_difference() {
652  std::transform(sst1, sst1 + npix, sst2, difference.begin(), std::minus<float>());
653  }
654  };
655 
660  struct Stats {
661  int current_scan = 0;
662  size_t npix;
663  size_t nscan;
664  size_t rad_x, rad_y, center;
665  size_t i_s, i_e;
666  l1qstr *l1qrec;
667  std::vector<float> var_box;
668  std::vector<int> mask_box;
669  std::vector<float> var_max, var_min, var_std, var_minmax;
670  std::vector<float> var_max_box, var_min_box;
671  std::unordered_map<std::string, float *> return_vals;
672 
676  void set_std() {
677  if (var_std.empty()) {
678  var_std = std::vector<float>(npix, BAD_FLT);
679  cldmsk::get_std_box(var_box.data(), mask_box.data(), npix, nscan, rad_x, rad_y, var_std.data(), center,
680  l1qrec);
681  return_vals["STD"] = var_std.data();
682  }
683  }
684 
688  void set_min() {
689  if (var_min.empty()) {
690  const float max_possible_val = std::abs(BAD_FLT);
691  var_min = std::vector<float>(npix, max_possible_val);
692  var_min_box = std::vector<float>(npix * nscan, max_possible_val);
693  return_vals["MIN"] = var_min.data();
694  for (size_t i_q = i_s; i_q <= i_e; i_q++) {
695  const size_t index = npix * i_q;
697  mask_box.data() + index, npix, rad_x);
698  }
700  }
701  }
702 
706  void set_max() {
707  if (var_max.empty()) {
708  var_max = std::vector<float>(npix, BAD_FLT);
709  var_max_box = std::vector<float>(npix * nscan, BAD_FLT);
710  return_vals["MAX"] = var_max.data();
711  for (size_t i_q = i_s; i_q <= i_e; i_q++) {
712  const size_t index = npix * i_q;
714  mask_box.data() + index, npix, rad_x);
715  }
717  }
718  }
719 
723  void set_maxmin() {
724  if (var_minmax.empty()) {
725  var_minmax = std::vector<float>(npix, BAD_FLT);
726  return_vals["MAXMIN"] = var_minmax.data();
727  set_max();
728  set_min();
729  std::transform(var_max.begin(), var_max.end(), var_min.begin(), var_minmax.begin(),
730  std::minus<float>());
731  }
732  }
733 
739  if (!var_std.empty()) {
740  cldmsk::get_std_box(var_box.data(), mask_box.data(), npix, nscan, rad_x, rad_y, var_std.data(), center,
741  l1qrec);
742  };
743  if (!var_max.empty()) {
744  for (size_t i_q = i_s; i_q < i_e; i_q++) {
745  const size_t index = (i_q + 1) * npix;
746  std::copy(var_max_box.begin() + index, var_max_box.begin() + index + npix,
747  var_max_box.begin() + index - npix);
748  }
749  const size_t index = npix * i_e;
751  npix, rad_x);
753  };
754  if (!var_min.empty()) {
755  for (size_t i_q = i_s; i_q < i_e; i_q++) {
756  const size_t index = (i_q + 1) * npix;
757  std::copy(var_min_box.begin() + index, var_min_box.begin() + index + npix,
758  var_min_box.begin() + index - npix);
759  }
760  const size_t index = npix * i_e;
762  npix, rad_x);
764  };
765  if (!var_minmax.empty()) {
766  std::transform(var_max.begin(), var_max.end(), var_min.begin(), var_minmax.begin(),
767  std::minus<float>());
768  };
769  };
770 
771  Stats() {};
772 
773  Stats(size_t npix, size_t nscan, size_t rad_x, size_t rad_y, size_t center, l1qstr *l1qrec) : npix(npix),
774  nscan(nscan),
775  rad_x(rad_x),
776  rad_y(rad_y),
777  center(center),
778  l1qrec(l1qrec) {
779  i_s = std::min(std::max(0, (int) center - (int) rad_y), (int) nscan - 1);
780  i_e = std::max(std::min(nscan - 1, center + rad_y), 0ul);
781  }
782 
783  float *operator()(const std::string &key) {
784  return return_vals.at(key);
785  }
786 
787  int *get_valid_mask() {
788  return mask_box.data() + npix * center;
789  }
790  };
791 
795  struct StatsSST : public Stats {
796  float *sst;
797 
798  StatsSST() {};
799 
800  StatsSST(size_t npix, size_t nscan, size_t rad_x, size_t rad_y, size_t center, l1qstr *l1qrec, float *sst,
801  bool std_requested = false, bool max_req = false, bool min_req = false, bool min_max_req = false)
803  return_vals["VAR"] = sst;
804  if (std_requested || max_req || min_req || min_max_req) {
805  var_box = std::vector<float>(npix * nscan, BAD_FLT);
806  mask_box = std::vector<int>(npix * nscan, 0);
807  for (size_t index = 0; index < npix * nscan; index++) {
808  var_box.at(index) = *(sst - npix * center + index);
810  }
811  if (std_requested)
812  set_std();
813  if (min_max_req)
814  set_maxmin();
815  else if (min_req)
816  set_min();
817  else if (max_req)
818  set_max();
819  }
820  }
821 
822  void rearrange() {
823  if (!var_box.empty()) {
824  for (size_t index = 0; index < npix * nscan; index++) {
825  var_box.at(index) = *(sst - npix * center + index);
827  }
828  }
829  rearrange_stats();
830  }
831  };
832 
836  struct StatsVarBand : public Stats {
837  int ib;
838  size_t nbands;
841 
843 
844  StatsVarBand(size_t npix, size_t nscan, size_t rad_x, size_t rad_y, size_t center, l1qstr *l1qrec, int ib,
845  size_t nbands, cldmsk::get_valid get_mask, cldmsk::get_value get_val, bool std_requested = false,
846  bool max_requested = false, bool min_requested = false, bool min_max_requested = false) : Stats(
848  get_val(get_val) {
849  var_box = std::vector<float>(npix * nscan, BAD_FLT);
850  mask_box = std::vector<int>(npix * nscan, 0);
851  for (size_t i_q = i_s; i_q <= i_e; i_q++) {
852  const l1str &q_ref = l1qrec->r[i_q];
853  cldmsk::get_var_vals(var_box.data() + npix * i_q, q_ref, npix, nbands, ib, get_val);
854  cldmsk::get_var_mask(mask_box.data() + npix * i_q, q_ref, npix, nbands, ib, get_mask);
855  }
856  return_vals["VAR"] = var_box.data() + npix * center;
857  if (std_requested) {
858  set_std();
859  }
860  if (min_max_requested) {
861  max_requested = true;
862  min_requested = true;
863  set_maxmin();
864  }
865  if (max_requested) {
866  set_max();
867  }
868  if (min_requested) {
869  set_min();
870  }
871  }
872 
873  void rearrange() {
874  for (size_t i_q = i_s; i_q < i_e; i_q++) {
875  const size_t index = (i_q + 1) * npix;
876  std::copy(var_box.begin() + index, var_box.begin() + index + npix, var_box.begin() + index - npix);
877  std::copy(mask_box.begin() + index, mask_box.begin() + index + npix, mask_box.begin() + index - npix);
878  }
879  {
880  const l1str &q_ref = l1qrec->r[i_e];
881  cldmsk::get_var_vals(var_box.data() + npix * i_e, q_ref, npix, nbands, ib, get_val);
882  cldmsk::get_var_mask(mask_box.data() + npix * i_e, q_ref, npix, nbands, ib, get_mask);
883  }
884  rearrange_stats();
885  }
886  };
887 }
888 #endif
const float SST3diff2
const float SSTvdiff
float tg_sen[NBANDS]
Definition: atrem_corl1.h:127
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
const get_valid rho_mask
const float Btmax40
const float Bt12unif2
std::vector< float > stdv
std::vector< float > var_min_box
Stats(size_t npix, size_t nscan, size_t rad_x, size_t rad_y, size_t center, l1qstr *l1qrec)
int32_t day
std::vector< int16_t > counts_out
const float dBtmax
bool(* get_valid)(const l1str &, int, int, int)
void get_ratio(l1str &l1str)
const float cldthresh
#define NBANDSIR
Definition: filehandle.h:23
std::vector< float > senz
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
const float SSTmaxn
const float dBt4min
const std::unordered_map< int, std::string > platforms
void get_window_1D_min(float *mins_1d, const float *inp1d, const int *mask1d, size_t npix, size_t radius)
the same as get_window_1D_max but for min
const std::unordered_map< std::string, const std::unordered_map< std::string, int > > bands_set
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in out
Definition: HISTORY.txt:422
std::vector< float > qual
float(* get_value)(const l1str &, int, int, int)
const float equatorialSouth
const float equatorialWest
read l1rec
#define VIIRSN
Definition: sensorDefs.h:23
PairOfSST(size_t npix, float *sst1, float *sst2, bool has_difference)
const std::unordered_set< int > viirs_sensors
void read_sst_bands(const std::string &sat, std::unordered_map< std::string, int > &bands)
StatsSST(size_t npix, size_t nscan, size_t rad_x, size_t rad_y, size_t center, l1qstr *l1qrec, float *sst, bool std_requested=false, bool max_req=false, bool min_req=false, bool min_max_req=false)
#define MODIST
Definition: sensorDefs.h:18
std::ostream & operator<<(std::ostream &stream, const std::vector< T > &data)
prints a vector
const float vhisenzv2
const float vhisenz
const float Bt11unif2
const float SST4diff2
std::vector< float > difference
const float equatorialEast
StatsVarBand(size_t npix, size_t nscan, size_t rad_x, size_t rad_y, size_t center, l1qstr *l1qrec, int ib, size_t nbands, cldmsk::get_valid get_mask, cldmsk::get_value get_val, bool std_requested=false, bool max_requested=false, bool min_requested=false, bool min_max_requested=false)
std::string get_ocdata_root()
Returns OCDATAROOT.
std::vector< float > var_max_box
std::vector< float > lat
@ string
std::vector< int16_t > counts
const get_value cirrus_value
Contains the data/stats for variables (BTs and SSTs) needed for cloud masks tests / ADT.
instr * input
#define M_PI
Definition: dtranbrdf.cpp:19
const float SST3diff1
#define BT_HI
Definition: l1.h:56
std::vector< float > bias_mean_out
const float Bt39unif2
std::string get_sensor(int key)
Get the sensor string.
const std::unordered_map< std::string, std::string(*)(const instr *)> get_sses_coeffs_path
returns the SSES LUTs paths to netCDF files depending on requested product
const float SSTdiff
void set_std()
Init and compute STD.
float * operator()(const std::string &key)
SSES data structure. Contains the LUTs and a look up routine. Reads NetCDF file LUTs file.
const get_valid cldrh_mask
A derived class to compute stats for an SST.
const float SSTmax
const float dBtmin
const double scan_time_modis_t_day2
const std::unordered_map< std::string, float > cldthresh_list
std::vector< float > bias_mean
Contains supportive functions and constants.
const float min_bt
const float Bt39unif1
void calculate_sses_bias_stats(const float *diff, const float *sst, const float *solz, const float *senz, size_t day, const float *lat, const int8_t *qual, const int16_t *flags, const float *glint)
Calculates the SSES bias vars.
std::vector< float > var_max
#define SSTF_BTNONUNIF
Definition: flags_sst.h:19
#define SSTF_HISENZ
Definition: flags_sst.h:24
const float el_corr_modis_t_1
void get_window_1D_max(float *maxs_1d, const float *inp1d, const int *mask1d, size_t npix, size_t radius)
Compute max values within a 1D window in one line in a queue.
const get_value cldrh_value
A derived class for a pair of BTs.
A derived class to compute stats for a BT.
float tg_sol[NBANDS]
Definition: atrem_corl1.h:127
std::vector< float > stdv_out
std::vector< float > var_std
const float invalid_val
const get_value rho_value
const float solznight
const float Btmin
std::vector< float > sst
const get_valid bt_mask
a context in which it is NOT documented to do so subscript which cannot be easily calculated when extracting TONS attitude data from the Terra L0 files Corrected several defects in extraction of entrained ephemeris and and as HDF file for both the L1A and Geolocation enabling retrieval of South Polar DEM data Resolved Bug by changing to opent the geolocation file only after a successful read of the L1A and also by checking for fatal errors from not restoring C5 and to report how many of those high resolution values were water in the new WaterPresent SDS Added valid_range attribute to Land SeaMask Changed to bilinearly interpolate the geoid_height to remove artifacts at one degree lines Made corrections to const qualification of pointers allowed by new version of M API library Removed casts that are no longer for same not the geoid Corrected off by one error in calculation of high resolution offsets Corrected parsing of maneuver list configuration parameter Corrected to set Height SDS to fill values when geolocation when for elevation and land water mask
Definition: HISTORY.txt:114
const std::unordered_map< std::string, size_t > bt_box_sizes
const float max_lt
void get_difference(l1str &l1str)
void get_total_2d_min(float *min_global, const float *inp2d, size_t npix, size_t nscan, size_t center, size_t radius)
the same as get_total_2d_max but for min
const float max_bt
A base class for a variable. Holds values for STD, MAX, MIN within a sliding windows.
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 pixel
Definition: HISTORY.txt:192
PairOfVarsBT(size_t npix, int ib1, int ib2, l1str &l1str, bool has_difference=false, bool has_ratio=false)
const float equatorialNorth
const float SSTmin
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
const float dBtrefmin
const float Bt11unif1
std::vector< float > var_min
std::vector< float > & month_data()
month variable for the desicion tree traversal
const float SST4diff1
SSESData(const std::string &nc_path_LUTs, size_t npix)
Construct a new SSESData object.
const get_valid cirrus_mask
#define BT
Definition: st_lt.h:10
def center
Definition: sort_gring.py:68
#define BAD_FLT
Definition: jplaeriallib.h:19
const float dBtrefmax
const std::unordered_map< std::string, std::string(*)(const instr *)> get_sst_coeffs_path
returns the sst coefs paths to netCDF files depending on requested product
flags
Definition: DDAlgorithm.h:22
void rearrange_stats()
when a l1 queue is updated, we need to updated the supporting arrays/queus as well
int32_t nbands
size_t operator()(size_t iqual, size_t ilat, size_t idiff, size_t isenz, size_t iquar, size_t iday, size_t isst)
const get_value bt_value
std::vector< float > var_minmax
void get_std_box(const float *inp2d, const int *mask2d, size_t npix, size_t nscan, size_t radius_x, size_t radius_y, float *out, size_t center, l1qstr *l1qrec)
Compute standard deviation in a window.
void get_total_2d_max(float *max_global, const float *inp2d, size_t npix, size_t nscan, size_t center, size_t radius)
Get the max value of center line with 2D sliding window using a 1D max queue computed with get_window...
std::vector< float > bias
void copy(double **aout, double **ain, int n)
std::vector< float > bias_out
const float Bt12unif1
const float Bt37unif2
const float Bt40unif1
float btrefdiffv6(int32_t ip, float BT39, float BT40, const l1str *l1rec, int fullscanpix)
Legacy BT diff test with interpolation adjutment.
const std::unordered_map< std::string, float *(*)(const l1str &)> get_non_BT_vars
const float Bt37unif1
#define BT_LO
Definition: l1.h:55
void set_maxmin()
Init and compute difference between MAX and MIN.
#define bool
Definition: usrmac.h:88
#define VIIRSJ2
Definition: sensorDefs.h:44
const double scan_time_modis_t_day1
A base class for a pair of vars. Holds the variables difference and ratio.
const float Btmax
void set_min()
Init and compute MIN.
void get_var_mask(int *mask, const l1str &l1str, size_t npix, int nbands, int ib, get_valid get_mask)
Get the mask (line) of a BT.
std::vector< int > mask_box
const float min_lt
const float Bt40unif2
void radius(double A)
Definition: proj_report.c:132
std::vector< float > ratio
#define abs(a)
Definition: misc.h:90
int i
Definition: decode_rs.h:71
const std::unordered_map< std::string, std::unordered_set< int > > all_sensors
const float dBt4max
void get_var_vals(float *BT, const l1str &l1str, size_t npix, int nbands, int ib, get_value get_val)
the same as get_var_mask but gets values
const float el_corr_modis_t_2
#define MODISA
Definition: sensorDefs.h:19
#define VIIRSJ1
Definition: sensorDefs.h:37
A derived class for a pair of SSTs.
el
Definition: decode_rs.h:168
int npix
Definition: get_cmp.c:28
const float glintmax
std::vector< float > var_box
const std::unordered_set< int > modis_sensors
std::vector< float > diff
std::unordered_map< std::string, float * > return_vals
void set_max()
Init and compute MAX.
l2prod max
const float hisenz