OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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