OB.DAAC Logo
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) {
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("Cannot apply GMAO total extinction absorbing aerosol test without anc_aerosol inputs defined\n");
111  exit(FATAL_ERROR);
112  }
113  break;
114  }
115  }
116  }
117 }
118 
119 
120 #define Between(a,x,b)(x >= a && x <= b)
121 
122 char isCoccolith(l2str *l2rec, int32_t ip) {
123  static float firstCall = 1;
124  static float c1, c2, c3, c4, c5, c6, c7, c8;
125  static int ib443;
126  static int ib510;
127  static int ib555;
128  static int ib670;
129 
130  int32_t nbands = l2rec->l1rec->l1file->nbands;
131 
132  if (firstCall) {
133  firstCall = 0;
134  float* wave = l2rec->l1rec->l1file->fwave;
135  ib443 = windex(443., wave, nbands);
136  ib510 = windex(510., wave, nbands);
137  ib555 = windex(550., wave, nbands);
138  ib670 = windex(670., wave, nbands);
139  c1 = input->coccolith[0];
140  c2 = input->coccolith[1];
141  c3 = input->coccolith[2];
142  c4 = input->coccolith[3];
143  c5 = input->coccolith[4];
144  c6 = input->coccolith[5];
145  c7 = input->coccolith[6];
146  c8 = input->coccolith[7];
147  }
148 
149  float nLw443 = l2rec->nLw[ip * nbands + ib443];
150  float nLw510 = l2rec->nLw[ip * nbands + ib510];
151  float nLw555 = l2rec->nLw[ip * nbands + ib555];
152  float La670 = l2rec->La [ip * nbands + ib670];
153 
154  if (nLw443 >= c1 &&
155  nLw510 >= 0.0 &&
156  nLw555 >= c2 &&
157  La670 <= 1.1 &&
158  Between(c3, nLw443 / nLw555, c4) &&
159  Between(c5, nLw510 / nLw555, c6) &&
160  Between(c7, nLw443 / nLw510, c8))
161 
162  return (1);
163  else
164  return (0);
165 }
166 
167 char isTurbid(l2str *l2rec, int32_t ip) {
168 
169  static float Fo;
170  static int firstCall = 1;
171  static int ib670;
172 
173  if (firstCall) {
174  firstCall = 0;
175  ib670 = windex(670., l2rec->l1rec->l1file->fwave, l2rec->l1rec->l1file->nbands);
176  if (l1_input->outband_opt >= 2)
177  Fo = l2rec->l1rec->l1file->Fonom[ib670];
178  else
179  Fo = l2rec->l1rec->l1file->Fobar[ib670];
180  }
181 
182  if (l2rec->nLw[ip * l2rec->l1rec->l1file->nbands + ib670] / Fo > 0.0012)
183  return (1);
184  else
185  return (0);
186 }
187 
188 char isOptShallow(l2str *l2rec, int32_t ip) {
189  static float firstCall = 1;
190  static int ib443;
191  static int ib555;
192  static int ib670;
193 
194  float* wave = l2rec->l1rec->l1file->fwave;
195 
196  if (firstCall) {
197  firstCall = 0;
198  ib443 = windex(443., wave, l2rec->l1rec->l1file->nbands);
199  ib555 = windex(550., wave, l2rec->l1rec->l1file->nbands);
200  ib670 = windex(670., wave, l2rec->l1rec->l1file->nbands);
201  }
202 
203  //McKinna, L. & P.J. Werdell (2018). Approach for identifying optically
204  //shallow pixels when processing ocean-color imagery, Opt. Express, 26(22),
205  //A915-A928, doi: 10.1364/OE.26.00A915
206 
207  float a670, a555, u670, u555, bbp670, bbp555, quasiC555, od555;
208 
209  float depth = 0. - l2rec->l1rec->dem[ip];
210 
211  float Rrs443 = l2rec->Rrs[ip * l2rec->l1rec->l1file->nbands + ib443];
212  float Rrs555 = l2rec->Rrs[ip * l2rec->l1rec->l1file->nbands + ib555];
213  float Rrs670 = l2rec->Rrs[ip * l2rec->l1rec->l1file->nbands + ib670];
214 
215  float rrsSub443 = Rrs443 / (0.52 + 1.7 * Rrs443);
216  float rrsSub555 = Rrs555 / (0.52 + 1.7 * Rrs555);
217  float rrsSub670 = Rrs670 / (0.52 + 1.7 * Rrs670);
218 
219  float aw670 = l2rec->l1rec->sw_a[ib670];
220  float bbw670 = l2rec->l1rec->sw_bb[ib670];
221  float bbw555 = l2rec->l1rec->sw_bb[ib555];
222 
223  //Check valid Rrs values
224  if (Rrs443 < 0.0 || Rrs555 < 0.0 || Rrs670 < 0.0 ) {
225  return 0;
226  }
227 
228  //1. if water column greater than 40m, not flagged.
229  if (depth > 40.0) {
230  return 0;
231  }
232 
233  //2. if water column shallower than 5m, flagged.
234  if ( depth <= 5.0) {
235  return 1;
236  }
237 
238  /*Use QAA methodology adopted by Barnes et al (2013) with reference*/
239  /* wavelength set to 667nm*/
240 
241  a670 = aw670 + 0.07 * pow((rrsSub670 / rrsSub443), 1.1);
242  u670 = (-0.089 + pow(0.089*0.089 + 4.0 * 0.125 * rrsSub670, 0.5)) / (2.0 * 0.125);
243  u555 = (-0.089 + pow(0.089*0.089 + 4.0 * 0.125 * rrsSub555, 0.5)) / (2.0 * 0.125);
244 
245  bbp670 = ((u670 * a670) / (1.0 - u670)) - bbw670;
246 
247  bbp555 = bbp670 * pow((wave[ib670] / wave[ib555]), 1.03373);
248  a555 = ((1.0 - u555)* (bbp555 + bbw555)) / u555;
249 
250  /*Estimate c555, the beam attenuation coefficient*/
251  /*Assume the mean particulate backscatter ratio is 0.015*/
252  quasiC555 = a555 + (bbp555/0.015) + (bbw555/0.5);
253 
254  /*Estimate the water column's optical depth at 547 nm*/
255  od555 = quasiC555 * depth;
256 
257  /*Flag as optically shallow pixel if od(547) less than or equal to 20 */
258  /*At this threshold, we assume the seafloor has a non-negligible effect*/
259  /* on the the short wave spectral water-leaving reflectance signal*/
260 
261  if (od555 <= 20.0) {
262  return 1;
263  } else {
264 
265  return 0;
266  }
267 
268 }
269 
270 /* ---------------------------------------------------------------------------------------- */
271 /* aerindx() - compute aerosol index for absorbing aerosol test */
272 
273 /* ---------------------------------------------------------------------------------------- */
274 float aerindex(l2str *l2rec, int32_t ip) {
275  static int firstCall = 1;
276 
277  static int ib412;
278  static int ib510;
279 
280  static int32_t mask = ATMFAIL | LAND | CHLFAIL;
281  static int i510;
282  static double poly_coeff[7];
283  static double poly_coeff_412[7] = {0.73172938, -0.54565918, 0.20022312, -0.12036241, 0.11687968, -0.018732825, -0.0095574674};
284 
285  double poly_coeff_510[7] = {0.58543905, -0.013125745, -0.059568208, -0.016141849, 0.0035106655, 0.0012957265, 1.4235445e-05};
286  double poly_coeff_531[7] = {0.51483158, 0.15893415, -0.051696975, -0.055474007, -0.0029635773, 0.0053882411, 0.0010106600};
287  double poly_coeff_551[7] = {0.47507366, 0.25216739, -0.0096908094, -0.070882408, -0.012501495, 0.0061436085, 0.0015798359};
288  double poly_coeff_555[7] = {0.43681192, 0.26663018, 0.016592559, -0.068132662, -0.015470602, 0.0051694309, 0.0015132129};
289  //double poly_coeff_412_S[7] = {0.72884183, -0.54380414, 0.20225533, -0.12180914, 0.11508257, -0.017784535, -0.0095387145};
290  //double poly_coeff_412_M[7] = {0.73461692, -0.54751422, 0.19819091, -0.11891567, 0.11867679, -0.019681115, -0.0095762203};
291  float tLw_pred;
292  float Lt_pred;
293  float Lt_meas;
294  float mu0;
295  float index = 100.0, index510, index412;
296  int32_t ipb, i;
297  double logchl, lchl, nLw_510, nLw_412;
298 
299 
300  /* load parameters */
301  if (firstCall) {
302 
303  /* determine index of the bands for the absorbing aerosol analysis */
304  ib412 = windex(412., l2rec->l1rec->l1file->fwave, l2rec->l1rec->l1file->nbands);
305  ib510 = windex(510., l2rec->l1rec->l1file->fwave, l2rec->l1rec->l1file->nbands);
306  i510 = l2rec->l1rec->l1file->iwave[ib510];
307 
308  switch (i510) {
309  case 510:
310  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_510[i];
311  break;
312  case 531:
313  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_531[i];
314  break;
315  case 551:
316  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_551[i];
317  break;
318  case 555:
319  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_555[i];
320  break;
321  default:
322  for (i = 0; i < 7; i++) poly_coeff[i] = poly_coeff_510[i];
323  break;
324  }
325  /*
326  switch (l2rec->sensorID) {
327  case SEAWIFS:
328  for (i=0; i<7; i++) poly_coeff_412[i] = poly_coeff_412_S[i];
329  break;
330  case MODISA:
331  case HMODISA:
332  for (i=0; i<7; i++) poly_coeff_412[i] = poly_coeff_412_M[i];
333  break;
334  default:
335  break;
336  }
337  */
338  firstCall = 0;
339 
340  }
341 
342 
343  if ((l2rec->l1rec->flags[ip] & mask) != 0)
344  return (index);
345 
346  if (l2rec->chl[ip] < 0.1)
347  return (index);
348 
349 
350  logchl = lchl = log10(l2rec->chl[ip]);
351  nLw_510 = poly_coeff[0] + logchl * poly_coeff[1];
352  nLw_412 = poly_coeff_412[0] + logchl * poly_coeff_412[1];
353  for (i = 0; i < 5; i++) {
354  logchl *= lchl;
355  nLw_510 += poly_coeff[i + 2] * logchl;
356  nLw_412 += poly_coeff_412[i + 2] * logchl;
357  }
358 
359  /* Convert water-leaving radiances from band-centered to band-averaged. */
360  /*
361  if (input->outband_opt >= 2) {
362  nLw_510 /= l2rec->f_outband[ip*l2rec->nbands + ib510];
363  nLw_412 /= l2rec->f_outband[ip*l2rec->nbands + ib412];
364  }
365  */
366 
367 
368  /* cos of solz */
369  mu0 = cos(l2rec->l1rec->solz[ip] / RADEG);
370 
371 
372  /* bring water-leaving radiance at 510nm to the TOA */
373  ipb = ip * l2rec->l1rec->l1file->nbands + ib510;
374  tLw_pred = (float) nLw_510 / l2rec->brdf[ipb]
375  * l2rec->l1rec->t_sol[ipb]
376  * l2rec->l1rec->t_sen[ipb]
377  * l2rec->l1rec->tg_sol[ipb]
378  * l2rec->l1rec->tg_sen[ipb]
379  * l2rec->l1rec->polcor[ipb]
380  * l2rec->l1rec->t_o2[ipb]
381  * mu0
382  * l2rec->l1rec->fsol;
383 
384  /* calculate TOA predicted radiance */
385  Lt_pred = tLw_pred
386  + ((l2rec->l1rec->TLg[ipb]
387  + l2rec->La[ipb])
388  * l2rec->l1rec->t_o2[ipb]
389  + l2rec->l1rec->Lr[ipb]
390  + l2rec->l1rec->tLf[ipb])
391  * l2rec->l1rec->polcor[ipb]
392  * l2rec->l1rec->tg_sol[ipb]
393  * l2rec->l1rec->tg_sen[ipb];
394 
395  /* get measured TOA radiance */
396  Lt_meas = l2rec->l1rec->Lt[ipb];
397 
398  /* obtain percent difference between measured and predicted TOA radiance */
399  index510 = 100.0 * (Lt_meas - Lt_pred) / Lt_meas;
400 
401  if (index510 >= 0.0)
402  return (index);
403 
404 
405  /* bring water-leaving radiance to the TOA */
406  ipb = ip * l2rec->l1rec->l1file->nbands + ib412;
407  tLw_pred = (float) nLw_412 / l2rec->brdf[ipb]
408  * l2rec->l1rec->t_sol[ipb]
409  * l2rec->l1rec->t_sen[ipb]
410  * l2rec->l1rec->tg_sol[ipb]
411  * l2rec->l1rec->tg_sen[ipb]
412  * l2rec->l1rec->polcor[ipb]
413  * l2rec->l1rec->t_o2[ipb]
414  * mu0
415  * l2rec->l1rec->fsol;
416 
417  /* calculate predicted TOA radiance */
418  Lt_pred = tLw_pred
419  + ((l2rec->l1rec->TLg[ipb]
420  + l2rec->La[ipb])
421  * l2rec->l1rec->t_o2[ipb]
422  + l2rec->l1rec->Lr[ipb]
423  + l2rec->l1rec->tLf[ipb])
424  * l2rec->l1rec->polcor[ipb]
425  * l2rec->l1rec->tg_sol[ipb]
426  * l2rec->l1rec->tg_sen[ipb];
427 
428  /* get measured TOA radiance */
429  Lt_meas = l2rec->l1rec->Lt[ipb];
430 
431  /* obtain percent difference between measured and predicted TOA radiance */
432  index412 = 100.0 * (Lt_meas - Lt_pred) / Lt_meas + 2.0;
433 
434  return (index412);
435 
436 }
437 
438 
char isCoccolith(l2str *l2rec, int32_t ip)
Definition: setflags_l2.c:122
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:274
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:167
char isOptShallow(l2str *l2rec, int32_t ip)
Definition: setflags_l2.c:188
#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:120
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:27
#define ABSAER
Definition: l2_flags.h:28