ocssw V2020
bioOptBandShift.c
Go to the documentation of this file.
1 /*
2  * @file bioOptBandShift.c
3  * @Brief Bio-optical band shift based on Melin & Sclep (2015).
4  * @author Erdem Karakoylu
5  * @author Ocean Biology Processing Group, Nasa Goddard Space Flight Center
6  *
7  * Version History:
8  * V. 0.1: - hardwired iop algorithm and LUT to get aw/bbw, adg(443),aph(443),bbp(443)
9  * - aph(443) computed w/ Bricaud 95 Coeffs.
10  * V. 0.2: - Removed hardwired qaa & and aw, bbw code
11  * - Added interface (qaa_iops_4_bshift) to get_qaa.c; returns adg(443),aph(443),bbp(443).
12  * - Added access to aph_bricaud_1995 in aph.c and get_aw_bbw
13  * V. 0.3: - Replaced aph_bricaud_1995 by aph_bricaud_1998
14  */
15 #include <stdio.h>
16 #include <stdlib.h>
17 #include <string.h>
18 #include <math.h>
19 #include "l12_proto.h"
20 #define REFWVL 443.0
21 
22 typedef struct _context {
23  int *ibvc_bin;
24  float greenband;
25  float *wvl_i, *wvl_o;
26  float *aw_i, *bbw_i;
27  float *aw_o, *bbw_o;
30 } ccstr;
31 
32 float BandShift(ccstr*, l2str*, int32_t, float);
33 void CheckNULL(void*, char*, int, char*);
34 float ToAboveWater(float);
35 float ToBelowWater(float);
36 float Get_bb(l2str*, ccstr*, int32_t, int, char);
37 float Get_a(l2str*, ccstr*, int32_t, int, char);
38 float BandInterp(ccstr*, float*, int, int);
39 void Setup(l2str*, ccstr*);
40 
41 void qaa_iops_4_bshift(l2str*, float*, float*);
42 //void pml_iops_4_bshift(l2str*,float*,float*);
43 // void gsm_iops_4_bshift(l2str*,float*,float*);
44 float aph_bricaud_1998(float, float);
45 
46 static int nband_vc = 7;
47 static float vc_bands[7] = {412., 443., 490., 510., 531., 555., 670.};
48 static int current_iscan = -1;
49 static int current_npix = -1;
50 static float *sh_prod_arr, *adg_ref, *bbp_ref;
51 
52 void bioOptBandShift(l2str *l2rec, l2prodstr *p, float prod[]) {
53 
54  int32_t ibvc = p->prod_ix;
55  int32_t ib, ip, ipb;
56  static ccstr vcc;
57  static int firstRun = 1;
58  if (firstRun) {
59  Setup(l2rec, &vcc);
60  firstRun = 0;
61  }
62  /* sh_prod_arr is a 2D array accessed as a simple pointer of dims npix rows, and nband_vc cols.
63  If band shift not necessary, corresponding cols filled w/ Rrs already in l2rec.
64  Occasional segment below -- happens only when switching to new row
65  avoids recalculation for multiple products otherwise*/
66  if (l2rec->l1rec->iscan != current_iscan) {
67  current_iscan = l2rec->l1rec->iscan;
68  current_npix = l2rec->l1rec->npix;
69  if (adg_ref != NULL)
70  free(adg_ref);
71  CheckNULL(adg_ref = (float*) malloc(current_npix * sizeof (float)),
72  "adg_ref", __LINE__, __FILE__);
73  if (bbp_ref != NULL)
74  free(bbp_ref);
75  CheckNULL(bbp_ref = (float*) malloc(current_npix * sizeof (float)),
76  "bbp_ref", __LINE__, __FILE__);
77  if (sh_prod_arr != NULL)
78  free(sh_prod_arr);
79  CheckNULL(sh_prod_arr = (float*) malloc(current_npix * nband_vc * sizeof (float))
80  , "sh_prod_arr", __LINE__, __FILE__);
81  qaa_iops_4_bshift(l2rec, adg_ref, bbp_ref);
82  //pml_iops_4_bshift(l2rec,adg_ref,bbp_ref);
83  //gsm_iops_4_bshift(l2rec,adg_ref,bbp_ref);
84  for (ip = 0; ip < current_npix; ip++) {
85  for (ib = 0; ib < nband_vc; ib++) {
86  // fill sh_prod_arr, by bandshifting if necessary, by copying l2rec rrs otherwise.
87  if (vcc.ibvc_bin[ib]) {
88  *(sh_prod_arr + ip * nband_vc + ib) = BandShift(&vcc, l2rec, ip, vc_bands[ib]);
89  } else {
90 
91  ipb = ip * l2rec->l1rec->l1file->nbands + bindex_get((int32_t) vc_bands[ib]);
92  *(sh_prod_arr + ip * nband_vc + ib) = l2rec->Rrs[ipb];
93  }
94  }
95  }
96  }
97 
98  for (ip = 0; ip < l2rec->l1rec->npix; ip++)
99  prod[ip] = *(sh_prod_arr + ip * nband_vc + ibvc);
100 
101 }
102 
103 void Setup(l2str *l2rec, ccstr *vccp) {
104  int kb, iflag;
105 
106  float lambda_i[4][5] = {
107  {510., 555., -1., -1., -1.}, //SWF
108  {488., 488., 531., 547., 667.}, //MODA
109  {510., 560., 560., 665., -1.}, //MER
110  {488., 488., 555., 555., 670.} //VIIRS
111  };
112 
113  float lambda_o[4][5] = {
114  {531., 531., -1., -1., -1.}, //SWF
115  {490., 510., 510., 555., 670.}, //MODA
116  {531., 531., 555., 670., -1.}, //MER
117  {490., 510., 510., 531., 670.} //VIIRS
118  };
119 
120  int ibvc_bin[4][7] = {//index of band needing bandshift for each sensor
121  {0, 0, 0, 0, 1, 0, 0}, //SWF
122  {0, 0, 1, 1, 0, 1, 1}, //MODA
123  {0, 0, 0, 0, 1, 1, 1}, //MER
124  {0, 0, 1, 1, 1, 0, 1} //VIIRS
125  };
126  vccp->wvl_ref = REFWVL;
127  switch (l2rec->l1rec->l1file->sensorID) {
128  case SEAWIFS:
129  {
130  vccp->nBandPairs = 2; // 2 input bands will be used to compute 531.
131  vccp->greenband = 555;
132  iflag = 0;
133  break;
134  }
135 
136  case MODISA:
137  {
138  vccp->nBandPairs = 5; //inp bands: 488,488,531,547,667, resp
139  vccp->greenband = 547;
140  iflag = 1;
141  break;
142  }
143 
144  case MERIS:
145  {
146  vccp->nBandPairs = 5;
147  vccp->greenband = 560;
148  iflag = 2;
149  break;
150  }
151  case VIIRSN:
152  {
153  vccp->nBandPairs = 5;
154  vccp->greenband = 555;
155  iflag = 3;
156  break;
157  }
158  }
159 
160  CheckNULL(vccp->ibvc_bin = (int*) malloc(nband_vc * sizeof (int)), "vccp->ibvc_bin", __LINE__, __FILE__);
161  CheckNULL(vccp->wvl_i = (float*) malloc(vccp->nBandPairs * sizeof (float)), "vccp->wvl_i", __LINE__, __FILE__);
162  CheckNULL(vccp->wvl_o = (float*) malloc(vccp->nBandPairs * sizeof (float)), "vccp->wvl_o", __LINE__, __FILE__);
163  CheckNULL(vccp->aw_i = (float*) malloc(vccp->nBandPairs * sizeof (float)), "vccp->aw_i", __LINE__, __FILE__);
164  CheckNULL(vccp->aw_o = (float*) malloc(vccp->nBandPairs * sizeof (float)), "vccp->aw_o", __LINE__, __FILE__);
165  CheckNULL(vccp->bbw_i = (float*) malloc(vccp->nBandPairs * sizeof (float)), "vccp->aw_all", __LINE__, __FILE__);
166  CheckNULL(vccp->bbw_o = (float*) malloc(vccp->nBandPairs * sizeof (float)), "vccp->bbw_o", __LINE__, __FILE__);
167 
168  for (kb = 0; kb < vccp->nBandPairs; kb++) {
169  vccp->wvl_i[kb] = lambda_i[iflag][kb];
170  vccp->wvl_o[kb] = lambda_o[iflag][kb];
171 
172  }
173  for (kb = 0; kb < nband_vc; kb++)
174  vccp->ibvc_bin[kb] = ibvc_bin[iflag][kb];
175 
176  get_aw_bbw(l2rec, vccp->wvl_i, vccp->nBandPairs, vccp->aw_i, vccp->bbw_i);
177  get_aw_bbw(l2rec, vccp->wvl_o, vccp->nBandPairs, vccp->aw_o, vccp->bbw_o);
178  get_aw_bbw(l2rec, &vccp->wvl_ref, 1, &vccp->aw_ref, &vccp->bbw_ref);
179 }
180 
181 float BandInterp(ccstr* vcc, float* rrs_sh, int idx_0, int num_in) {
182  /* Inverse distance interpolation
183  * E.g. for bandshift from SeaWiFS, rrs(510) and rrs(555) are used to calculate
184  * two separate rrs(531). These are used here to calculate a final rrs(531).
185  * Max of 2 input bands is assumed.
186  */
187 
188  float *wts, wts_sum = 0;
189  float rrs_sh_interp = 0;
190  int idx, ib;
191  CheckNULL(wts = (float*) malloc(num_in * sizeof (float)), "wts", __LINE__, __FILE__);
192 
193  for (idx = 0; idx < num_in; idx++) {
194  ib = idx_0 + idx;
195  wts[idx] = 1.0 / fabs(vcc->wvl_o[ib] - vcc->wvl_i[ib]);
196  wts_sum += wts[idx];
197  }
198 
199  for (idx = 0; idx < num_in; idx++) {
200  rrs_sh_interp += (wts[idx] * rrs_sh[idx]) / (wts_sum);
201 
202  }
203  free(wts);
204  return rrs_sh_interp;
205 }
206 
207 float BandShift(ccstr* vccp, l2str *l2rec, int32_t ip, float target_band) {
208  float g0 = 0.089, g1 = 0.1245;
209  float rrs_in, *rrs_e_f, rrs_f_i, rrs_f_f, rrs_e_temp, Rrs_e_f;
210  int ipb, ib, idx, t_num = 0, t_start;
211  float bb_i, bb_f;
212  float a_i, a_f, u_i, u_f;
213 
214  // locate target_band
215  for (ib = 0; ib < vccp->nBandPairs; ib++) {
216  if (vccp->wvl_o[ib] == target_band) {
217  if (t_num == 0)
218  t_start = ib;
219  t_num++;
220  }
221  }
222 
223  if (t_num == 0) {
224  printf("\n -E- line %d Unable to find target band,%f => t_num = 0 :( \n", __LINE__, target_band);
225  exit(FATAL_ERROR);
226  }
227 
228  CheckNULL(rrs_e_f = (float*) malloc(t_num * sizeof (float)), "rrs_e_f", __LINE__, __FILE__);
229 
230  for (idx = 0; idx < t_num; idx++) {
231  ib = t_start + idx;
232  ipb = ip * l2rec->l1rec->l1file->nbands + bindex_get((int32_t) vccp->wvl_i[ib]);
233  rrs_in = ToBelowWater(l2rec->Rrs[ipb]);
234  // get IOPs
235  bb_i = Get_bb(l2rec, vccp, ip, ib, 'i');
236  bb_f = Get_bb(l2rec, vccp, ip, ib, 'f');
237  a_i = Get_a(l2rec, vccp, ip, ib, 'i');
238  a_f = Get_a(l2rec, vccp, ip, ib, 'f');
239  // get u
240  u_i = bb_i / (a_i + bb_i);
241  u_f = bb_f / (a_f + bb_f);
242  // rrs (input & output) from forward model
243  rrs_f_i = g0 * u_i + g1 * pow(u_i, 2);
244  rrs_f_f = g0 * u_f + g1 * pow(u_f, 2);
245  // estimate shifted rrs
246  rrs_e_f[idx] = rrs_in * rrs_f_f / rrs_f_i;
247 
248  }
249  if (t_num > 1) {//run inverse distance weighted interpolation if > 1 input band
250  rrs_e_temp = BandInterp(vccp, rrs_e_f, t_start, t_num);
251  Rrs_e_f = ToAboveWater(rrs_e_temp);
252  } else {
253  Rrs_e_f = ToAboveWater(rrs_e_f[0]);
254  }
255  free(rrs_e_f);
256  return Rrs_e_f;
257 }
258 
259 void CheckNULL(void *ptr, char* ptrName, int lineNum, char* filename) {
260  if (ptr == NULL) {
261  printf("\n --E-- pointer %s at line %d is %s not allocated :(\n", ptrName, lineNum, filename);
262  exit(FATAL_ERROR);
263  }
264 }
265 
266 float Get_bb(l2str *l2rec, ccstr *vccp, int32_t ip, int vccidx, char key) {
267  int refIdx, greenIdx, ipb;
268  float rrsRef, rrsGreen, eta, wavl;
269  float bbw, bbp, bb;
270 
271  refIdx = bindex_get((int32_t) REFWVL);
272  ipb = ip * l2rec->l1rec->l1file->nbands + refIdx;
273  rrsRef = ToBelowWater(l2rec->Rrs[ipb]);
274 
275  greenIdx = bindex_get(vccp->greenband);
276  ipb = ip * l2rec->l1rec->l1file->nbands + greenIdx;
277  rrsGreen = ToBelowWater(l2rec->Rrs[ipb]);
278 
279  if (key == 'i') {
280  wavl = vccp->wvl_i[vccidx];
281  bbw = vccp->bbw_i[vccidx];
282  } else {
283  wavl = vccp->wvl_o[vccidx];
284  bbw = vccp->bbw_o[vccidx];
285  }
286  eta = 2. * (1 - 1.2 * exp(-0.9 * rrsRef / rrsGreen));
287  bbp = bbp_ref[ip] * pow((REFWVL / wavl), eta);
288  bb = bbp + bbw;
289  return (bb);
290 }
291 
292 float Get_a(l2str *l2rec, ccstr *vccp, int32_t ip, int vccidx, char key) {
293 
294  int refIdx, greenIdx, ipb;
295  float rrsRef, rrsGreen, wavl, s;
296  float aw, aph, adg, a;
297 
298  refIdx = bindex_get((int32_t) REFWVL);
299  ipb = ip * l2rec->l1rec->l1file->nbands + refIdx;
300  rrsRef = ToBelowWater(l2rec->Rrs[ipb]);
301 
302  greenIdx = bindex_get((int32_t) vccp->greenband);
303  ipb = ip * l2rec->l1rec->l1file->nbands + greenIdx;
304  rrsGreen = ToBelowWater(l2rec->Rrs[ipb]);
305 
306  if (key == 'i') {
307  wavl = vccp->wvl_i[vccidx];
308  aw = vccp->aw_i[vccidx];
309  } else {
310  wavl = vccp->wvl_o[vccidx];
311  aw = vccp->aw_o[vccidx];
312  }
313  s = 0.015 + 0.002 / (0.6 + (rrsRef / rrsGreen));
314  adg = adg_ref[ip] * exp(-s * (wavl - REFWVL));
315  aph = aph_bricaud_1998(wavl, l2rec->chl[ip]);
316  a = aph + adg + aw;
317  return (a);
318 }
319 
320 float ToBelowWater(float rrsAbove) {
321  float rrsBelow;
322  rrsBelow = rrsAbove / (0.52 + 1.7 * rrsAbove);
323  return rrsBelow;
324 }
325 
326 float ToAboveWater(float rrsBelow) {
327  float rrsAbove;
328  rrsAbove = 0.52 * rrsBelow / (1 - 1.7 * rrsBelow);
329  return rrsAbove;
330 }
float BandInterp(ccstr *, float *, int, int)
#define REFWVL
float * bbw_o
int16 iflag[18338]
float bbw_ref
#define NULL
Definition: decode_rs.h:63
void Setup(l2str *, ccstr *)
#define VIIRSN
Definition: sensorDefs.h:23
#define MERIS
Definition: sensorDefs.h:22
float aph_bricaud_1998(float, float)
Definition: aph.c:116
int bindex_get(int32_t wave)
Definition: windex.c:43
float * wvl_i
void qaa_iops_4_bshift(l2str *, float *, float *)
Definition: get_qaa.c:351
float ToAboveWater(float)
float BandShift(ccstr *, l2str *, int32_t, float)
float ToBelowWater(float)
float greenband
float wvl_ref
float * aw_o
float * bbw_i
float * aw_i
#define FATAL_ERROR
Definition: swl0_parms.h:5
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
void get_aw_bbw(l2str *l2rec, float wave[], int nwave, float *aw, float *bbw)
void bioOptBandShift(l2str *l2rec, l2prodstr *p, float prod[])
void CheckNULL(void *, char *, int, char *)
float Get_bb(l2str *, ccstr *, int32_t, int, char)
#define fabs(a)
Definition: misc.h:93
float * wvl_o
float Get_a(l2str *, ccstr *, int32_t, int, char)
data_t s[NROOTS]
Definition: decode_rs.h:75
#define SEAWIFS
Definition: sensorDefs.h:12
float aw_ref
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
#define MODISA
Definition: sensorDefs.h:19
float p[MODELMAX]
Definition: atrem_corl1.h:131
int * ibvc_bin