OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
gac_st.c
Go to the documentation of this file.
1 /*-----------------------------------------------------------------------------
2  Function: stray_light_gac
3 
4  Returns: int32 (status)
5  Returns as 0, "not done", to indicate that the routine did not
6  return any useful information for the last call, but was filling
7  its buffers with scan line data; and 1, "done", to indicate that
8  the l1b_data and sl_flag arrays for scan line sl_scan have been
9  set.
10 
11  Description:
12  The function stray_light_gac processes GAC scan lines by detecting
13  bright targets (BTs) and labelling pixels as being within the BTs
14  or adjacent to them where they are likely to be contaminated by
15  stray light (ST) from the BTs. The Level-1B data for the ST pixels
16  in the along-scan direction are corrected for ST contamination
17  according to the method in Barnes et al. (1995). The routine
18  processes one scan line at a time using a rolling window for along-
19  track processing of three scan lines for GAC data. If along-track
20  processing is indicated (nscan < 3), each scan line is processed
21  individually without the rolling window.
22 
23  Parameters: (in calling order)
24  Type Name I/O Description
25  ---- ---- --- -----------
26  int32 *initial I Flag, should be set to "true" (a value
27  of 1) if it is the fst call for a scene
28  float32 Ltyp_frac I fraction of Ltypical for band 8
29  float32 Styp_frac I fraction value that will be applied
30  to knee value to calculate straylight
31  threshold
32  int32 nscans I the number of scan lines in the scene
33  (to be processed). Must be set to <5
34  to process only in the along scan
35  int32 nsamples I number of pixel data values in
36  l1b_data and sl_flag
37  int32 scan_no I scan line number of l1b_data input
38  int16 gn I Band 8 gain of scan line scan_no
39  float32 *rads I Radiances of calibration knees[8,4,5]
40  float32 *l1b_data I/O array of nsamples long containing
41  l1b_data of scan line scan_no
42  int32 *sl_scan O scan line number for which the
43  returned l1b_data apply
44  int32 *sl_flag O array of nsamples long containing
45  codes identifying BT pixels and the
46  proximity of others to BTs.
47  Notes:
48 
49  Modification history:
50  Programmer Organization Date Description of change
51  -------------- ------------ -------- ---------------------
52  W. Robinson SAIC 01/28/02 enable correction for
53  2nd pixel from BT and
54  fix problem in flagging
55  and correcting
56  Joel Gales Futuretech 10/16/00 1) Set AT in flag_buf
57  when initially setting
58  BT.
59  2) Set left edge found
60  if first pixel in scan
61  is above knee_val
62  Lakshmi Kumar Hughes STX 02/12/97 Added Styp_frac to input
63  parameter list
64  Removed PROTOTYPE defns.
65  corrected algorithm that
66  finds edges
67  Lakshmi Kumar Hughes STX 01/26/96
68  Lakshmi Kumar Hughes STX 10/12/95 Original development
69 ------------------------------------------------------------------------------*/
70 
71 #include "st_lt.h"
72 #include <mfhdf.h>
73 
74 int32_t stray_light_gac(int32_t *initial, float Ltyp_frac, float Styp_frac,
75  int32_t nscans, int32_t nsamples, int32_t scan_no, int16_t gn,
76  float *rads, float *l1b_data, int32_t *sl_scan, int32_t *sl_flag) {
77 
78  /*-----------------------------------------------------------------------------
79  ** i -- dummy index
80  ** n -- pixel pointer on the scan line being processed
81  ** p -- dummy index
82  ** band -- dummy index
83  ** status -- returns as 0 (NOTDONE) to indicate that the routine did not
84  return any useful information for the last call,
85  but was filling its buffers with scan line data;
86  and 1 (DONE) to indicate that the l1b_data and sl_flag
87  arrays for scan line sl_scan have been set
88  ** edge -- flag that indicates whether an edge was found or not
89  ** scans_done -- number of scans processed (if nscans < 5, scans_done = 2
90  ** scan -- pointer to the scan line being processed in the buffers--
91  the central scan line in a rolling window of five
92  consecutive scan lines
93  ** scan_p1 -- pointer to the scan line following scan in the buffers
94  (scan plus 1)
95  ** prev_scans -- the number of scans preceding the one represented by scan
96  ** next_scans -- the number of scans following the one represented by scan
97  ** line_no -- a 5-integer array of the scan line numbers of the scan
98  lines in the buffers
99  ** flag_buf -- a 2-dimensional buffer (5 x nsamples) for marking BT and
100  SL codes for five consecutive scan lines
101  ** scan_m1 -- pointer to the scan line preceding scan lin the buffers
102  (scan minus 1)
103  ** Ltypical_8 -- Ltypical value for band 8; = 1.09
104  ** Ltyp_thresh -- value used to test if the band-8 radiance difference from
105  one pixel to the next indicates the edge
106  ** l1b_buf -- a 3-dimensional buffer (8 x 5 x nsamples) for storing five
107  consecutive scan lines of Level-1B data for 8 bands
108  ** GAC_RRANGE -- Maximum number of pixels to the right of a BT that are
109  affected by ST from that BT; = 3
110  ** GAC_LRANGE -- Maximum number of pixels to the left (along scan) of a
111  BT that are affected by ST from that BT; = 3
112  ** K -- BT response constants for along-scan stray light correction
113  array of two dimensions: 8 bands by GAC_LRANGE+GAC_RRANGE+1;
114  values are from Table 8 of Barnes et al. (1995),
115  except for value of central pixel which equals 1 minus
116  table value
117  ----------------------------------------------------------------------------*/
118 
119  int32 i, n, p, farthest, status, pix;
120  float32 delta, Ltyp_thresh, Ltypical_8 = 1.09;
121  float32 knee_val = 0, stray_thresh = 0;
122  int32 edge = 0, right_edge = -1, left_edge = -1;
123  static int32 scans_done, scan, scan_p1, prev_scans, next_scans;
124  static int32 scan_m1, rt_edge_found = FALSE;
125  static int32 line_no[GACLINES], flag_buf[GACLINES][GACSAMPS];
126  static float32 l1b_buf[MAXBANDS][GACLINES][GACSAMPS];
127  float32 Ctyp_frac = 1.25, Ctyp_thresh;
128 
129  /* If first call, initialize the counters and pointers and fill buffer
130  ** with first scan line. If along-track processing is requested (nscans > 2)
131  ** then return for more scans to fill the buffers; otherwise, set
132  ** scans_done = 2.
133  */
134 
135  if (*initial == TRUE) {
136  scans_done = 0;
137  scan = 1;
138  scan_p1 = 0;
139  prev_scans = 0;
140  next_scans = 0;
141  line_no[scan - 1] = scan_no;
142  for (i = 0; i < MAXBANDS; i++)
143  memcpy(l1b_buf[i][scan - 1], &l1b_data[i * nsamples],
144  (sizeof (float)*nsamples));
145  for (i = 0; i < nsamples; i++)
146  flag_buf[scan - 1][i] = BLANK;
147  status = NOTDONE;
148  if (nscans > GACLINES - 1) {
149  *initial = FALSE;
150  return status;
151  }
152  }
153 
154  if (*initial == FALSE) {
155  status = NOTDONE;
156  /* If only 1 scan has been obtained, get the 2nd and return for more */
157  if (scan_p1 == 0) {
158  scan_p1 = 2;
159  next_scans = 1;
160  line_no[scan_p1 - 1] = scan_no;
161  for (i = 0; i < MAXBANDS; i++)
162  memcpy((float *) l1b_buf[i][scan_p1 - 1],
163  &l1b_data[i * nsamples], sizeof (float)*nsamples);
164  for (i = 0; i < nsamples; i++)
165  flag_buf[scan_p1 - 1][i] = BLANK;
166  }
167 
168  /* If only 2 scans have been obtained, rotate the array pointers
169  ** (since the 1st scan line pointed to by scan has been processed),
170  ** get the 3rd scan line, and process scan = 2.
171  */
172  if (scans_done == 1) {
173  scan_m1 = 1;
174  scan = 2;
175  scan_p1 = 3;
176  prev_scans = 1;
177  next_scans = 1;
178  line_no[scan_p1 - 1] = scan_no;
179  for (i = 0; i < MAXBANDS; i++)
180  memcpy((float *) l1b_buf[i][scan_p1 - 1],
181  &l1b_data[i * nsamples], sizeof (float)*nsamples);
182  for (i = 0; i < nsamples; i++)
183  flag_buf[scan_p1 - 1][i] = BLANK;
184  }
185 
186  /* If the buffers are full (as happens when scans_done > 1), rotate
187  ** the array pointers */
188  if (scans_done > 1) {
189  scan_m1 = (scan_m1) % 3 + 1;
190  scan = (scan) % 3 + 1;
191  scan_p1 = (scan_p1) % 3 + 1;
192 
193  /* If not all scan lines have been obtained, get the next one into
194  ** scan_p1 buffers */
195  if (scans_done < nscans - 1) {
196  line_no[scan_p1 - 1] = scan_no;
197  for (i = 0; i < MAXBANDS; i++)
198  memcpy((float *) l1b_buf[i][scan_p1 - 1],
199  &l1b_data[i * nsamples], sizeof (float)*nsamples);
200  for (i = 0; i < nsamples; i++)
201  flag_buf[scan_p1 - 1][i] = BLANK;
202  }
203  /* Otherwise, finish ones that haven't been processed. If all have
204  ** been processed (next_scans < 0), flush out those remaining */
205  else {
206  next_scans = next_scans - 1;
207  if (next_scans < 0) {
208  *sl_scan = line_no[scan_m1 - 1];
209  for (i = 0; i < MAXBANDS; i++)
210  memcpy((float *) &l1b_data[i * nsamples],
211  l1b_buf[i][scan_m1 - 1], sizeof (float)*nsamples);
212  for (i = 0; i < nsamples; i++)
213  sl_flag[i] = flag_buf[scan_m1 - 1][i];
214  status = DONE;
215  return status;
216  } /* end if */
217  } /* end if-else */
218  } /* end if */
219  } /* end if */
220 
221  if (*initial == TRUE) {
222  scans_done = 2;
223  scan_m1 = scan;
224  }
225 
226 
227  /* setup clear water threshold value Ctyp_thresh */
228  Ctyp_thresh = Ctyp_frac * Ltypical_8;
229 
230  /* if input radiance is greater than its knee, mark that pix as BT */
231  knee_val = rads[7 * GAINS * KNEES + gn * KNEES + 1];
232  stray_thresh = Styp_frac * knee_val;
233 
234  for (pix = 0; pix < nsamples; pix++) {
235  if (l1b_buf[7][scan - 1][pix] > knee_val) {
236  flag_buf[scan - 1][pix] = BT;
237  if ((prev_scans >= 1) && (flag_buf[scan_m1 - 1][pix]) != BT)
238  flag_buf[scan_m1 - 1][pix] = AT;
239  if ((next_scans >= 1) && (flag_buf[scan_p1 - 1][pix] != BT))
240  flag_buf[scan_p1 - 1][pix] = AT;
241  }
242  } /* end for */
243 
244 
245 
246  /* Start loop to process the scan line pointed to by scan */
247  rt_edge_found = FALSE;
248  n = 0;
249  while (n < nsamples - 1) {
250  edge = NO;
251 
252  if (n == 0 && flag_buf[scan - 1][n] == BT) {
253  edge = LEFT;
254  rt_edge_found = FALSE;
255  }
256 
257 
258  while (edge == NO && (n < nsamples - 1)) {
259 
260  if (flag_buf[scan - 1][n] != BT || flag_buf[scan - 1][n + 1] != BT) {
261  delta = l1b_buf[7][scan - 1][n + 1] - l1b_buf[7][scan - 1][n];
262  if (delta > 0)
263  Ltyp_thresh =
264  max((Ltyp_frac * (l1b_buf[7][scan - 1][n + 1] - Ltypical_8)),
265  Ltyp_frac * Ltypical_8);
266  else
267  Ltyp_thresh =
268  max((Ltyp_frac * (l1b_buf[7][scan - 1][n] - Ltypical_8)),
269  Ltyp_frac * Ltypical_8);
270  if (delta > Ltyp_thresh && l1b_buf[7][scan - 1][n + 1] > stray_thresh) {
271  edge = LEFT;
272  rt_edge_found = FALSE;
273  }
274  if (-1 * delta > Ltyp_thresh && l1b_buf[7][scan - 1][n] > stray_thresh) {
275  if (!(rt_edge_found)) {
276  edge = RIGHT;
277  rt_edge_found = TRUE;
278  }
279  }
280 
281  }
282  n++;
283  }
284 
285  if (edge != NO) {
286  n--;
287 
288  /* If right edge is found first, set the scan-line start as the left
289  ** edge, and correct the ST pixels which are 1 and two pixels away
290  ** from the right edge.
291  ** CF_right = along scan stray-light GAC correctionn factors for
292  ** pixels to the right of a BT; array of two dimensions: 8 bands by
293  ** GAC_RRANGE; values are from Table 15 of Barnes et al. (1995).
294  ** The above table "CF_right" is defined in cf.h file */
295 
296  if (edge == RIGHT) {
297  right_edge = n;
298  for (i = n; i >= 0; i--) if (l1b_buf[7][scan - 1][i] < Ctyp_thresh) break;
299  if (i < 0) i = 0;
300  left_edge = i;
301  p = right_edge + 1;
302  farthest = min((right_edge + GAC_RRANGE), (nsamples - 1));
303  while (p <= farthest) {
304  if (flag_buf[scan - 1][p] == BT)
305  p = farthest + 1;
306  else {
307  if (flag_buf[scan - 1][p] == BLANK) {
308  flag_buf[scan - 1][p] = p - right_edge; /*flag*/
309  if (p == (right_edge + GAC_RRANGE)) {
310  for (i = 0; i < MAXBANDS; i++)
311  l1b_buf[i][scan - 1][p] += /*correction*/
312  CF_right[i][2] * l1b_buf[i][scan - 1][right_edge];
313  } else if (p == (right_edge + GAC_RRANGE - 1)) {
314  for (i = 0; i < MAXBANDS; i++)
315  l1b_buf[i][scan - 1][p] += /*correction*/
316  CF_right[i][1] * l1b_buf[i][scan - 1][right_edge];
317  }
318  }
319  } /* end else */
320  p += 1;
321  } /* end while */
322  }/* end if */
323  else {
324  /* if a left edge is found, correct the 2nd and third pixels
325  ** prior to the left edge.
326  ** CF_left = along-scan stray-light GAC correction factors for
327  ** pixels to the left of a BT; array of two dimensions: 8 bands by
328  ** Lrange; values are from Table 15 of Barnes et al. (1995).
329  ** The table "CF_left" is defined in cf.h file */
330 
331  left_edge = n + 1;
332  p = n;
333  farthest = max(left_edge - GAC_LRANGE, 0);
334  while (p >= farthest) {
335  if (flag_buf[scan - 1][p] == BT)
336  p = farthest - 1;
337  else {
338  if (flag_buf[scan - 1][p] == BLANK) {
339  flag_buf[scan - 1][p] = left_edge - p; /* flag */
340  if (p == (left_edge - GAC_LRANGE)) {
341  for (i = 0; i < MAXBANDS; i++)
342  l1b_buf[i][scan - 1][p] += /* correction*/
343  (CF_left[i][2] * l1b_buf[i][scan - 1][left_edge]);
344  } else if (p == (left_edge - GAC_LRANGE + 1)) {
345  for (i = 0; i < MAXBANDS; i++)
346  l1b_buf[i][scan - 1][p] += /* correction*/
347  (CF_left[i][1] * l1b_buf[i][scan - 1][left_edge]);
348  }
349  } else if (flag_buf[scan - 1][p] > 0) {
350  flag_buf[scan - 1][p] += 1000 * (left_edge - p); /* flag */
351  if (p == (left_edge - GAC_LRANGE)) {
352  for (i = 0; i < MAXBANDS; i++)
353  l1b_buf[i][scan - 1][p] += /* correction*/
354  (CF_left[i][2] * l1b_buf[i][scan - 1][left_edge]);
355  } else if (p == (left_edge - GAC_LRANGE + 1)) {
356  for (i = 0; i < MAXBANDS; i++)
357  l1b_buf[i][scan - 1][p] += /* correction*/
358  (CF_left[i][1] * l1b_buf[i][scan - 1][left_edge]);
359  }
360  }
361  } /* end else */
362  p = p - 1;
363  } /* end while */
364 
365  /* search for the associated right edge */
366  right_edge = nsamples - 1;
367  n = n + 1;
368  while ((n < nsamples - 1) && (edge == LEFT)) {
369  if (flag_buf[scan - 1][n] != BT || flag_buf[scan - 1][n + 1] != BT) {
370  if (l1b_buf[7][scan - 1][n + 1] < Ctyp_thresh) goto R_found;
371  delta = l1b_buf[7][scan - 1][n + 1] - l1b_buf[7][scan - 1][n];
372  if (delta > 0)
373  Ltyp_thresh =
374  max((Ltyp_frac * (l1b_buf[7][scan - 1][n + 1] - Ltypical_8)),
375  Ltyp_frac * Ltypical_8);
376  else
377  Ltyp_thresh =
378  max((Ltyp_frac * (l1b_buf[7][scan - 1][n] - Ltypical_8)),
379  Ltyp_frac * Ltypical_8);
380  if (-1 * delta > Ltyp_thresh &&
381  l1b_buf[7][scan - 1][n] > stray_thresh) {
382  R_found:
383  edge = RIGHT;
384  rt_edge_found = TRUE;
385  }
386  } /*if BT */
387  n++;
388  } /*while ((n < nsamples - 1) && (edge == LEFT) */
389 
390  if (edge == LEFT) /* Cound not find right_edge */
391  rt_edge_found = TRUE; /* setting the last pix as rt edge*/
392  /* When the right edge is found, apply the correction to the
393  ** 2nd, third pixels down from the right edge */
394  /* WDR fix flag_buf check to account for corrections from
395  previous right edges */
396  if (edge == RIGHT) {
397  right_edge = n - 1;
398  p = right_edge + 1;
399  farthest = min((right_edge + GAC_RRANGE), (nsamples - 1));
400  while (p <= farthest) {
401  if (flag_buf[scan - 1][p] == BT)
402  p = farthest + 1;
403  else {
404  if ((flag_buf[scan - 1][p] == BLANK) ||
405  (flag_buf[scan - 1][p] > 0)) {
406  flag_buf[scan - 1][p] = p - right_edge; /* flag */
407  if (p == (right_edge + GAC_RRANGE)) {
408  for (i = 0; i < MAXBANDS; i++) /* correction */
409  l1b_buf[i][scan - 1][p] +=
410  (CF_right[i][2] * l1b_buf[i][scan - 1][right_edge]);
411  } else if (p == (right_edge + GAC_RRANGE - 1)) {
412  for (i = 0; i < MAXBANDS; i++) /* correction */
413  l1b_buf[i][scan - 1][p] +=
414  (CF_right[i][1] * l1b_buf[i][scan - 1][right_edge]);
415  }
416  }
417  } /* end else */
418  p += 1;
419  } /* end while */
420  } /* end if */
421  } /* end else */
422  /* Label all pixels between both edges as BT pixels. Some of the
423  ** pixels might already be BT, will get set again to avoid another
424  ** if check */
425 
426  for (p = left_edge; p <= right_edge; p++) {
427  flag_buf[scan - 1][p] = BT;
428  if ((prev_scans >= 1) && (flag_buf[scan_m1 - 1][p]) != BT)
429  flag_buf[scan_m1 - 1][p] = AT;
430  if ((next_scans >= 1) && (flag_buf[scan_p1 - 1][p] != BT))
431  flag_buf[scan_p1 - 1][p] = AT;
432  } /* end for */
433  } /* end if */
434 
435  if (rt_edge_found && n < nsamples - 1)
436  n = right_edge + 1;
437 
438  } /* end while */
439 
440 
441  /* If three scan lines have been obtained (and one has been processed),
442  ** return the earliest scan line (scan_m1) in the window. */
443 
444  scans_done = scans_done + 1;
445  if (scans_done > 1) {
446  *sl_scan = line_no[scan_m1 - 1];
447 
448  for (i = 0; i < MAXBANDS; i++)
449  memcpy(&l1b_data[i * nsamples],
450  l1b_buf[i][scan_m1 - 1], (sizeof (float)*nsamples));
451  for (i = 0; i < nsamples; i++)
452  sl_flag[i] = flag_buf[scan_m1 - 1][i];
453 
454  status = DONE;
455  } /* end if */
456  return status;
457 }
458 
float rads[BANDS_DIMS_1A][GAINS_DIMS_1A][KNEES_DIMS_1A]
Definition: l1a_seawifs.c:47
#define BLANK
Definition: st_lt.h:13
int status
Definition: l1_czcs_hdf.c:32
#define NO
Definition: l1.h:44
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
#define GAINS
Definition: calib_get_cal.h:10
#define FALSE
Definition: rice.h:164
MOD_PR01 Production producing one five minute granule of output data in each run It can be configured to produce as many as three five minute granules per run Each execution with one construction record and one date file for each dataset In normal these are created by which splits them out of the hour datasets For LANCE they are created by which merges all session MODIS L0 datasets overlapping the requested time and extracts from the merged data those packets which fall within that time period Each scan of data is stored in the L1A granule that covers the start time of that scan
Definition: MOD_PR01_pr.txt:19
#define TRUE
Definition: rice.h:165
#define GACSAMPS
Definition: st_lt.h:29
#define NOTDONE
Definition: l1a_seawifs.c:19
int32_t stray_light_gac(int32_t *initial, float Ltyp_frac, float Styp_frac, int32_t nscans, int32_t nsamples, int32_t scan_no, int16_t gn, float *rads, float *l1b_data, int32_t *sl_scan, int32_t *sl_flag)
Definition: gac_st.c:74
#define GAC_RRANGE
Definition: cf.h:8
#define AT
Definition: st_lt.h:11
#define KNEES
Definition: st_lt.h:32
float Ltyp_frac
Definition: l1a_seawifs.c:101
#define GACLINES
Definition: st_lt.h:27
const double delta
#define RIGHT
Definition: st_lt.h:15
#define MAXBANDS
Definition: l1_viirs_h5.c:23
#define BT
Definition: st_lt.h:10
int i
Definition: decode_rs.h:71
#define DONE
Definition: st_lt.h:24
#define GAC_LRANGE
Definition: cf.h:7
#define LEFT
Definition: st_lt.h:16
float p[MODELMAX]
Definition: atrem_corl1.h:131
l2prod max