NASA Logo
Ocean Color Science Software

ocssw V2022
setflags_l2.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 
3 char isCoccolith(l2str *l2rec, int32_t ip);
4 char isTurbid(l2str *l2rec, int32_t ip);
5 float aerindex(l2str *l2rec, int32_t ip);
6 char isOptShallow(l2str *l2rec, int32_t ip);
7 
9  RHOWINDEX = 1,
11 };
12 
13 
14 void setflagbits_l2(l2str *l2rec, int32_t ipix) {
15  static int firstCall = 1;
16  static int ib490;
17  static int ib510;
18  static int ib555;
19 
20  int32_t npix, spix, epix;
21  int32_t ip;
22  int32_t nwave;
23 
24  l1str *l1rec = l2rec->l1rec;
25 
26  if (l1rec == NULL) {
27  printf("-E- %s line %d: attempt to set flags from NULL L1 record.\n",
28  __FILE__, __LINE__);
29  exit(FATAL_ERROR);
30  }
31 
32  npix = l1rec->npix;
33  nwave = l1rec->l1file->nbands;
34 
35  if (ipix < 0) {
36  spix = 0;
37  epix = npix - 1;
38  } else {
39  spix = ipix;
40  epix = ipix;
41  }
42 
43  if (l2rec == NULL) {
44  printf("-E- %s line %d: attempt to set flags from NULL L2 record.\n",
45  __FILE__, __LINE__);
46  exit(FATAL_ERROR);
47  }
48 
49  if (firstCall) {
50  firstCall = 0;
51  float* wave = l2rec->l1rec->l1file->fwave;
52  ib490 = windex(490., wave, nwave);
53  ib510 = windex(510., wave, nwave);
54  ib555 = windex(550., wave, nwave);
55  }
56 
57  for (ip = spix; ip <= epix; ip++) {
58  if (l2rec->eps[ip] < input->epsmin ||
59  l2rec->eps[ip] > input->epsmax || l2rec->chi2[ip]>5.) {
60  l1rec->flags[ip] |= ATMWARN;
61  }
62 
63  if (!l1rec->land[ip] && !l1rec->cloud[ip] && !l1rec->ice[ip]){
64  if (isOptShallow(l2rec,ip))
65  l1rec->flags[ip] |= OPSHAL;
66  }
67 
68  if ((l2rec->Lw[ip * nwave + ib490] < 0.0) ||
69  (l2rec->Lw[ip * nwave + ib510] < 0.0) ||
70  (l2rec->Lw[ip * nwave + ib555] < 0.0)) {
71  l1rec->flags[ip] |= ATMWARN;
72  }
73 
74  if (l2rec->nLw[ip * nwave + ib555] < input->nlwmin)
75  l1rec->flags[ip] |= LOWLW;
76 
77  if (isCoccolith(l2rec, ip))
78  l1rec->flags[ip] |= COCCOLITH;
79 
80  if (input->aer_opt != AERWANGSWIR && input->aer_opt != AERRHSWIR) {
81  if (isTurbid(l2rec, ip))
82  l1rec->flags[ip] |= TURBIDW;
83  }
84 
85  if (l2rec->chl[ip] == BAD_FLT)
86  l1rec->flags[ip] |= CHLFAIL;
87 
88  else if (l2rec->chl[ip] > CHL_MAX || l2rec->chl[ip] < CHL_MIN)
89  l1rec->flags[ip] |= CHLWARN;
90 
91  if (l2rec->num_iter[ip] > input->aer_iter_max) {
92  l1rec->flags[ip] |= MAXAERITER;
93  l1rec->flags[ip] |= ATMWARN;
94  }
95 
96  if (input->absaer_opt > 0) {
97  switch (input->absaer_opt) {
98  case RHOWINDEX:
99  l2rec->aerindex[ip] = aerindex(l2rec, ip);
100  if (l2rec->aerindex[ip] < input->absaer)
101  l1rec->flags[ip] |= ABSAER;
102  break;
103  case TOTEXINCTION:
104  if (l1rec->anc_aerosol) {
105  float aaod = l1rec->anc_aerosol->total_aerosol_ext[ip]
106  - l1rec->anc_aerosol->total_aerosol_scat[ip];
107  if (aaod > 0.01)
108  l1rec->flags[ip] |= ABSAER;
109  } else {
110  printf("Warning: Cannot apply GMAO total extinction absorbing aerosol test without\n");
111  printf(" anc_aerosol inputs defined\n");
112  printf(" Setting absaer_opt=0\n\n");
113  input->absaer_opt = 0;
114  }
115  break;
116  }
117  }
118  }
119 }
120 
121 
122 #define Between(a,x,b)(x >= a && x <= b)
123 
124 char isCoccolith(l2str *l2rec, int32_t ip) {
125  static float firstCall = 1;
126  static float c1, c2, c3, c4, c5, c6, c7, c8;
127  static int ib443;
128  static int ib510;
129  static int ib555;
130  static int ib670;
131 
132  int32_t nbands = l2rec->l1rec->l1file->nbands;
133 
134  if (firstCall) {
135  firstCall = 0;
136  float* wave = l2rec->l1rec->l1file->fwave;
137  ib443 = windex(443., wave, nbands);
138  ib510 = windex(510., wave, nbands);
139  ib555 = windex(550., wave, nbands);
140  ib670 = windex(670., wave, nbands);
141  c1 = input->coccolith[0];
142  c2 = input->coccolith[1];
143  c3 = input->coccolith[2];
144  c4 = input->coccolith[3];
145  c5 = input->coccolith[4];
146  c6 = input->coccolith[5];
147  c7 = input->coccolith[6];
148  c8 = input->coccolith[7];
149  }
150 
151  float nLw443 = l2rec->nLw[ip * nbands + ib443];
152  float nLw510 = l2rec->nLw[ip * nbands + ib510];
153  float nLw555 = l2rec->nLw[ip * nbands + ib555];
154  float La670 = l2rec->La [ip * nbands + ib670];
155 
156  if (nLw443 >= c1 &&
157  nLw510 >= 0.0 &&
158  nLw555 >= c2 &&
159  La670 <= 1.1 &&
160  Between(c3, nLw443 / nLw555, c4) &&
161  Between(c5, nLw510 / nLw555, c6) &&
162  Between(c7, nLw443 / nLw510, c8))
163 
164  return (1);
165  else
166  return (0);
167 }
168 
169 char isTurbid(l2str *l2rec, int32_t ip) {
170 
171  static float Fo;
172  static int firstCall = 1;
173  static int ib670;
174 
175  if (firstCall) {
176  firstCall = 0;
177  ib670 = windex(670., l2rec->l1rec->l1file->fwave, l2rec->l1rec->l1file->nbands);
178  if (l1_input->outband_opt >= 2)
179  Fo = l2rec->l1rec->l1file->Fonom[ib670];
180  else
181  Fo = l2rec->l1rec->l1file->Fobar[ib670];
182  }
183 
184  if (l2rec->nLw[ip * l2rec->l1rec->l1file->nbands + ib670] / Fo > 0.0012)
185  return (1);
186  else
187  return (0);
188 }
189 
190 char isOptShallow(l2str *l2rec, int32_t ip) {
191  static float firstCall = 1;
192  static int ib443;
193  static int ib555;
194  static int ib670;
195 
196  float* wave = l2rec->l1rec->l1file->fwave;
197 
198  if (firstCall) {
199  firstCall = 0;
200  ib443 = windex(443., wave, l2rec->l1rec->l1file->nbands);
201  ib555 = windex(550., wave, l2rec->l1rec->l1file->nbands);
202  ib670 = windex(670., wave, l2rec->l1rec->l1file->nbands);
203  }
204 
205  //McKinna, L. & P.J. Werdell (2018). Approach for identifying optically
206  //shallow pixels when processing ocean-color imagery, Opt. Express, 26(22),
207  //A915-A928, doi: 10.1364/OE.26.00A915
208 
209  float a670, a555, u670, u555, bbp670, bbp555, quasiC555, od555;
210 
211  float depth = 0. - l2rec->l1rec->dem[ip];
212 
213  float Rrs443 = l2rec->Rrs[ip * l2rec->l1rec->l1file->nbands + ib443];
214  float Rrs555 = l2rec->Rrs[ip * l2rec->l1rec->l1file->nbands + ib555];
215  float Rrs670 = l2rec->Rrs[ip * l2rec->l1rec->l1file->nbands + ib670];
216 
217  float rrsSub443 = Rrs443 / (0.52 + 1.7 * Rrs443);
218  float rrsSub555 = Rrs555 / (0.52 + 1.7 * Rrs555);
219  float rrsSub670 = Rrs670 / (0.52 + 1.7 * Rrs670);
220 
221  float aw670 = l2rec->l1rec->sw_a[ib670];
222  float bbw670 = l2rec->l1rec->sw_bb[ib670];
223  float bbw555 = l2rec->l1rec->sw_bb[ib555];
224 
225  //Check valid Rrs values
226  if (Rrs443 < 0.0 || Rrs555 < 0.0 || Rrs670 < 0.0 ) {
227  return 0;
228  }
229 
230  //1. if water column greater than 40m, not flagged.
231  if (depth > 40.0) {
232  return 0;
233  }
234 
235  //2. if water column shallower than 5m, flagged.
236  if ( depth <= 5.0) {
237  return 1;
238  }
239 
240  /*Use QAA methodology adopted by Barnes et al (2013) with reference*/
241  /* wavelength set to 667nm*/
242 
243  a670 = aw670 + 0.07 * pow((rrsSub670 / rrsSub443), 1.1);
244  u670 = (-0.089 + pow(0.089*0.089 + 4.0 * 0.125 * rrsSub670, 0.5)) / (2.0 * 0.125);
245  u555 = (-0.089 + pow(0.089*0.089 + 4.0 * 0.125 * rrsSub555, 0.5)) / (2.0 * 0.125);
246 
247  bbp670 = ((u670 * a670) / (1.0 - u670)) - bbw670;
248 
249  bbp555 = bbp670 * pow((wave[ib670] / wave[ib555]), 1.03373);
250  a555 = ((1.0 - u555)* (bbp555 + bbw555)) / u555;
251 
252  /*Estimate c555, the beam attenuation coefficient*/
253  /*Assume the mean particulate backscatter ratio is 0.015*/
254  quasiC555 = a555 + (bbp555/0.015) + (bbw555/0.5);
255 
256  /*Estimate the water column's optical depth at 547 nm*/
257  od555 = quasiC555 * depth;
258 
259  /*Flag as optically shallow pixel if od(547) less than or equal to 20 */
260  /*At this threshold, we assume the seafloor has a non-negligible effect*/
261  /* on the the short wave spectral water-leaving reflectance signal*/
262 
263  if (od555 <= 20.0) {
264  return 1;
265  } else {
266 
267  return 0;
268  }
269 
270 }
271 
272 /* ---------------------------------------------------------------------------------------- */
273 /* aerindx() - compute aerosol index for absorbing aerosol test */
274 
275 /* ---------------------------------------------------------------------------------------- */
276 float aerindex(l2str *l2rec, int32_t ip) {
277  static int firstCall = 1;
278 
279  static int ib412;
280  static int ib510;
281 
282  static int32_t mask = ATMFAIL | LAND | CHLFAIL;
283  static int i510;
284  static double poly_coeff[7];
285  static double poly_coeff_412[7] = {0.73172938, -0.54565918, 0.20022312, -0.12036241, 0.11687968, -0.018732825, -0.0095574674};
286 
287  double poly_coeff_510[7] = {0.58543905, -0.013125745, -0.059568208, -0.016141849, 0.0035106655, 0.0012957265, 1.4235445e-05};
288  double poly_coeff_531[7] = {0.51483158, 0.15893415, -0.051696975, -0.055474007, -0.0029635773, 0.0053882411, 0.0010106600};
289  double poly_coeff_551[7] = {0.47507366, 0.25216739, -0.0096908094, -0.070882408, -0.012501495, 0.0061436085, 0.0015798359};
290  double poly_coeff_555[7] = {0.43681192, 0.26663018, 0.016592559, -0.068132662, -0.015470602, 0.0051694309, 0.0015132129};
291  //double poly_coeff_412_S[7] = {0.72884183, -0.54380414, 0.20225533, -0.12180914, 0.11508257, -0.017784535, -0.0095387145};
292  //double poly_coeff_412_M[7] = {0.73461692, -0.54751422, 0.19819091, -0.11891567, 0.11867679, -0.019681115, -0.0095762203};
293  float tLw_pred;
294  float Lt_pred;
295  float Lt_meas;
296  float mu0;
297  float index = 100.0, index510, index412;
298  int32_t ipb, i;
299  double logchl, lchl, nLw_510, nLw_412;
300 
301 
302  /* load parameters */
303  if (firstCall) {
304 
305  /* determine index of the bands for the absorbing aerosol analysis */
306  ib412 = windex(412., l2rec->l1rec->l1file->fwave, l2rec->l1rec->l1file->nbands);
307  ib510 = windex(510., l2rec->l1rec->l1file->fwave, l2rec->l1rec->l1file->nbands);
308  i510 = l2rec->l1rec->l1file->iwave[ib510];
309 
310  switch (i510) {
311  case 510:
312  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_510[i];
313  break;
314  case 531:
315  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_531[i];
316  break;
317  case 551:
318  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_551[i];
319  break;
320  case 555:
321  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_555[i];
322  break;
323  default:
324  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_510[i];
325  break;
326  }
327  /*
328  switch (l2rec->sensorID) {
329  case SEAWIFS:
330  for (i=0; i<7; i++) poly_coeff_412[i] = poly_coeff_412_S[i];
331  break;
332  case MODISA:
333  case HMODISA:
334  for (i=0; i<7; i++) poly_coeff_412[i] = poly_coeff_412_M[i];
335  break;
336  default:
337  break;
338  }
339  */
340  firstCall = 0;
341 
342  }
343 
344 
345  if ((l2rec->l1rec->flags[ip] & mask) != 0)
346  return (index);
347 
348  if (l2rec->chl[ip] < 0.1)
349  return (index);
350 
351 
352  logchl = lchl = log10(l2rec->chl[ip]);
353  nLw_510 = poly_coeff[0] + logchl * poly_coeff[1];
354  nLw_412 = poly_coeff_412[0] + logchl * poly_coeff_412[1];
355  for (i = 0; i < 5; i++) {
356  logchl *= lchl;
357  nLw_510 += poly_coeff[i + 2] * logchl;
358  nLw_412 += poly_coeff_412[i + 2] * logchl;
359  }
360 
361  /* Convert water-leaving radiances from band-centered to band-averaged. */
362  /*
363  if (input->outband_opt >= 2) {
364  nLw_510 /= l2rec->f_outband[ip*l2rec->nbands + ib510];
365  nLw_412 /= l2rec->f_outband[ip*l2rec->nbands + ib412];
366  }
367  */
368 
369 
370  /* cos of solz */
371  mu0 = cos(l2rec->l1rec->solz[ip] / RADEG);
372 
373 
374  /* bring water-leaving radiance at 510nm to the TOA */
375  ipb = ip * l2rec->l1rec->l1file->nbands + ib510;
376  tLw_pred = (float) nLw_510 / l2rec->brdf[ipb]
377  * l2rec->l1rec->t_sol[ipb]
378  * l2rec->l1rec->t_sen[ipb]
379  * l2rec->l1rec->tg_sol[ipb]
380  * l2rec->l1rec->tg_sen[ipb]
381  * l2rec->l1rec->polcor[ipb]
382  * l2rec->l1rec->t_o2[ipb]
383  * mu0
384  * l2rec->l1rec->fsol;
385 
386  /* calculate TOA predicted radiance */
387  Lt_pred = tLw_pred
388  + ((l2rec->l1rec->TLg[ipb]
389  + l2rec->La[ipb])
390  * l2rec->l1rec->t_o2[ipb]
391  + l2rec->l1rec->Lr[ipb]
392  + l2rec->l1rec->tLf[ipb])
393  * l2rec->l1rec->polcor[ipb]
394  * l2rec->l1rec->tg_sol[ipb]
395  * l2rec->l1rec->tg_sen[ipb];
396 
397  /* get measured TOA radiance */
398  Lt_meas = l2rec->l1rec->Lt[ipb];
399 
400  /* obtain percent difference between measured and predicted TOA radiance */
401  index510 = 100.0 * (Lt_meas - Lt_pred) / Lt_meas;
402 
403  if (index510 >= 0.0)
404  return (index);
405 
406 
407  /* bring water-leaving radiance to the TOA */
408  ipb = ip * l2rec->l1rec->l1file->nbands + ib412;
409  tLw_pred = (float) nLw_412 / l2rec->brdf[ipb]
410  * l2rec->l1rec->t_sol[ipb]
411  * l2rec->l1rec->t_sen[ipb]
412  * l2rec->l1rec->tg_sol[ipb]
413  * l2rec->l1rec->tg_sen[ipb]
414  * l2rec->l1rec->polcor[ipb]
415  * l2rec->l1rec->t_o2[ipb]
416  * mu0
417  * l2rec->l1rec->fsol;
418 
419  /* calculate predicted TOA radiance */
420  Lt_pred = tLw_pred
421  + ((l2rec->l1rec->TLg[ipb]
422  + l2rec->La[ipb])
423  * l2rec->l1rec->t_o2[ipb]
424  + l2rec->l1rec->Lr[ipb]
425  + l2rec->l1rec->tLf[ipb])
426  * l2rec->l1rec->polcor[ipb]
427  * l2rec->l1rec->tg_sol[ipb]
428  * l2rec->l1rec->tg_sen[ipb];
429 
430  /* get measured TOA radiance */
431  Lt_meas = l2rec->l1rec->Lt[ipb];
432 
433  /* obtain percent difference between measured and predicted TOA radiance */
434  index412 = 100.0 * (Lt_meas - Lt_pred) / Lt_meas + 2.0;
435 
436  return (index412);
437 
438 }
439 
440 
char isCoccolith(l2str *l2rec, int32_t ip)
Definition: setflags_l2.c:124
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
#define COCCOLITH
Definition: l2_flags.h:21
#define MAXAERITER
Definition: l2_flags.h:30
#define CHLFAIL
Definition: l2_flags.h:26
#define NULL
Definition: decode_rs.h:63
read l1rec
abs_aer_option
Definition: setflags_l2.c:8
instr * input
float aerindex(l2str *l2rec, int32_t ip)
Definition: setflags_l2.c:276
void setflagbits_l2(l2str *l2rec, int32_t ipix)
Definition: setflags_l2.c:14
#define TURBIDW
Definition: l2_flags.h:22
#define AERRHSWIR
Definition: l12_parms.h:33
@ RHOWINDEX
Definition: setflags_l2.c:9
l1_input_t * l1_input
Definition: l1_options.c:9
#define FATAL_ERROR
Definition: swl0_parms.h:5
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
#define ATMWARN
Definition: l2_flags.h:33
#define RADEG
Definition: czcs_ctl_pt.c:5
#define LAND
Definition: l2_flags.h:12
#define ATMFAIL
Definition: l2_flags.h:11
char isTurbid(l2str *l2rec, int32_t ip)
Definition: setflags_l2.c:169
char isOptShallow(l2str *l2rec, int32_t ip)
Definition: setflags_l2.c:190
#define BAD_FLT
Definition: jplaeriallib.h:19
#define CHLWARN
Definition: l2_flags.h:32
int32_t nbands
int32 spix
Definition: l1_czcs_hdf.c:21
#define Between(a, x, b)
Definition: setflags_l2.c:122
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
#define OPSHAL
Definition: l2_flags.h:34
float mu0
int32 epix
Definition: l1_czcs_hdf.c:23
#define CHL_MIN
Definition: l12_parms.h:42
@ TOTEXINCTION
Definition: setflags_l2.c:10
#define CHL_MAX
Definition: l12_parms.h:43
#define LOWLW
Definition: l2_flags.h:25
int i
Definition: decode_rs.h:71
#define AERWANGSWIR
Definition: l12_parms.h:31
int npix
Definition: get_cmp.c:28
#define ABSAER
Definition: l2_flags.h:28