A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
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 /* -----------------------------------------------------------------------*/
488 void raman_and_sensor_iops() {
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.
729  raman_and_sensor_iops();
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.
790  raman_and_sensor_iops();
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 /* ---------------------------------------------------------------------------*/
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
sensorId2SensorDir
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:315
OLCIS3A
#define OLCIS3A
Definition: sensorDefs.h:32
SGLI
#define SGLI
Definition: sensorDefs.h:33
denom
int16_t * denom[MAXNFILES]
Definition: l2bin.cpp:93
OCI
#define OCI
Definition: sensorDefs.h:42
allocateMemory
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
NULL
#define NULL
Definition: decode_rs.h:63
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
VIIRSN
#define VIIRSN
Definition: sensorDefs.h:23
raman_pixel_alloc
void raman_pixel_alloc(l2str *l2rec)
Definition: raman.c:192
sumsq
unsigned long long sumsq(signed short *in, int cnt)
set_raman_aph_uv
void set_raman_aph_uv(l2str *l2rec, int ip)
Definition: raman.c:328
MERIS
#define MERIS
Definition: sensorDefs.h:22
get_raman_coeffs
void get_raman_coeffs(l2str *l2rec)
Definition: raman.c:246
MODIST
#define MODIST
Definition: sensorDefs.h:18
lambda
data_t lambda[NROOTS+1]
Definition: decode_rs.h:75
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
input
instr * input
Definition: main_l2binmatch.cpp:27
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
l12_proto.h
raman_cor_lee2
void raman_cor_lee2(l2str *l2rec, int ip)
Definition: raman.c:771
WAVEMAX
#define WAVEMAX
Definition: raman.c:54
raman_cor_lee1
void raman_cor_lee1(l2str *l2rec, int ip)
Definition: raman.c:662
WAVEMIN
#define WAVEMIN
Definition: raman.c:53
linterp
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
WESTBERRY2013
#define WESTBERRY2013
Definition: raman.c:51
run_raman_cor
void run_raman_cor(l2str *l2rec, int ip)
Definition: raman.c:817
aw_spectra
float aw_spectra(int32_t wl, int32_t width)
Definition: water_spectra.c:105
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
OCIS
#define OCIS
Definition: sensorDefs.h:43
raman_cor_westberry
void raman_cor_westberry(l2str *l2rec, int ip)
Definition: raman.c:706
bbw_spectra
float bbw_spectra(int32_t wl, int32_t width)
Definition: water_spectra.c:142
NORAMAN
#define NORAMAN
Definition: raman.c:49
filename
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
LEE1994
#define LEE1994
Definition: raman.c:52
OCTS
#define OCTS
Definition: sensorDefs.h:14
raman_k_func
void raman_k_func(l2str *l2rec, int ip)
Definition: raman.c:521
LEE2013
#define LEE2013
Definition: raman.c:50
getRamanCoeff
void getRamanCoeff(int grpid, char *varName, float *variable)
Definition: raman.c:145
numer
int16_t * numer[MAXNFILES]
Definition: l2bin.cpp:93
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
u
data_t u
Definition: decode_rs.h:74
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
VIIRSJ2
#define VIIRSJ2
Definition: sensorDefs.h:44
raman_qaa
void raman_qaa(l2str *l2rec, int ip)
Definition: raman.c:395
rdsensorinfo
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
color_dtdb.alpha
alpha
Definition: color_dtdb.py:352
subsensorId2SubsensorDir
const char * subsensorId2SubsensorDir(int subsensorId)
Definition: sensorInfo.c:329
OLCIS3B
#define OLCIS3B
Definition: sensorDefs.h:41
raman_radtran_ed
void raman_radtran_ed(l2str *l2rec, int ip)
Definition: raman.c:550
fit_raman_tsol
void fit_raman_tsol(l2str *l2rec, int ip)
Definition: raman.c:365
SEAWIFS
#define SEAWIFS
Definition: sensorDefs.h:12
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
MODISA
#define MODISA
Definition: sensorDefs.h:19
VIIRSJ1
#define VIIRSJ1
Definition: sensorDefs.h:37
SEAWIFS_LAC
#define SEAWIFS_LAC
Definition: sensorDefs.h:56
raman_and_sensor_iops
void raman_and_sensor_iops()
Definition: raman.c:488


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster