NASA Logo
Ocean Color Science Software

ocssw V2022
get_rhown_nir.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 
3 
4 /* ---------------------------------------------------------------------- */
5 /* Convert Rrs[0+] to Rrs[0-] */
6 
7 /* ---------------------------------------------------------------------- */
8 float above_to_below(float Rrs) {
9  return (Rrs / (0.52 + 1.7 * Rrs));
10 }
11 
12 /* ---------------------------------------------------------------------- */
13 /* Convert Rrs[0-] to Rrs[0+] */
14 
15 /* ---------------------------------------------------------------------- */
16 float below_to_above(float Rrs) {
17  return (Rrs * 0.52 / (1 - 1.7 * Rrs));
18 }
19 
20 /* ------------------------------------------------------------------- */
21 /* Description: */
22 /* This computes the normalized water-leaving reflectances */
23 /* at the CZCS 670 channel. */
24 /* */
25 /* Outputs: */
26 /* rhown(670) */
27 /* */
28 /* Algorithm Author: */
29 /* Sean Bailey, Futuretech Corporation */
30 /* based largely on the work of R. Arnone and R. Stumpf */
31 /* using eta relationship from Z. Lee */
32 
33 /* ------------------------------------------------------------------- */
34 
35 void rhown_red(char *fqfile, float chl, float aw[], float bbw[], float Rrs[],
36  float wave[], int32_t nwave, int32_t ib_red, float rhown[]) {
37  static int firstCall = 1;
38  static int32_t ib5;
39  static float eta = 0.5;
40  float static chl_min = 0.2;
41  float static chl_max = 30.;
42 
43  float aw_red, apg_red, bb5;
44  float a, bb, chl_in;
45  float salbedo;
46  float Rrs_red;
47  float Rrs5;
48  float foq[nwave];
49 
50  if (firstCall) {
51  ib5 = windex(555, wave, nwave);
52  if (fabs(555 - wave[ib5]) > 15) {
53  printf("%s line %d: can't find reasonable green band\n", __FILE__, __LINE__);
54  printf("looking for 555, found %f\n", wave[ib5]);
55  exit(EXIT_FAILURE);
56  }
57 
58  firstCall = 0;
59  }
60 
61  chl_in = MAX(MIN(chl, chl_max), chl_min);
62 
63  aw_red = aw [ib_red];
64 
65  Rrs5 = Rrs[ib5];
66 
67  if (Rrs5 <= 0.0)
68  Rrs5 = 0.001;
69 
70  foqint_morel(fqfile, wave, nwave, 0.0, 0.0, 0.0, chl_in, foq);
71 
72  /* Compute total absorption at 670 */
73  // NOMAD fit of apg670 to chl
74  apg_red = exp(log(chl) * 0.9389 - 3.7589);
75  apg_red = MIN(MAX(apg_red, 0.0), 0.5);
76  a = aw_red + apg_red;
77 
78  /* Compute backscatter at 550 from Carder/Lee */
79  //bb5 = (-0.00182 + 2.058*Rrs5 + bbw5);
80  bb5 = (-0.00182 + 2.058 * Rrs5);
81 
82  /* Translate bb to NIR wavelength */
83  bb = bb5 * pow((wave[ib5] / wave[ib_red]), eta) + bbw[ib_red];
84 
85  /* Remote-sensing reflectance */
86  salbedo = bb / (a + bb);
87  Rrs_red = foq[ib_red] * salbedo;
88  /* Normalized water-leaving reflectance */
89  Rrs_red = below_to_above(Rrs_red);
90  rhown[ib_red] = PI*Rrs_red;
91 }
92 
93 
94 
95 /* ------------------------------------------------------------------- */
96 /* Description: */
97 /* This computes the normalized water-leaving reflectances */
98 /* for NIR bands using a bio-optical model and assuming that the */
99 /* NIR absorption is due to the water only. */
100 /* */
101 /* Algorithm Author: */
102 /* Sean Bailey, Futuretech Corporation */
103 /* based largely on the work of R. Arnone and R. Stumpf */
104 /* using eta relationship from Z. Lee */
105 
106 /* ------------------------------------------------------------------- */
107 
108 void rhown_nir(char *fqfile, float chl, float aw[], float bbw[], float Rrs[], float wave[],
109  int32_t nwave, float solz, float senz, float phi,
110  int32_t nir_s, int32_t nir_l, float rhown[],l2str *l2rec, int32_t ip) {
111  static int firstCall = 1;
112  static int32_t ib2;
113  static int32_t ib5;
114  static int32_t ib6;
115 
116  uncertainty_t *uncertainty=l2rec->l1rec->uncertainty;
117 
118  float *dchl=NULL,*drhown=NULL,*covariance;
119 
120  if(uncertainty){
121  dchl=&uncertainty->dchl;
122  drhown=uncertainty->drhown_nir;
123  if(input->proc_uncertainty==2)
124  covariance=uncertainty->covaraince_matrix;
125  else
126  covariance=uncertainty->pixel_covariance;
127  }
128 
129  float a6, aw6, apg6, bbp6;
130  float a, bb, eta;
131  float static chl_min = 0.2;
132  float static chl_max = 30.;
133  float foq[nwave];
134  float salbedo;
135  float Rrs_nir;
136  float Rrs2, Rrs5, Rrs6, Rrs6_star;
137  int32_t ib;
138  int32_t dnir = MAX(nir_l - nir_s, 1);
139  if(input->aer_opt==AERRHSM)
140  dnir=1;
141 
142 
143  float da6;
144  float dRrs2, dRrs5, dRrs6, dRrs6_star;
145  float da,dbb,deta;
146  float dbbp6;
147  float dsalbedo;
148  float dRrs_nir, Rrs_below;
149 
150  if (firstCall) {
151  ib6 = windex(670, wave, nwave);
152  if (fabs(670 - wave[ib6]) > 50) {
153  printf("%s line %d: can't find reasonable red band\n", __FILE__, __LINE__);
154  printf("looking for 670, found %f\n", wave[ib6]);
155  exit(EXIT_FAILURE);
156  }
157 
158  ib5 = windex(555, wave, nwave);
159  if (fabs(555 - wave[ib5]) > 15) {
160  printf("%s line %d: can't find reasonable green band\n", __FILE__, __LINE__);
161  printf("looking for 555, found %f\n", wave[ib5]);
162  exit(EXIT_FAILURE);
163  }
164  ib2 = windex(443, wave, nwave);
165  if (fabs(443 - wave[ib2]) > 5) {
166  printf("%s line %d: can't find reasonable blue band\n", __FILE__, __LINE__);
167  printf("looking for 443, found %f\n", wave[ib2]);
168  exit(EXIT_FAILURE);
169  }
170  firstCall = 0;
171  }
172 
173  Rrs2 = Rrs[ib2];
174  Rrs5 = Rrs[ib5];
175  Rrs6 = Rrs[ib6];
176 
177  if (Rrs6 <= 0.0) {
178  for (ib = nir_s; ib <= nir_l; ib += dnir){
179  rhown[ib] = 0.0;
180  if(uncertainty)
181  drhown[ib]=0.0;
182  }
183  return;
184  }
185 
186  chl = MAX(MIN(chl, chl_max), chl_min);
187 
188  if(uncertainty && (fabs(chl-chl_max)<0.00001 || fabs(chl-chl_min)<0.00001))
189  *dchl=0.0;
190 
191  // NOMAD fit of apg670 to chl
192  apg6 = exp(log(chl) * 0.9389 - 3.7589);
193  apg6 = MIN(MAX(apg6, 0.0), 0.5);
194 
195  /* Compute total absorption at 670 */
196  aw6 = aw [ib6];
197  a6 = aw6 + apg6;
198 
199  if(uncertainty)
200  da6=apg6*0.9389*1.0/chl* (*dchl);
201 
202  /* Go below... */
203  Rrs2 = above_to_below(Rrs2);
204  Rrs5 = above_to_below(Rrs5);
205  Rrs6 = above_to_below(Rrs6);
206 
207  if(uncertainty){
208  dRrs2 = 0.52 * pow(Rrs2 / Rrs[ib2], 2) * sqrt(covariance[ib2*nwave+ib2]);
209  dRrs5 = 0.52 * pow(Rrs5 / Rrs[ib5], 2) * sqrt(covariance[ib5*nwave+ib5]);
210  dRrs6 = 0.52 * pow(Rrs6 / Rrs[ib6], 2) * sqrt(covariance[ib6*nwave+ib6]);
211  }
212  //foqint_morel(wave,nwave,0.0,0.0,0.0,chl_in,foq);
213  foqint_morel(fqfile, wave, nwave, solz, senz, phi, chl, foq);
214 
215  /* Compute the backscatter slope ala Lee */
216  eta = 0.0;
217  deta=0.0;
218  if (Rrs5 > 0.0 && Rrs2 > 0.0) {
219  eta = 2.0 * (1. - 1.2 * exp(-0.9 * (Rrs2 / Rrs5)));
220  eta = MIN(MAX(eta, 0.0), 1.0);
221 
222  if(uncertainty){
223  if(eta>1.|| eta<0.)
224  deta=0.0;
225  else
226  deta=2.4*0.9*exp(-0.9 * (Rrs2 / Rrs5))* sqrt( pow(dRrs2/Rrs5,2)+ pow(Rrs2*dRrs5/(Rrs5*Rrs5),2)-2*Rrs2/pow(Rrs5,3)*covariance[ib2*nwave+ib5]);
227  }
228  }
229 
230  /* Compute total backscatter at 670 */
231  Rrs6_star = Rrs6 / foq[ib6];
232  if(uncertainty)
233  dRrs6_star=dRrs6/foq[ib6];
234 
235  bbp6 = (Rrs6_star * a6 / (1. - Rrs6_star)) - bbw[ib6];
236  if(uncertainty)
237  dbbp6= sqrt( pow(Rrs6_star*da6/(1-Rrs6_star),2) + pow(dRrs6_star*a6/pow(1-Rrs6_star,2) ,2) );
238 
239  /* Compute normalized water-leaving reflectance at each NIR wavelength */
240  for (ib = nir_s; ib <= nir_l; ib += dnir) {
241 
242  if (ib == ib6) {
243  a = a6;
244  if(uncertainty)
245  da = da6;
246  } else {
247  a = aw[ib];
248  if(uncertainty)
249  da = 0;
250  }
251 
252  /* Translate bb to NIR wavelength */
253  bb = bbp6 * pow((wave[ib6] / wave[ib]), eta) + bbw[ib];
254 
255  if(uncertainty)
256  dbb=pow((wave[ib6] / wave[ib]), eta)* sqrt( dbbp6*dbbp6 + pow(bbp6*log(wave[ib6]/wave[ib])*deta,2) );
257 
258  /* Remote-sensing reflectance */
259  salbedo = bb / (a + bb);
260  if(uncertainty)
261  dsalbedo=salbedo*salbedo* sqrt( pow(da/bb,2) + pow(a*dbb/(bb*bb),2) );
262  //dsalbedo=sqrt( pow(bb*da,2)+ pow(a*dbb,2) )/pow(a+bb,2);
263 
264  Rrs_nir = foq[ib6] * salbedo;
265  if(uncertainty) {
266  dRrs_nir=foq[ib6] *dsalbedo;
267  Rrs_below=Rrs_nir;
268  }
269  /* Normalized water-leaving reflectance */
270  Rrs_nir = below_to_above(Rrs_nir);
271  rhown[ib] = PI*Rrs_nir;
272  if(uncertainty)
273  drhown[ib]=PI*1/0.52* Rrs_nir*Rrs_nir/(Rrs_below*Rrs_below)*dRrs_nir;
274  }
275 }
276 
277 
278 /* ------------------------------------------------------------------- */
279 /* get_rhown_nir(): calls the appropriate function to compute the */
280 /* normalized water-leaving reflectance contribution in the NIR */
281 /* */
282 /* B.A. Franz, OBPG 25 January 2006 */
283 
284 /* ------------------------------------------------------------------- */
285 
286 void get_rhown_eval(char *fqfile, float Rrs[], float wave[], int32_t nir_s, int32_t nir_l,
287  int32_t nwave, float aw[], float bbw[], float chl,
288  float solz, float senz, float phi, float rhown[],l2str *l2rec, int32_t ip) {
289 
290  if (wave[nir_l] < MAXWAVE_VIS) {
291  rhown_red(fqfile, chl, aw, bbw, Rrs, wave, nwave, nir_l, rhown);
292  } else {
293  rhown_nir(fqfile, chl, aw, bbw, Rrs, wave, nwave, solz, senz, phi, nir_s, nir_l, rhown, l2rec,ip);
294  }
295 
296  return;
297 }
298 
#define MAX(A, B)
Definition: swl0_utils.h:25
#define MIN(x, y)
Definition: rice.h:169
#define AERRHSM
Definition: l12_parms.h:37
#define NULL
Definition: decode_rs.h:63
void get_rhown_eval(char *fqfile, float Rrs[], float wave[], int32_t nir_s, int32_t nir_l, int32_t nwave, float aw[], float bbw[], float chl, float solz, float senz, float phi, float rhown[], l2str *l2rec, int32_t ip)
instr * input
#define PI
Definition: l3_get_org.c:6
#define MAXWAVE_VIS
Definition: sensorDefs.h:51
void rhown_red(char *fqfile, float chl, float aw[], float bbw[], float Rrs[], float wave[], int32_t nwave, int32_t ib_red, float rhown[])
Definition: get_rhown_nir.c:35
float above_to_below(float Rrs)
Definition: get_rhown_nir.c:8
void foqint_morel(char *file, float wave[], int32_t nwave, float solz, float senzp, float phi, float chl, float brdf[])
Definition: brdf.c:314
#define fabs(a)
Definition: misc.h:93
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
void rhown_nir(char *fqfile, float chl, float aw[], float bbw[], float Rrs[], float wave[], int32_t nwave, float solz, float senz, float phi, int32_t nir_s, int32_t nir_l, float rhown[], l2str *l2rec, int32_t ip)
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
float below_to_above(float Rrs)
Definition: get_rhown_nir.c:16