OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
viirs_oxt.c
Go to the documentation of this file.
1 #include "viirs_sim_sdr.h"
2 #include <stdio.h>
3 #include <string.h>
4 #include <stdlib.h>
5 #include <ctype.h>
6 #define N_INT_RNG 16 /* # of inter-band ranges */
7 static float *ib_field[N_INT_RNG];
8 
9 int viirs_oxt(ctl_struc *ctl, in_rec_struc *in_rec)
10 /*-----------------------------------------------------------------------------
11  Program: viirs_oxt.c
12 
13  Description: Add the VIIRS optical crosstalk artifact to the simulated
14  radiance field
15 
16  Arguments:
17  Type Name I/O Description
18  ---- ---- --- -----------
19  ctl_struc * ctl I input controls
20  in_rec_struc * in_rec I/O controls for input record reading
21 
22  Modification history:
23 
24  W. Robinson, SAIC 26 May 2010 Original development
25 
26 ----------------------------------------------------------------------------*/ {
27  static int oxt_init = 0, n_vis_asg = 7;
28  static int vis_asg[] = {1, 3, 5, 7, 10, 12, 14}, ib_status[N_INT_RNG];
29  static int int_low_loc[N_INT_RNG], int_hi_loc[N_INT_RNG];
30  static float cw_wav[N_INT_RNG], av_rad[N_INT_RNG];
31  int irng, last_ib2, last_ib, iend, ib_targ, ib_ct, ib_pair[2];
32  int in_loc, nval;
33  int loc1, loc2, ilin, ipix;
34  float rad0, rad1, rad2, lam0, lam1, lam2, r_bnd1, r_bnd2, rfill;
35 
36  /*
37  * do initialization
38  */
39  if (oxt_init == 0) {
40  oxt_init = 1;
41 
42  printf("OXT parms:\n");
43  printf("oxt_mode = %d\n", ctl->oxt_mode);
44  printf("oxt_coef = %s\n", ctl->oxt_coef);
45  printf("inter_band = %s\n", ctl->inter_band);
46  /*
47  * read in the inter-band TOA radiance information
48  */
49  if (viirs_oxt_ib_read(ctl->inter_band, cw_wav, av_rad) != 0)
50  return 1;
51  /*
52  * set up the inter-band storage area
53  */
54  for (irng = 0; irng < N_INT_RNG; irng++) {
55  if ((ib_field[irng] = (float *)
56  malloc(in_rec->npix * NDET * sizeof (float))) == NULL) {
57  printf("%s, %d: unable to allocate interband field storage\n",
58  __FILE__, __LINE__);
59  return 1;
60  }
61  }
62  /*
63  * set up the inter-band status (0-unfilled, 1 filled) and the low and
64  * high bands to use in the interpolation
65  */
66  for (irng = 0; irng < N_INT_RNG; irng++)
67  ib_status[irng] = 0;
68  for (irng = 0; irng < n_vis_asg; irng++)
69  ib_status[ vis_asg[irng] ] = 1;
70  /*
71  * last_ib and last_ib2 will retain locations of up to 2 of the previous
72  * inter-bands that are filled and are updated as loop is run. ib_pair
73  * will contain up to 2 next inter-bands that are filled. These are used
74  * to set the best inter-bands to interpolate with in filling the
75  * empty inter-bands
76  */
77  last_ib2 = -1;
78  last_ib = -1;
79  for (irng = 0; irng < N_INT_RNG; irng++) {
80  if (ib_status[irng] == 0) {
81  /*
82  * get the next filled bands
83  */
84  iend = 0;
85  ib_targ = irng;
86  ib_ct = 0;
87  while (iend != 1) {
88  ib_targ++;
89  if (ib_targ == N_INT_RNG)
90  iend = 1;
91  else
92  if (ib_status[ib_targ] == 1) ib_pair[ib_ct++] = ib_targ;
93  if (ib_ct == 2) iend = 1;
94  }
95  /*
96  * use the previously and next filled bands to get best
97  * interpolation pair
98  */
99  if ((last_ib != -1) && (ib_ct >= 1)) { /* 2 surrounding locations, next lower and higher */
100  ib_pair[1] = ib_pair[0];
101  ib_pair[0] = last_ib;
102  } else if ((ib_ct == 0) && (last_ib2 != -1)) { /* previous 2 locs */
103  ib_pair[0] = last_ib2;
104  ib_pair[1] = last_ib;
105  } else if ((ib_ct == 2) && (last_ib == -1)) { /* next 2 locations they are already in ib_pair */
106  ;
107  } else { /* insufficient points to do the fill - should not happen */
108  printf("%s, %d: insufficient data to interpolate\n",
109  __FILE__, __LINE__);
110  return 1;
111  }
112  int_low_loc[irng] = ib_pair[0];
113  int_hi_loc[irng] = ib_pair[1];
114  } else { /* filled band found, maintain previous filled band locations */
115  last_ib2 = last_ib;
116  last_ib = irng;
117  }
118  }
119  }
120  /*
121  * insert the band radiances into the inter-band storage
122  * in_loc is any offset due to a line margin which te storage dosen't need
123  */
124  in_loc = in_rec->margin[0] * in_rec->npix;
125  nval = sizeof (float) * in_rec->npix * NDET;
126  for (irng = 0; irng < n_vis_asg; irng++) {
127  memcpy(ib_field[vis_asg[irng]], in_rec->bnd_lt[irng] + in_loc, nval);
128  }
129  /*
130  * fill the empty inter-band values with an interpolation of surrounding
131  * radiances, conditioned by the expected TOA spectrum
132  */
133  for (ilin = 0; ilin < NDET; ilin++)
134  for (ipix = 0; ipix < in_rec->npix; ipix++)
135  for (irng = 0; irng < N_INT_RNG; irng++) {
136  if (ib_status[irng] == 0) {
137  loc1 = int_low_loc[irng];
138  loc2 = int_hi_loc[irng];
139  lam1 = cw_wav[loc1];
140  lam2 = cw_wav[loc2];
141  lam0 = cw_wav[irng];
142  rad1 = av_rad[loc1];
143  rad2 = av_rad[loc2];
144  rad0 = av_rad[irng];
145  r_bnd1 = *(ib_field[loc1] + ipix + ilin * in_rec->npix);
146  r_bnd2 = *(ib_field[loc2] + ipix + ilin * in_rec->npix);
147  rfill = rad0 / (lam2 - lam1) *
148  (r_bnd1 * (lam2 - lam0) / rad1 +
149  r_bnd2 * (lam0 - lam1) / rad2);
150  *(ib_field[irng] + ipix + ilin * in_rec->npix) = rfill;
151  }
152  }
153  /*
154  * Get and apply the coefficients to proper detector in the inter-band field
155  * and modify the radiances with it
156  */
157  if (viirs_oxt_comp(ctl, in_rec) != 0)
158  return 1;
159  /*
160  printf( "PREMATURE END\n" );
161  return 1;
162  */
163  return 0;
164 }
165 
166 int viirs_oxt_ib_read(char *ib_file, float *cw_wav, float *av_rad)
167 /*-----------------------------------------------------------------------------
168  routine: viirs_oxt_read
169 
170  Description: set up the information needed for the inter-band radiance
171  determination
172 
173  Arguments:
174  Type Name I/O Description
175  ---- ---- --- -----------
176  char * ib_file I file containing the inter-band radiance
177  information for typical TOA radiances
178  in VIIRS iinterband wavelength ranges
179  float * cw_wav O radiance-weighted center wavelengths
180  of the TOA spectrum
181  float * av_rad O average radiance in the inter-bands
182 
183  Modification history:
184 
185  W. Robinson, SAIC 26 May 2010 Original development
186 
187 ----------------------------------------------------------------------------*/ {
188  FILE *fp;
189  int items_set[N_INT_RNG * 2], ipar, ip1, ip2;
190  char line[80], name[80], value[80], parm[80];
191  char *p, *p1, *p2;
192  char *ib_par_nm[] = {"CW_WAVE", "AV_RAD"};
193  char *ib_par_rng[] = {"[IB01]", "[M1]", "[IB12]", "[M2]", "[IB23]", "[M3]",
194  "[IB34]", "[M4]", "[IB4I1]", "[I1]", "[M5]", "[IB56]", "[M6]", "[IB67]",
195  "[M7]", "[IB70]"};
196  float fval;
197  for (ipar = 0; ipar < N_INT_RNG * 2; ipar++)
198  items_set[ipar] = 0;
199  /*
200  * open the file
201  */
202  if ((fp = fopen(ib_file, "r")) == NULL) {
203  printf("%s, line %d: unable to open inter-band file: %s\n",
204  __FILE__, __LINE__, ib_file);
205  return 1;
206  }
207  /*
208  * Loop through parameters to define
209  */
210  while (fgets(line, 80, fp)) {
211  memset(name, '\0', sizeof ( name));
212  memset(value, '\0', sizeof ( value));
213  /*
214  * Skip comment lines, empty lines, and lines without a
215  * name = value pair.
216  */
217  if (line[0] == '#' || line[0] == '\n')
218  continue;
219  if (!(p = strchr(line, '=')))
220  continue;
221 
222  /*
223  * Parse parameter name string
224  */
225  p1 = line;
226  while (isspace(*p1))
227  p1++;
228  p2 = p - 1;
229  while (isspace(*p2))
230  p2--;
231  strncpy(name, p1, p2 - p1 + 1);
232 
233  /*
234  * Parse parameter value string
235  */
236  p1 = p + 1;
237  while (isspace(*p1))
238  p1++;
239  p2 = p1;
240  while (!isspace(*p2))
241  p2++;
242  strncpy(value, p1, p2 - p1);
243  printf("int band values are: %s = %s\n", name, value);
244 
245  /*
246  * for each parameter, set the array value
247  */
248  for (ipar = 0; ipar < N_INT_RNG * 2; ipar++) {
249  ip1 = 1;
250  if (ipar < N_INT_RNG)
251  ip1 = 0;
252  ip2 = ipar - ip1 * N_INT_RNG;
253  sprintf(parm, "%s%s", ib_par_nm[ip1], ib_par_rng[ip2]);
254  /*
255  *
256  */
257  if (strcmp(name, parm) == 0) {
258  items_set[ipar] = 1;
259  fval = (float) atof(value);
260  if (ip1 == 0)
261  cw_wav[ip2] = fval;
262  else
263  av_rad[ip2] = fval;
264  break;
265  }
266  }
267  }
268  /*
269  * make sure all values were read from the file
270  */
271  for (ipar = 0; ipar < N_INT_RNG * 2; ipar++)
272  if (items_set[ipar] == 0) {
273  printf("%s, %d: An item in the file %s was not set\n",
274  __FILE__, __LINE__, ib_file);
275  return 1;
276  }
277  fclose(fp);
278  return 0;
279 }
280 
281 int viirs_oxt_comp(ctl_struc *ctl, in_rec_struc *in_rec)
282 /*-----------------------------------------------------------------------------
283  routine: viirs_oxt_comp
284 
285  Description: compute the optical crosstalk artifact and modify the
286  radiance field with it
287 
288  Returns type: int - 0 if good
289 
290  Arguments:
291  Type Name I/O Description
292  ---- ---- --- -----------
293  ctl_struc * ctl I input controls
294  in_rec_struc * in_rec I/O controls for input record reading
295 
296  Note that ib_field - the field of inter-band radiances, is common
297  with this routine and the viirs_oxt routine
298 
299  Modification history:
300 
301  W. Robinson, SAIC 1 Jun 2010 Original development
302 
303 ----------------------------------------------------------------------------*/ {
304  static int ifirst = 0, n_rec_bnd, n_rec_det, n_int_rng, n_snd_rng, n_snd_det;
305  /* registration offset for the re-ordered sender bands:
306  (where is band M? relative to M1?)
307  M1 M2 M3 M4 M5 M6 M7 I1 I2 */
308  static int snd_off[] = {0, -3, -9, -6, -18, -21, -15, -12, -14};
309  static float *oxt_coeff;
310  h5io_str oxt_fid;
311  int irlin, irpix, irbnd, ilam, isbnd, isdet, reg_off, isoff;
312  float mod_rad;
313  /*
314  * initialization will read in the crosstalk coefficients
315  */
316  if (ifirst == 0) {
317  ifirst = 1;
318  /*
319  * open the optical crosstalk coefficient file
320  */
321  if (h5io_openr(ctl->oxt_coef, 0, &oxt_fid) != 0) {
322  printf("%s, %d - could not open HDF 5 Optical Xtalk file: %s\n",
323  __FILE__, __LINE__, ctl->oxt_coef);
324  return 1;
325  }
326  /*
327  * read array sizes:
328  */
329  if (h5io_rd_attr(&oxt_fid, "number of receiver bands", &n_rec_bnd) != 0) {
330  printf("%s, %d - failed to read Optical Xtalk attrib n_rec_bnd\n",
331  __FILE__, __LINE__);
332  return 1;
333  }
334  if (h5io_rd_attr(&oxt_fid, "number of receiver detectors", &n_rec_det)
335  != 0) {
336  printf("%s, %d - failed to read Optical Xtalk attrib n_rec_det\n",
337  __FILE__, __LINE__);
338  return 1;
339  }
340  if (h5io_rd_attr(&oxt_fid, "number of inter-band ranges", &n_int_rng)
341  != 0) {
342  printf("%s, %d - failed to read Optical Xtalk attrib n_int_rng\n",
343  __FILE__, __LINE__);
344  return 1;
345  }
346  if (h5io_rd_attr(&oxt_fid, "number of sender bands", &n_snd_rng) != 0) {
347  printf("%s, %d - failed to read Optical Xtalk attrib n_snd_rng\n",
348  __FILE__, __LINE__);
349  return 1;
350  }
351  if (h5io_rd_attr(&oxt_fid, "number of sender detectors", &n_snd_det)
352  != 0) {
353  printf("%s, %d - failed to read Optical Xtalk attrib n_snd_det\n",
354  __FILE__, __LINE__);
355  return 1;
356  }
357  printf("%s, %d: optical crosstalk array dimensions\n",
358  __FILE__, __LINE__);
359  printf(" n_rec_bnd = %d, n_rec_det = %d, n_int_rng = %d\n",
360  n_rec_bnd, n_rec_det, n_int_rng);
361  printf(" n_snd_rng = %d, n_snd_det = %d\n", n_snd_rng, n_snd_det);
362  /*
363  * make storage for the coefficients and read them in
364  */
365  if ((oxt_coeff = malloc(n_rec_bnd * n_rec_det * n_int_rng *
366  n_snd_rng * n_snd_det * sizeof ( float))) == NULL) {
367  printf("%s, %d: Unable to allocate optical crosstalk coeff storage\n",
368  __FILE__, __LINE__);
369  return 1;
370  }
371 
372  if (h5io_grab_ds(&oxt_fid, "optical crosstalk coefficients",
373  (void *) oxt_coeff) != 0) {
374  printf("%s, %d: Unable to read crosstalk coefficients\n",
375  __FILE__, __LINE__);
376  return 1;
377  }
378  /*
379  * and finish with the file
380  */
381  if (h5io_close(&oxt_fid) != 0) {
382  printf("%s, %d: Unable to close crosstalk coefficient file\n",
383  __FILE__, __LINE__);
384  return 1;
385  }
386  }
387  /*
388  * Apply the coefficients to the inter-band field and modify the radiances
389  * with the artifact
390  */
391  for (irlin = 0; irlin < NDET; irlin++) {
392  for (irpix = 0; irpix < in_rec->npix; irpix++) {
393  for (irbnd = 0; irbnd < N_VNIR_BND; irbnd++) {
394  mod_rad = 0.;
395  for (isdet = 0; isdet < NDET; isdet++) {
396  for (isbnd = 0; isbnd < 9; isbnd++) {
397  /*
398  * isoff is the offset to the sender detector, including the
399  * registration correction required
400  * the registration offset of the sender relative to M1, snd_off
401  * can be shown below
402  * sender M1 M2 M3 M4 M5 M6 M7 I1 I2
403  * order 0 1 2 3 4 5 6 7 8
404  * snd_off 0, -3, -9, -6, -18, -21, -15, -12, -14
405  */
406  reg_off = snd_off[ isbnd ] - snd_off[ irbnd ];
407  /*
408  * note that if data beyond margin is needed, skip including it
409  */
410  if ((reg_off + irpix >= 0) && (reg_off + irpix < in_rec->npix)) {
411  isoff = isdet * in_rec->npix + irpix + reg_off;
412  for (ilam = 0; ilam < 16; ilam++) {
413  /*
414  * oxt_coeff[ irbnd, irlin, ilam, isbnd, islin ]
415  */
416  mod_rad +=
417  *(oxt_coeff + irbnd + N_VNIR_BND * (irlin + NDET *
418  (ilam + 16 * (isbnd + 9 * isdet)))) *
419  *(ib_field[ilam] + isoff);
420  }
421  }
422  }
423  }
424  /*
425  * add the artifact to the radiance
426  */
427  *(in_rec->bnd_lt[irbnd] + irpix +
428  in_rec->npix * (irlin + in_rec->margin[0])) += mod_rad;
429  }
430  }
431  }
432  /*
433  * scan modified, return
434  */
435  return 0;
436 }
int32 value
Definition: Granule.c:1235
int h5io_openr(char *file, int opt, h5io_str *id)
Definition: h5io.c:4
#define NULL
Definition: decode_rs.h:63
int h5io_close(h5io_str *id)
Definition: h5io.c:115
int h5io_rd_attr(h5io_str *id, char *attr_name, void *data)
Definition: h5io.c:412
int viirs_oxt_ib_read(char *ib_file, float *cw_wav, float *av_rad)
Definition: viirs_oxt.c:166
int viirs_oxt_comp(ctl_struc *ctl, in_rec_struc *in_rec)
Definition: viirs_oxt.c:281
#define isspace(c)
int h5io_grab_ds(h5io_str *id, char *path_name, void *data)
Definition: h5io.c:1347
#define N_INT_RNG
Definition: viirs_oxt.c:6
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second and prod_ix set to PARAM_TYPE_INT name_prefix is compared with the beginning of the product name If name_suffix is not empty the it must match the end of the product name The characters right after the prefix are read as an integer and prod_ix is set to that number strncpy(l2prod->name_prefix, "myprod", UNITLEN)
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
#define NDET
Definition: polcor.c:13
#define N_VNIR_BND
Definition: viirs_sim_sdr.h:31
int viirs_oxt(ctl_struc *ctl, in_rec_struc *in_rec)
Definition: viirs_oxt.c:9
int npix
Definition: get_cmp.c:27
float p[MODELMAX]
Definition: atrem_corl1.h:131