ocssw V2020
get_habs.c
Go to the documentation of this file.
1 /*
2  * get_hab.c
3  *
4  * Harmful Algal Blooms
5  * Created on: Aug 31, 2015
6  * Author: Rick Healy (richard.healy@nasa.gov)
7  */
8 #include <stdlib.h>
9 #include <math.h>
10 #include "l12_proto.h"
11 #include <sensorInfo.h>
12 #include "mph_flags.h"
13 /*
14  * Harmful Algal Bloom Indexes and related functions
15  * for flagging products and cloud masking
16  * R. Healy 9/1/2015 (richard.healy@nasa.gov)
17  *
18  * Cyanobacteria Index
19  * Wynne,T.T.; Stumpf, R.P.; 2013. Spatial and Temporal Patterns in the Seasonal Distribution of
20  Toxic Cyanobacteria in Western Lake Erie from 2002–2014,Toxins 2015, 7, 1649-1663; doi:10.3390/toxins7051649
21 
22  Maximum Chlorophyll Index
23  C.E. Binding ⁎, T.A. Greenberg, R.P. Bukata, The MERIS Maximum Chlorophyll Index; its merits and limitations for inland water
24  algal bloom monitoring, JGLR-00579
25  *
26  */
27 
28 static int numFlagMPHPixels = -1;
29 static int numFlagHABSPixels = -1;
30 static uint8_t *flags_mph = NULL;
31 static uint8_t *flags_habs = NULL;
32 static int32_t mask = LAND;
33 
34 void allocateFlagsMPH(int numPix) {
35  if((numFlagMPHPixels != numPix) || !flags_mph) {
36  numFlagMPHPixels = numPix;
37  if(flags_mph)
38  free(flags_mph);
39  flags_mph = (uint8_t*) malloc(numPix);
40  }
41 }
42 void allocateFlagsHABS(int numPix) {
43  if((numFlagHABSPixels != numPix) || !flags_habs) {
44  numFlagHABSPixels = numPix;
45  if(flags_habs)
46  free(flags_habs);
47  flags_habs = (uint8_t*) malloc(numPix);
48  }
49 }
50 
51 void habs_meris_ci_corr(float rhos443, float rhos490, float rhos560,
52  float rhos620, float rhos665, float rhos709, float rhos865, float *ci) {
53 
54  float kd, kd_709, ss_560;
55 
56  kd = (0.7 * ((((rhos620 + rhos665) / 2.) - rhos865) / (((rhos443 + rhos490) / 2.) - rhos865)));
57 
58  //calculate Kd using rhos_709 instead of rhos_620 & rhos_665
59  kd_709 = (0.7 * ((rhos709 - rhos865) / (((rhos443 + rhos490) / 2.) - rhos865)));
60 
61  //709 switching modification for scum conditions
62  if (kd_709 > kd) kd = kd_709;
63 
64  //expect a 560 peak with cyano blooms
65  ss_560 = rhos560 - rhos443 + (rhos443 - rhos620)*(560 - 443) / (620 - 443);
66 
67  //identify over-corrected pixels that would result in negative Kd's
68  if (((((rhos620 + rhos665) / 2.) - rhos865) > 0) &&
69  ((((rhos443 + rhos490) / 2.) - rhos865) > 0)) {
70  if ((kd < 0.25) &&
71  ((rhos865 <= rhos490) ||
72  (rhos865 <= rhos665) ||
73  (rhos865 <= rhos709)) &&
74  (ss_560 < 0.01)) {
75  *ci = 0;
76  }
77  }
78 }
79 
80 void get_habs_ci(l2str *l2rec, l2prodstr *p, float ci[]) {
81 
82  int ib0, ib1, ib2, ib3, ib4, ib5, ib6, ib7;
83  float wav0, wav1, wav2, wav3, wav4, wav5, wav6, wav7, nonci, fac;
84 
85  int ip, ipb;
86  static productInfo_t* ci_product_info;
87 
88  switch (l2rec->l1rec->l1file->sensorID) {
89  case MODISA:
90  case MODIST:
91  fac = 1.3;
92  break;
93  default:
94  fac = 1.0;
95  }
96  switch (p->cat_ix) {
97 
98  case CAT_CI_stumpf:
99  case CAT_CI_cyano:
100  case CAT_CI_noncyano:
101  // Cyanobacteria Index
102  // Wynne, Stumpf algorithm 2013
103  switch (l2rec->l1rec->l1file->sensorID) {
104  case MODISA:
105  case MODIST:
106  wav0 = 547;
107  wav1 = 667;
108  wav2 = 678;
109  wav3 = 748;
110  break;
111  default:
112  wav0 = 620;
113  wav1 = 665;
114  wav2 = 681;
115  wav3 = 709;
116  wav4 = 443;
117  wav5 = 490;
118  wav6 = 560;
119  wav7 = 865;
120  //wav8 = 754;
121  ib4 = bindex_get(wav4);
122  ib5 = bindex_get(wav5);
123  ib6 = bindex_get(wav6);
124  ib7 = bindex_get(wav7);
125  //ib8 = bindex_get(wav8);
126  }
127  ib0 = bindex_get(wav0);
128  break;
129 
130  case CAT_MCI_stumpf:
131  if (l2rec->l1rec->l1file->sensorID != MERIS &&
132  l2rec->l1rec->l1file->sensorID != OLCIS3A &&
133  l2rec->l1rec->l1file->sensorID != OLCIS3B) {
134  printf("MCI not supported for this sensor (%s).\n",
135  sensorId2SensorName(l2rec->l1rec->l1file->sensorID));
136  exit(EXIT_FAILURE);
137  }
138  wav1 = 681;
139  wav2 = 709;
140  wav3 = 754;
141  break;
142 
143  default:
144  printf("HABS_CI: Hmm, something's really messed up.\n");
145  exit(EXIT_FAILURE);
146  }
147 
148  if (ci_product_info == NULL) {
149  ci_product_info = allocateProductInfo();
150  findProductInfo("CI_cyano", l2rec->l1rec->l1file->sensorID, ci_product_info);
151  }
152 
153  ib1 = bindex_get(wav1);
154  ib2 = bindex_get(wav2);
155  ib3 = bindex_get(wav3);
156 
157  if (ib1 < 0 || ib2 < 0 || ib3 < 0) {
158  printf("(M)CI_stumpf: incompatible sensor wavelengths for this algorithm\n");
159  exit(EXIT_FAILURE);
160  }
161 
162  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
163 
164  ipb = l2rec->l1rec->l1file->nbands*ip;
165  // check for near coastal and inland waters (height == 0 elev < -10)
166  // and valid data for processing CI
167  if ((l2rec->l1rec->height[ip] == 0 && l2rec->l1rec->elev[ip] < -10) ||
168  (l2rec->l1rec->flags[ip] & mask) != 0 ||
169  l2rec->l1rec->rhos[ipb + ib1] <= 0.0 ||
170  l2rec->l1rec->rhos[ipb + ib2] <= 0.0 ||
171  l2rec->l1rec->rhos[ipb + ib3] <= 0.0) {
172  ci[ip] = BAD_FLT;
173  l2rec->l1rec->flags[ip] |= PRODFAIL;
174  } else {
175  switch (p->cat_ix) {
176 
177  case CAT_CI_stumpf:
178  case CAT_CI_cyano:
179  case CAT_CI_noncyano:
180  ci[ip] = fac * ((l2rec->l1rec->rhos[ipb + ib3] - l2rec->l1rec->rhos[ipb + ib1])*(wav2 - wav1) / (wav3 - wav1)
181  - (l2rec->l1rec->rhos[ipb + ib2] - l2rec->l1rec->rhos[ipb + ib1]));
182 
183  //following corrections currently only applicable for MERIS/OLCI
184  if (l2rec->l1rec->l1file->sensorID == MERIS || l2rec->l1rec->l1file->sensorID == OLCIS3A
185  || l2rec->l1rec->l1file->sensorID == OLCIS3B) {
186 
187  habs_meris_ci_corr(l2rec->l1rec->rhos[ipb + ib4], l2rec->l1rec->rhos[ipb + ib5],
188  l2rec->l1rec->rhos[ipb + ib6], l2rec->l1rec->rhos[ipb + ib0],
189  l2rec->l1rec->rhos[ipb + ib1], l2rec->l1rec->rhos[ipb + ib3],
190  l2rec->l1rec->rhos[ipb + ib7], &ci[ip]);
191 
192  if (p->cat_ix == CAT_CI_cyano || p->cat_ix == CAT_CI_noncyano) {
193  nonci = l2rec->l1rec->rhos[ipb + ib1]
194  - l2rec->l1rec->rhos[ipb + ib0]
195  + (l2rec->l1rec->rhos[ipb + ib0]
196  - l2rec->l1rec->rhos[ipb + ib2])
197  *(wav1 - wav0) / (wav2 - wav0);
198 
199  if (l2rec->l1rec->rhos[ipb + ib0] <= 0 && p->cat_ix == CAT_CI_noncyano) {
200  ci[ip] = BAD_FLT;
201  l2rec->l1rec->flags[ip] |= PRODFAIL;
202  } else {
203  if (p->cat_ix == CAT_CI_noncyano) {
204 
205  if (nonci >= 0) ci[ip] = 0;
206  } else {
207  if (nonci < 0) ci[ip] = 0;
208  }
209  }
210  }
211  }
212  break;
213  case CAT_MCI_stumpf:
214  ci[ip] = fac * (l2rec->l1rec->rhos[ipb + ib2] - l2rec->l1rec->rhos[ipb + ib1]
215  - (l2rec->l1rec->rhos[ipb + ib3] - l2rec->l1rec->rhos[ipb + ib1])*(wav2 - wav1) / (wav3 - wav1));
216  break;
217  default:
218  ci[ip] = BAD_FLT;
219  break;
220  }
221 
222  // set non-detect levels of CI to the minimum valid value in product.xml definition
223  if (ci[ip] < ci_product_info->validMin) {
224  ci[ip] = ci_product_info->validMin;
225  }
226  }
227  }
228 
229  flags_habs = get_flags_habs(l2rec);
230  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
231  if (flags_habs[ip] != 0){
232  ci[ip] = BAD_FLT;
233  l2rec->l1rec->flags[ip] |= PRODFAIL;
234  }
235  }
236  if (l2rec->l1rec->iscan == l2rec->l1rec->l1file->nscan) {
237  freeProductInfo(ci_product_info);
238  }
239 }
240 /*
241  * Maximum Peak Height of chlorophyll for MERIS
242  *
243  * Updated: J. Scott (6/1/2018) joel.scott@nasa.gov; removed Rrs check, fixed mistyped coefficient, removed negative chl filter
244  *
245  * R. Healy (9/1/2015) richard.healy@nasa.gov
246  *
247  * Mark William Matthews , Daniel Odermatt
248  * Remote Sensing of Environment 156 (2015) 374–382
249  *
250  *
251  */
252 void get_habs_mph(l2str *l2rec, l2prodstr *p, float chl_mph[]) {
253  float wav6, wav7, wav8, wav9, wav10, wav14;
254  int ip, ipb, ib6, ib7, ib8, ib9, ib10, ib14;
255  float Rmax0, Rmax1, wavmax0, wavmax1, ndvi;
256  float sipf, sicf, bair, mph0, mph1;
257  float *rhos = l2rec->l1rec->rhos;
258 
259  if ((l2rec->l1rec->l1file->sensorID != MERIS) &&
260  (l2rec->l1rec->l1file->sensorID != OLCIS3A) &&
261  (l2rec->l1rec->l1file->sensorID != OLCIS3B)) {
262  printf("MPH not supported for this sensor (%s).\n",
263  sensorId2SensorName(l2rec->l1rec->l1file->sensorID));
264  exit(EXIT_FAILURE);
265  }
266 
267  wav6 = 620;
268  wav7 = 664;
269  wav8 = 681;
270  wav9 = 709;
271  wav10 = 753;
272  wav14 = 885;
273 
274  ib6 = bindex_get(wav6);
275  ib7 = bindex_get(wav7);
276  ib8 = bindex_get(wav8);
277  ib9 = bindex_get(wav9);
278  ib10 = bindex_get(wav10);
279  ib14 = bindex_get(wav14);
280 
281  if (ib6 < 0 || ib7 < 0 || ib8 < 0 || ib9 < 0 || ib10 < 0 || ib14 < 0) {
282  printf("MPH_stumpf: incompatible sensor wavelengths for this algorithm\n");
283  exit(EXIT_FAILURE);
284  }
285 
286  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
287  ipb = l2rec->l1rec->l1file->nbands*ip;
288 
289  //if (l2rec->Rrs[ipb + ib6] <= 0.0) {
290  // if(l2rec->Rrs[ipb+ib6] <= 0.0 || l2rec->Rrs[ipb+ib7] <= 0.0 || l2rec->Rrs[ipb+ib8] <= 0.0
291  // || l2rec->Rrs[ipb+ib9] <= 0.0) {
292  //chl_mph[ip] = BAD_FLT;
293  //l2rec->l1rec->flags[ip] |= PRODFAIL;
294  //} else {
295  if (rhos[ipb + ib8] > rhos[ipb + ib9]) {
296  wavmax0 = wav8;
297  Rmax0 = rhos[ipb + ib8];
298  } else {
299  wavmax0 = wav9;
300  Rmax0 = rhos[ipb + ib9];
301  }
302  if (Rmax0 > rhos[ipb + ib10]) {
303  wavmax1 = wavmax0;
304  Rmax1 = Rmax0;
305  } else {
306  wavmax1 = wav10;
307  Rmax1 = rhos[ipb + ib10];
308  }
309 
310  //sun-induced phycocyanin absorption fluorescence
311  sipf = rhos[ipb + ib7] - rhos[ipb + ib6] - (rhos[ipb + ib8] - rhos[ipb + ib6]) * (664 - 619) / (681 - 619);
312  //sun induced chlorophyll fluorescence
313  sicf = rhos[ipb + ib8] - rhos[ipb + ib7] - (rhos[ipb + ib9] - rhos[ipb + ib7]) * (681 - 664) / (709 - 664);
314  //normalised difference vegetation index
315  ndvi = (rhos[ipb + ib14] - rhos[ipb + ib7]) / (rhos[ipb + ib14] + rhos[ipb + ib7]);
316  //backscatter and absorption induced reflectance
317  bair = rhos[ipb + ib9] - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (709 - 664) / (885 - 664);
318  mph0 = Rmax0 - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (wavmax0 - 664) / (885 - 664);
319  mph1 = Rmax1 - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (wavmax1 - 664) / (885 - 664);
320 
321  if (wavmax1 != wav10) {
322 
323  if (sicf >= 0 || sipf <= 0 || bair <= 0.002) {
324  chl_mph[ip] = 5.24e9 * pow(mph0, 4) - 1.95e8 * pow(mph0, 3) + 2.46e6 * pow(mph0, 2) + 4.02e3 * mph0 + 1.97;
325  } else {
326  chl_mph[ip] = 22.44 * exp(35.79 * mph1);
327  }
328  } else {
329  if (mph1 >= 0.02 || ndvi >= 0.2) {
330 
331  if (sicf < 0 && sipf > 0) {
332  chl_mph[ip] = 22.44 * exp(35.79 * mph1);
333 
334  } else {
335  chl_mph[ip] = BAD_FLT;
336  }
337  } else {
338  chl_mph[ip] = 5.24e9 * pow(mph0, 4) - 1.95e8 * pow(mph0, 3) + 2.46e6 * pow(mph0, 2) + 4.02e3 * mph0 + 1.97;
339  }
340  }
341  //}
342  //if (chl_mph[ip] < 0.) chl_mph[ip] = 0.;
343  }
344 
345 }
346 
347 uint8_t* get_flags_habs_mph(l2str *l2rec) {
348  //MPH - Maximum Peak Height
349  // from "Improved algorithm for routine monitoring of cyanobacteria and
350  // eutrophication in inland and near-coastal waters"
351  // Matthews and Odermatt, Remote Sensing of Environment (doi:10.1016/j.rse.2014.10.010)
352  //
353 
354  float wav6, wav7, wav8, wav9, wav10, wav14;
355  int ip, ipb, ib6, ib7, ib8, ib9, ib10, ib14;
356  float Rmax0, Rmax1, wavmax0, wavmax1, ndvi;
357  float sipf, sicf, bair, mph1, chl_mph;
358  float *rhos = l2rec->l1rec->rhos;
359  static float thresh = 350;
360 
361  if ((l2rec->l1rec->l1file->sensorID != MERIS) &&
362  (l2rec->l1rec->l1file->sensorID != OLCIS3A) &&
363  (l2rec->l1rec->l1file->sensorID != OLCIS3B)) {
364  printf("MPH not supported for this sensor (%s).\n",
365  sensorId2SensorName(l2rec->l1rec->l1file->sensorID));
366  exit(EXIT_FAILURE);
367  }
368 
369  allocateFlagsMPH(l2rec->l1rec->npix);
370 
371  wav6 = 620;
372  wav7 = 664;
373  wav8 = 681;
374  wav9 = 709;
375  wav10 = 753;
376  wav14 = 885;
377 
378  ib6 = bindex_get(wav6);
379  ib7 = bindex_get(wav7);
380  ib8 = bindex_get(wav8);
381  ib9 = bindex_get(wav9);
382  ib10 = bindex_get(wav10);
383  ib14 = bindex_get(wav14);
384 
385  if (ib6 < 0 || ib7 < 0 || ib8 < 0 || ib9 < 0 || ib10 < 0 || ib14 < 0) {
386  printf("MPH_stumpf: incompatible sensor wavelengths for this algorithm\n");
387  exit(EXIT_FAILURE);
388  }
389 
390  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
391 
392  flags_mph[ip] = 0;
393  ipb = l2rec->l1rec->l1file->nbands*ip;
394 
395  // if(l2rec->rhos[ipb+ib6] <= 0.0 ) {
396  if (l2rec->l1rec->rhos[ipb + ib6] <= 0.0 ||
397  l2rec->l1rec->rhos[ipb + ib7] <= 0.0 ||
398  l2rec->l1rec->rhos[ipb + ib8] <= 0.0 ||
399  l2rec->l1rec->rhos[ipb + ib9] <= 0.0 ||
400  (l2rec->l1rec->flags[ip] & LAND) != 0 ||
401  (l2rec->l1rec->flags[ip] & NAVFAIL) != 0) {
402  l2rec->l1rec->flags[ip] |= PRODFAIL;
403  flags_mph[ip] |= MPH_BADINPUT;
404  } else {
405  if (rhos[ipb + ib8] > rhos[ipb + ib9]) {
406  wavmax0 = wav8;
407  Rmax0 = rhos[ipb + ib8];
408  } else {
409  wavmax0 = wav9;
410  Rmax0 = rhos[ipb + ib9];
411  }
412  if (Rmax0 > rhos[ipb + ib10]) {
413  wavmax1 = wavmax0;
414  Rmax1 = Rmax0;
415  } else {
416  wavmax1 = wav10;
417  Rmax1 = rhos[ipb + ib10];
418  }
419 
420  //sun-induced phycocyanin absorption fluorescence
421  sipf = rhos[ipb + ib7] - rhos[ipb + ib6] - (rhos[ipb + ib8] - rhos[ipb + ib6]) * (664 - 619) / (681 - 619);
422  //sun induced chlorophyll fluorescence
423  sicf = rhos[ipb + ib8] - rhos[ipb + ib7] - (rhos[ipb + ib9] - rhos[ipb + ib7]) * (681 - 664) / (709 - 664);
424  //normalised difference vegetation index
425  ndvi = (rhos[ipb + ib14] - rhos[ipb + ib7]) / (rhos[ipb + ib14] + rhos[ipb + ib7]);
426  //backscatter and absorption induced reflectance
427  bair = rhos[ipb + ib9] - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (709 - 664) / (885 - 664);
428  //mph0 = Rmax0 - rhos[ipb+ib7] - ( rhos[ipb+ib14] - rhos[ipb+ib7]) * (wavmax0 - 664)/(885 - 664);
429  mph1 = Rmax1 - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (wavmax1 - 664) / (885 - 664);
430 
431  if (wavmax1 != wav10) {
432 
433  if (sicf < 0 && sipf > 0 && bair > 0.002) {
434  flags_mph[ip] |= MPH_CYANO;
435  chl_mph = 22.44 * exp(35.79 * mph1);
436  if (chl_mph > thresh)
437  flags_mph[ip] |= MPH_FLOAT;
438  }
439  } else {
440  if (mph1 >= 0.02 || ndvi >= 0.2) {
441  flags_mph[ip] |= MPH_FLOAT;
442 
443  if (sicf < 0 && sipf > 0) {
444  flags_mph[ip] |= MPH_CYANO;
445  chl_mph = 22.44 * exp(35.79 * mph1);
446  if (chl_mph > thresh)
447  flags_mph[ip] |= MPH_FLOAT;
448 
449  }
450  } else {
451  flags_mph[ip] |= MPH_ADJ;
452  }
453  }
454 
455  }
456 
457  }
458  return flags_mph;
459 }
460 
461 uint8_t* get_flags_habs_meris(l2str *l2rec) {
462  // Cloud Masking for MERIS
463 
464  int ib443, ib490, ib510, ib560, ib620, ib665, ib681, ib709, ib754, ib865, ib885;
465  int ip, ipb;
466  float *rhos = l2rec->l1rec->rhos;
467  float mdsi, cv, mean, sum, sdev;
468  int i, n=7;
469 
470  allocateFlagsHABS(l2rec->l1rec->npix);
471 
472  if (l2rec->l1rec->l1file->sensorID != MERIS &&
473  l2rec->l1rec->l1file->sensorID != OLCIS3A &&
474  l2rec->l1rec->l1file->sensorID != OLCIS3B) {
475  printf("HABS flags not supported for this sensor (%s).\n",
476  sensorId2SensorName(l2rec->l1rec->l1file->sensorID));
477  exit(EXIT_FAILURE);
478  }
479 
480  ib443 = bindex_get(443);
481  ib490 = bindex_get(490);
482  ib510 = bindex_get(510);
483  ib560 = bindex_get(560);
484  ib620 = bindex_get(620);
485  ib665 = bindex_get(665);
486  ib681 = bindex_get(681);
487  ib709 = bindex_get(709);
488  ib754 = bindex_get(754);
489  ib865 = bindex_get(865);
490  ib885 = bindex_get(885);
491 
492  if (ib443 < 0 || ib490 < 0 || ib510 < 0 || ib560 < 0 || ib620 < 0 || ib665 < 0 || ib681 < 0 || ib709 < 0 || ib754 < 0 || ib865 < 0 || ib885 < 0) {
493  printf("get_flags_habs: incompatible sensor wavelengths for this algorithm\n");
494  exit(EXIT_FAILURE);
495  }
496 
497  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
498  flags_habs[ip] = 0;
499  ipb = l2rec->l1rec->l1file->nbands*ip;
500 
501  // if the eval cloud masking is used, just use the cloud mask
502  // otherwise, run the get_cldmask function
503  if ((input->evalmask | 2) == 2){
504  flags_habs[ip] = (uint8_t) l2rec->l1rec->cloud[ip];
505  } else {
506  if (get_cldmask(l2rec->l1rec, ip)) {
507  flags_habs[ip] |= HABS_CLOUD;
508  }
509  }
510 
511  if (rhos[ipb + ib443] >= 0.0 &&
512  rhos[ipb + ib490] >= 0.0 &&
513  rhos[ipb + ib510] >= 0.0 &&
514  rhos[ipb + ib560] >= 0.0 &&
515  rhos[ipb + ib620] >= 0.0 &&
516  rhos[ipb + ib665] >= 0.0 &&
517  rhos[ipb + ib681] >= 0.0 &&
518  rhos[ipb + ib709] >= 0.0 &&
519  rhos[ipb + ib754] >= 0.0 &&
520  rhos[ipb + ib885] >= 0.0) {
521 
522  if (rhos[ipb + ib885] > rhos[ipb + ib620] &&
523  rhos[ipb + ib885] > rhos[ipb + ib709] &&
524  rhos[ipb + ib885] > rhos[ipb + ib754] &&
525  rhos[ipb + ib885] > 0.01) {
526  flags_habs[ip] |= HABS_NONWTR;
527  }
528 
529  // test dry lake condition: (rhos_620 > rhos_560)
530  if ((rhos[ipb + ib620] > rhos[ipb + ib560]) &&
531  (rhos[ipb + ib560] > 0.15) &&
532  (rhos[ipb + ib885] > 0.15)){
533  flags_habs[ip] |= HABS_NONWTR;
534  }
535 
536  // test snow/ice condition - meris differential snow index
537  mdsi = (rhos[ipb + ib865] - rhos[ipb + ib885]) / (rhos[ipb + ib865] + rhos[ipb + ib885]);
538 
539  // test snow/ice condition - exclude potential bloom conditions (based on coefficient of variation in visible bands)
540  int b[7] = {ib443, ib490, ib510, ib560, ib620, ib665, ib681};
541  sum = 0;
542  for(i = 0; i < n; ++i) {
543  sum += rhos[ipb + b[i]];
544  }
545  mean = sum / n;
546 
547  sdev = 0;
548  for(i = 0; i < n; ++i) {
549  sdev += pow((rhos[ipb + b[i]] - mean),2);
550  }
551  cv = sqrt(sdev / (n - 1)) / mean;
552 
553  if ((mdsi > 0.01) && (rhos[ipb + ib885] > 0.15) && (cv < 0.1)) {
554  flags_habs[ip] |= HABS_SNOWICE;
555  }
556  //adjacency flagging
557  float mci = 1.0 * (rhos[ipb + ib709] - rhos[ipb + ib681]
558  + (rhos[ipb + ib681] - rhos[ipb + ib754])*(709 - 681) / (754 - 681));
559  float ci = 1.0 * ((rhos[ipb + ib709] - rhos[ipb + ib665])*(681 - 665) / (709 - 665)
560  - (rhos[ipb + ib681] - rhos[ipb + ib665]));
561 
562  habs_meris_ci_corr(rhos[ipb + ib443], rhos[ipb + ib490],
563  rhos[ipb + ib560], rhos[ipb + ib620],
564  rhos[ipb + ib665], rhos[ipb + ib709],
565  rhos[ipb + ib865], &ci);
566 
567  if ((mci < 0) && (ci > 0)) {
568  flags_habs[ip] |= HABS_ADJ;
569  }
570  } else {
571  if (rhos[ipb + ib620] < 0.0 ||
572  rhos[ipb + ib665] < 0.0 ||
573  rhos[ipb + ib681] < 0.0 ||
574  rhos[ipb + ib709] < 0.0) {
575  flags_habs[ip] |= HABS_BADINPUT;
576  }
577  }
578  }
579  return flags_habs;
580 }
581 
582 uint8_t* get_flags_habs_modis(l2str *l2rec) {
583  // Cloud Masking for MERIS
584 
585  int ib469, ib555, ib645, ib667, ib859, ib1240, ib2130;
586  int ip, ipb;
587  float *rhos = l2rec->l1rec->rhos, *Rrs = l2rec->Rrs, *cloud_albedo = l2rec->l1rec->cloud_albedo;
588  float ftemp, ftemp2, ftemp3;
589  float cloudthr = 0.027;
590 
591  allocateFlagsHABS(l2rec->l1rec->npix);
592 
593  if (l2rec->l1rec->l1file->sensorID != MODISA && l2rec->l1rec->l1file->sensorID != MODIST) {
594  printf("HABS flags not supported for this sensor (%s).\n",
595  sensorId2SensorName(l2rec->l1rec->l1file->sensorID));
596  exit(EXIT_FAILURE);
597  }
598 
599  ib469 = bindex_get(469);
600  ib555 = bindex_get(555);
601  ib645 = bindex_get(645);
602  ib667 = bindex_get(667);
603  ib859 = bindex_get(859);
604  ib1240 = bindex_get(1240);
605  ib2130 = bindex_get(2130);
606 
607  if (ib469 < 0 || ib555 < 0 || ib645 < 0 || ib667 < 0 || ib859 < 0 || ib1240 < 0 || ib2130 < 0) {
608  printf("get_flags_habs: incompatible sensor wavelengths for this algorithm\n");
609  exit(EXIT_FAILURE);
610  }
611 
612  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
613  flags_habs[ip] = 0;
614  ipb = l2rec->l1rec->l1file->nbands*ip;
615  // TODO: make consistent with get_cloudmask_modis and use it instead of
616  // this duplication of code...
617 
618  if (l2rec->chl[ip] < 0)
619  ftemp = Rrs[ipb + ib667];
620  else
621  ftemp = Rrs[ipb + ib667]*(0.45 + l2rec->chl[ip]*0.005) / 4.3;
622  // first correct for turbid water
623  if (Rrs[ipb + ib667] < 0.0) ftemp = 0.0;
624  ftemp2 = cloud_albedo[ip] - ftemp;
625 
626  if (ftemp2 > 0.027) flags_habs[ip] |= HABS_CLOUD;
627  // non-water check 1240 is bright relative to 859 and the combination is bright
628  // this may hit glint by accident, need to be checked.
629 
630  if (rhos[ipb + ib1240] / rhos[ipb + ib859] > 0.5 && (rhos[ipb + ib1240] + rhos[ipb + ib2130]) > 0.10) flags_habs[ip] |= HABS_CLOUD;
631 
632  // now try to correct for glint
633  // region check was thrown out {IF (region = "OM") cloudthr = 0.04} rjh 11/2/2015
634 
635  ftemp = rhos[ipb + ib645] - rhos[ipb + ib555] + (rhos[ipb + ib555] - rhos[ipb + ib859])*(645.0 - 555.0) / (859.0 - 555.0);
636  ftemp2 = cloud_albedo[ip] + ftemp;
637  if (ftemp2 < cloudthr) flags_habs[ip] = 0;
638  if (rhos[ipb + ib859] / rhos[ipb + ib1240] > 4.0) flags_habs[ip] = 0;
639 
640  // scum areas
641 
642  if ((rhos[ipb + ib859] - rhos[ipb + ib469]) > 0.01 && cloud_albedo[ip] < 0.30) flags_habs[ip] = 0;
643  if ((rhos[ipb + ib859] - rhos[ipb + ib645]) > 0.01 && cloud_albedo[ip] < 0.15) flags_habs[ip] = 0;
644  if (rhos[ipb + ib1240] < 0.2)
645  ftemp2 = ftemp2 - (rhos[ipb + ib859] - rhos[ipb + ib1240]) * fabs(rhos[ipb + ib859] - rhos[ipb + ib1240]) / cloudthr;
646 
647  ftemp3 = ftemp2;
648  if (ftemp2 < cloudthr * 2) {
649  if ((rhos[ipb + ib555] - rhos[ipb + ib1240]) > (rhos[ipb + ib469] - rhos[ipb + ib1240])) {
650  ftemp3 = ftemp2 - (rhos[ipb + ib555] - rhos[ipb + ib1240]);
651  } else {
652  ftemp3 = ftemp2 - (rhos[ipb + ib469] - rhos[ipb + ib1240]);
653  }
654  }
655 
656  if (ftemp3 < cloudthr) flags_habs[ip] = 0;
657 
658  if (rhos[ipb + ib555] >= 0 && rhos[ipb + ib1240] >= 0.0 && rhos[ipb + ib1240] > rhos[ipb + ib555])
659  flags_habs[ip] |= HABS_NONWTR;
660  }
661  return flags_habs;
662 }
663 
664 uint8_t* get_flags_habs(l2str *l2rec) {
665  static int32_t currentLine = -1;
666 
667  if (currentLine == l2rec->l1rec->iscan )
668  return flags_habs;
669  currentLine = l2rec->l1rec->iscan;
670  switch (l2rec->l1rec->l1file->sensorID) {
671  case MERIS:
672  case OLCIS3A:
673  case OLCIS3B:
674  return get_flags_habs_meris(l2rec);
675  break;
676  case MODISA:
677  case MODIST:
678  return get_flags_habs_modis(l2rec);
679  break;
680  default:
681  printf("HABS flags not supported for this sensor (%s).\n",
682  sensorId2SensorName(l2rec->l1rec->l1file->sensorID));
683  exit(EXIT_FAILURE);
684  }
685  return NULL;
686 }
get_cldmask
char get_cldmask(l1str *l1rec, int32_t ip)
Definition: cloud_flag.c:345
get_habs_mph
void get_habs_mph(l2str *l2rec, l2prodstr *p, float chl_mph[])
Definition: get_habs.c:252
OLCIS3A
#define OLCIS3A
Definition: sensorDefs.h:32
freeProductInfo
void freeProductInfo(productInfo_t *info)
Definition: productInfo.cpp:130
HABS_ADJ
#define HABS_ADJ
Definition: mph_flags.h:21
HABS_CLOUD
#define HABS_CLOUD
Definition: mph_flags.h:18
mean
float mean(float *xs, int sample_size)
Definition: numerical.c:81
NULL
#define NULL
Definition: decode_rs.h:63
CAT_CI_noncyano
#define CAT_CI_noncyano
Definition: l2prod.h:371
MERIS
#define MERIS
Definition: sensorDefs.h:22
HABS_NONWTR
#define HABS_NONWTR
Definition: mph_flags.h:19
MODIST
#define MODIST
Definition: sensorDefs.h:18
get_habs_ci
void get_habs_ci(l2str *l2rec, l2prodstr *p, float ci[])
Definition: get_habs.c:80
MPH_BADINPUT
#define MPH_BADINPUT
Definition: mph_flags.h:15
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:39
get_flags_habs_meris
uint8_t * get_flags_habs_meris(l2str *l2rec)
Definition: get_habs.c:461
bindex_get
int bindex_get(int32_t wave)
Definition: windex.c:43
input
instr * input
Definition: main_l2binmatch.cpp:27
habs_meris_ci_corr
void habs_meris_ci_corr(float rhos443, float rhos490, float rhos560, float rhos620, float rhos665, float rhos709, float rhos865, float *ci)
Definition: get_habs.c:51
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
l12_proto.h
fac
#define fac
Definition: atmocor1_land.c:111
get_flags_habs
uint8_t * get_flags_habs(l2str *l2rec)
Definition: get_habs.c:664
HABS_BADINPUT
#define HABS_BADINPUT
Definition: mph_flags.h:22
allocateProductInfo
productInfo_t * allocateProductInfo()
Definition: productInfo.cpp:118
MPH_FLOAT
#define MPH_FLOAT
Definition: mph_flags.h:12
MPH_ADJ
#define MPH_ADJ
Definition: mph_flags.h:14
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
allocateFlagsHABS
void allocateFlagsHABS(int numPix)
Definition: get_habs.c:42
findProductInfo
int findProductInfo(const char *productName, int sensorId, productInfo_t *info)
Definition: productInfo.cpp:871
get_flags_habs_modis
uint8_t * get_flags_habs_modis(l2str *l2rec)
Definition: get_habs.c:582
b
data_t b[NROOTS+1]
Definition: decode_rs.h:77
allocateFlagsMPH
void allocateFlagsMPH(int numPix)
Definition: get_habs.c:34
MPH_CYANO
#define MPH_CYANO
Definition: mph_flags.h:13
sensorId2SensorName
const char * sensorId2SensorName(int sensorId)
Definition: sensorInfo.c:175
CAT_MCI_stumpf
#define CAT_MCI_stumpf
Definition: l2prod.h:327
fabs
#define fabs(a)
Definition: misc.h:93
HABS_SNOWICE
#define HABS_SNOWICE
Definition: mph_flags.h:20
OLCIS3B
#define OLCIS3B
Definition: sensorDefs.h:41
i
int i
Definition: decode_rs.h:71
mph_flags.h
get_flags_habs_mph
uint8_t * get_flags_habs_mph(l2str *l2rec)
Definition: get_habs.c:347
MODISA
#define MODISA
Definition: sensorDefs.h:19
CAT_CI_stumpf
#define CAT_CI_stumpf
Definition: l2prod.h:326
LAND
@ LAND
Definition: Granule.h:41
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
CAT_CI_cyano
#define CAT_CI_cyano
Definition: l2prod.h:370
NAVFAIL
#define NAVFAIL
Definition: l2_flags.h:34