NASA Logo
Ocean Color Science Software

ocssw V2022
raman.c
Go to the documentation of this file.
1 /* =================================================================== */
2 /* module raman.c - Raman scattering correction for Rrs */
3 /* */
4 /* This module contains the functions to correct sensor-observed */
5 /* above-water remote sensing reflectances, Rrs. for Raman Scattering */
6 /* effects. */
7 /* */
8 /* References: */
9 /* Lee et al (1994) Model for the interpretation of hyperspectral */
10 /* remote-sensing reflectance, Applied Optics, 33(24), 5721-5732 */
11 /* doi:10.1364/AO.33.005721 */
12 /* */
13 /* Lee et al (2013) Penetration of UV-visible solar radiation in the */
14 /* global oceans: Insights from ocean color remote sensing, Journal of */
15 /* Geophysical Research Oceans, 118, 4241-4255, doi:10.1002/jgrc.20308 */
16 /* */
17 /* McKinna et al. (2016) Implementation of an analytical Raman */
18 /* scattering correction for satellite ocean-color processing, */
19 /* Opt. Express, 24(14), A1123-A1137, doi: 10.1364/OE.24.0A1123 */
20 /* */
21 /* Mobley. C.D. (2010) Hydrolight Technical Note 10: Interpretation */
22 /* of Raman Scattering Computations, Sequoia Scientific. */
23 /* */
24 /* Sathyendranath and Platt (1998) Ocean-color model incorporating */
25 /* transspectral processes, Applied Optics, 37(12), 2216-2227, */
26 /* doi:10.1364/AO.37.002216 */
27 /* */
28 /* Walrafen (1967) Raman spectral studies of the effects of temperature*/
29 /* on water structure, Journal of Chemical Physics, 47, 118-121, */
30 /* doi:10.1063/1.711834 */
31 /* */
32 /* Westberry et al. (2013) Influence if Raman scattering on ocean color*/
33 /* inversion models, Applied Optics, 52(22), 5552-5561, */
34 /* doi:10.1364/AO.52.005552 */
35 /* */
36 /* */
37 /* Implementation: */
38 /* L. McKinna NASA/OBPG/SAIC, April 2015 */
39 /* */
40 /* Notes: */
41 /* Feb 2019: Downwelling irradiance model updated to be consistent */
42 /* with atmospheric correction model. Added support for OLCIA, OLCIB, */
43 /* SGLI, and VIIRS-J1. */
44 /* */
45 /* May 2021: replaced hard-coded coefficients with netCDF files. */
46 /* Added support for PACE OCI. */
47 /* =================================================================== */
48 
49 #define NORAMAN 0
50 #define LEE2013 1
51 #define WESTBERRY2013 2
52 #define LEE1994 3
53 #define WAVEMIN 380
54 #define WAVEMAX 700
55 
56 #include <stdio.h>
57 #include <math.h>
58 #include <netcdf.h>
59 #include <gsl/gsl_fit.h>
60 #include "l12_proto.h"
61 #include "sensorDefs.h"
62 
63 static int32_t nbandVis; //Number of visible bands
64 static float *lambda; //Sensor wavelengths
65 
66 //Wavelength indices
67 static int idx412;
68 static int idx443;
69 static int idx488;
70 static int idx550;
71 static int idx670;
72 static int bandNum550; //Number of bands shorter than 500 nm;
73 static int limitVisFlag;
74 
75 //static float angstEst; //Estimate Angstrom exponent (Taua power law slope)
76 
77 static float nWater = 1.344; //Refractive index of seawater
78 
79 //Bio-optical model parameters
80 static float *aw; //water absorption coefficient at sensor bands
81 static float *bbw; //water backscattering coefficient at sensor bands
82 static float *awR; //water absorption coeficient at Raman excitation bands
83 static float *bbwR; //water backscattering coefficient at Raman excitation bandss
84 static float *bbtot; //total backscattering modelled at sensor bands
85 static float *atot; //total absorption modelled at sensor bands
86 static float *atotR; //total absorption coefficient at Raman excitation bands
87 static float *bbtotR; //total backscattering coefficient at Raman excitation bands
88 static float *aphStar; //Normalised phytoplankton absorption spectral shape at sensor bands
89 static float *aphStarR; //Normalised phytoplankton absorption spectral shape at Raman bands
90 static float *kdSen; //Diffuse attenuation coefficient at sensor bands
91 static float *kdRam; //Diffuse attenuation coefficient at raman bands
92 static float *kappaSen;
93 static float *kappaRam;
94 // at sensor bands
95 static float aph443; //QAA-derived phytoplankton absorption at 443 nm
96 static float adg443; //QAA-derived Absorption of CDOM and colored detrital matter at 443 nm
97 static float bbp443; //QAA-derived backscattering coefficient at 443 nm
98 static float Ybbp; //Power exponent of bbp
99 static float Sdg; //Exponential slope of adg
100 
101 static float *aph1Sen; //Bricaud aph0 coefficient at sensor bands
102 static float *aph2Sen; //Bricaud aph1 coefficient at sensor bands
103 static float *aph1Ram; //Bricaud aph0 coefficient at Raman bands
104 static float *aph2Ram; //Bricaud aph1 coefficient at Raman bands
105 
106 
107 //Radtran atmospheric model parameters
108 static float *eDRam; //Down-welling irradiance at Raman excitation bands
109 static float *eDSen; //Down-welling irradiance at sensor bands
110 
111 static float *t_sol_ram; //Down-welling irradiance at Raman excitation bands
112 
113 //static float *tauaRam; //Aerosol optical thickness at Raman bands
114 
115 static float *ramLam; //Raman excitation wavelengths
116 static float *ramBandW; //Raman band widths
117 static float *bRex; //Raman scattering coefficient
118 
119 static float *f0BarRam; //TOA solar irradiance at Raman bands
120 static float *rayRam; //Rayleigh transmittance at Raman bands
121 static float *ozRam; //Ozone absorption at Raman bands
122 static float *wvRam; //Water vapor absorption at Raman bands
123 static float *oxyRam; //Oxygen absorption at Raman bands
124 static float *no2Ram; //Oxygen absorption at Raman bands
125 
126 static float *Rrs_ram; //Rrs due to Raman scattering
127 
128 //Lee et al. 2013 Raman correction 1 coefficients
129 //Pointers for interpolation/extrapolation
130 static float *alphaCor;
131 static float *betaCor0;
132 static float *betaCor1;
133 
134 //Lee et al 2013 Empirical coefficients
135 static float lamCor1[6] = {412.0, 443.0, 488.0, 531.0, 551.0, 667.0};
136 static float alpha[6] = {0.003, 0.004, 0.011, 0.015, 0.017, 0.018};
137 static float beta0[6] = {0.014, 0.015, 0.010, 0.010, 0.010, 0.010};
138 static float beta1[6] = {-0.022, -0.023, -0.051, -0.070, -0.080, -0.081};
139 
140 
141 /* ------------------------------------------------------------------- */
142 /* getRamanCoeffs() get spectral coefficients from netCDF file */
143 /* ------------------------------------------------------------------- */
144 
145 void getRamanCoeff(int grpid, char* varName, float* variable) {
146 
147  int varid;
148  int dimid;
149  size_t varLen;
150 
151  //get the variable id
152  if (nc_inq_varid(grpid, varName, &varid) != NC_NOERR) {
153  fprintf(stderr, "-E- %s line %d: could not find netCDF variable \"%s\" in netCDF File.\n",
154  __FILE__, __LINE__, varName);
155  exit(1);
156  }
157 
158  //get the dimension id for the variable
159  if (nc_inq_dimid(grpid,varName, &dimid) != NC_NOERR ) {
160  fprintf(stderr, "-E- %s line %d: could not find %s in netCDF File.\n",
161  __FILE__, __LINE__, varName);
162  }
163 
164  //get the dimension length of the variable
165  if (nc_inq_dimlen(grpid, dimid, &varLen) != NC_NOERR ) {
166  fprintf(stderr, "-E- %s line %d: could not find netCDF variable dimensions in netCDF File.\n",
167  __FILE__, __LINE__);
168  exit(1);
169  }
170 
171  //Check to see if the dimensions of the coefficient arrays are consistent with
172  //visible bands of the sensor
173  if (varLen != nbandVis && limitVisFlag!=1) {
174  fprintf(stderr, "-E- %s line %d: number of Raman coefficients %ld does not equal sensors visible"
175  " bands %d.\n",
176  __FILE__, __LINE__,varLen,nbandVis);
177  exit(1);
178  }
179 
180  //get the spectral coefficient from the file
181  if (nc_get_var_float (grpid, varid, variable) != NC_NOERR ) {
182  fprintf(stderr, "-E- %s line %d: could not read netCDF variable %s in netCDF File.\n",
183  __FILE__, __LINE__,varName);
184  exit(1);
185  }
186 
187 }
188 
189 /* ---------------------------------------------------------------------------*/
190 /* raman_pixel_alloc() Allocate pointer memory for Raman computations */
191 /* ---------------------------------------------------------------------------*/
192 void raman_pixel_alloc(l2str *l2rec) {
193 
194  //Allocate static pointer memory
195  lambda = (float*) allocateMemory(nbandVis * sizeof (float), "lambda");
196 
197  aw = (float*) allocateMemory(nbandVis * sizeof (float), "aw");
198  bbw = (float*) allocateMemory(nbandVis * sizeof (float), "bbw");
199  awR = (float*) allocateMemory(nbandVis * sizeof (float), "awR");
200  bbwR = (float*) allocateMemory(nbandVis * sizeof (float), "bbwR");
201  aphStar = (float*) allocateMemory(nbandVis * sizeof (float), "aphStar");
202  aphStarR = (float*) allocateMemory(nbandVis * sizeof (float), "aphStarR");
203  bbtot = (float*) allocateMemory(nbandVis * sizeof (float), "bbtot");
204  atot = (float*) allocateMemory(nbandVis * sizeof (float), "atot");
205  bbtotR = (float*) allocateMemory(nbandVis * sizeof (float), "bbtotR");
206  atotR = (float*) allocateMemory(nbandVis * sizeof (float), "atotR");
207 
208  kdSen = (float*) allocateMemory(nbandVis * sizeof (float), "kdSen");
209  kdRam = (float*) allocateMemory(nbandVis * sizeof (float), "kdRam");
210  kappaSen = (float*) allocateMemory(nbandVis * sizeof (float), "kappaSen");
211  kappaRam = (float*) allocateMemory(nbandVis * sizeof (float), "kappaRam");
212 
213  betaCor0 = (float*) allocateMemory(nbandVis * sizeof (float), "betaCor0");
214  betaCor1 = (float*) allocateMemory(nbandVis * sizeof (float), "betaCor1");
215  alphaCor = (float*) allocateMemory(nbandVis * sizeof (float), "alphaCor");
216 
217  eDRam = (float*) allocateMemory(nbandVis * sizeof (float), "eDRam");
218  eDSen = (float*) allocateMemory(nbandVis * sizeof (float), "eDSen");
219 
220  t_sol_ram = (float*) allocateMemory(nbandVis * sizeof (float), "t_sol_ram");
221 
222  Rrs_ram = (float*) allocateMemory(l2rec->l1rec->npix *
223  l2rec->l1rec->l1file->nbands * sizeof (float), "Rrs_ram");
224 
225  ramLam = (float*) allocateMemory(nbandVis * sizeof (float), "ramLam");
226  ramBandW = (float*) allocateMemory(nbandVis * sizeof (float), "ramBandW");
227  bRex = (float*) allocateMemory(nbandVis * sizeof (float), "bRex");
228  f0BarRam = (float*) allocateMemory(nbandVis * sizeof (float), "foBarRam");
229  rayRam = (float*) allocateMemory(nbandVis * sizeof (float), "rayRam");
230  ozRam = (float*) allocateMemory(nbandVis * sizeof (float), "ozRam");
231  wvRam = (float*) allocateMemory(nbandVis * sizeof (float), "wvRam");
232  oxyRam = (float*) allocateMemory(nbandVis * sizeof (float), "oxyRam");
233  no2Ram = (float*) allocateMemory(nbandVis * sizeof (float), "no2Ram");
234  aph1Ram = (float*) allocateMemory(nbandVis * sizeof (float), "aph1Ram");
235  aph2Ram = (float*) allocateMemory(nbandVis * sizeof (float), "aph2Ram");
236 
237  aph1Sen = (float*) allocateMemory(nbandVis * sizeof (float), "aph1Sen");
238  aph2Sen = (float*) allocateMemory(nbandVis * sizeof (float), "aph2Sen");
239 
240 }
241 
242 /* ---------------------------------------------------------------------------*/
243 /* get_raman_coeffs() Allocate pointer memory for Raman computations */
244 
245 /* ---------------------------------------------------------------------------*/
246 void get_raman_coeffs(l2str *l2rec) {
247 
248  char filename[FILENAME_MAX];
249  char *filedir;
250 
251  if ((filedir = getenv("OCDATAROOT")) == NULL) {
252  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
253  exit(1);
254  }
255 
256  filehandle *l1file = l2rec->l1rec->l1file;
257  strcpy(filename, filedir);
258  strcat(filename, "/");
259  strcat(filename, sensorId2SensorDir(l1file->sensorID));
260  if(l1file->subsensorID > SEAWIFS_LAC) { // ignore seawifs subsensor ID
261  strcat(filename, "/");
262  strcat(filename, subsensorId2SubsensorDir(l1file->subsensorID));
263  }
264  strcat(filename, "/raman.nc");
265 
266  if(l1file->sensorID == OCIS) {
267  limitVisFlag = 1;
268  } else if(l1file->sensorID == OCI) {
269  limitVisFlag = 1;
270  } else {
271  limitVisFlag = 0;
272  }
273 
274  printf("Loading Raman coefficients from: %s.\n", filename);
275 
276  int ncid;
277  int grpid;
278  //Open netCDF file
279  if (nc_open(filename, NC_NOWRITE, &ncid) != NC_NOERR) {
280  fprintf(stderr, "-E- %s line %d: could not open netCDF File \"%s\".\n",
281  __FILE__, __LINE__, filename);
282  exit(1);
283  }
284 
285  //get the group id for raman_coeffs
286  if (nc_inq_grp_ncid(ncid, "raman_coeffs", &grpid) != NC_NOERR) {
287  fprintf(stderr, "-E- %s line %d: could not find netCDF group \"raman_coeffs\".\n",
288  __FILE__, __LINE__);
289  exit(1);
290  }
291 
292  getRamanCoeff(grpid, "ramLam",ramLam);
293  getRamanCoeff(grpid, "ramBandW",ramBandW);
294  getRamanCoeff(grpid, "bRex",bRex);
295  getRamanCoeff(grpid, "f0BarRam",f0BarRam);
296  getRamanCoeff(grpid, "rayRam",rayRam);
297  getRamanCoeff(grpid, "ozRam",ozRam);
298  getRamanCoeff(grpid, "wvRam",wvRam);
299  getRamanCoeff(grpid, "oxyRam",oxyRam);
300  getRamanCoeff(grpid, "no2Ram",no2Ram);
301  getRamanCoeff(grpid, "aph1Ram",aph1Ram);
302  getRamanCoeff(grpid, "aph2Ram",aph2Ram);
303 
304  //get the group id for sensor_coeffs
305  if (nc_inq_grp_ncid(ncid, "sensor_coeffs", &grpid) != NC_NOERR) {
306  fprintf(stderr, "-E- %s line %d: could not find netCDF group \"sensor_coeffs\".\n",
307  __FILE__, __LINE__);
308  exit(1);
309  }
310 
311  getRamanCoeff(grpid, "aph1Sen",aph1Sen);
312  getRamanCoeff(grpid, "aph2Sen",aph2Sen);
313 
314  //close the netcdf
315  if (nc_close(ncid) != NC_NOERR){
316  fprintf(stderr, "-E- %s line %d: could not close netCDF \"%s\".\n",
317  __FILE__, __LINE__,filename);
318  exit(1);
319  }
320 
321 }
322 /* -----------------------------------------------------------------------*/
323 /* set_raman_aph_uv() - Extrapolates aph* into the UV (below 400 nm). */
324 /* This methods assumes that aph can be linearly */
325 /* extrapolated into the UV */
326 
327 /* -----------------------------------------------------------------------*/
328 void set_raman_aph_uv(l2str *l2rec, int ip) {
329 
330  int iw;
331  float aphM, aphC;
332  float aphStar443;
333 
334  //Calculate aphStar at 443 nm using Bricaud et al 1998
335  aphStar443 = aph1Sen[idx443] * pow(l2rec->chl[ip], (aph2Sen[idx443]));
336 
337  //Using sensor-specific coefficients, calculate aphstar at all bands and
338  //normalize to 1.0 at 443 nm
339  for (iw = 0; iw < nbandVis; iw++) {
340  aphStar[iw] = (aph1Sen[iw] * pow(l2rec->chl[ip], (aph2Sen[iw])))
341  / aphStar443;
342  aphStarR[iw] = (aph1Ram[iw] * pow(l2rec->chl[ip], (aph2Ram[iw])))
343  / aphStar443;
344 
345  }
346 
347  //Linear model of aph betweeen 412 and 443 nm
348  aphM = (aphStar[idx443] - aphStar[idx412]) / (lambda[idx443] - lambda[idx412]);
349  aphC = aphStar[idx443] - aphM * lambda[idx443];
350 
351  //If outside the Bricaud in UV, extrapolate as linear function.
352  for (iw = 0; iw < nbandVis; iw++) {
353  if (ramLam[iw] <= 400) {
354  aphStarR[iw] = (aphM * ramLam[iw] + aphC);
355  }
356  }
357 
358 }
359 
360 /* -----------------------------------------------------------------------*/
361 /* fit_raman_tsol() - Interpolates/extrapolates tramsittanc.e from sun to */
362 /* surface (t_sol) at Raman excitation bands. The */
363 /* assumption is that t_sol follows a power law. */
364 /* -----------------------------------------------------------------------*/
365 void fit_raman_tsol(l2str *l2rec, int ip) {
366 
367  int iw;
368 
369  double logTsol[nbandVis];
370  double waveRatio[nbandVis];
371 
372  double c0, c1, cov00, cov01, cov11, sumsq;
373  int nbands = l2rec->l1rec->l1file->nbands;
374 
375  //Log transform data so function has a linear form
376  for (iw = 0; iw < nbandVis; iw++) {
377  logTsol[iw] = log(l2rec->l1rec->t_sol[ip * nbands + iw] / l2rec->l1rec->t_sol[ip * nbands + idx488]);
378  waveRatio[iw] = log(lambda[iw] / lambda[idx488]);
379  }
380 
381  //Use a linear fit on long-transformed data to compute power law exponent (range 412-490 nm)
382  gsl_fit_linear(waveRatio, 1, logTsol, 1, nbandVis, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
383 
384 
385  //Extrapolate/interpolate using power law model
386  for (iw = 0; iw < nbandVis; iw++) {
387  t_sol_ram[iw] = l2rec->l1rec->t_sol[ip * nbands + idx488] * pow((ramLam[iw] / lambda[idx488]), c1);
388  }
389 
390 }
391 
392 /* -----------------------------------------------------------------------*/
393 /* raman_qaa() - computes a first order estimate of the constituent IOPs */
394 /* -----------------------------------------------------------------------*/
395 void raman_qaa(l2str *l2rec, int ip) {
396 
397  int iw, idxref;
398 
399  float aref, bbpref, numer, denom, rho;
400  float rat, zeta, xi, term1, term2;
401 
402  float rrs_a[nbandVis];
403  float rrs_s[nbandVis];
404  float u[nbandVis];
405  float at[nbandVis];
406  float bbt[nbandVis];
407 
408  float g0 = 0.08945;
409  float g1 = 0.1247;
410  float acoefs[3] = {-1.146, -1.366, -0.469};
411 
412  //Calculate the sub-surface remote sensing reflectance
413  for (iw = 0; iw < nbandVis; iw++) {
414  rrs_a[iw] = l2rec->Rrs[ip * l2rec->l1rec->l1file->nbands + iw];
415  rrs_s[iw] = rrs_a[iw] / (0.52 + 1.7 * rrs_a[iw]);
416  }
417 
418  /*pre-test of Rrs550*/
419  if (rrs_a[idx550] <= 0.0)
420  rrs_a[idx550] = 0.001;
421 
422  /* pre-test 670 */
423  if ((rrs_a[idx670] > 20.0 * pow(rrs_a[idx550], 1.5)) || (rrs_a[idx670] < 0.9 * pow(rrs_a[idx550], 1.7))) {
424  rrs_a[idx670] = 1.27 * pow(rrs_a[idx550], 1.47) + 0.00018 * pow(rrs_a[idx488] / rrs_a[idx550], -3.19);
425  }
426 
427  /*Quadratic formula to get the quantity u*/
428  for (iw = 0; iw < nbandVis; iw++) {
429  u[iw] = (sqrt(g0 * g0 + 4.0 * g1 * rrs_s[iw]) - g0) / (2.0 * g1);
430  }
431 
432  /*Determine whether to use 550 or 670 nm as reference then compute total*/
433  /*absorption at the reference wavelength, aref*/
434  if (rrs_s[idx670] >= 0.0015) {
435  aref = aw[idx670] + 0.07 * pow(rrs_a[idx670] / rrs_s[idx443], 1.1);
436  idxref = idx670;
437 
438  } else {
439  numer = rrs_a[idx443] + rrs_a[idx488];
440  denom = rrs_a[idx550] + 5 * rrs_a[idx670]*(rrs_a[idx670] / rrs_a[idx488]);
441  rho = log10(numer / denom);
442  rho = acoefs[0] + acoefs[1] * rho + acoefs[2] * rho*rho;
443  aref = aw[idx550] + pow(10.0, rho);
444  idxref = idx550;
445  }
446 
447  /*Calculate the backscattering coefficient at the reference wavelength*/
448  bbpref = ((u[idxref] * aref) / (1.0 - u[idxref])) - bbw[idxref];
449 
450  /*Steps to compute power exponent of bbp coefficient*/
451  rat = rrs_s[idx443] / rrs_s[idx550];
452  Ybbp = 2.0 * (1.0 - 1.2 * exp(-0.9 * rat));
453 
454  /*Compute backscattering coefficient at 443nm, this is needed later for */
455  /*calculating bbp at the Raman excitation bands.s */
456  bbp443 = bbpref * pow((lambda[idxref] / lambda[idx443]), Ybbp);
457 
458  /*Compute the total backscattering coefficient for sensor bands*/
459  for (iw = 0; iw < nbandVis; iw++) {
460  bbt[iw] = bbpref * pow((lambda[idxref] / lambda[iw]), Ybbp) + bbw[iw];
461  }
462 
463  /*Calculate the total absorption coefficient at sensor bands*/
464  for (iw = 0; iw < nbandVis; iw++) {
465  at[iw] = ((1.0 - u[iw]) * bbt[iw]) / u[iw];
466  }
467 
468  /*Calculate the exponential slope coefficient of absorption by colored */
469  /*dissolved and detrital matter */
470  Sdg = 0.015 + 0.002 / (0.6 + rrs_s[idx443] / rrs_s[idxref]);
471 
472  /*Compute the absorption coefficient of colored dissolved and detrital */
473  /*matter at 443 nm, adg443*/
474  zeta = 0.74 + 0.2 / (0.8 + rrs_s[idx443] / rrs_s[idxref]);
475  xi = exp(Sdg * (lambda[idx443] - lambda[idx412]));
476  term1 = (at[idx412] - zeta * at[idx443]) - (aw[idx412] - zeta * aw[idx443]);
477  term2 = xi - zeta;
478  adg443 = term1 / term2;
479 
480  /*Calculate the absorption of phytoplankton at 443 nm*/
481  aph443 = at[idx443] - adg443 - aw[idx443];
482 
483 }
484 
485 /* -----------------------------------------------------------------------*/
486 /* raman_iops() - computes IOPs at raman excitation bands using a */
487 /* -----------------------------------------------------------------------*/
489 
490  int iw;
491  float bbpR, adgR, aphR;
492  float bbp, adg, aph;
493 
494  /*loop over sensor and raman bands*/
495  for (iw = 0; iw < nbandVis; iw++) {
496 
497  /*Compute IOPs at Raman bands*/
498  bbpR = bbp443 * pow((lambda[idx443] / ramLam[iw]), Ybbp);
499  adgR = adg443 * exp(-Sdg * (ramLam[iw] - lambda[idx443]));
500  aphR = aph443 * aphStarR[iw];
501 
502  /*Compute IOPs at sensor bands*/
503  bbp = bbp443 * pow((lambda[idx443] / lambda[iw]), Ybbp);
504  adg = adg443 * exp(-Sdg * (lambda[iw] - lambda[idx443]));
505  aph = aph443 * aphStar[iw];
506 
507  /*total absorption and backscattering coefficients at Raman bands*/
508  atotR[iw] = awR[iw] + adgR + aphR;
509  bbtotR[iw] = bbwR[iw] + bbpR;
510 
511  /*total absorption and backscattering coefficients at sensor bands*/
512  atot[iw] = aw[iw] + adg + aph;
513  bbtot[iw] = bbw[iw] + bbp;
514  }
515 
516 }
517 
518 /* -----------------------------------------------------------------------*/
519 /* k_functions() - computes k functions from IOP bio-optical models */
520 /* -----------------------------------------------------------------------*/
521 void raman_k_func(l2str *l2rec, int ip) {
522 
523  int iw;
524  float solzRad, solzRadSw, muD, muU;
525 
526  //Convert solar and sensor zenith angles from degrees to radianss
527  solzRad = l2rec->l1rec->solz[ip]*(M_PI / 180.0);
528  solzRadSw = asin(sin(solzRad) / nWater); //subsurface solar zenith angle
529 
530  //Means cosines
531  muD = cos(solzRadSw);
532  muU = 0.5;
533 
534  for (iw = 0; iw < nbandVis; iw++) {
535 
536  //Diffuse attenuation coefficient
537  kdSen[iw] = (atot[iw] + bbtot[iw]) / muD;
538  kdRam[iw] = (atotR[iw] + bbtotR[iw]) / muD;
539 
540  kappaSen[iw] = (atot[iw] + bbtot[iw]) / muU;
541  kappaRam[iw] = (atotR[iw] + bbtotR[iw]) / muU;
542 
543  }
544 
545 }
546 
547 /* -----------------------------------------------------------------------*/
548 /* radtran_raman_ed() - calculate total downwelling irradiance */
549 /* -----------------------------------------------------------------------*/
550 void raman_radtran_ed(l2str *l2rec, int ip) {
551 
552  int32_t iw, ipb;
553 
554  l1str *l1rec = l2rec->l1rec;
555 
556  //Constants
557  float p0 = 29.92 * 33.8639; //standard pressure (convert mmHg -> HPa);
558  float radCon = M_PI / 180.0; //degrees -> radians conversion
559  float ws = l2rec->l1rec->ws[ip]; //wind speed
560  float nw2= pow(nWater,2.); //refractive index of water squared
561 
562  //Define model variables
563  float sin2Solz, rf0, a0,a1,a2,a3, rho_fres;
564  float solz, solzRad, cosSolz, mPres, mAtm, f0Ex;
565  float th2oEx, to2Ex, tno2Ex, to3Ex, tgEx;
566  float a_285R, a_225R,tau_to200R;
567  float no2_tr200R, eDRam_above;
568 
569  //Load per-pixel variables from L2 record
570  solz = l1rec->solz[ip]; //solar zenith
571 
572  //Estimate the transmittance from sun to surface at Raman excitation bands
573  fit_raman_tsol(l2rec, ip);
574 
575  //Convert solar zenith in degrees to radians
576  solzRad = l1rec->solz[ip]*radCon;
577 
578  //Calculate cosine of the solar zenith angle
579  cosSolz = cos(solzRad);
580 
581  //Claculate sin2(solz)
582  sin2Solz = pow(sin(solzRad),2.);
583 
584  //Compute Fresnel reflection coefficient for a wind-blown surface following
585  //Haltrin 2002
586  //Fresnel specular reflection coefficient
587  rf0 = 0.5*( pow( (cosSolz - pow(nw2 - sin2Solz ,0.5))/(cosSolz + pow(nw2 + sin2Solz ,0.5)) ,2.)
588  + pow( (nw2*cosSolz - pow(nw2 - sin2Solz ,0.5))/(nw2*cosSolz + pow(nw2 + sin2Solz ,0.5)) ,2.) );
589 
590  //Wind-specific coefficients
591  a0 = 0.001*(6.944831 - 1.912076*ws + 0.03654833*ws*ws);
592  a1 = 0.7431368 + 0.0679787*ws - 0.0007171*ws*ws;
593  a2 = 0.5650262 + 0.0061502*ws - 0.0239810*ws*ws + 0.0010695*ws*ws*ws;
594  a3 = -0.4128083 - 0.1271037*ws + 0.0283907*ws*ws - 0.0011706*ws*ws*ws;
595 
596  //Fresnel reflectance as a function of wind speed
597  rho_fres = a0 + rf0*(a1 + rf0*(a2 + a3*rf0));
598 
599 
600  //---------------------
601  //Calculate the atmospheric pathlengths
602  //Use the Kasten abd Young (1989) coefficients
603  // mAtm = 1.0 / ( cosSolz + 0.50572*(96.07995 - solz)**(-1.6364) );
604  mAtm = 1.0 / (cosSolz + 0.50572 * pow((86.07995 - solz), -1.6364));
605 
606  //Greg and Carder (1990) values **NOT USED**
607  //mAtm = 1.0 / ( np.cos(solzRad) + 0.15* pow((93.885 - solz),-1.253) );
608 
609  //Pathlength corrected for non-standard atmoshperic pressure
610  mPres = mAtm * (l1rec->pr[ip] / p0);
611 
612 
613  //compute no2 above 200m per get_trans.c
614  if (l2rec->l1rec->no2_tropo[ip] > 0.0)
615  //compute tropo no2 above 200m (Z.Ahmad)
616  no2_tr200R = l1rec->no2_frac[ip] *l1rec->no2_tropo[ip];
617  else
618  no2_tr200R = 0.0;
619 
620  //---------------------//
621  //Spectral calculations - loop over sensor
622  for (iw = 0; iw < nbandVis; iw++) {
623 
624  ipb = ip*l1rec->l1file->nbands + iw;
625 
626  //Correct for sun-earth distance
627  f0Ex = f0BarRam[iw] * l1rec->fsol; //TOA solar irradiance at sensor bands
628 
629  //Compute ozone transmittance
630  to3Ex = exp(-(ozRam[iw] * l1rec->oz[ip] / l1rec->csolz[ip]) );
631 
632  //compute no2 absorption transmittance
633  a_285R = no2Ram[iw] * (1.0 - 0.003 * (285.0 - 294.0));
634  a_225R = no2Ram[iw] * (1.0 - 0.003 * (225.0 - 294.0));
635  tau_to200R = a_285R * no2_tr200R + a_225R * l1rec->no2_strat[ip];
636  tno2Ex = exp(-(tau_to200R / l1rec->csolz[ip]) );
637 
638  //Compute water vapour transmittance
639  th2oEx = exp((-0.2385 * wvRam[iw] * l1rec->wv[ip] * mAtm) /
640  (pow((1.0 + 20.07 * wvRam[iw] * l1rec->wv[ip] * mAtm), 0.45)));
641 
642 
643  //Gas Transmittance due to oxygen/gases
644  to2Ex = exp((-1.41 * oxyRam[iw] * mPres) / (pow((1.0 + 118.3 * oxyRam[iw] * mPres), 0.45)));
645 
646  tgEx = to3Ex*tno2Ex*th2oEx*to2Ex;
647 
648  //
649  eDRam_above = f0Ex * tgEx * t_sol_ram[iw]* cos(solzRad) * 10.0;
650  eDSen[iw] = l1rec->Fo[iw] * l1rec->tg_sol[ipb] * l1rec->t_sol[ipb]* cos(solzRad) * 10.0;
651 
652  eDRam[iw] = 1.03*(1.0 - rho_fres)*eDRam_above;
653 
654  }
655 
656 }
657 
658 /* -------------------------------------------------------------------------*/
659 /* raman_cor_1() - computes the Rrs corrected for Raman scattering reported */
660 /* by Lee et al, 2013 */
661 /* -------------------------------------------------------------------------*/
662 void raman_cor_lee1(l2str *l2rec, int ip) {
663 
664  int iw;
665  float Rrs443, Rrs550;
666  float rFactor;
667  float rrs_a;
668  int nbands = l2rec->l1rec->l1file->nbands;
669  //If pixel is already masked, return and do not attempt the correction.
670  if (l2rec->l1rec->mask[ip]) {
671  return;
672  }
673 
674  //Get Rrs443 and Rrs550 from the L2 record
675  Rrs443 = l2rec->Rrs[ip * nbands + idx443];
676  Rrs550 = l2rec->Rrs[ip * nbands + idx550];
677 
678  //Loop over vis bands
679  for (iw = 0; iw < nbandVis; iw++) {
680 
681  if (lambda[iw] >= WAVEMIN && lambda[iw] <= WAVEMAX) {
682 
683  //Raman correction factor - Eq. 11 in Lee et al 2013
684  rFactor = alphaCor[iw] * Rrs443 / Rrs550 + betaCor0[iw] * pow(Rrs550, betaCor1[iw]);
685 
686  rrs_a = l2rec->Rrs[ip * nbands + iw];
687 
688  Rrs_ram[ip * nbands + iw] = rrs_a - rrs_a / (1.0 + rFactor);
689 
690  } else {
691  Rrs_ram[ip * nbands + iw] = 0.0;
692  }
693 
694  }
695 
696 }
697 
698 /* -----------------------------------------------------------------------*/
699 /* raman_cor_westberry() - computes the Rrs corrected for Raman scattering*/
700 /* following the method of Westberry et al 2013 */
701 /* */
702 /* Note: only the first-order scattering term is */
703 /* used as the higher order term often to */
704 /* resulted in: Rrs_raman > Rrs. */
705 /* -----------------------------------------------------------------------*/
706 void raman_cor_westberry(l2str *l2rec, int ip) {
707 
708  //We are only considering the first-order Raman scattering
709 
710  int iw;
711  int nbands = l2rec->l1rec->l1file->nbands;
712  float term0, term1, numer1, denom1;
713  float Rrs_rc;
714 
715 
716  //If pixel is already masked, return and do not attempt the correction.
717  if (l2rec->l1rec->mask[ip]) {
718  return;
719  }
720 
721  //Extrapolate aphStar into the UV Raman excitation bands
722  set_raman_aph_uv(l2rec, ip);
723 
724  //Run QAA to estimate IOP magnitudes and spectral slopes
725  raman_qaa(l2rec, ip);
726 
727  //Get total absorption and scattering coefficients at both sensor
728  //and Raman excitation bands.
730 
731  //Get the k-functions
732  raman_k_func(l2rec, ip);
733 
734  //Get the down-welling irradiances
735  raman_radtran_ed(l2rec, ip);
736 
737  //From Mobley 2010 - assume isotropic emission
738  term0 = 1.0 / (4.0 * M_PI * nWater * nWater);
739 
740 
741  //Loop over visible bands
742  for (iw = 0; iw < nbandVis; iw++) {
743 
744  //Only calculate correction for visible bands
745  if (lambda[iw] >= WAVEMIN && lambda[iw] <= WAVEMAX) {
746 
747  //First term in in Westberry et al. 2013, Eq 7.
748  numer1 = bRex[iw] * eDRam[iw];
749  denom1 = (kdRam[iw]+kappaSen[iw])*eDSen[iw];
750  term1 = numer1 / denom1;
751 
752  //Raman corrected Rrs (first order scattering term only)
753  Rrs_rc = term0 * term1;
754 
755  //Fill Rrs_raman array
756  Rrs_ram[ip * nbands + iw] = Rrs_rc;
757 
758  } else {
759  Rrs_ram[ip * nbands + iw] = 0.0;
760  }
761 
762  }
763 
764 }
765 
766 /* -----------------------------------------------------------------------*/
767 /* raman_cor_lee2() - computes the Rrs corrected for Raman scattering */
768 /* following the method of Lee et al. 1994 */
769 
770 /* -----------------------------------------------------------------------*/
771 void raman_cor_lee2(l2str *l2rec, int ip) {
772 
773  int iw;
774  int nbands = l2rec->l1rec->l1file->nbands;
775  float Rrs_rc;
776 
777  //If pixel is already masked, return and do not attempt the correction.
778  if (l2rec->l1rec->mask[ip]) {
779  return;
780  }
781 
782  //Extrapolate aphStar into the UV Raman excitation bands
783  set_raman_aph_uv(l2rec, ip);
784 
785  //Run QAA to estimate IOP magnitudes and spectral slopes
786  raman_qaa(l2rec, ip);
787 
788  //Get total absorption and scattering coefficients at both sensor
789  //and Raman excitation bands.
791 
792  //Get the down-welling irradiances
793  raman_radtran_ed(l2rec, ip);
794 
795  //Loop over vis bands
796  for (iw = 0; iw < nbandVis; iw++) {
797 
798  if (lambda[iw] >= WAVEMIN && lambda[iw] <= WAVEMAX) {
799  //Raman-corrected Rrs
800  Rrs_rc = 0.072 * (bRex[iw] * eDRam[iw]) / ((2.0 * atot[iw] + atotR[iw]) * eDSen[iw]);
801 
802  //Fill Rrs_raman array
803  Rrs_ram[ip * nbands + iw] = Rrs_rc;
804 
805  } else {
806  Rrs_ram[ip * nbands + iw] = 0.0;
807  }
808 
809  }
810 
811 }
812 
813 /* -----------------------------------------------------------------------*/
814 /* run_raman_cor() - run the Raman scattering correction algorithm */
815 
816 /* -----------------------------------------------------------------------*/
817 void run_raman_cor(l2str *l2rec, int ip) {
818 
819  int iw;
820  static int firstCall = 1;
821  static int sensorSupported;
822  //int32_t ocSensorID;
823 
824  //Number of visible bands
825  nbandVis = rdsensorinfo(l2rec->l1rec->l1file->sensorID, l1_input->evalmask, "NbandsVIS", NULL);
826 
827  //First call
828  if (firstCall) {
829 
830  //Supported sensors list
831  switch (l2rec->l1rec->l1file->sensorID){
832  case MODISA:
833  case MODIST:
834  case SEAWIFS:
835  case VIIRSN:
836  case VIIRSJ1:
837  case VIIRSJ2:
838  case MERIS:
839  case OCTS:
840  case OLCIS3A:
841  case OLCIS3B:
842  case SGLI:
843  case OCIS:
844  case OCI:
845  sensorSupported = 1;
846  break;
847  default:
848  //Force l2gen input to raman_opt=0
849  input->raman_opt = NORAMAN;
850  //Set flag to 1 - supported sensor
851  sensorSupported = 0;
852  }
853 
854 
855  //No Raman correction selected
856  if (input->raman_opt == 0) {
857 
858  printf("\n");
859  if (!sensorSupported) {
860  printf("Raman correction unsupported for this sensor. \n");
861  }
862 
863  printf("No Raman scattering correction calculated for Rrs. \n");
864  printf("\n");
865 
866  //Raman correction selected
867  } else if (input->raman_opt > 0 && input->raman_opt <= 3) {
868 
869  printf("\n");
870  printf("Compute Raman scattering correction #%d. \n", input->raman_opt);
871  printf("\n");
872 
873  //Invalid Raman correction selected
874  } else {
875  printf("-E- %s line %d : '%d' is not a valid Raman correction option.\n",
876  __FILE__, __LINE__, input->raman_opt);
877  exit(1);
878  }
879 
880 
881  /*Allocate memory according to number of sensor bands*/
882  raman_pixel_alloc(l2rec);
883 
884  //If sensor is supported for Raman correction, coefficients are allocated
885  if (sensorSupported) {
886 
887  /*Get the sensor-specific coefficients*/
888  get_raman_coeffs(l2rec);
889 
890  /*Loop over visible bands*/
891  for (iw = 0; iw < nbandVis; iw++) {
892 
893  //Get the visible band centers
894  lambda[iw] = l2rec->l1rec->l1file->fwave[iw];
895 
896  //Interpolate the Lee et al. (2013) correction coefficients to sensor resolution
897  alphaCor[iw] = linterp(lamCor1, alpha, 6, lambda[iw]);
898  betaCor0[iw] = linterp(lamCor1, beta0, 6, lambda[iw]);
899  betaCor1[iw] = linterp(lamCor1, beta1, 6, lambda[iw]);
900 
901  //Populate the pure water absorption data//
902  awR[iw] = aw_spectra(ramLam[iw], ramBandW[iw]);
903  bbwR[iw] = bbw_spectra(ramLam[iw], ramBandW[iw]);
904  aw[iw] = l2rec->l1rec->l1file->aw[iw];
905  bbw[iw] = l2rec->l1rec->l1file->bbw[iw];
906 
907 
908  }
909 
910  //Find the sensor wavelength indices, do this once.
911  idx412 = windex(412.0, lambda, nbandVis);
912  idx443 = windex(443.0, lambda, nbandVis);
913  idx488 = windex(488.0, lambda, nbandVis);
914  idx550 = windex(547.0, lambda, nbandVis);
915  idx670 = windex(667.0, lambda, nbandVis);
916 
917 
918  //How many bands between 412 and 490?
919  for (iw = 0; iw < nbandVis; iw++) {
920  if (lambda[iw] > 550.0) {
921  bandNum550 = iw - 1;
922  break;
923  }
924  }
925 
926  }
927  firstCall = 0;
928  }
929 
930 
931  /*No raman correction is applied, raman_opt=0*/
932  if (input->raman_opt == NORAMAN) {
933  /*Fill Rrs_raman with zeros*/
934  for (iw = 0; iw < l2rec->l1rec->l1file->nbands; iw++) {
935  Rrs_ram[ip * l2rec->l1rec->l1file->nbands + iw] = 0.0;
936  }
937 
938  /*Apply the Lee et al 2013 correction, raman_opt=1*/
939  } else if (input->raman_opt == LEE2013) {
940  raman_cor_lee1(l2rec, ip);
941 
942  /*Apply the Westberry et al 2013 correction, raman_opt=2*/
943  } else if (input->raman_opt == WESTBERRY2013) {
944  raman_cor_westberry(l2rec, ip);
945 
946  /*Apply the Lee et al 1994 Raman correction, raman_opt=3*/
947  } else if (input->raman_opt == LEE1994) {
948  raman_cor_lee2(l2rec, ip);
949  }
950 
951 
952  l2rec->Rrs_raman = Rrs_ram;
953 
954 }
955 
956 /* ---------------------------------------------------------------------------*/
957 /* ---------------------------------------------------------------------------*/
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:315
#define OLCIS3A
Definition: sensorDefs.h:32
#define SGLI
Definition: sensorDefs.h:33
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:93
#define OCI
Definition: sensorDefs.h:42
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
#define NULL
Definition: decode_rs.h:63
read l1rec
#define VIIRSN
Definition: sensorDefs.h:23
void raman_pixel_alloc(l2str *l2rec)
Definition: raman.c:192
unsigned long long sumsq(signed short *in, int cnt)
void set_raman_aph_uv(l2str *l2rec, int ip)
Definition: raman.c:328
#define MERIS
Definition: sensorDefs.h:22
void get_raman_coeffs(l2str *l2rec)
Definition: raman.c:246
#define MODIST
Definition: sensorDefs.h:18
data_t lambda[NROOTS+1]
Definition: decode_rs.h:75
instr * input
#define M_PI
Definition: dtranbrdf.cpp:19
void raman_cor_lee2(l2str *l2rec, int ip)
Definition: raman.c:771
#define WAVEMAX
Definition: raman.c:54
void raman_cor_lee1(l2str *l2rec, int ip)
Definition: raman.c:662
#define WAVEMIN
Definition: raman.c:53
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
#define WESTBERRY2013
Definition: raman.c:51
void run_raman_cor(l2str *l2rec, int ip)
Definition: raman.c:817
float aw_spectra(int32_t wl, int32_t width)
l1_input_t * l1_input
Definition: l1_options.c:9
#define OCIS
Definition: sensorDefs.h:43
void raman_cor_westberry(l2str *l2rec, int ip)
Definition: raman.c:706
float bbw_spectra(int32_t wl, int32_t width)
#define NORAMAN
Definition: raman.c:49
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
#define LEE1994
Definition: raman.c:52
#define OCTS
Definition: sensorDefs.h:14
void raman_k_func(l2str *l2rec, int ip)
Definition: raman.c:521
#define LEE2013
Definition: raman.c:50
void getRamanCoeff(int grpid, char *varName, float *variable)
Definition: raman.c:145
int16_t * numer[MAXNFILES]
Definition: l2bin.cpp:93
int32_t nbands
data_t u
Definition: decode_rs.h:74
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
#define VIIRSJ2
Definition: sensorDefs.h:44
void raman_qaa(l2str *l2rec, int ip)
Definition: raman.c:395
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
const char * subsensorId2SubsensorDir(int subsensorId)
Definition: sensorInfo.c:329
#define OLCIS3B
Definition: sensorDefs.h:41
void raman_radtran_ed(l2str *l2rec, int ip)
Definition: raman.c:550
void fit_raman_tsol(l2str *l2rec, int ip)
Definition: raman.c:365
#define SEAWIFS
Definition: sensorDefs.h:12
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
#define MODISA
Definition: sensorDefs.h:19
#define VIIRSJ1
Definition: sensorDefs.h:37
#define SEAWIFS_LAC
Definition: sensorDefs.h:56
void raman_and_sensor_iops()
Definition: raman.c:488