NASA Logo
Ocean Color Science Software

ocssw V2022
sst_cloud_mask.cpp
Go to the documentation of this file.
1 #include "sst_adt.hpp"
3 #include "sst.h"
4 #include "timeutils.h"
9 namespace {
10  const double CtoK = 273.15e0; // Celcius to Kelvin
19  template<class T>
20  std::ostream &operator<<(std::ostream &stream, const std::vector<T> &data) {
21  for (const auto &el: data) {
22  stream << el << " ";
23  }
24  stream << "\n";
25  return stream;
26  }
27 
37  std::vector<std::string> parse_test_name(const std::string &test_name) {
38  std::vector<std::string> out;
39  boost::split(out, test_name, boost::is_any_of("_"));
40  return out;
41  }
42 
52  template<class T>
53  void binary_search(const std::vector<T> &data, const T &val, size_t &start,
54  size_t &end) {
55  while (end - start > 1) {
56  size_t middle = (end + start) / 2;
57  if (data.at(middle) >= val) {
58  end = middle;
59  binary_search(data, val, start, end);
60  } else {
61  start = middle;
62  binary_search(data, val, start, end);
63  }
64  }
65  }
66 
67 // bands
68  std::unordered_map<std::string, int> ibands;
69  size_t n_bands;
70 // variables needed for quality flags setting
71  const boost::unordered_map<std::pair<std::string, std::string>,
72  std::set<std::string>>
73  needed_stats_for_cloud_mask = {
74  {{"viirs", "SST3"},
75  {"BT37_MAXMIN", "BT12_MAXMIN", "BT11_MAXMIN", "BT40", "d_BT37_BT12",
76  "d_BT11_BT12"}},
77  {{"viirs", "SST"},
78  {"BT12_MAXMIN", "BT11_MAXMIN", "BT40", "CLDRHred_MAXMIN"}},
79  {{"modis", "SST"},
80  {"BT12_MAXMIN", "BT11_MAXMIN", "BT40_MAXMIN", "BT39_MAXMIN",
81  "CLDRHred_MAXMIN"}},
82  {{"modis", "SST4"},
83  {"BT12_MAXMIN", "BT11_MAXMIN", "BT40_MAXMIN", "BT39_MAXMIN",
84  "d_BT39_BT40"}}};
85 
86 } // namespace
87 
130 private:
131  std::string this_product;
132  std::vector<float> sst;
133  std::vector<int> valid_mask;
134  std::vector<int16_t> flags_sst;
135  std::vector<int8_t> qual_sst;
136  std::vector<float> treesum;
137  double scan_time;
138  int16_t day;
139  size_t npix;
140  size_t n_q_size;
141  size_t fullscanpix = 1354;
142  const float latwin = 2.5f;
143  size_t bt_box = 5; // default is 5
144  size_t i_center, i_s, i_e;
145  int32_t processed_scan_sst = -1;
146  int32_t processed_scan_flags = -1;
147  int32_t processed_scan_sses = -1;
148  std::string sensor;
149  std::string platform;
150  std::vector<float> sst_regression_coefficients;
151  std::vector<float> bounds_lat, bounds_search;
152  std::unique_ptr<cldmsk::SSESData> sses_bias_data;
153  // regression model name
154  std::string regression_model;
155  size_t nlatbands, ncoeffs;
156  int ib11, ib12, ib37, ib40, ib39, ib85;
157  std::unordered_map<std::string, SeaSurfaceTemperatureCalcuations *>
158  external_products;
159  adt::Treenode *desicion_tree = nullptr;
160  std::set<std::string> test_names;
161  std::unordered_map<std::string, std::vector<std::string>>
162  test_classifications;
163  std::unordered_map<std::string, float *> test_data;
164  std::unordered_map<std::string, int *> test_masks;
165  // maps
166  boost::unordered_map<std::string, bstats::StatsVarBand> StatsBTs;
167  boost::unordered_map<std::pair<std::string, std::string>,
169  PairOfVarsBTs;
170  boost::unordered_map<std::string, bstats::StatsSST> StatsSSTs;
171  boost::unordered_map<std::pair<std::string, std::string>, bstats::PairOfSST>
172  PairOfSSTs;
173  std::string path_to_cloud_mask_config;
174  // modis reference
175  std::vector<float> modis_ref;
176  // modis dust correction
177  std::vector<float> dust_correction_data;
178 
179  // modis A BT40 data for detector #0
180  std::vector<float> modis_a_bt40;
181  l1qstr *l1qstr_modis_a = nullptr;
182  // modis T electrionic correction
183  float elecor_terra = 0.0f;
184  bool is_aqua = false;
185  bool is_terra = false;
186 
194  void set_modis_a_bt_40_data(l1str *l1rec, l1qstr *l1qrec, instr *input,
195  const std::string &product_type) {
196  if (!is_aqua) return;
197  std::vector<float> &var_box = StatsBTs.at("40").var_box;
198  const size_t center = StatsBTs.at("40").center;
199  if (l1rec->detnum == 0 && l1qrec->r[i_center].detnum == 0 && i_center > 0) {
200  modis_a_bt40 = std::vector<float>(npix, BAD_FLT);
201  const size_t index_prev = (i_center - 1) * npix;
202  const size_t index_next = (i_center + 1) * npix;
203  for (size_t i_p = 0; i_p < npix; i_p++) {
204  float sum = 0.0f;
205  int count = 0;
206  const float prev_val = var_box.at(index_prev + i_p);
207  const float next_val = var_box.at(index_next + i_p);
208  const int mask_prev = cldmsk::bt_mask(l1qstr_modis_a->r[i_center - 1], i_p, NBANDSIR, ib40);
209  const int mask_next = cldmsk::bt_mask(l1qstr_modis_a->r[i_center + 1], i_p, NBANDSIR, ib40);
210  if (mask_prev) {
211  count++;
212  sum += prev_val;
213  }
214  if (mask_next) {
215  count++;
216  sum += next_val;
217  }
218  if (count > 0) {
219  modis_a_bt40.at(i_p) = sum / count;
220  } else {
221  const int mask_cur = cldmsk::bt_mask(l1qstr_modis_a->r[i_center], i_p, NBANDSIR, ib40);
222  if (mask_cur) {
223  modis_a_bt40.at(i_p) =
224  cldmsk::bt_mask(l1qstr_modis_a->r[i_center], i_p, NBANDSIR, ib40);
225  }
226  }
227  }
228  StatsBTs.at("40").return_vals.at("VAR") = modis_a_bt40.data();
229  for (auto &pair_bt : PairOfVarsBTs) {
230  if (pair_bt.first.first == "40") {
231  const auto ib_pair = ibands.at("ib" + pair_bt.first.second);
232  for (size_t i_p = 0; i_p < npix; i_p++) {
233  pair_bt.second.difference.at(i_p) =
234  modis_a_bt40.at(i_p) - cldmsk::bt_value(*l1rec, i_p, NBANDSIR, ib_pair);
235  }
236  }
237  if (pair_bt.first.second == "40") {
238  const auto ib_pair = ibands.at("ib" + pair_bt.first.first);
239  for (size_t i_p = 0; i_p < npix; i_p++) {
240  pair_bt.second.difference.at(i_p) =
241  cldmsk::bt_value(*l1rec, i_p, NBANDSIR, ib_pair) - modis_a_bt40.at(i_p);
242  }
243  }
244  }
245  } else {
246  StatsBTs.at("40").return_vals.at("VAR") =
247  var_box.data() + npix * center;
248  }
249  };
250 
251  void modis_dust_correction(const l1str &q_ref, size_t i_q) {
252  // get dust correction
253  if (!dust_correction_data.empty() && this_product != "SST4") {
254  const int sensor_id = q_ref.l1file->sensorID;
255  for (size_t i_p = 0; i_p < npix; i_p++)
256  if (q_ref.solz[i_p] >= cldmsk::solznight) {
257  const float dustExtinction =
258  q_ref.anc_aerosol->dust_ext[i_p];
259  float csenz = q_ref.csenz[i_p];
260  const float Bt39 =
261  cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib39);
262  const float Bt85 =
263  cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib85);
264  const float Bt11 =
265  cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib11);
266  const float Bt12 =
267  cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib12);
268  dust_correction_data.at(i_p) =
269  dust_correction(dustExtinction, csenz, Bt39, Bt85, Bt11,
270  Bt12, sensor_id);
271  sst.at(i_p + i_q * npix) += dust_correction_data.at(i_p);
272  }
273  }
274  }
275 
284  float *get_reference(l1str *l1rec, l1qstr *l1qrec, instr *input,
285  const std::string &product_type, size_t is) //
286  {
287  if (sensor != "modis" || this_product == "SST4")
288  return l1qrec->r[is].sstref;
289  else {
290  if (modis_ref.empty())
291  modis_ref = std::vector<float>(npix, BAD_FLT);
292  for (size_t i_p = 0; i_p < npix; i_p++) {
293  const float solz = l1qrec->r[is].solz[i_p];
294  if (solz >= cldmsk::solznight) {
295  modis_ref.at(i_p) = external_products.at("SST4")->get_sst(
296  l1rec, l1qrec, input,
297  "SST4")[i_p + (is - i_center) * npix];
298  } else {
299  modis_ref.at(i_p) = l1qrec->r[is].sstref[i_p];
300  }
301  }
302  return modis_ref.data();
303  }
304  }
305 
312  void get_valid_mask_sst(const l1str &q_ref, int *mask,
313  const float *sst_ref) {
314  for (size_t i_p = 0; i_p < npix; i_p++) {
315  mask[i_p] = cldmsk::bt_mask(q_ref, i_p, NBANDSIR, ib11) &&
316  cldmsk::bt_mask(q_ref, i_p, NBANDSIR, ib12) &&
317  ((q_ref.flags[i_p] & LAND) == 0) &&
318  ((q_ref.flags[i_p] & NAVFAIL) == 0) &&
319  (sst_ref[i_p] >= BAD_FLT + 1.0f);
320  }
321  }
322 
326  void get_valid_mask_sst3(const l1str &q_ref, int *mask,
327  const float *sst_ref) {
328  for (size_t i_p = 0; i_p < npix; i_p++) {
329  mask[i_p] = cldmsk::bt_mask(q_ref, i_p, NBANDSIR, ib37) &&
330  cldmsk::bt_mask(q_ref, i_p, NBANDSIR, ib11) &&
331  cldmsk::bt_mask(q_ref, i_p, NBANDSIR, ib12) &&
332  ((q_ref.flags[i_p] & LAND) == 0) &&
333  ((q_ref.flags[i_p] & NAVFAIL) == 0) &&
334  (sst_ref[i_p] >= BAD_FLT + 1.0f);
335  }
336  }
337 
341  void get_valid_mask_sst4(const l1str &q_ref, int *mask,
342  const float *sst_ref) {
343  for (size_t i_p = 0; i_p < npix; i_p++) {
344  mask[i_p] = cldmsk::bt_mask(q_ref, i_p, NBANDSIR, ib39) &&
345  cldmsk::bt_mask(q_ref, i_p, NBANDSIR, ib40) &&
346  ((q_ref.flags[i_p] & LAND) == 0) &&
347  ((q_ref.flags[i_p] & NAVFAIL) == 0) &&
348  (sst_ref[i_p] >= BAD_FLT + 1.0f);
349  }
350  }
351 
359  size_t find_lat_band_interpolation_point(const l1str &q_ref, size_t i_p) {
360  size_t ic;
361  {
362  const float lat = q_ref.lat[i_p];
363  size_t start = 0;
364  size_t end = bounds_search.size() - 1;
365  if (lat - latwin > bounds_search[0])
366  binary_search(bounds_search, lat - latwin, start, end);
367  ic = start;
368  }
369  return ic;
370  }
371 
375  void get_nlsst_v6_cmc(const l1str &q_ref, float *sst_ptr, const int *mask,
376  const float *ref) {
377  for (size_t i_p = 0; i_p < npix; i_p++) {
378  sst_ptr[i_p] = BAD_FLT;
379  if (mask[i_p]) {
380  const size_t interp_pnt =
381  find_lat_band_interpolation_point(q_ref, i_p);
382  const float lat = q_ref.lat[i_p];
383  const float mu = q_ref.csenz[i_p];
384  const float satzen = q_ref.senz[i_p];
385  const float Bt11 =
386  cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib11) + CtoK;
387  const float dBT =
388  Bt11 - cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib12) - CtoK;
389  const float cmc = std::max(ref[i_p], 0.0f);
390  const float satzdir =
391  (q_ref.pixnum[i_p] < (int32_t)fullscanpix / 2) ? -1.0 : 1.0;
392  float *coeffptr =
393  sst_regression_coefficients.data() + ncoeffs * interp_pnt;
394  const float lsst = coeffptr[0] + coeffptr[1] * Bt11 +
395  coeffptr[2] * dBT * cmc +
396  coeffptr[3] * dBT * ((1.0 / mu) - 1.0) +
397  coeffptr[4] * satzen * satzdir +
398  coeffptr[5] * satzen * satzen;
399  if (lat < bounds_search[interp_pnt + 1] - latwin ||
400  interp_pnt == (nlatbands - 1)) {
401  sst_ptr[i_p] = lsst - CtoK;
402  } else {
403  coeffptr += ncoeffs; // ncoeffs must be 6 for viirs
404  const float dBtlo = bounds_search[interp_pnt + 1] - latwin;
405  const float hsst = coeffptr[0] + coeffptr[1] * Bt11 +
406  coeffptr[2] * dBT * cmc +
407  coeffptr[3] * dBT * ((1.0 / mu) - 1.0) +
408  coeffptr[4] * satzen * satzdir +
409  coeffptr[5] * satzen * satzen;
410  sst_ptr[i_p] =
411  lsst + ((lat - dBtlo) / latwin / 2.0) * (hsst - lsst) -
412  CtoK;
413  }
414  }
415  }
416  }
417 
421  void get_nlsst_v7_cmc(const l1str &q_ref, float *sst_ptr, const int *mask,
422  const float *ref) {
423  const int mside = q_ref.mside;
424  for (size_t i_p = 0; i_p < npix; i_p++) {
425  sst_ptr[i_p] = BAD_FLT;
426  if (mask[i_p]) {
427  const size_t interp_pnt =
428  find_lat_band_interpolation_point(q_ref, i_p);
429  const float lat = q_ref.lat[i_p];
430  const float mu = q_ref.csenz[i_p];
431  const float satzen = q_ref.senz[i_p];
432  const float Bt11 = cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib11);
433  const float dBT =
434  Bt11 - cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib12);
435  //const float solz = q_ref.solz[i_p];
436  const float cmc = std::max(ref[i_p], 0.0f);
437  const float satzdir =
438  (q_ref.pixnum[i_p] < (int32_t)fullscanpix / 2) ? -1.0 : 1.0;
439  float *coeffptr =
440  sst_regression_coefficients.data() + ncoeffs * interp_pnt;
441  const float lsst =
442  coeffptr[0] + coeffptr[1] * Bt11 + coeffptr[2] * dBT * cmc +
443  coeffptr[3] * dBT * ((1.0 / mu) - 1.0) +
444  coeffptr[4] * mside + coeffptr[5] * satzen * satzdir +
445  coeffptr[6] * satzen * satzen;
446  if (lat < bounds_search[interp_pnt + 1] - latwin ||
447  interp_pnt == (nlatbands - 1)) {
448  sst_ptr[i_p] = lsst;
449  } else {
450  coeffptr += ncoeffs; // ncoeffs must be 7 for modis
451  const float dBtlo = bounds_search[interp_pnt + 1] - latwin;
452  const float hsst = coeffptr[0] + coeffptr[1] * Bt11 +
453  coeffptr[2] * dBT * cmc +
454  coeffptr[3] * dBT * ((1.0 / mu) - 1.0) +
455  coeffptr[4] * mside +
456  coeffptr[5] * satzen * satzdir +
457  coeffptr[6] * satzen * satzen;
458  sst_ptr[i_p] =
459  lsst + ((lat - dBtlo) / latwin / 2.0) * (hsst - lsst);
460  }
461  }
462  }
463  }
467  void get_mcsst3_v7(const l1str &q_ref, float *sst_ptr, const int *mask, const float *ref) {
468  const int mside = q_ref.mside;
469  for (size_t i_p = 0; i_p < npix; i_p++) {
470  sst_ptr[i_p] = BAD_FLT;
471  if (mask[i_p]) {
472  const size_t interp_pnt = find_lat_band_interpolation_point(q_ref, i_p);
473  const float lat = q_ref.lat[i_p];
474  const float mu = q_ref.csenz[i_p];
475  const float satzen = q_ref.senz[i_p];
476  const float Bt37 = cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib37);
477  const float Bt11 = cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib11);
478  const float dBT = Bt37 - cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib12);
479  float *coeffptr = sst_regression_coefficients.data() + ncoeffs * interp_pnt;
480  const float satzdir =
481  (q_ref.pixnum[i_p] < (int32_t)fullscanpix / 2) ? -1.0 : 1.0;
482  const float lsst = coeffptr[0] + coeffptr[1] * Bt11 + coeffptr[2] * dBT +
483  coeffptr[3] * ((1.0 / mu) - 1.0) + coeffptr[5] * satzen * satzdir +
484  coeffptr[6] * satzen * satzen + coeffptr[4] * mside;
485  if (lat < bounds_search[interp_pnt + 1] - latwin || interp_pnt == (nlatbands - 1)) {
486  sst_ptr[i_p] = lsst;
487  } else {
488  coeffptr += ncoeffs; // ncoeffs must be 7 for modis
489  const float dBtlo = bounds_search[interp_pnt + 1] - latwin;
490  const float hsst = coeffptr[0] + coeffptr[1] * Bt11 + coeffptr[2] * dBT +
491  coeffptr[3] * ((1.0 / mu) - 1.0) + coeffptr[5] * satzen * satzdir+
492  coeffptr[6] * satzen * satzen + coeffptr[4] * mside;
493  sst_ptr[i_p] = lsst + ((lat - dBtlo) / latwin / 2.0) * (hsst - lsst);
494  }
495  }
496  }
497  }
501  void get_nlsst_modis4(const l1str &q_ref, float *sst_ptr, const int *mask,
502  const float *ref) {
503  const int mside = q_ref.mside;
504  const int iscan = q_ref.iscan;
505  int center = i_center;
506  if (iscan > 0) center = i_e;
507  bool modis_a_correction = false;
508  if (is_aqua)
509  modis_a_correction =
510  (l1qstr_modis_a->r[center].detnum == 0) && q_ref.detnum == 0 && center > 0;
511  for (size_t i_p = 0; i_p < npix; i_p++) {
512  sst_ptr[i_p] = BAD_FLT;
513  if (mask[i_p]) {
514  const size_t interp_pnt =
515  find_lat_band_interpolation_point(q_ref, i_p);
516  const float lat = q_ref.lat[i_p];
517  const float mu = q_ref.csenz[i_p];
518  const float satzen = q_ref.senz[i_p];
519  const float Bt39 = cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib39);
520  float Bt40 = cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib40);
521  if (modis_a_correction) {
522  float sum = 0.0f;
523  int count = 0;
524  const float Bt40_prev = cldmsk::bt_value(
525  l1qstr_modis_a->r[center - 1], i_p, NBANDSIR, ib40);
526  const float Bt40_next = cldmsk::bt_value(
527  l1qstr_modis_a->r[center + 1], i_p, NBANDSIR, ib40);
528  const int mask_prev = cldmsk::bt_mask(
529  l1qstr_modis_a->r[center - 1], i_p, NBANDSIR, ib40);
530  const int mask_next = cldmsk::bt_mask(
531  l1qstr_modis_a->r[center + 1], i_p, NBANDSIR, ib40);
532  if (mask_prev) {
533  sum += Bt40_prev;
534  count++;
535  }
536  if (mask_next) {
537  sum += Bt40_next;
538  count++;
539  }
540  if (count > 0) {
541  Bt40 = sum / count;
542  } else
543  continue;
544  }
545 
546  const float dBT = Bt39 - Bt40;
547  const float satzdir =
548  (q_ref.pixnum[i_p] < (int32_t)fullscanpix / 2) ? -1.0 : 1.0;
549  float *coeffptr =
550  sst_regression_coefficients.data() + ncoeffs * interp_pnt;
551  const float lsst =
552  coeffptr[0] + coeffptr[1] * Bt39 + coeffptr[2] * dBT +
553  coeffptr[3] * ((1.0 / mu) - 1.0) + coeffptr[4] * mside +
554  coeffptr[5] * satzen * satzdir +
555  coeffptr[6] * satzen * satzen;
556  if (lat < bounds_search[interp_pnt + 1] - latwin ||
557  interp_pnt == (nlatbands - 1)) {
558  sst_ptr[i_p] = lsst;
559  } else {
560  coeffptr += ncoeffs; // ncoeffs must be 7 for modis
561  const float dBtlo = bounds_search[interp_pnt + 1] - latwin;
562  const float hsst =
563  coeffptr[0] + coeffptr[1] * Bt39 + coeffptr[2] * dBT +
564  coeffptr[3] * ((1.0 / mu) - 1.0) + coeffptr[4] * mside +
565  coeffptr[5] * satzen * satzdir +
566  coeffptr[6] * satzen * satzen;
567  sst_ptr[i_p] =
568  lsst + ((lat - dBtlo) / latwin / 2.0) * (hsst - lsst);
569  }
570  }
571  if (is_terra) {
572  sst_ptr[i_p] += elecor_terra;
573  }
574  }
575  }
576 
580  void get_nlsst_viirs3(const l1str &q_ref, float *sst_ptr, const int *mask,
581  const float *ref) {
582  for (size_t i_p = 0; i_p < npix; i_p++) {
583  sst_ptr[i_p] = BAD_FLT;
584  if (mask[i_p]) {
585  const size_t interp_pnt =
586  find_lat_band_interpolation_point(q_ref, i_p);
587  const float lat = q_ref.lat[i_p];
588  const float mu = q_ref.csenz[i_p];
589  const float satzen = q_ref.senz[i_p];
590  const float Bt11 =
591  cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib11) + CtoK;
592  const float dBT = cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib37) -
593  cldmsk::bt_value(q_ref, i_p, NBANDSIR, ib12);
594  const float cmc = std::max(ref[i_p], 0.0f);
595  const float satzdir =
596  (q_ref.pixnum[i_p] < (int32_t)fullscanpix / 2) ? -1.0 : 1.0;
597  float *coeffptr =
598  sst_regression_coefficients.data() + ncoeffs * interp_pnt;
599  const float lsst = coeffptr[0] + coeffptr[1] * Bt11 +
600  coeffptr[2] * dBT * cmc +
601  coeffptr[3] * ((1.0 / mu) - 1.0) +
602  coeffptr[4] * satzen * satzdir +
603  coeffptr[5] * satzen * satzen;
604  if (lat < bounds_search[interp_pnt + 1] - latwin ||
605  interp_pnt == (nlatbands - 1)) {
606  sst_ptr[i_p] = lsst - CtoK;
607  } else {
608  coeffptr += ncoeffs; // ncoeffs must be 6 for viirs
609  const float dBtlo = bounds_search[interp_pnt + 1] - latwin;
610  const float hsst = coeffptr[0] + coeffptr[1] * Bt11 +
611  coeffptr[2] * dBT * cmc +
612  coeffptr[3] * ((1.0 / mu) - 1.0) +
613  coeffptr[4] * satzen * satzdir +
614  coeffptr[5] * satzen * satzen;
615  sst_ptr[i_p] =
616  lsst + ((lat - dBtlo) / latwin / 2.0) * (hsst - lsst) -
617  CtoK;
618  }
619  }
620  }
621  }
622 
628  void cloud_mask(const l1str &q_ref, instr *input) {
629  cldmsk::product_types case_switcher;
630  if (this_product == "SST")
631  case_switcher = cldmsk::SST;
632  else if (this_product == "SST3")
633  case_switcher = cldmsk::SST3;
634  else
635  case_switcher = cldmsk::SST4;
636  // only VIIRS and ONLY SST
637  const size_t day2 = day;
638  const float xdoy = day2 + 284.0;
639  const float xrad = (360.0 / 365.0) * xdoy / RADEG;
640  const float subsolar = 23.45 * std::sin(xrad);
641  const bool is_viirs =
642  sensor ==
643  "viirs"; // technically cthe compiler should optimize it out
644  const bool is_modis = sensor == "modis";
645  for (size_t i_p = 0; i_p < npix; i_p++) {
646  if (((q_ref.flags[i_p] & LAND) > 0) ||
647  ((q_ref.flags[i_p] & NAVFAIL) > 0)) {
648  flags_sst[i_p] |= SSTF_ISMASKED;
649  qual_sst[i_p] = 4;
650  continue;
651  }
652  // bt bad
653  {
654  bool good_bt = true;
655  switch (case_switcher) {
656  case cldmsk::SST3:
657  good_bt &= StatsBTs.at("37").get_valid_mask()[i_p];
658  case cldmsk::SST:
659  good_bt &= (StatsBTs.at("11").get_valid_mask()[i_p] &&
660  StatsBTs.at("12").get_valid_mask()[i_p]);
661  break;
662  case cldmsk::SST4:
663  good_bt &= (StatsBTs.at("39").get_valid_mask()[i_p] &&
664  StatsBTs.at("40").get_valid_mask()[i_p]);
665  break;
666  default:
667  break;
668  }
669  if (!good_bt) {
670  flags_sst[i_p] |= SSTF_BTBAD;
671  qual_sst[i_p] = 4;
672  continue;
673  }
674  }
675  bool night = q_ref.solz[i_p] >= cldmsk::solznight;
676  bool glint = q_ref.glint_coef[i_p] <= cldmsk::glintmax && !night;
677  {
678  if (!night) {
679  switch (case_switcher) {
680  case cldmsk::SST4:
681  case cldmsk::SST3:
682  qual_sst[i_p] = 3;
683  break;
684  case cldmsk::SST:
685  if (q_ref.glint_coef[i_p] > cldmsk::glintmax)
686  qual_sst[i_p] = 1;
687  default:
688  break;
689  }
690  }
691  }
692  // bt range
693  {
694  bool bt_range = true;
695  bool bt_range_max_sst4 = true;
696  bool bt_range_min_sst4 = true;
697  bool bt_range_max_sst3 = true;
698  bool bt_range_min_sst3 = true;
699  bool bt_range_sst = true;
700  if (case_switcher != cldmsk::SST4) {
701  float bt11 = StatsBTs.at("11")("VAR")[i_p];
702  float bt12 = StatsBTs.at("12")("VAR")[i_p];
703  bt_range_sst &=
704  (bt11 <= cldmsk::Btmax) && (bt12 <= cldmsk::Btmax) &&
705  (bt11 >= cldmsk::Btmin) && (bt12 >= cldmsk::Btmin);
706  }
707  bt_range &= bt_range_sst;
708  if (is_modis) {
709  float bt39 = StatsBTs.at("39")("VAR")[i_p];
710  float bt40 = StatsBTs.at("40")("VAR")[i_p];
711  bt_range_min_sst4 &=
712  (bt39 >= cldmsk::Btmin) && (bt40 >= cldmsk::Btmin);
713  bt_range_max_sst4 &=
714  (bt39 <= cldmsk::Btmax) && (bt40 <= cldmsk::Btmax40);
715  }
716 
717  if (is_viirs) {
718  float bt37 = StatsBTs.at("37")("VAR")[i_p];
719  float bt40 = StatsBTs.at("40")("VAR")[i_p];
720  bt_range_min_sst3 &=
721  (bt37 >= cldmsk::Btmin) && (bt40 >= cldmsk::Btmin);
722  bt_range_max_sst3 &=
723  (bt37 <= cldmsk::Btmax) && (bt40 <= cldmsk::Btmax40);
724  }
725  switch (case_switcher) {
726  case cldmsk::SST4:
727  bt_range &= bt_range_min_sst4 && bt_range_max_sst4;
728  break;
729  case cldmsk::SST3:
730  bt_range &= bt_range_min_sst3 && bt_range_max_sst3;
731  break;
732  case cldmsk::SST:
733  if (night) {
734  bt_range &= bt_range_min_sst3 && bt_range_max_sst3;
735  } else if (glint) {
736  bt_range &= bt_range_min_sst3;
737  }
738  if (is_modis) {
739  if (night) {
740  bt_range &=
741  bt_range_min_sst4 && bt_range_max_sst4;
742  } else if (glint) {
743  bt_range &= bt_range_min_sst4;
744  }
745  }
746  break;
747  default:
748  break;
749  }
750  if (!bt_range) {
751  flags_sst[i_p] |= SSTF_BTRANGE;
752  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
753  }
754  }
755  // BT diff test for SST/SST3
756  {
757  bool bt_diff = true;
758  switch (case_switcher) {
759  case cldmsk::SST3:
760  case cldmsk::SST: {
761  float diffBT = PairOfVarsBTs.at({"11", "12"})()[i_p];
762  if (diffBT < cldmsk::dBtmin || diffBT > cldmsk::dBtmax)
763  bt_diff = false;
764  }
765  break;
766  case cldmsk::SST4: {
767  float diffBT = PairOfVarsBTs.at({"39", "40"})()[i_p];
768  if (diffBT < cldmsk::dBt4min ||
769  diffBT > cldmsk::dBt4max)
770  bt_diff = false;
771  }
772  break;
773  default:
774  break;
775  }
776  if (!bt_diff) {
777  flags_sst[i_p] |= SSTF_BTDIFF;
778  }
779  }
780  // BT diff for SST4/SST
781  {
782  switch (case_switcher) {
783  case cldmsk::SST3:
784  break;
785  case cldmsk::SST:
786  case cldmsk::SST4:
787  if (is_modis) {
788  const float BT39 = StatsBTs.at("39")("VAR")[i_p];
789  const float BT40 = StatsBTs.at("40")("VAR")[i_p];
790  const float d3940ref = cldmsk::btrefdiffv6(
791  i_p, BT39, BT40, &q_ref, fullscanpix);
792  if ((d3940ref < cldmsk::dBtrefmin ||
793  d3940ref > cldmsk::dBtrefmax) &&
794  night) {
795  flags_sst[i_p] |= SSTF_BT4REFDIFF;
796  qual_sst[i_p] =
797  std::max(qual_sst[i_p], (int8_t) 3);
798  }
799  }
800  break;
801  default:
802  break;
803  }
804  }
805  // sst range
806  {
807  bool sst_range = true;
808  float sst_val = sst.at(i_center * npix + i_p);
809  bool sst_min = sst_val >= cldmsk::SSTmin;
810  bool sst_max_night = sst_val <= cldmsk::SSTmaxn;
811  bool sst_max_day = sst_val <= cldmsk::SSTmax;
812 
813  switch (case_switcher) {
814  case cldmsk::SST4:
815  case cldmsk::SST3:
816  sst_range &= sst_min && sst_max_night;
817  break;
818  case cldmsk::SST:
819  sst_range &= sst_min && ((sst_max_night && night) ||
820  (sst_max_day && !night));
821  break;
822  default:
823  break;
824  }
825  if (!sst_range) {
826  flags_sst[i_p] |= SSTF_SSTRANGE;
827  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
828  }
829  }
830  // sst diff continue
831  {
832  const float dt_analysis =
833  sst.at(i_center * npix + i_p) - q_ref.sstref[i_p];
834  if ((l1_input->evalmask & SSTMODS) == 0) {
835  float dsst = cldmsk::SSTdiff;
836  dsst = input->sstrefdif;
837  if ((dt_analysis < -dsst) &&
838  q_ref.lat[i_p] >= cldmsk::equatorialSouth &&
839  q_ref.lat[i_p] <= cldmsk::equatorialNorth &&
840  q_ref.lon[i_p] >= cldmsk::equatorialWest &&
841  q_ref.lon[i_p] <= cldmsk::equatorialEast) {
842  flags_sst[i_p] |= SSTF_SSTREFDIFF;
843  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 2);
844  }
845  };
846  if (night) {
847  if (std::abs(dt_analysis) > cldmsk::SSTdiff) {
848  flags_sst[i_p] |= SSTF_SSTREFDIFF;
849  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 2);
850  }
851  if (std::abs(dt_analysis) > cldmsk::SSTvdiff) {
852  flags_sst[i_p] |= SSTF_SSTREFVDIFF;
853  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
854  }
855  } else {
856  if (dt_analysis < -cldmsk::SSTdiff ||
857  q_ref.sstref[i_p] < cldmsk::invalid_val) {
858  flags_sst[i_p] |= SSTF_SSTREFDIFF;
859  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 2);
860  }
861  if (dt_analysis < -cldmsk::SSTvdiff ||
862  q_ref.sstref[i_p] < cldmsk::invalid_val) {
863  flags_sst[i_p] |= SSTF_SSTREFVDIFF;
864  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
865  }
866  }
867  }
868  // difference between sst
869  {
870  float dsst, lim1, lim2;
871  bool flag1(false), flag2(false);
872  if (is_viirs) {
873  dsst = PairOfSSTs.at({"SST", "SST3"})()[i_p];
874  lim1 = cldmsk::SST3diff1;
875  lim2 = cldmsk::SST3diff2;
876  flag1 = std::abs(dsst) > lim1;
877  flag2 = std::abs(dsst) > lim2;
878  }
879  if (is_modis) {
880  dsst = PairOfSSTs.at({"SST", "SST4"})()[i_p];
881  // this is actually a bug
882  // dsst = std::abs(dsst);
883  lim1 = cldmsk::SST4diff1;
884  lim2 = cldmsk::SST4diff2;
885  flag1 = dsst < lim1;
886  flag2 = dsst < lim2;
887  }
888  switch (case_switcher) {
889  case cldmsk::SST4:
890  if (flag1) flags_sst[i_p] |= SSTF_SST4DIFF;
891  if (flag2) flags_sst[i_p] |= SSTF_SST4VDIFF;
892  break;
893  case cldmsk::SST3:
894  if (flag1) flags_sst[i_p] |= SSTF_SST3DIFF;
895  if (flag2) flags_sst[i_p] |= SSTF_SST3VDIFF;
896  break;
897  case cldmsk::SST:
898  if (night) {
899  if (is_modis) {
900  if (flag1) flags_sst[i_p] |= SSTF_SST4DIFF;
901  if (flag2) flags_sst[i_p] |= SSTF_SST4VDIFF;
902  }
903  if (is_viirs) {
904  if (flag1) flags_sst[i_p] |= SSTF_SST3DIFF;
905  if (flag2) flags_sst[i_p] |= SSTF_SST3VDIFF;
906  }
907  }
908  break;
909  default:
910  break;
911  }
912  if (night) {
913  if (flag1) {
914  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 2);
915  }
916  if (flag2) {
917  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 1);
918  }
919  }
920  }
921  // some vza tests
922  {
923  if (q_ref.senz[i_p] > cldmsk::hisenz) {
924  flags_sst[i_p] |= SSTF_HISENZ;
925  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 1);
926  }
927  if (q_ref.pixnum[i_p] < 2 ||
928  q_ref.pixnum[i_p] > ((int32_t)fullscanpix - 3))
929  flags_sst[i_p] |= SSTF_VHISENZ;
930  float vhisenz = cldmsk::vhisenz;
931  if (is_viirs) {
933  }
934  if (is_terra && is_modis) {
935  if (q_ref.pixnum[i_p] > 1349)
936  flags_sst[i_p] |= SSTF_VHISENZ;
937  }
938  switch (case_switcher) {
939  case cldmsk::SST4:
940  case cldmsk::SST3:
941  if (q_ref.senz[i_p] > cldmsk::vhisenz) {
942  flags_sst[i_p] |= SSTF_VHISENZ;
943  }
944  break;
945  case cldmsk::SST:
946  if (q_ref.senz[i_p] > vhisenz) {
947  flags_sst[i_p] |= SSTF_VHISENZ;
948  }
949  break;
950  default:
951  break;
952  }
953  {
954  if ((flags_sst[i_p] & SSTF_VHISENZ) > 0) {
955  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
956  }
957  }
958  }
959  {
960  // normalize rho test
961  if (case_switcher == cldmsk::SST) {
962  {
963  const float cldthreshold =
964  cldmsk::cldthresh_list.at(sensor);
965  const float val_to_threshold =
966  StatsBTs.at("RHred")("MAXMIN")[i_p];
967  if (glint) {
968  const float dt_analysis =
969  sst.at(i_center * npix + i_p) -
970  q_ref.sstref[i_p];
971  const float cos_solz = q_ref.csolz[i_p];
972  if (dt_analysis < cldmsk::cldthresh &&
973  (val_to_threshold == 0.0 ||
974  val_to_threshold / cos_solz > cldthreshold ||
975  val_to_threshold < cldmsk::invalid_val ||
976  val_to_threshold > cldmsk::max_bt)) {
977  flags_sst[i_p] |= SSTF_REDNONUNIF;
978  qual_sst[i_p] =
979  std::max(qual_sst[i_p], (int8_t) 2);
980  }
981  }
982  }
983  if (is_viirs) {
984  const float rhotRED = StatsBTs.at("red")("VAR")[i_p];
985  const float rhot16 = StatsBTs.at("16")("VAR")[i_p];
986  if (rhotRED > 0.3 && rhot16 >= 0.006 && rhot16 < 0.1 &&
987  ((q_ref.lat[i_p] > (subsolar + 30.0)) ||
988  (q_ref.lat[i_p] < (subsolar - 30.0))) &&
989  glint) {
990  flags_sst[i_p] |= SSTF_CLOUD;
991  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
992  }
993  }
994  }
995  }
996  // uniformity tests
997  {
998  switch (case_switcher) {
999  case cldmsk::SST3:
1000  if (StatsBTs.at("37")("MAXMIN")[i_p] >
1002  flags_sst[i_p] |= SSTF_BTNONUNIF;
1003  if (StatsBTs.at("37")("MAXMIN")[i_p] >
1005  flags_sst[i_p] |= SSTF_BTVNONUNIF;
1006  case cldmsk::SST:
1007  if (StatsBTs.at("11")("MAXMIN")[i_p] >
1009  flags_sst[i_p] |= SSTF_BTNONUNIF;
1010  if (StatsBTs.at("11")("MAXMIN")[i_p] >
1012  flags_sst[i_p] |= SSTF_BTVNONUNIF;
1013  if (StatsBTs.at("12")("MAXMIN")[i_p] >
1015  flags_sst[i_p] |= SSTF_BTNONUNIF;
1016  if (StatsBTs.at("12")("MAXMIN")[i_p] >
1018  flags_sst[i_p] |= SSTF_BTVNONUNIF;
1019  break;
1020  case cldmsk::SST4:
1021  if (StatsBTs.at("39")("MAXMIN")[i_p] >
1023  flags_sst[i_p] |= SSTF_BTNONUNIF;
1024  if (StatsBTs.at("39")("MAXMIN")[i_p] >
1026  flags_sst[i_p] |= SSTF_BTVNONUNIF;
1027  if (StatsBTs.at("40")("MAXMIN")[i_p] >
1029  flags_sst[i_p] |= SSTF_BTNONUNIF;
1030  if (StatsBTs.at("40")("MAXMIN")[i_p] >
1032  flags_sst[i_p] |= SSTF_BTVNONUNIF;
1033  break;
1034  default:
1035  break;
1036  }
1037  {
1038  if ((flags_sst[i_p] & SSTF_BTVNONUNIF) > 0) {
1039  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
1040  }
1041  if ((flags_sst[i_p] & SSTF_BTNONUNIF) > 0) {
1042  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 1);
1043  }
1044  }
1045  }
1046  {
1047  if (treesum.at(i_p) <= 0.0) {
1048  flags_sst[i_p] |= SSTF_CLOUD;
1049  qual_sst[i_p] = std::max(qual_sst[i_p], (int8_t) 3);
1050  }
1051  }
1052  {
1053  if (case_switcher == cldmsk::SST) {
1054  if (night && is_modis) {
1055  bool decrease_night_qual_non_uniform =
1056  (StatsBTs.at("39")("MAXMIN")[i_p] >
1057  cldmsk::Bt39unif1) ||
1058  (StatsBTs.at("40")("MAXMIN")[i_p] >
1060  if (decrease_night_qual_non_uniform) {
1061  qual_sst[i_p] =
1062  std::min((int8_t) 3, ++qual_sst[i_p]);
1063  }
1064  }
1065  }
1066  }
1067  }
1068  }
1069 
1074  void set_treesum(int iscan) {
1075  adt::VarsAtPixel vars;
1076  {
1077  for (const auto &test_data_: test_data) {
1078  vars.vars[test_data_.first] = BAD_FLT;
1079  vars.masks[test_data_.first] = true;
1080  }
1081  }
1082  std::fill(treesum.begin(), treesum.end(), 0.0);
1083  std::fill(flags_sst.begin(), flags_sst.end(), 0);
1084  std::fill(qual_sst.begin(), qual_sst.end(), 0);
1085  for (size_t i_p = 0; i_p < npix; i_p++) {
1086  if (valid_mask.at(
1087  i_p +
1088  i_center * npix)) //(valid_mask.at(i_p + i_center * npix))
1089  {
1090  for (const auto &test_var: test_data) {
1091  vars.vars.at(test_var.first) = test_var.second[i_p];
1092  }
1093  for (const auto &test_mask: test_masks) {
1094  vars.masks.at(test_mask.first) = test_mask.second[i_p];
1095  }
1096  vars.i_p = i_p;
1097  vars.i_scan = iscan;
1098  adt::tree_traversal(treesum.data() + i_p,
1099  flags_sst.data() + i_p, vars,
1100  desicion_tree);
1101  adt::tree_traversal(desicion_tree);
1102  }
1103  }
1104  };
1105 
1110  void desicion_tree_ini(const std::string &product_type) {
1111  treesum = std::vector<float>(npix, 0.0f);
1112  if (desicion_tree == nullptr) {
1113  desicion_tree = new adt::Treenode();
1114  const std::string ocdata_root = envset::get_ocdata_root();
1115  if (product_type == "SST3")
1116  path_to_cloud_mask_config = ocdata_root + "/" + sensor + "/" +
1117  platform +
1118  "/cal/cloud_mask_sst3.json";
1119  if (product_type == "SST")
1120  path_to_cloud_mask_config = ocdata_root + "/" + sensor + "/" +
1121  platform +
1122  "/cal/cloud_mask_sst.json";
1123  if (product_type == "SST4")
1124  path_to_cloud_mask_config = ocdata_root + "/" + sensor + "/" +
1125  platform +
1126  "/cal/cloud_mask_sst4.json";
1127  adt::build_tree(path_to_cloud_mask_config, desicion_tree,
1128  test_names);
1129  }
1130  for (const auto &test_name: test_names) {
1131  test_classifications[test_name] = parse_test_name(test_name);
1132  }
1133  tree_traversal(desicion_tree);
1134  };
1135 
1147  void set_data_cldmsk(
1148  l1str *l1rec, l1qstr *l1qrec, instr *input,
1149  const std::string &product_type,
1150  const std::unordered_map<std::string, std::vector<std::string>>
1151  &tests_classifiers,
1152  bool set_adt_data) {
1153  for (const auto &test_pair: tests_classifiers) {
1154  const std::string &name = test_pair.first;
1155  const std::vector<std::string> &classifier = test_pair.second;
1156  const std::vector<std::string> key_vars = {"BT", "RHO", "CLD"};
1157  std::string pos;
1158  size_t pos_c = 0;
1159  bool key_ided = false;
1160  for (const auto &key_: key_vars) {
1161  if (boost::algorithm::contains(classifier.at(0), key_)) {
1162  pos = key_;
1163  key_ided = true;
1164  pos_c = classifier.at(0).find(pos);
1165  { auto bt = classifier.at(0).substr(pos_c + pos.size()); }
1166  break;
1167  }
1168  }
1169  const size_t size_of_classifier = classifier.size();
1170  if (size_of_classifier == 3) // difference
1171  {
1172  if (boost::algorithm::contains(classifier.at(1), "SST") &&
1173  boost::algorithm::contains(classifier.at(2), "SST")) {
1174  auto product_name_1 = classifier.at(1);
1175  auto sst1 =
1176  external_products.at(product_name_1)
1177  ->get_sst(l1rec, l1qrec, input, product_name_1);
1178  auto product_name_2 = classifier.at(2);
1179  auto sst2 =
1180  external_products.at(product_name_2)
1181  ->get_sst(l1rec, l1qrec, input, product_name_2);
1182  if (PairOfSSTs.count({product_name_1, product_name_2}) ==
1183  0) {
1184  PairOfSSTs[{product_name_1, product_name_2}] =
1185  bstats::PairOfSST(npix, sst1, sst2, true);
1186  PairOfSSTs.at({product_name_1, product_name_2})
1187  .get_difference();
1188  }
1189  if (set_adt_data)
1190  test_data[name] =
1191  PairOfSSTs.at({product_name_1, product_name_2})();
1192  } else if (boost::algorithm::contains(classifier.at(1), "BT") &&
1193  boost::algorithm::contains(classifier.at(2), "BT")) {
1194  assert((classifier.at(1).size() == 4) == 1);
1195  assert((classifier.at(2).size() == 4) == 1);
1196  auto bt1 = classifier.at(1).substr(2, 2);
1197  auto bt2 = classifier.at(2).substr(2, 2);
1198  auto ib1 = ibands.at("ib" + bt1);
1199  auto ib2 = ibands.at("ib" + bt2);
1200  if (PairOfVarsBTs.count({bt1, bt2}) == 0) {
1201  PairOfVarsBTs[{bt1, bt2}] =
1202  bstats::PairOfVarsBT(npix, ib1, ib2, *l1rec, true);
1203  PairOfVarsBTs.at({bt1, bt2}).get_difference(*l1rec);
1204  }
1205  if (set_adt_data)
1206  test_data[name] = PairOfVarsBTs.at({bt1, bt2})();
1207  } else {
1209  path_to_cloud_mask_config, name);
1210  }
1211  } else if (size_of_classifier == 2) // stats
1212  {
1213  if (boost::algorithm::contains(classifier.at(0), "SST")) {
1214  auto product_name = classifier.at(0);
1215  auto sst =
1216  external_products.at(product_name)
1217  ->get_sst(l1rec, l1qrec, input, product_name);
1218  bool std = classifier.at(1) == "STD";
1219  bool maxmin = classifier.at(1) == "MAXMIN";
1220  bool max_ = classifier.at(1) == "MAX";
1221  bool min_ = classifier.at(1) == "MIN";
1222  if (StatsSSTs.count(product_name) == 0) {
1223  StatsSSTs[product_name] = bstats::StatsSST(
1224  npix, n_q_size, bt_box / 2, bt_box / 2, i_center,
1225  l1qrec, sst, std, max_, min_, maxmin);
1226  }
1227  test_data[name] =
1228  StatsSSTs.at(product_name)(classifier.at(1));
1229  } else if (key_ided) {
1230  auto bt = classifier.at(0).substr(pos_c + pos.size());
1231  int ib = -1;
1232  if (!boost::algorithm::contains(classifier.at(0),
1233  "RHOCIRRUS") &&
1234  !boost::algorithm::contains(classifier.at(0), "CLD"))
1235  ib = ibands.at("ib" + bt);
1236  if (boost::algorithm::contains(classifier.at(0), "CLD"))
1237  ib = ibands.at("ibred");
1238  bool std = classifier.at(1) == "STD";
1239  bool maxmin = classifier.at(1) == "MAXMIN";
1240  bool max_ = classifier.at(1) == "MAX";
1241  bool min_ = classifier.at(1) == "MIN";
1242  if (StatsBTs.count(bt) == 0) {
1243  if (boost::algorithm::contains(classifier.at(0), "BT"))
1244  StatsBTs[bt] = bstats::StatsVarBand(
1245  npix, n_q_size, bt_box / 2, bt_box / 2,
1246  i_center, l1qrec, ib, NBANDSIR, cldmsk::bt_mask,
1247  cldmsk::bt_value, std, max_, min_, maxmin);
1248  if (boost::algorithm::contains(classifier.at(0), "RHO"))
1249  StatsBTs[bt] = bstats::StatsVarBand(
1250  npix, n_q_size, bt_box / 2, bt_box / 2,
1251  i_center, l1qrec, ib, n_bands, cldmsk::rho_mask,
1252  cldmsk::rho_value, std, max_, min_, maxmin);
1253  if (boost::algorithm::contains(classifier.at(0),
1254  "RHOCIRRUS"))
1255  StatsBTs[bt] = bstats::StatsVarBand(
1256  npix, n_q_size, bt_box / 2, bt_box / 2,
1257  i_center, l1qrec, ib, n_bands,
1259  max_, min_, maxmin);
1260  if (boost::algorithm::contains(classifier.at(0), "CLD"))
1261  StatsBTs[bt] = bstats::StatsVarBand(
1262  npix, n_q_size, bt_box / 2, bt_box / 2,
1263  i_center, l1qrec, ib, n_bands,
1265  max_, min_, maxmin);
1266  } else {
1267  if (std) StatsBTs.at(bt).set_std();
1268  if (maxmin) StatsBTs.at(bt).set_maxmin();
1269  if (min_) StatsBTs.at(bt).set_min();
1270  if (max_) StatsBTs.at(bt).set_max();
1271  }
1272  if (set_adt_data) {
1273  test_data[name] = StatsBTs.at(bt)(classifier.at(1));
1274  test_masks[name] = StatsBTs.at(bt).get_valid_mask();
1275  }
1276  } else {
1278  path_to_cloud_mask_config, name);
1279  }
1280  } else if (size_of_classifier == 1) // variable
1281  {
1282  if (boost::algorithm::contains(classifier.at(0), "SST")) {
1283  auto product_name = classifier.at(0);
1284  auto sst =
1285  external_products.at(product_name)
1286  ->get_sst(l1rec, l1qrec, input, product_name);
1287  if (StatsSSTs.count(product_name) == 0) {
1288  StatsSSTs[product_name] =
1289  bstats::StatsSST(npix, n_q_size, bt_box / 2,
1290  bt_box / 2, i_center, l1qrec, sst);
1291  }
1292  test_data[name] = StatsSSTs.at(product_name)("VAR");
1293  continue;
1294  } else if (key_ided) {
1295  auto bt = classifier.at(0).substr(pos_c + pos.size());
1296  int ib = -1;
1297  if (!boost::algorithm::contains(classifier.at(0),
1298  "RHOCIRRUS"))
1299  ib = ibands.at("ib" + bt);
1300  if (StatsBTs.count(bt) == 0) {
1301  if (boost::algorithm::contains(classifier.at(0), "BT"))
1302  StatsBTs[bt] = bstats::StatsVarBand(
1303  npix, n_q_size, bt_box / 2, bt_box / 2,
1304  i_center, l1qrec, ib, NBANDSIR, cldmsk::bt_mask,
1306  if (boost::algorithm::contains(classifier.at(0), "RHO"))
1307  StatsBTs[bt] = bstats::StatsVarBand(
1308  npix, n_q_size, bt_box / 2, bt_box / 2,
1309  i_center, l1qrec, ib, n_bands, cldmsk::rho_mask,
1311  if (boost::algorithm::contains(classifier.at(0),
1312  "RHOCIRRUS"))
1313  StatsBTs[bt] = bstats::StatsVarBand(
1314  npix, n_q_size, bt_box / 2, bt_box / 2,
1315  i_center, l1qrec, ib, n_bands,
1317  }
1318  if (set_adt_data) {
1319  test_data[name] = StatsBTs.at(bt)("VAR");
1320  test_masks[name] = StatsBTs.at(bt).get_valid_mask();
1321  }
1322  continue;
1323  } else if (boost::algorithm::contains(classifier.at(0),
1324  "Tdeflong")) {
1325  std::string bt1 = "11";
1326  std::string bt2 = "12";
1327  auto ib1 = ibands.at("ib" + bt1);
1328  auto ib2 = ibands.at("ib" + bt2);
1329  if (PairOfVarsBTs.count({bt1, bt2}) == 0) {
1330  PairOfVarsBTs[{bt1, bt2}] = bstats::PairOfVarsBT(
1331  npix, ib1, ib2, *l1rec, true, true);
1332  PairOfVarsBTs.at({bt1, bt2}).get_difference(*l1rec);
1333  PairOfVarsBTs.at({bt1, bt2}).get_ratio(*l1rec);
1334  } else {
1335  PairOfVarsBTs.at({bt1, bt2}).get_ratio(*l1rec);
1336  }
1337  if (set_adt_data)
1338  test_data[name] =
1339  PairOfVarsBTs.at({bt1, bt2}).get_Tdeflong();
1340  continue;
1341  } else {
1344  path_to_cloud_mask_config, name);
1345  }
1346  if (set_adt_data)
1347  test_data[name] = cldmsk::get_non_BT_vars.at(name)(*l1rec);
1348  }
1349  }
1350  }
1351 
1352  typedef void (SeaSurfaceTemperatureCalcuations::*nlsst_ptr)(const l1str &,
1353  float *,
1354  const int *,
1355  const float *);
1356 
1357  typedef void (SeaSurfaceTemperatureCalcuations::*valid_mask_ptr)(
1358  const l1str &, int *, const float *);
1359 
1360  typedef float *(
1361  SeaSurfaceTemperatureCalcuations::*get_desicion_tree_value_ptr)(
1362  l1str *, l1qstr *, instr *, const std::string &, size_t);
1363 
1364  std::unordered_map<std::string, nlsst_ptr> sst_calls = {
1365  {"SST_v6_cmc", &SeaSurfaceTemperatureCalcuations::get_nlsst_v6_cmc},
1366  {"SST3", &SeaSurfaceTemperatureCalcuations::get_nlsst_viirs3},
1367  {"SST_v7_cmc", &SeaSurfaceTemperatureCalcuations::get_nlsst_v7_cmc},
1368  {"SST3_v7", &SeaSurfaceTemperatureCalcuations::get_mcsst3_v7},
1369  {"SST4",
1370  &SeaSurfaceTemperatureCalcuations::get_nlsst_modis4}}; // ulitmately need to put a model name
1371  std::unordered_map<std::string, valid_mask_ptr> valid_mask_calls = {
1372  {"SST", &SeaSurfaceTemperatureCalcuations::get_valid_mask_sst},
1373  {"SST3", &SeaSurfaceTemperatureCalcuations::get_valid_mask_sst3},
1374  {"SST4", &SeaSurfaceTemperatureCalcuations::get_valid_mask_sst4}};
1375 
1376  public:
1377  SeaSurfaceTemperatureCalcuations(l1str *l1rec, l1qstr *l1qrec, instr *input,
1378  const char *product_type) {
1379  this_product = product_type;
1380  { std::cout << "Product type to ini " << this_product << std::endl; }
1381  npix = l1rec->npix;
1382  const int sensor_id_int = l1rec->l1file->sensorID;
1383  auto it_s = cldmsk::platforms.find(sensor_id_int);
1384  if (it_s != cldmsk::platforms.end())
1385  platform = it_s->second;
1386  else {
1387  std::cerr << "Platform is not found. Exiting ... " << sensor_id_int
1388  << std::endl;
1389  exit(EXIT_FAILURE);
1390  }
1391  if (platform == "aqua") {
1392  l1qstr_modis_a = l1qrec;
1393  is_aqua = true;
1394  {
1395  std::cout << "MODIS AQUA correction to BT 4 micron for "
1396  "detector #0 will be applied\n";
1397  }
1398  } else if (platform == "terra") {
1399  is_terra = true;
1400  std::cout << "MODIS TERRA specific corrections will be applied\n";
1401  }
1402  {
1403  // set sensor
1404  if (cldmsk::modis_sensors.count(sensor_id_int) > 0) {
1405  sensor = "modis";
1406  bt_box = cldmsk::bt_box_sizes.at(sensor);
1407  if (l1_input->resolution == 250) {
1408  fullscanpix = 5416;
1409  } else if (l1_input->resolution == 500) {
1410  fullscanpix = 2708;
1411  }
1412  if (this_product == "SST") {
1413  regression_model = "SST_v7_cmc";
1414  } else {
1415  regression_model = "SST4";
1416  }
1417  // dust correction
1418  {
1419  if (l1rec->anc_aerosol && strlen(input->dsdicoeffile)) {
1420  dust_correction_data = std::vector<float>(npix, 0.0);
1421  {
1422  std::cout << "Applying dust correction "
1423  << std::endl;
1424  }
1425  }
1426  }
1427  } else if (cldmsk::viirs_sensors.count(sensor_id_int) > 0) {
1428  sensor = "viirs";
1429  bt_box = cldmsk::bt_box_sizes.at(sensor);
1430  fullscanpix = 3200;
1431  if (this_product == "SST") {
1432  if (platform == "npp") {
1433  regression_model = "SST_v6_cmc";
1434  } else {
1435  regression_model = "SST_v7_cmc";
1436  }
1437  } else {
1438  regression_model = "SST3";
1439  }
1440  } else {
1441  std::cerr << "Sensor is not found " << sensor_id_int
1442  << std::endl;
1443  exit(EXIT_FAILURE);
1444  }
1445  }
1446  if (ibands.empty()) {
1447  cldmsk::read_sst_bands(sensor, ibands);
1448  n_bands = l1rec->l1file->nbands;
1449  const std::set<std::string> not_adj_bands = {"ib16", "ibred",
1450  "ib07", "ib08"};
1451  for (auto &pair: ibands) {
1452  if (not_adj_bands.find(pair.first) == not_adj_bands.end())
1453  pair.second -= n_bands;
1454  }
1455  }
1456  {
1457  if (ibands.count("ib11") > 0) {
1458  ib11 = ibands.at("ib11");
1459  }
1460  if (ibands.count("ib12") > 0) {
1461  ib12 = ibands.at("ib12");
1462  }
1463  if (ibands.count("ib37") > 0) {
1464  ib37 = ibands.at("ib37");
1465  }
1466  if (ibands.count("ib39") > 0) {
1467  ib39 = ibands.at("ib39");
1468  }
1469  if (ibands.count("ib40") > 0) {
1470  ib40 = ibands.at("ib40");
1471  }
1472  if (ibands.count("ib85") > 0) {
1473  ib85 = ibands.at("ib85");
1474  }
1475  }
1476  { desicion_tree_ini(product_type); }
1477  {
1478  const std::string sst_coef_file =
1480  {
1481  std::cout << "SST Coefficients file for product "
1482  << product_type << " " << sst_coef_file << std::endl;
1483  }
1484  {
1485  netCDF::NcFile sst_reg_coefficients_file(sst_coef_file,
1486  netCDF::NcFile::read);
1487  netCDF::NcGroupAtt reg_model =
1488  sst_reg_coefficients_file.getAtt("regression_model");
1489  if (!reg_model.isNull()) {
1490  size_t attr_len = reg_model.getAttLength();
1491  std::string model_name;
1492  model_name.reserve(attr_len);
1493  reg_model.getValues(model_name);
1494  regression_model = model_name;
1495  std::cout
1496  << "The regression model from the SST coefficient file "
1497  << regression_model << std::endl;
1498  }
1499  std::vector<float> coeffs;
1500  netCDF::NcVar coeffs_nc =
1501  sst_reg_coefficients_file.getVar("coefficients");
1502  // bounds
1503  netCDF::NcVar bounds_nc =
1504  sst_reg_coefficients_file.getVar("latbands");
1505  auto dims = coeffs_nc.getDims();
1506  auto lbdims = bounds_nc.getDims();
1507  size_t reserved_size = 1;
1508  size_t reserved_size_ld = 1;
1509  std::vector<size_t> dims_all;
1510  std::vector<size_t> lddims_all;
1511  for (const auto &dim: dims) {
1512  reserved_size *= dim.getSize();
1513  dims_all.push_back(dim.getSize());
1514  }
1515  for (const auto &dim: lbdims) {
1516  reserved_size_ld *= dim.getSize();
1517  lddims_all.push_back(dim.getSize());
1518  }
1519 
1520  coeffs.resize(reserved_size);
1521  bounds_lat.resize(reserved_size_ld);
1522  coeffs_nc.getVar(coeffs.data());
1523  bounds_nc.getVar(bounds_lat.data());
1524  double scantime = l1rec->scantime;
1525  int16_t year, month;
1526  double secs;
1527  unix2ymds(scantime, &year, & month, & day,& secs);
1528  {
1529  std::cout << "Scantime - unix time: " << scantime << "; ";
1530  std::cout << "scan day : "
1531  << day << "; ";
1532  std::cout
1533  << "scan month : "
1534  << month
1535  << "; ";
1536  cldmsk::month_data() = std::vector<float>(npix, -1.0f);
1537  for (size_t i_p = 0; i_p < npix; i_p++) {
1538  cldmsk::month_data().at(i_p) =
1539  month;
1540  }
1541  std::cout << "scan year : "
1542  << year << "\n";
1543  }
1544  const auto month_counter =
1545  month - 1;
1546 
1547  if (is_terra) {
1548  if (scan_time < cldmsk::scan_time_modis_t_day1) {
1549  std::cout << "Applying electronic correction for MODIS "
1550  "TERRA, 2000 October 29"
1551  << std::endl;
1552  elecor_terra = cldmsk::el_corr_modis_t_1;
1553  } else if (scan_time < cldmsk::scan_time_modis_t_day2) {
1554  std::cout << "Applying electronic correction for MODIS "
1555  "TERRA, 2001 July 1"
1556  << std::endl;
1557  elecor_terra = cldmsk::el_corr_modis_t_2;
1558  }
1559  }
1560  {
1561  std::cout << "Targeted month " << month_counter + 1
1562  << std::endl;
1563  sst_regression_coefficients = std::vector<float>(
1564  coeffs.data() +
1565  dims_all.at(1) * dims_all.at(2) * month_counter,
1566  coeffs.data() + dims_all.at(1) * dims_all.at(2) *
1567  (month_counter + 1));
1568  }
1569  // sst regression coefficients bounds
1570  ncoeffs = dims_all.at(2);
1571  { std::cout << "The regression coefficients used: \n"; }
1572  for (size_t j = 0; j < dims_all.at(1); j++) {
1573  {
1574  const size_t i_lat_st = j * lddims_all.at(1);
1575  const size_t i_lat_e = j * lddims_all.at(1) + 1;
1576  const float lat_st = bounds_lat.at(i_lat_st);
1577  const float lat_e = bounds_lat.at(i_lat_e);
1578  std::cout << "Lat_start : " << lat_st
1579  << "; lat_end : " << lat_e
1580  << "; Coefficients : ";
1581  }
1582  for (size_t k = 0; k < dims_all.at(2); k++) {
1583  const size_t index = j * dims_all.at(2) + k;
1584  std::cout << "C" << k << " = "
1585  << sst_regression_coefficients.at(index)
1586  << "; ";
1587  }
1588  std::cout << "\n";
1589  }
1590  bounds_search.push_back(bounds_lat.at(0));
1591  nlatbands = lddims_all.at(0);
1592  for (size_t j = 0; j < lddims_all.at(0); j++) {
1593  for (size_t k = 0; k < lddims_all.at(1); k++) {
1594  const size_t index = j * lddims_all.at(1) + k;
1595  if (k == 1) {
1596  bounds_search.push_back(bounds_lat.at(index));
1597  }
1598  }
1599  std::cout << "\n";
1600  }
1601  // reading the SSES bias
1602  {
1603  const std::string sses_luts_path =
1605  sses_bias_data = std::make_unique<cldmsk::SSESData>(
1606  sses_luts_path, npix);
1607  }
1608  }
1609  }
1610  };
1611 
1612  float *get_sst(l1str *l1rec, l1qstr *l1qrec, instr *input,
1613  const std::string &product_type) {
1614  const int iscan = l1rec->iscan;
1615  if (sst.empty()) {
1616  {
1617  n_q_size = l1qrec->nq;
1618  i_center = n_q_size / 2;
1619  i_s = std::min(std::max(0, (int) i_center - (int) bt_box / 2),
1620  (int) n_q_size - 1);
1621  i_e = std::max(std::min((int) n_q_size - 1, (int)i_center + (int) bt_box / 2),
1622  (int)0);
1623  }
1624  sst = std::vector<float>(n_q_size * npix, 0.0f);
1625  valid_mask = std::vector<int>(n_q_size * npix, 0);
1626  for (size_t i_q = i_s; i_q <= i_e; i_q++) {
1627  l1str &q_ref = l1qrec->r[i_q];
1628  float *sst_ref =
1629  get_reference(l1rec, l1qrec, input, product_type, i_q);
1630  (this->*valid_mask_calls.at(product_type))(
1631  q_ref, valid_mask.data() + i_q * npix, sst_ref);
1632  (this->*sst_calls.at(regression_model))(
1633  q_ref, sst.data() + i_q * npix,
1634  valid_mask.data() + i_q * npix, sst_ref);
1635  modis_dust_correction(q_ref, i_q);
1636  }
1637  processed_scan_sst = iscan;
1638  } else {
1639  if (processed_scan_sst == iscan ||
1640  processed_scan_sst != iscan - 1) {
1641  return sst.data() + i_center * npix;
1642  }
1643  for (size_t i_q = i_s; i_q < i_e; i_q++) {
1644  const size_t index = (i_q + 1) * npix;
1645  std::copy(valid_mask.begin() + index,
1646  valid_mask.begin() + index + npix,
1647  valid_mask.begin() + index - npix);
1648  }
1649  for (size_t i_q = i_s; i_q < i_e; i_q++) {
1650  const size_t index = (i_q + 1) * npix;
1651  std::copy(sst.begin() + index, sst.begin() + index + npix,
1652  sst.begin() + index - npix);
1653  }
1654  {
1655  size_t i_q = i_e;
1656  l1str &q_ref = l1qrec->r[i_q];
1657  float *sst_ref =
1658  get_reference(l1rec, l1qrec, input, product_type, i_q);
1659  (this->*valid_mask_calls.at(product_type))(
1660  q_ref, valid_mask.data() + i_q * npix, sst_ref);
1661  (this->*sst_calls.at(regression_model))(
1662  q_ref, sst.data() + i_q * npix,
1663  valid_mask.data() + i_q * npix, sst_ref);
1664  modis_dust_correction(q_ref, i_q);
1665  }
1666  processed_scan_sst = iscan;
1667  }
1668  return sst.data() + i_center * npix;
1669 
1670  }
1671 
1672  float *get_std(l1str *l1rec, l1qstr *l1qrec, instr *input,
1673  const std::string &product_type) {
1674  get_bias(l1rec, l1qrec, input, product_type);
1675  return sses_bias_data->stdv_out.data();
1676  }
1677 
1678  float *get_bias(l1str *l1rec, l1qstr *l1qrec, instr *input,
1679  const std::string &product_type) {
1680  const int iscan = l1rec->iscan;
1681  if (processed_scan_sses == iscan || processed_scan_sses != iscan - 1) {
1682  return sses_bias_data->bias_out.data();
1683  } else {
1684  float *diff;
1685  const int16_t *flags_ptr = get_flags_sst(l1rec, l1qrec, input, product_type);
1686  if (this_product == "SST")
1687  diff = PairOfVarsBTs.at({"11", "12"})();
1688  else if (this_product == "SST3")
1689  diff = PairOfVarsBTs.at({"37", "12"})();
1690  else if (this_product == "SST4")
1691  diff = PairOfVarsBTs.at({"39", "40"})();
1692  sses_bias_data->calculate_sses_bias_stats(
1693  diff, sst.data() + npix * i_center, l1rec->solz, l1rec->senz,
1694  day, l1rec->lat, qual_sst.data(), flags_ptr,
1695  l1rec->glint_coef);
1696  }
1697  processed_scan_sses = iscan;
1698  return sses_bias_data->bias_out.data();
1699  }
1700 
1701  float *get_bias_mean(l1str *l1rec, l1qstr *l1qrec, instr *input,
1702  const std::string &product_type) {
1703  get_bias(l1rec, l1qrec, input, product_type);
1704  return sses_bias_data->bias_mean_out.data();
1705  }
1706 
1707  int16_t *get_flags_sst(l1str *l1rec, l1qstr *l1qrec, instr *input,
1708  const std::string &product_type) {
1709  const int iscan = l1rec->iscan;
1710  if (flags_sst.empty()) {
1711  flags_sst = std::vector<int16_t>(npix, 0);
1712  qual_sst = std::vector<int8_t>(npix, 0);
1713  external_products[this_product] =
1714  this; // cloud mask may need other's SSTs
1715  {
1716  std::cout << "Setting up flags_sst for " << product_type
1717  << std::endl;
1718  }
1719  {
1720  set_data_cldmsk(l1rec, l1qrec, input, product_type,
1721  test_classifications, true);
1722  {
1723  const auto &vars =
1724  needed_stats_for_cloud_mask.at({sensor, product_type});
1725  std::unordered_map<std::string, std::vector<std::string>>
1726  non_tree_test_classifications;
1727  for (const auto &test_name: vars) {
1728  non_tree_test_classifications[test_name] =
1729  parse_test_name(test_name);
1730  }
1731  set_data_cldmsk(l1rec, l1qrec, input, product_type,
1732  non_tree_test_classifications, false);
1733  {
1734  std::cout << "The cloud mask has been set for "
1735  << product_type << std::endl;
1736  }
1737  }
1738  {
1739  set_modis_a_bt_40_data(l1rec, l1qrec, input, product_type);
1740  set_treesum(iscan);
1741  cloud_mask(*l1rec, input);
1742  }
1743  }
1744  processed_scan_flags = iscan;
1745  } else {
1746  if (processed_scan_flags == iscan ||
1747  iscan - 1 != processed_scan_flags) {
1748  return flags_sst.data();
1749  }
1750  for (auto &pairSST: PairOfSSTs) {
1751  pairSST.second.sst1 =
1752  external_products.at(pairSST.first.first)
1753  ->get_sst(l1rec, l1qrec, input, pairSST.first.first);
1754  pairSST.second.sst2 =
1755  external_products.at(pairSST.first.second)
1756  ->get_sst(l1rec, l1qrec, input, pairSST.first.second);
1757  pairSST.second.get_difference();
1758  }
1759  for (auto &pairBT: PairOfVarsBTs) {
1760  pairBT.second.get_difference(*l1rec);
1761  }
1762  for (auto &statSST: StatsSSTs) {
1763  statSST.second.sst =
1764  external_products.at(statSST.first)
1765  ->get_sst(l1rec, l1qrec, input, statSST.first);
1766  statSST.second.rearrange();
1767  }
1768  for (auto &statBT: StatsBTs) {
1769  statBT.second.rearrange();
1770  }
1771  {
1772  set_modis_a_bt_40_data(l1rec, l1qrec, input, product_type);
1773  set_treesum(iscan);
1774  cloud_mask(*l1rec, input);
1775  }
1776  processed_scan_flags = iscan;
1777  }
1778  return flags_sst.data();
1779  }
1780 
1782  std::unique_ptr<SeaSurfaceTemperatureCalcuations> &ext_product,
1783  const std::string &ext_product_name) {
1784  external_products[ext_product_name] = ext_product.get();
1785  {
1786  std::cout << "External product reference added : "
1787  << ext_product_name << std::endl;
1788  }
1789  };
1790 
1792  return external_products.size();
1793  }
1794 
1795  int8_t *get_quality_flags(l1str *l1rec, l1qstr *l1qrec, instr *input,
1796  const std::string &product_type) {
1798  return qual_sst.data();
1799  }
1800 
1801  int16_t *get_counts(l1str *l1rec, l1qstr *l1qrec, instr *input,
1802  const std::string &product_type) {
1803  get_bias(l1rec, l1qrec, input, product_type);
1804  return sses_bias_data->counts_out.data();
1805  }
1806 
1807  float *get_dust_correction(l1str *l1rec, l1qstr *l1qrec, instr *input,
1808  const std::string &product_type) {
1809  if (dust_correction_data.empty()) {
1810  std::cerr << "Dust correction is not availible. Exiting ... "
1811  << std::endl;
1812  exit(EXIT_FAILURE);
1813  }
1814  return dust_correction_data.data();
1815  }
1816 
1818 };
1819 
1826 class SST {
1827 public:
1828  std::unordered_map<std::string,
1829  std::unique_ptr<SeaSurfaceTemperatureCalcuations>>
1831  std::set<std::string> initilized_products;
1832  std::unordered_map<
1833  std::string, std::unordered_map<std::string, std::vector<std::string>>>
1834  products_needed = {{"SST", {{"viirs", {"SST3"}}, {"modis", {"SST4"}}}},
1835  {"SST3", {{"viirs", {"SST"}}}},
1836  {"SST4", {{"modis", {"SST"}}}}};
1837 
1838  void ini_product(l1str *l1rec, l1qstr *l1qrec, instr *input,
1839  const char *product_type) {
1840  if (initilized_products.count(product_type) > 0) return;
1841  std::vector<std::string> needed_products;
1842  // setting up needed products
1843  {
1844  const auto sensorId = l1rec->l1file->sensorID;
1845  auto it = products_needed.find(product_type);
1846  if (it != products_needed.end()) {
1847  const auto &entry = it->second;
1848  auto it_e = entry.find(cldmsk::get_sensor(sensorId));
1849  if (it_e != entry.end()) {
1850  needed_products = it_e->second;
1851  }
1852  }
1853  }
1854  if (sst.count(product_type) == 0) {
1855  sst[product_type] =
1856  std::make_unique<SeaSurfaceTemperatureCalcuations>(
1857  l1rec, l1qrec, input, product_type);
1858  for (const auto &needed_product: needed_products) {
1859  ini_product(l1rec, l1qrec, input, needed_product.c_str());
1860  }
1861  }
1862 
1863  if (sst.at(product_type)->get_number_of_external_products() <
1864  needed_products.size()) {
1865  for (const auto &needed_product: needed_products) {
1866  sst.at(product_type)
1867  ->set_ext_products(sst.at(needed_product), needed_product);
1868  }
1869  }
1871  }
1872 
1873  float *get_sst(l1str *l1rec, l1qstr *l1qrec, instr *input,
1874  const char *product_type) {
1875  ini_product(l1rec, l1qrec, input, product_type);
1876  return sst.at(product_type)
1877  ->get_sst(l1rec, l1qrec, input, std::string(product_type));
1878  }
1879 
1880  float *get_sst_stdev(l1str *l1rec, l1qstr *l1qrec, instr *input,
1881  const char *product_type) {
1882  ini_product(l1rec, l1qrec, input, product_type);
1883  return sst.at(product_type)
1884  ->get_std(l1rec, l1qrec, input, std::string(product_type));
1885  }
1886 
1887  float *get_sst_bias(l1str *l1rec, l1qstr *l1qrec, instr *input,
1888  const char *product_type) {
1889  ini_product(l1rec, l1qrec, input, product_type);
1890  return sst.at(product_type)
1891  ->get_bias(l1rec, l1qrec, input, std::string(product_type));
1892  }
1893 
1894  float *get_sst_mean_bias(l1str *l1rec, l1qstr *l1qrec, instr *input,
1895  const char *product_type) {
1896  ini_product(l1rec, l1qrec, input, product_type);
1897  return sst.at(product_type)
1898  ->get_bias_mean(l1rec, l1qrec, input, std::string(product_type));
1899  }
1900 
1901  float *get_dust_correction(l1str *l1rec, l1qstr *l1qrec, instr *input,
1902  const char *product_type) {
1903  ini_product(l1rec, l1qrec, input, product_type);
1904  return sst.at(product_type)
1905  ->get_dust_correction(l1rec, l1qrec, input,
1907  }
1908 
1909  int16_t *get_flags_sst(l1str *l1rec, l1qstr *l1qrec, instr *input,
1910  const char *product_type) {
1911  ini_product(l1rec, l1qrec, input, product_type);
1912  return sst.at(product_type)
1913  ->get_flags_sst(l1rec, l1qrec, input, std::string(product_type));
1914  }
1915 
1916  int8_t *get_qual_flags(l1str *l1rec, l1qstr *l1qrec, instr *input,
1917  const char *product_type) {
1918  ini_product(l1rec, l1qrec, input, product_type);
1919  return sst.at(product_type)
1920  ->get_quality_flags(l1rec, l1qrec, input,
1922  }
1923 
1924  int16_t *get_counts(l1str *l1rec, l1qstr *l1qrec, instr *input,
1925  const char *product_type) {
1926  ini_product(l1rec, l1qrec, input, product_type);
1927  return sst.at(product_type)
1928  ->get_counts(l1rec, l1qrec, input, std::string(product_type));
1929  }
1930 };
1931 namespace cld_mask_sst {
1933 }
1934 
1935 void call_sst(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, float **sst_out_2) {
1936  *sst_out_2 = cld_mask_sst::sst_obj.get_sst(l1rec, l1qrec, input, product_type);
1937 }
1938 void call_sses_bias(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type,
1939  float **sses_bias_out) {
1940  *sses_bias_out = cld_mask_sst::sst_obj.get_sst_bias(l1rec, l1qrec, input, product_type);
1941 }
1942 void call_sses_std(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type,
1943  float **sses_std_out) {
1944  *sses_std_out = cld_mask_sst::sst_obj.get_sst_stdev(l1rec, l1qrec, input, product_type);
1945 }
1946 void call_sst_flags(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type,
1947  int16_t **sst_flags) {
1949 }
1950 void call_qual_flags(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type,
1951  int8_t **qual_flags) {
1952  *qual_flags = cld_mask_sst::sst_obj.get_qual_flags(l1rec, l1qrec, input, product_type);
1953 }
1954 void call_sses_bias_mean(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type,
1955  float **sses_bias_mean_out) {
1956  *sses_bias_mean_out = cld_mask_sst::sst_obj.get_sst_mean_bias(l1rec, l1qrec, input, product_type);
1957 }
1958 void call_dust_correction(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type,
1959  float **dust_correction) {
1961 }
1962 void call_sses_counts(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type,
1963  int16_t **sses_counts) {
1964  *sses_counts = cld_mask_sst::sst_obj.get_counts(l1rec, l1qrec, input, product_type);
1965 }
const float SST3diff2
const float SSTvdiff
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
#define SSTF_SST3VDIFF
Definition: flags_sst.h:17
float * get_bias(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
float * get_sst(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
int j
Definition: decode_rs.h:73
#define SSTF_SST3DIFF
Definition: flags_sst.h:14
int32_t day
const float dBtmax
void call_dust_correction(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, float **dust_correction)
std::ostream & operator<<(std::ostream &out, focs::DataProvider &in)
const float cldthresh
#define SSTF_SST4DIFF
Definition: flags_sst.h:13
#define SSTF_SSTREFDIFF
Definition: flags_sst.h:12
#define NBANDSIR
Definition: filehandle.h:23
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 set_ext_products(std::unique_ptr< SeaSurfaceTemperatureCalcuations > &ext_product, const std::string &ext_product_name)
int n_bands
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
#define SSTF_VHISENZ
Definition: flags_sst.h:25
const float equatorialSouth
float * get_bias_mean(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
const float equatorialWest
read l1rec
int16_t * get_flags_sst(l2str *l2rec)
Retrive Cloud Mask Flags for SST.
Definition: sst.c:179
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
#define SSTF_REDNONUNIF
Definition: flags_sst.h:23
const std::unordered_set< int > viirs_sensors
void read_sst_bands(const std::string &sat, std::unordered_map< std::string, int > &bands)
float mu
float32 * pos
Definition: l1_czcs_hdf.c:35
const float vhisenzv2
const float vhisenz
const float Bt11unif2
#define SSTF_BTBAD
Definition: flags_sst.h:8
void ini_product(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
float * get_dust_correction(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
std::unordered_map< std::string, float > vars
Definition: sst_adt.hpp:44
const float SST4diff2
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 thereby resolving MODur00108 Changed header guard macro names to avoid reserved name resolving MODur00104 Maneuver list file for Terra satellite was updated to include data from Jecue DuChateu Maneuver list files for both Terra and Aqua were updated to include two maneuvers from recent IOT weekly reports The limits for Z component of angular momentum was and to set the fourth GEO scan quality flag for a scan depending upon whether or not it occurred during one of them Added _FillValue attributes to many and changed the fill value for the sector start times from resolving MODur00072 Writes boundingcoordinate metadata to L1A archived metadata For PERCENT *ECS change to treat rather resolving GSFcd01518 Added a LogReport Changed to create the Average Temperatures vdata even if the number of scans is
Definition: HISTORY.txt:301
const float equatorialEast
std::string get_ocdata_root()
Returns OCDATAROOT.
@ string
#define SSTF_BT4REFDIFF
Definition: flags_sst.h:21
const get_value cirrus_value
instr * input
#define SSTF_SSTREFVDIFF
Definition: flags_sst.h:26
void tree_traversal(float *inp_treesum, int16_t *inp_flags_sst, const VarsAtPixel &vars, Treenode *node)
Definition: sst_adt.cpp:81
const float SST3diff1
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.
void call_sst(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, float **sst_out_2)
void print_error_message_for_adt(const std::string &file_path, const std::string &key_word, const std::string &message)
Prints error message for unrecognized keywords in the ADT config file.
Definition: sst_adt.cpp:249
const get_valid cldrh_mask
A derived class to compute stats for an SST.
#define SSTF_ISMASKED
Definition: flags_sst.h:7
const float SSTmax
double precision function f(R1)
Definition: tmd.lp.f:1454
const double scan_time_modis_t_day2
const std::unordered_map< std::string, float > cldthresh_list
int32_t mside
void call_sses_bias_mean(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, float **sses_bias_mean_out)
const float Bt39unif1
#define SSTF_BTNONUNIF
Definition: flags_sst.h:19
SeaSurfaceTemperatureCalcuations(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
#define SSTF_HISENZ
Definition: flags_sst.h:24
const float el_corr_modis_t_1
const get_value cldrh_value
void build_tree(const std::string &config_path, Treenode *decision_tree, std::set< std::string > &tests_name)
Building the from a json file *.
Definition: sst_adt.cpp:198
void call_qual_flags(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, int8_t **qual_flags)
A derived class for a pair of BTs.
A derived class to compute stats for a BT.
#define SSTF_SSTRANGE
Definition: flags_sst.h:11
#define SSTF_BTDIFF
Definition: flags_sst.h:10
#define CtoK
Definition: sstref.c:34
const float invalid_val
const get_value rho_value
const float solznight
#define SSTF_SST4VDIFF
Definition: flags_sst.h:16
l1_input_t * l1_input
Definition: l1_options.c:9
void call_sst_flags(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, int16_t **sst_flags)
int16_t * get_counts(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
const float Btmin
float * get_dust_correction(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
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
void deallocate_tree(Treenode *root)
Definition: sst_adt.cpp:72
#define SSTF_BTVNONUNIF
Definition: flags_sst.h:20
const std::unordered_map< std::string, size_t > bt_box_sizes
int16_t * get_flags_sst(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
subroutine diff(x, conec, n, dconecno, dn, dconecmk, units, u, inno, i, outno, o, input, deriv)
Definition: ffnet.f:205
const float max_bt
float * get_sst_bias(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
const float equatorialNorth
#define RADEG
Definition: czcs_ctl_pt.c:5
const float SSTmin
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
float * get_sst_mean_bias(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
const float dBtrefmin
const float Bt11unif1
int8_t * get_quality_flags(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
std::vector< float > & month_data()
month variable for the desicion tree traversal
float dust_correction(float dustExtinction, float csenz, float BT39, float BT85, float BT11, float BT12, int32_t sensorID)
Definition: sst_dsdi.c:9
int8_t * get_qual_flags(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
std::unordered_map< std::string, std::unique_ptr< SeaSurfaceTemperatureCalcuations > > sst
const float SST4diff1
The decision tree.
Definition: sst_adt.hpp:57
const get_valid cirrus_mask
float * get_sst(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
def center
Definition: sort_gring.py:68
#define BAD_FLT
Definition: jplaeriallib.h:19
const float dBtrefmax
float * get_std(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
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
#define SSTF_BTRANGE
Definition: flags_sst.h:9
void call_sses_counts(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, int16_t **sses_counts)
int16_t * get_counts(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
const get_value bt_value
int32_t iscan
int16_t * get_flags_sst(l1str *l1rec, l1qstr *l1qrec, instr *input, const std::string &product_type)
#define SSTMODS
Definition: l1.h:82
void copy(double **aout, double **ain, int n)
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
float * get_sst_stdev(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type)
const float Bt37unif1
std::set< std::string > initilized_products
const double scan_time_modis_t_day1
const float Btmax
The driver class to interface with the pure legacy C code.
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
const float Bt40unif2
#define SSTF_CLOUD
Definition: flags_sst.h:27
std::unordered_map< std::string, bool > masks
Definition: sst_adt.hpp:45
#define abs(a)
Definition: misc.h:90
const float dBt4max
const float el_corr_modis_t_2
A derived class for a pair of SSTs.
std::unordered_map< std::string, std::unordered_map< std::string, std::vector< std::string > > > products_needed
int k
Definition: decode_rs.h:73
el
Definition: decode_rs.h:168
int npix
Definition: get_cmp.c:28
const float glintmax
const std::unordered_set< int > modis_sensors
l2prod max
void call_sses_bias(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, float **sses_bias_out)
const float hisenz
Class to compute the SST products (the sst, cloud mask etc)
void call_sses_std(l1str *l1rec, l1qstr *l1qrec, instr *input, const char *product_type, float **sses_std_out)
int count
Definition: decode_rs.h:79