ocssw V2020
convert_band.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <math.h>
3 #include "l12_proto.h"
4 
5 #define MINWVDIFF 10
6 #define REFWVL 443
7 
8 typedef struct _context {
9  float *wvl_i;
10  float *tar_rrs;
11  float *aw_i, *bbw_i, *aBric_i, *bBric_i;
12  float *aw_o, *bbw_o, *aBric_o, *bBric_o;
13  float *aphRef, *adgRef, *bbpRef;
15 } ccstr;
16 
17 float conv_rrs_to_555(float Rrs, float wave) {
18  float sw, a1, b1, a2, b2;
19 
20  if (fabs(wave - 555) > 2) {
21  if (fabs(wave - 547) <= 2) {
22  sw = 0.001723;
23  a1 = 0.986;
24  b1 = 0.081495;
25  a2 = 1.031;
26  b2 = 0.000216;
27  } else if (fabs(wave - 550) <= 2) {
28  sw = 0.001597;
29  a1 = 0.988;
30  b1 = 0.062195;
31  a2 = 1.014;
32  b2 = 0.000128;
33  } else if (fabs(wave - 560) <= 2) {
34  sw = 0.001148;
35  a1 = 1.023;
36  b1 = -0.103624;
37  a2 = 0.979;
38  b2 = -0.000121;
39  } else if (fabs(wave - 565) <= 2) {
40  sw = 0.000891;
41  a1 = 1.039;
42  b1 = -0.183044;
43  a2 = 0.971;
44  b2 = -0.000170;
45  } else {
46  printf("-E- %s line %d: Unable to convert Rrs at %f to 555nm.\n", __FILE__, __LINE__, wave);
47  exit(1);
48  }
49 
50  if (Rrs < sw)
51  Rrs = pow(10.0, a1 * log10(Rrs) - b1);
52  else
53  Rrs = a2 * Rrs - b2;
54  }
55 
56  return (Rrs);
57 }
58 
59 /*
60 int windex(float wave, float twave[], int ntwave)
61 {
62  int iw, index;
63  float wdiff;
64  float wdiffmin = 99999.;
65 
66  for (iw=0; iw<ntwave; iw++) {
67 
68  // break on exact match
69  if (twave[iw] == wave) {
70  index = iw;
71  break;
72  }
73 
74  // look for closest
75  wdiff = fabs(twave[iw]-wave);
76  if (wdiff < wdiffmin) {
77  wdiffmin = wdiff;
78  index = iw;
79  }
80  }
81 
82  return(index);
83 }
84  */
85 void grabBricaud(float wvl[], int nwvl, float* a_bricaud, float* b_bricaud) {
86 
87  float bricaud_L_A_B[][17] ={
88  {410, 412, 413, 440, 443, 488, 490, 510, 530, 531, 547, 550, 555, 560, 665, 667, 670},
89  {0.0313, 0.0323, 0.032775, 0.0403, 0.0394, 0.0279, 0.0274, 0.018, 0.0117, 0.0115, 0.00845,
90  0.008, 0.007, 0.0062, 0.0152, 0.01685, 0.0189},
91  {0.283, 0.286, 0.28775, 0.332, 0.3435, 0.369, 0.361, 0.260, 0.139, 0.134, 0.0625,
92  0.052, 0.0315, 0.016, 0.134, 0.140, 0.149}
93  };
94 
95  int wvlToSrch = 17;
96  int kb = 0, ko = 0;
97  int srchCounter = 0;
98 
99  /* if (wvl[nwvl-1] > bricaud_L_A_B[0][wvlToSrch - 1]){
100  printf ("-E- %s line %d: Bricaud coefficient not available for band %f :(\n",
101  __FILE__,__LINE__,wvl[ko]);
102  exit(1);
103  }*/
104  while (ko < nwvl && srchCounter < wvlToSrch) {
105 
106  if (wvl[ko] == bricaud_L_A_B[0][kb]) {
107  a_bricaud[ko] = bricaud_L_A_B[1][kb];
108  b_bricaud[ko] = bricaud_L_A_B[2][kb];
109  ko++;
110  srchCounter = 0;
111  } else if (wvl[ko] < bricaud_L_A_B[0][kb]) {
112  kb--;
113  srchCounter += 1;
114  } else {
115  kb++;
116  srchCounter += 1;
117  }
118  }
119 }
120 
121 void grabAwBw(float wvl[], int nwvl, float aw[], float bw[]) {
122  int i;
123  for (i = 0; i < nwvl; i++) {
124  aw[i] = aw_spectra(wvl[i], 11);
125  bw[i] = bbw_spectra(wvl[i], 11);
126  }
127  // float pwLAwBw[][17] =
128  // {
129  // {410,412,413,440,443,488,490,510,530,531,547,550,555,560,665,667,670},
130  // {0.00473,0.00455056,0.00449607,0.00635,0.00706914,0.0145167,0.015,0.0325,0.0434,
131  // 0.0439153,0.0531686,0.0565,0.0596,0.0619,0.429,0.434888,0.439},
132  // {0.0059145,0.00579201,0.00573196,0.00437075,0.00424592,0.00281659,0.00276835,0.00233870,
133  // 0.00198953,0.00197385,0.00174280,0.00170333,0.00163999,0.00157958,0.000772104,
134  // 0.000762543,0.000748479}
135  // };
136  //
137  // int wvlToSrch = 17;
138  // int kb=0,ko=0;
139  // int srchCounter = 0;
140  //
141  // /*if (wvl[nwvl-1] > pwLAwBw[0][wvlToSrch-1]){
142  // printf ("-E-: purewater coefficient not found for band %f >:(\n",wvl[ko]);
143  // exit(1);
144  // }
145  // */
146  // while (ko<nwvl && srchCounter<wvlToSrch){
147  // if (wvl[ko] == pwLAwBw[0][kb]){
148  // aw[ko] = pwLAwBw[1][kb];
149  // bw[ko] = pwLAwBw[2][kb];
150  // ko++;
151  // srchCounter = 0;
152  // }
153  //
154  // else if (wvl[ko] < pwLAwBw[0][kb]){
155  // kb--;
156  // srchCounter+=1;
157  // }
158  // else{
159  // kb++;
160  // srchCounter+=1;
161  // }
162  // }
163 }
164 
165 void calcQAA443(float ab_surf_rrs[], float wvl[], int nwvl, ccstr *ctxt) {
166  int ib;
167  float *blw_surf_rrs, *u_ptr;
168  float c1 = 0.52, c2 = 1.7;
169  float g0 = 0.089, g1 = 0.125;
170  float a0, bbp0, bbp0_exp, a0_exp;
171  float ratio_aph, slope_adg, ratio_adg;
172  static int firstCall = 1;
173  static int rIdx = -1, gIdx = -1, bgIdx = -1, refIdx = -1, purIdx = -1;
174 
175  blw_surf_rrs = (float*) malloc(nwvl * sizeof (float));
176  u_ptr = (float*) malloc(nwvl * sizeof (float));
177 
178  // get some indices:
179  if (firstCall) {
180  rIdx = windex(670, wvl, nwvl); // red (~670) band, idx=5 in Melin's code
181  gIdx = windex(550, wvl, nwvl); // green (~555) , idx=4 ---------------
182  bgIdx = windex(490, wvl, nwvl); // blue-green (~490),idx=2 ---------------
183  refIdx = windex(REFWVL, wvl, nwvl); //ref (~443), idx=1 ---------------
184  purIdx = windex(412, wvl, nwvl); // purple (~412), idx=0 ---------------
185  firstCall = 0;
186  }
187  float bbp[refIdx + 1], a[refIdx + 1];
188  // BEG ----- red band correction
189  float up_667 = 20 * pow(ab_surf_rrs[gIdx], 1.5);
190  float lw_667 = 0.9 * pow(ab_surf_rrs[gIdx], 1.7);
191  if (ab_surf_rrs[rIdx] > up_667 || ab_surf_rrs[rIdx] < lw_667) {
192  ab_surf_rrs[rIdx] = 1.27 * pow(ab_surf_rrs[gIdx], 1.47) +
193  0.00018 * pow(ab_surf_rrs[gIdx] / (ab_surf_rrs[bgIdx]), -3.19);
194 
195  }
196  // END ----- red band correciton
197 
198  for (ib = 0; ib < nwvl; ib++) {
199 
200  *(blw_surf_rrs + ib) = ab_surf_rrs[ib] / (c1 + c2 * ab_surf_rrs[ib]);
201  *(u_ptr + ib) = (-g0 + pow(pow(g0, 2) + 4 * g1 *
202  *(blw_surf_rrs + ib), 0.5)) / (2 * g1);
203  }
204  a0_exp = log10(
205  (*(blw_surf_rrs + refIdx) + *(blw_surf_rrs + bgIdx)) /
206  (*(blw_surf_rrs + gIdx) + 5 * (*(blw_surf_rrs + rIdx) /
207  *(blw_surf_rrs + bgIdx)) *
208  *(blw_surf_rrs + rIdx))
209  );
210  a0 = *(ctxt->aw_i + gIdx) + pow(10.0, (-1.146 - (1.366 * a0_exp) -(0.469 * pow(a0_exp, 2))));
211  bbp0 = *(u_ptr + gIdx) * a0 / (1 - *(u_ptr + gIdx)) - *(ctxt->bbw_i + gIdx);
212  bbp0_exp = 2 * (1 - 1.2 * exp(-0.9 * *(blw_surf_rrs + refIdx) / *(blw_surf_rrs + gIdx)));
213  ratio_aph = 0.74 + (0.2 / (0.8 + *(blw_surf_rrs + refIdx) / *(blw_surf_rrs + gIdx)));
214  slope_adg = 0.015 + 0.002 / (0.6 + *(blw_surf_rrs + refIdx) / *(blw_surf_rrs + gIdx));
215  ratio_adg = exp(slope_adg * (wvl[refIdx] - wvl[purIdx]));
216  for (ib = 0; ib <= refIdx; ib++) {
217 
218  bbp[ib] = bbp0 * pow((wvl[gIdx] / wvl[ib]), bbp0_exp);
219  a[ib] = (1 - *(u_ptr + ib)) * (bbp[ib] + *(ctxt->bbw_i + ib)) / *(u_ptr + ib);
220  }
221  *(ctxt->adgRef) = ((a[purIdx] - ratio_aph * a[refIdx]) - (*(ctxt->aw_i + purIdx) - ratio_aph * *(ctxt->aw_i + refIdx))) /
222  (ratio_adg - ratio_aph);
223  *(ctxt->aphRef) = a[refIdx] - *(ctxt->aw_i + refIdx) - *(ctxt->adgRef);
224  *(ctxt->bbpRef) = bbp[refIdx];
225 }
226 
227 /*
228  *Interpolation for cases with more than 2 input bands used
229  */
230 
231 float invDistInterp(ccstr *ctxt, float xout) {
232  /*
233  float invDistInterp(float xin[],float yin[],int ninout,float xout){*/
234  float wt, interpOut;
235  float wtSum = 0.0;
236  float prodSum = 0.0;
237  int in;
238  for (in = 0; in<*ctxt->num_wvl_i; in++) {
239  wt = 1.0 / fabs(ctxt->wvl_i[in] - xout);
240  prodSum += wt * ctxt->tar_rrs[in];
241  wtSum += wt;
242 
243  }
244  interpOut = prodSum / wtSum;
245  return interpOut;
246 }
247 
248 /*
249  * Function to identify band or bands (at most 2) to be used
250  */
251 
252 void idInputBand(float inputWvl[], float targetWvl, ccstr *ctxt, int nin) {
253  /* Function to match the target Wvl with one or two input Wvl.
254  * Two wavelengths are deemed necessary if the target wvl and the closest input wvl
255  * are more than 10nm apart. If this is the case and the input wvl is smaller than the
256  * target wvl, then the next higher wvl is also recruited. Similarly, if the input wvl is greater
257  * than the target wvl, the next smaller wvl is also recruited for shifting.
258  */
259  int closeIdx, kb;
260  int nout = 2;
261  int idx[2];
262  closeIdx = windex(targetWvl, inputWvl, nin);
263  if (((inputWvl[closeIdx] - targetWvl) > MINWVDIFF) && (closeIdx > 0)) {
264  idx[0] = closeIdx - 1;
265  idx[1] = closeIdx;
266  } else if (((inputWvl[closeIdx] - targetWvl) < -MINWVDIFF) && (closeIdx < (nin - 1))) {
267  idx[0] = closeIdx;
268  idx[1] = closeIdx + 1;
269  } else {
270  nout = 1;
271  idx[0] = closeIdx;
272  }
273  ctxt->sh_strt_idx = (int*) malloc(nout * sizeof (int));
274  ctxt->wvl_i = (float*) malloc(nout * sizeof (float));
275  ctxt->num_wvl_i = (int*) malloc(sizeof (int));
276  *(ctxt->num_wvl_i) = nout;
277  //*startWvlIdx = (int*)malloc(*nout*sizeof(int));
278  for (kb = 0; kb < nout; kb++) {
279  *(ctxt->sh_strt_idx + kb) = idx[kb];
280  *(ctxt->wvl_i + kb) = inputWvl[idx[kb]];
281  }
282  //*(*(startWvlIdx)+kb) = idx[kb];
283 
284 }
285 
286 /*
287  * Actual band shifting function
288  */
289 
290 void shiftBand(float inputWvl[], float inputRrs[], float chla, int numBands, float tarWvl, ccstr* ctxt)
291  {
292  //float* wvlIn,*rrsIn;
293  float wvlIn, rrsIn, refwvl;
294  float s = 0.015, g0 = 0.089, g1 = 0.125, c1 = -0.52, c2 = 1.7;
295  static int gIdx = -1, refIdx = -1; //,firstCall=1;
296  float b2g, sdg, yy; //,chla;
297  float ll_i, aph_i, adg_i, bbp_i;
298  float a_tot_i, bb_tot_i, qaa_fwd_i, qaa_rrs_bbw_i, qaa_rrs_aw_i;
299  float ll_o, aph_o, adg_o, bbp_o;
300  float a_tot_o, bb_tot_o, qaa_fwd_o, qaa_rrs_bbw_o, qaa_rrs_aw_o, correc_factor;
301  int kb, nin, idx;
302 
303 
304  //if(firstCall){
305  gIdx = windex(550, inputWvl, numBands);
306  refIdx = windex(REFWVL, inputWvl, numBands);
307  //firstCall=0;
308  refwvl = inputWvl[refIdx];
309  //}
310  nin = *(ctxt->num_wvl_i);
311 
312  b2g = inputRrs[refIdx] / inputRrs[gIdx];
313  sdg = s + 0.002 / (0.6 + b2g);
314  yy = 2 * (1 - 1.2 * exp(-0.9 * b2g));
315 
316  // chla = pow((*(ctxt->aphRef) / *(ctxt->aBric_i+refIdx)),1/(1-*(ctxt->bBric_i+refIdx)));
317  for (kb = 0; kb < nin; kb++) {
318  idx = *(ctxt->sh_strt_idx + kb);
319  wvlIn = inputWvl[idx];
320  rrsIn = inputRrs[idx];
321 
322  ll_i = wvlIn - refwvl;
323  aph_i = aph_bricaud(wvlIn, chla); //pow( (*(ctxt->aBric_i+idx) * chla),*(ctxt->bBric_i+idx));
324  adg_i = *(ctxt->adgRef) * exp(-sdg * ll_i);
325  bbp_i = *(ctxt->bbpRef) * pow((443 / wvlIn), yy);
326  a_tot_i = aph_i + adg_i + *(ctxt->aw_i + idx);
327  bb_tot_i = bbp_i + *(ctxt->bbw_i + idx);
328  qaa_fwd_i = bb_tot_i / (a_tot_i + bb_tot_i);
329  qaa_rrs_bbw_i = (g0 + g1 * qaa_fwd_i) * qaa_fwd_i;
330  qaa_rrs_aw_i = (c1 * qaa_rrs_bbw_i) / ((c2 * qaa_rrs_bbw_i) - 1);
331 
332  ll_o = tarWvl - REFWVL;
333  aph_o = aph_bricaud(tarWvl, chla); //aph_o = pow( ( *(ctxt->aBric_o) * chla),*(ctxt->bBric_o));
334  adg_o = *(ctxt->adgRef) * exp(-sdg * ll_o);
335  bbp_o = *(ctxt->bbpRef) * pow((443 / tarWvl), yy);
336  a_tot_o = aph_o + adg_o + *(ctxt->aw_o);
337  bb_tot_o = bbp_o + *(ctxt->bbw_o);
338  qaa_fwd_o = bb_tot_o / (a_tot_o + bb_tot_o);
339  qaa_rrs_bbw_o = (g0 + g1 * qaa_fwd_o) * qaa_fwd_o;
340  qaa_rrs_aw_o = (c1 * qaa_rrs_bbw_o) / ((c2 * qaa_rrs_bbw_o) - 1);
341  correc_factor = qaa_rrs_aw_o / qaa_rrs_aw_i;
342  *(ctxt->tar_rrs + kb) = rrsIn * correc_factor;
343  }
344 
345 }
346 
347 void deMalloc(ccstr* ctxt) {
348  free(ctxt->aw_i);
349  free(ctxt->bbw_i);
350  free(ctxt->aBric_i);
351  free(ctxt->bBric_i);
352  free(ctxt->aw_o);
353  free(ctxt->bbw_o);
354  free(ctxt->aBric_o);
355  free(ctxt->bBric_o);
356  free(ctxt->aphRef);
357  free(ctxt->adgRef);
358  free(ctxt->bbpRef);
359  free(ctxt->tar_rrs);
360  free(ctxt->num_wvl_i);
361  free(ctxt->wvl_i);
362  free(ctxt->sh_strt_idx);
363 
364 }
365 
366 float bioBandShift(float wvl[], float rrs[], float chla, int nw, float tarWvl) {
367  ccstr sh_ctxt; //band-shifting context structure
368  float result;
369  idInputBand(wvl, tarWvl, &sh_ctxt, nw);
370  sh_ctxt.aw_i = (float*) malloc(nw * sizeof (float));
371  sh_ctxt.bbw_i = (float*) malloc(nw * sizeof (float));
372  sh_ctxt.aBric_i = (float*) malloc(nw * sizeof (float));
373  sh_ctxt.bBric_i = (float*) malloc(nw * sizeof (float));
374  sh_ctxt.aw_o = (float*) malloc(sizeof (float));
375  sh_ctxt.bbw_o = (float*) malloc(sizeof (float));
376  sh_ctxt.aBric_o = (float*) malloc(sizeof (float));
377  sh_ctxt.bBric_o = (float*) malloc(sizeof (float));
378  sh_ctxt.aphRef = (float*) malloc(sizeof (float));
379  sh_ctxt.adgRef = (float*) malloc(sizeof (float));
380  sh_ctxt.bbpRef = (float*) malloc(sizeof (float));
381  sh_ctxt.tar_rrs = (float*) malloc(sizeof (float));
382 
383  grabAwBw(wvl, nw, sh_ctxt.aw_i, sh_ctxt.bbw_i);
384  grabAwBw(&tarWvl, 1, sh_ctxt.aw_o, sh_ctxt.bbw_o);
385 
386  calcQAA443(rrs, wvl, nw, &sh_ctxt);
387  shiftBand(wvl, rrs, chla, nw, tarWvl, &sh_ctxt);
388  if (*sh_ctxt.num_wvl_i > 1)
389  result = invDistInterp(&sh_ctxt, tarWvl);
390 
391  else
392  result = *(sh_ctxt.tar_rrs);
393  deMalloc(&sh_ctxt);
394  return result;
395 }
float * bbw_o
float * tar_rrs
Definition: convert_band.c:10
float * adgRef
Definition: convert_band.c:13
void grabBricaud(float wvl[], int nwvl, float *a_bricaud, float *b_bricaud)
Definition: convert_band.c:85
float * bbpRef
Definition: convert_band.c:13
float * aphRef
Definition: convert_band.c:13
float * bBric_i
Definition: convert_band.c:11
void grabAwBw(float wvl[], int nwvl, float aw[], float bw[])
Definition: convert_band.c:121
float * wvl_i
float * aBric_o
Definition: convert_band.c:12
#define REFWVL
Definition: convert_band.c:6
void deMalloc(ccstr *ctxt)
Definition: convert_band.c:347
float * aw_o
float * bbw_i
float invDistInterp(ccstr *ctxt, float xout)
Definition: convert_band.c:231
float * aw_i
int * num_wvl_i
Definition: convert_band.c:14
void calcQAA443(float ab_surf_rrs[], float wvl[], int nwvl, ccstr *ctxt)
Definition: convert_band.c:165
int * sh_strt_idx
Definition: convert_band.c:14
float aw_spectra(int32_t wl, int32_t width)
float * bBric_o
Definition: convert_band.c:12
float bbw_spectra(int32_t wl, int32_t width)
#define MINWVDIFF
Definition: convert_band.c:5
float * aBric_i
Definition: convert_band.c:11
instead the metadata field ProcessingEnvinronment is filled in from the output of a call to the POSIX compliant function uname from within the L1B code A small bug in L1B_Tables an incorrect comparison of RVS coefficients for TEBs to RVS coefficients for RSBs was being made This was replaced with a comparison between TEB coefficients This error never resulted in an incorrect RVS correction but did lead to recalculating the coefficients for each detector in a thermal band even if the coefficients were the same for all detectors To reduce to overall size of the reflective LUT HDF fill values were eliminated from all LUTs previously dimensioned where and where NUM_TIMES is the number of time dependent table pieces In Preprocess a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed This counter is internal only and is not yet used for any chiefly to casting of were added to make it LINUX compatible Output of code run on LINUX machines displays differences of at most scaled sector incalculable values of the Emissive calibration factor b1
Definition: HISTORY.txt:576
void shiftBand(float inputWvl[], float inputRrs[], float chla, int numBands, float tarWvl, ccstr *ctxt)
Definition: convert_band.c:290
#define fabs(a)
Definition: misc.h:93
data_t s[NROOTS]
Definition: decode_rs.h:75
float bioBandShift(float wvl[], float rrs[], float chla, int nw, float tarWvl)
Definition: convert_band.c:366
int i
Definition: decode_rs.h:71
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to a
Definition: HISTORY.txt:424
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:71
float aph_bricaud(float wave, float chl)
Definition: aph.c:181
void idInputBand(float inputWvl[], float targetWvl, ccstr *ctxt, int nin)
Definition: convert_band.c:252
float conv_rrs_to_555(float Rrs, float wave)
Definition: convert_band.c:17