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
get_sdp.c
Go to the documentation of this file.
1 /*=========================================================================*/
2 /* get_sdp.c */
3 /* */
4 /* Description: */
5 /* This algorithm calculates multiple phytoplankton pigment concentrations */
6 /* based on spectral derivative of hyperspectral radiometry */
7 /* */
8 /* */
9 /* Notes: */
10 /* */
11 /* References: */
12 /* Kramer, S., Siegel, D.A., Maritorena, S. and D. Catlett. (2022) Modeling*/
13 /* surface ocean phytoplankton pigments from hyperspectral remote sensing */
14 /* reflectance on global scales, Remote Sens. Environ., 112879 */
15 /* doi: 10.1016/j.rse.2021.112879. */
16 /* */
17 /* Implementation: */
18 /* L. McKinna (Go2Q), 19 Dec 2022 */
19 /*=========================================================================*/
20 
21 #include <math.h>
22 #include <gsl/gsl_interp.h>
23 #include <gsl/gsl_filter.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_statistics.h>
26 #include <levmar.h>
27 #include "l12_proto.h"
28 #include "get_sdp.h"
29 
30 static int32_t LastRecNum = -1;
31 
32 static int npig = 13; //Number of pigments in model
33 static int nAcoeff = 299; //Number of model coefficents
34 static int nreps = 100; //Number of replicates in model coefficients
35 static char *pignames[13] = {"tchla","zea","dvchla","butfuco","hexfuco","allo",
36  "mvchlb","neo","viola","fuco","chlc12","chlc3","perid"};
37 
38 
39 static int16 *sdpflags; //Binary quality control flags
40 static float *pigArr; //Pigment array pointer (this is where we'll store the solutions)
41 static float *pigArr_unc; //Pigment uncertinty array pointer (this is where we'll store the
42  //uncertainty in solutions)
43 static float *mRrs_out; //Array holding forward-modelled Rrs
44 static float *Rrs_diff_out; //Array holdig difference between modeled and observed Rrs.
45 static float *d2Rrs_out; //Array holding 2nd Derivative of RrsDiff
46 
47 
48 //static int nfitbands = 6;
49 //static float fitbands[6] = {412,443,490,530,555,670};
50 
51 //static int nfitbands = 12;
52 //static float fitbands[12] = {400, 425, 450, 475, 500, 525, 550, 575, 600, 625, 650, 675};
53 
54 static int nfitbands = 56;
55 static float fitbands[56] = {400,405,410,415,420,425,430,435,440,445,450,455,460,465,470,475,
56  480,485,490,495,500,505,510,515,520,525,530,535,540,545,550,555,
57  560,565,570,575,580,585,590,595,600,605,610,615,620,625,630,635,
58  640,645,650,655,660,665,670,675};
59 
60 
61 void sdp_qssa( double *p, double *mrrs_s, int m, int n, void *data);
62 
63 /* ------------------------------------------------------------------- */
64 /*sdp_ran() - record if run_sdp has been executed for scan line */
65 
66 /* ------------------------------------------------------------------- */
67 int sdp_ran(int recnum) {
68  if (recnum == LastRecNum)
69  return 1;
70  else
71  return 0;
72 }
73 
74 /* -----------------------------------------------------------------*/
75 /* aph_kramer() - get A and B model coefficients for modelling aph */
76 /* as a function of chl */
77 
78 /****************Note: may moved into aph.c after testing************/
79 
80 /* -----------------------------------------------------------------*/
81 int aph_kramer_2022(float wave, float *A, float *B) {
82 
83  static int firstCall = 1;
84 
85  static float *table;
86  static float *wtab;
87  static float *Atab, *Btab;
88  static float *dAtab, *dBtab;
89  static int ntab = 0;
90 
91 
92  if (firstCall) {
93 
94  char *filedir;
95  char filename[FILENAME_MAX];
96  int ncol, nrow;
97 
98  if ((filedir = getenv("OCDATAROOT")) == NULL) {
99  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
100  exit(1);
101  }
102  strcpy(filename, filedir);
103  strcat(filename, "/common/aph_kramer_2022.txt");
104  printf("Reading aph* for sdp model from %s.\n", filename);
105 
106  nrow = table_row_count(filename);
107 
108  if (nrow <= 0) {
109  printf("-E- %s line %d: error opening (%s) file", __FILE__, __LINE__, filename);
110  exit(1);
111  }
112 
113  ncol = table_column_count(filename);
114  table = table_read_r4(filename, ncol, nrow);
115  ntab = nrow;
116 
117  wtab = &table[nrow * 0]; //wavelength
118  Atab = &table[nrow * 1]; //Spectral A coefficient
119  Btab = &table[nrow * 2]; //Spectral B coefficent
120 
121  // precompute derivatives for spline interpolation
122  if ((dAtab = (float*) calloc(ntab, sizeof (float))) == NULL) {
123  printf("-E- %s line %d : error allocating memory for Kramer aph.\n",
124  __FILE__, __LINE__);
125  exit(1);
126  }
127  if ((dBtab = (float*) calloc(ntab, sizeof (float))) == NULL) {
128  printf("-E- %s line %d : error allocating memory for Kramer aph.\n",
129  __FILE__, __LINE__);
130  exit(1);
131  }
132  spline(wtab, Atab, ntab, 1e30, 1e30, dAtab);
133  spline(wtab, Btab, ntab, 1e30, 1e30, dBtab);
134 
135  firstCall = 0;
136  }
137 
138  if (wave > wtab[ntab - 1]){
139  wave = wtab[ntab - 1];
140  }
141 
142  // interpolate coefficients to wavelength to get A and B
143  splint(wtab, Atab, dAtab, ntab, wave, A);
144  splint(wtab, Btab, dBtab, ntab, wave, B);
145 
146  return 0;
147 
148 }
149 
150 /* ------------------------------------------------------------------*/
151 /* get_sdp_aw() - read netcf file and return purewater absorption */
152 /* coefficients */
153 
154 int get_sdp_aw(float wave, float *A) {
155 
156  static int firstCall = 1;
157 
158  int32_t status, netcdf_input;
159  size_t nlam;
160  int ncid, grpid, varid1, varid2, dimid ;
161 
162  static float *wtab;
163  static float *atab;
164  static float *dAtab;
165  static int ntab = 0;
166 
167  if (firstCall) {
168 
169  char *filedir;
170  char filename[FILENAME_MAX];
171 
172  if ((filedir = getenv("OCDATAROOT")) == NULL) {
173  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
174  exit(1);
175  }
176 
177  strcpy(filename, filedir);
178  strcat(filename, "/common/sdp_pure_water.nc");
179  printf("Reading SDP pure water coefficients from %s.\n", filename);
180 
181  /* Does the file exist? */
182  if (access(filename, F_OK) || access(filename, R_OK)) {
183  printf("-E- %s: SDP pure water coefficient file '%s' does not exist or cannot open.\n",
184  __FILE__, filename);
185  exit(EXIT_FAILURE);
186  }
187 
188  /* test for NetCDF input file */
189  status = nc_open(filename, NC_NOWRITE, &ncid);
190  netcdf_input = (status == NC_NOERR);
191 
192  if (netcdf_input) {
193 
194  /* Get the group ids of our two groups. */
195  if ((nc_inq_ncid(ncid, "water_absorption", &grpid)) == NC_NOERR) {
196 
197  //Get the dimensions of the variable
198  nc_inq_dimid(grpid,"n_wavelength", &dimid);
199  nc_inq_dimlen(grpid, dimid, &nlam);
200  ntab = (int) nlam;
201 
202  //Allocate memory
203  atab = (float*) allocateMemory(ntab * sizeof (float), "atab");
204  wtab = (float*) allocateMemory(ntab * sizeof (float), "wtab");
205  dAtab = (float*) allocateMemory(ntab * sizeof (float), "dAtab");
206 
207  nc_inq_varid(grpid, "aw", &varid1);
208  nc_get_var_float(grpid, varid1, atab);
209 
210  nc_inq_varid(grpid, "wavelength", &varid2);
211  nc_get_var_float(grpid, varid2, wtab);
212 
213  }
214 
215 
216  spline(wtab, atab, nlam, 1e30, 1e30, dAtab);
217 
218  firstCall = 0;
219 
220  }
221  }
222 
223  if (wave > wtab[ntab - 1]){
224  wave = wtab[ntab - 1];
225  }
226 
227  // interpolate coefficients to wavelength to get A
228  splint(wtab, atab, dAtab, ntab, wave, A);
229 
230  return 0;
231 
232 }
233 
234 
235 /* -----------------------------------------------------------------*/
236 /* get_sdp_coeff() - read netcf file and return the A and C model */
237 /* coefficients */
238 /* input: */
239 /* pigment name (char) */
240 /* output: */
241 /* A (float array) */
242 /* C (float) */
243 
244 /* -----------------------------------------------------------------*/
245 int get_sdp_coeff(double ** A, double **C, char **pigments, int n, int m) {
246 
247  static int firstCall = 1;
248 
249  int32_t i, status, netcdf_input;
250  int ncid, varidA, varidC, grpidA, grpidC;
251 
252  char varNameA[16];
253  char varNameC[16];
254 
255  if (firstCall) {
256 
257  char *filedir;
258  char filename[FILENAME_MAX];
259 
260  if ((filedir = getenv("OCDATAROOT")) == NULL) {
261  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
262  exit(1);
263  }
264 
265  strcpy(filename, filedir);
266  strcat(filename, "/common/sdp_AC_coeffs_all.nc");
267  printf("Reading SDP model coefficients from %s.\n", filename);
268 
269  /* Does the file exist? */
270  if (access(filename, F_OK) || access(filename, R_OK)) {
271  printf("-E- %s: SDP coefficient file '%s' does not exist or cannot open.\n",
272  __FILE__, filename);
273  exit(EXIT_FAILURE);
274  }
275 
276  /* test for NetCDF input file */
277  status = nc_open(filename, NC_NOWRITE, &ncid);
278  netcdf_input = (status == NC_NOERR);
279 
280  if (netcdf_input) {
281 
282  for (i=0; i<npig; i++) {
283 
284  strncpy(varNameA, "A",16);
285  strncat(varNameA,pigments[i],16);
286  strncpy(varNameC, "C",16);
287  strncat(varNameC,pigments[i],16);
288 
289  /* Get the group ids of our two groups. */
290  if ((nc_inq_ncid(ncid, "A_coeff", &grpidA)) == NC_NOERR) {
291  if ((nc_inq_varid(grpidA, varNameA, &varidA)) == NC_NOERR) {
292  /* Read the coefficeint data. */
293  nc_get_var_double(grpidA, varidA, A[i]);
294  }
295  }
296 
297  if (( nc_inq_ncid(ncid, "C_coeff", &grpidC)) == NC_NOERR) {
298  if ((nc_inq_varid(grpidC, varNameC, &varidC)) == NC_NOERR) {
299  nc_get_var_double(grpidC, varidC, C[i]);
300  }
301  }
302  }
303  }
304 
305  firstCall = 0;
306  }
307 
308  return 0;
309 }
310 
311 
312 /* ------------------------------------------------------------------- */
313 /* sdp_interp() - linear 1D interpolation function */
314 
315 /* Inputs */
316 /* x: initial x array */
317 /* y: intial y array */
318 /* xinterp: new x array */
319 /* n: length or x and y */
320 /* m: length of xinterp and yinterp */
321 /* */
322 /* Output */
323 /* yinterp: new y array */
324 
325 /* ------------------------------------------------------------------- */
326 
327 int sdp_interp(float *x, float *y, float *xinterp, float *yinterp,
328  int m, int n) {
329  int iw;
330  double *xd;
331  double *yd;
332 
333  xd = (double*)calloc(n, sizeof(double));
334  yd = (double*)calloc(n, sizeof(double));
335 
336  gsl_interp *interp_work;
337 
338  for (iw=0; iw < n; iw++) {
339  xd[iw] = x[iw];
340  yd[iw] = y[iw];
341  }
342 
343 
344  interp_work = gsl_interp_alloc(gsl_interp_linear, m);
345  gsl_interp_init(interp_work, xd, yd, m);
346 
347  for (iw=0; iw < n; iw++) {
348  yinterp[iw] = (float) gsl_interp_eval(interp_work, xd, yd, (double) xinterp[iw], NULL);
349  }
350 
351  gsl_interp_free(interp_work);
352  free(xd);
353  free(yd);
354  return 0;
355 }
356 
357 /* ------------------------------------------------------------------- */
358 /* sdp_mean_filter() - simple mean moving average filter */
359 
360 /* ------------------------------------------------------------------- */
361 
362 int sdp_mean_filter(float *y, float *yfilt, int n, int window) {
363 
364  int iw, ih, halfwin;
365  float temp;
366 
367  if(window % 2 == 0) {
368  printf("-E- %s line %d : mean moving average filter window must be odd value.\n",
369  __FILE__, __LINE__);
370  exit(1);
371  }
372 
373  halfwin = (window-1)/2;
374 
375  for (iw=0; iw< n; iw++) {
376 
377  if (iw <= halfwin) {
378  temp = y[iw];
379  } else if (iw >= halfwin) {
380  temp = y[iw];
381  } else {
382  for (ih=0;ih<window;ih++) {
383  temp += y[iw+ih];
384  }
385  temp /= window;
386  }
387  yfilt[iw] = temp;
388  }
389  return 0;
390 }
391 
392 /* ------------------------------------------------------------------- */
393 /* sdp_filter() - simple median filter */
394 
395 /* ------------------------------------------------------------------- */
396 int sdp_filter(float *y, float *yfilt, int n, int window) {
397 
398  int iw;
399  int success;
400  float temp;
401 
402  gsl_filter_median_workspace *medfilt;
403  medfilt = gsl_filter_median_alloc(window);
404 
405  gsl_vector *yin = gsl_vector_alloc(n);
406  gsl_vector *yout = gsl_vector_alloc(n);
407 
408  for (iw=0; iw<n; iw++) {
409  gsl_vector_set(yin, iw, y[iw]);
410  gsl_vector_set(yout, iw, 0.0);
411  }
412 
413  success= gsl_filter_median(GSL_FILTER_END_PADVALUE, yin, yout, medfilt);
414 
415  for (iw=0; iw<n; iw++) {
416  temp = gsl_vector_get(yout, iw) ;
417  yfilt[iw] = temp;
418  }
419 
420  //cleanup
421  gsl_filter_median_free(medfilt);
422  gsl_vector_free(yin);
423  gsl_vector_free(yout);
424  return success;
425 }
426 
427 /* ------------------------------------------------------------------- */
428 /* sdp_deriv2() - numerical second derivative */
429 
430 /* ------------------------------------------------------------------- */
431 
432 int sdp_deriv2(float *x, float *y, float *deriv2, int n) {
433 
434  int i;
435  float d1, d2;
436 
437  //zero pad first and last element of output
438  deriv2[0] =0.0;
439  deriv2[n] = 0.0;
440 
441  for (i=1; i<=n-1; i++) {
442 
443  //1st derative at two points
444  d1 = (y[i] - y[i-1]) /(x[i]-x[i-1]); //1st Deriv at point 1
445  d2 = (y[i+1] - y[i])/(x[i+1]-x[i]); //1st Deriv at point 2
446 
447  //Compute 2nd deriv from d1 and d2
448  deriv2[i] = (d2 - d1)/ ((x[i+1] - x[i-1])/2.);
449 
450 
451  }
452 
453  return 0;
454 }
455 
456 /* ------------------------------------------------------------------- */
457 
458 
459 /* ------------------------------------------------------------------- */
460 /* sdp_qssa() - Quasi-single scattering approximation (QSSA) forward */
461 /* reflectance model */
462 
463 /* ------------------------------------------------------------------- */
464 
465 
466 void sdp_qssa( double *p, double *mrrs_s, int m, int n, void *data){
467 
468  sdpstr *s = ((sdpstr *)data);
469 
470  int iw;
471  double atot, bbtot, u;
472  double aph, adg, bbp;
473 
474  for (iw = 0; iw < n; iw++) {
475 
476  //Calculate the total IOPs
477  aph = s->aphA[iw]*pow(p[0],s->aphB[iw]);
478  adg = p[1]*s->adgstar[iw];
479  bbp = p[2]*s->bbpstar[iw];
480  atot = s->aw[iw] + aph + adg;
481  bbtot = s->bbw[iw] + bbp;
482 
483 
484  u = bbtot/(atot + bbtot);
485  mrrs_s[iw] = (s->g[0] + s->g[1]*u)*u;
486 
487  }
488 }
489 
490 /* ------------------------------------------------------------------- */
491 /* sdp_qssaf() - Quasi-single scattering approximation (QSSA) forward */
492 /* reflectance model at fine resolution 1nm */
493 
494 /* ------------------------------------------------------------------- */
495 
496 
497 void sdp_qssaf( double *p, double *mrrs_s, int m, int n, void *data){
498 
499  sdpstr *s = ((sdpstr *)data);
500 
501  int iw;
502  double atot, bbtot, u;
503  double aph, adg, bbp;
504 
505  for (iw = 0; iw < n; iw++) {
506 
507  //Calculate the total IOPs
508  aph = s->faphA[iw]*pow(p[0],s->faphB[iw]);
509  adg = p[1]*s->fadgstar[iw];
510  bbp = p[2]*s->fbbpstar[iw];
511  atot = s->faw[iw] + aph + adg;
512  bbtot = s->fbbw[iw] + bbp;
513 
514  u = bbtot/(atot + bbtot);
515  mrrs_s[iw] = (s->g[0] + s->g[1]*u)*u;
516  }
517 }
518 
519 /* ------------------------------------------------------------------- */
520 /* cal_adstar() - compute shape of adgstar */
521 
522 
523 /*two options, 0 for spectral fit model, 1 for 1nm forward model */
524 /* ------------------------------------------------------------------- */
525 void calc_adgstar(sdpstr *s, int opt) {
526 
527  int iw;
528 
529  switch (opt) {
530  case 0:
531  for (iw = 0; iw < s->nfitbands; iw++) {
532  /*Colored dissolved and detrital matter absorption shape*/
533  s->adgstar[iw] = exp((s->bands[iw] - s->bands[s->i440])*s->adg_s);
534  }
535  break;
536  case 1:
537  for (iw = 0; iw < s->nfbands; iw++) {
538  /*Colored dissolved and detrital matter absorption shape*/
539  s->fadgstar[iw] = exp((s->fbands[iw] - s->fbands[s->i440f])*s->adg_s);
540  }
541  break;
542  default:
543  printf(
544  "-E- %s, %d: Spectral Derivative Pigment adgstar "
545  "model invalid input option value\n",
546  __FILE__, __LINE__);
547  exit(1);
548  break;
549  }
550 }
551 
552 
553 /* ------------------------------------------------------------------- */
554 /* cal_bbpstar() - compute shape of adgstar */
555 
556 /*two options, 0 for spectral fit model, 1 for 1nm forward model */
557 /* ------------------------------------------------------------------- */
558 void calc_bbpstar(sdpstr *s, int opt) {
559 
560  int iw;
561 
562  switch (opt) {
563  case 0:
564  for (iw = 0; iw < s->nfitbands; iw++) {
565  /*Particle backscattering shape*/
566  s->bbpstar[iw] = pow((s->bands[iw]/s->bands[s->i440] ),s->bbp_s);
567  }
568  break;
569  case 1:
570  for (iw = 0; iw < s->nfbands; iw++) {
571  /*Particle backscattering shape*/
572  s->fbbpstar[iw] = pow(( s->fbands[iw]/s->fbands[s->i440f] ),s->bbp_s);
573  }
574  break;
575  default:
576  printf(
577  "-E- %s, %d: Spectral Derivative Pigment bbpstar "
578  "model invalid input option value\n",
579  __FILE__, __LINE__);
580  exit(1);
581  break;
582  }
583 }
584 
585 
586 /* ------------------------------------------------------------------- */
587 /* sdp_jacqssa() - Compute Jacobian of cost function */
588 /* e = y - y(x) */
589 /* dedx = - dy(x)/dx */
590 
591 /* dimension of jac is: npar x nbands */
592 
593 /* ------------------------------------------------------------------- */
594 void sdp_jacqssa(double *p, double *jac, int n, int m, void* data) {
595 
596  sdpstr *s = ((sdpstr *)data);
597 
598  int i, j;
599  double atot, bbtot, aph, adg, bbp, u;
600  double A, B;
601  double drdchl;
602  double drdadg;
603  double drdbbp;
604 
605  for (i=j=0; i< m; i++) {
606 
607  aph = s->aphA[i]*pow(p[0], s->aphB[i]);
608  adg = p[1] * s->adgstar[i];
609  bbp = p[2] * s->bbpstar[i];
610  atot = s->aw[i] + aph + adg;
611  bbtot = s->bbw[i] + bbp;
612 
613  u = bbtot/(atot + bbtot);
614  A = s->aphA[i];
615  B = s->aphB[i];
616 
617  drdchl = -(s->g[0] + 2.*s->g[1]*u)*(bbtot*A*B*pow(p[0],(B-1.))) / pow( bbtot + atot,2.);
618  drdadg = -(s->g[0] + 2.*s->g[1]*u)*(bbtot*s->adgstar[i])/ pow( bbtot + atot, 2.);
619  drdbbp = (s->g[0] + 2.*s->g[1]*u)*(atot*s->bbpstar[i]) / pow( bbtot + atot, 2.);
620 
621  jac[j++] = drdchl;
622  jac[j++] = drdadg;
623  jac[j++] = drdbbp;
624 
625  }
626 }
627 
628 
629 /* ------------------------------------------------------------------- */
630 /* rund sdp_saa() - runs sdp semianalytical algorithm to get to mRrs */
631 
632 /* ------------------------------------------------------------------- */
633 int run_sdp_saa(sdpstr *s) {
634 
635  int iw;
636  int statusLM;
637  int maxit = 500;
638 
639  double p[3] = {0.01, 0.1, 0.0005};
640  double info[LM_INFO_SZ];
641  double opts[LM_OPTS_SZ];
642 
643  double *rrs_s;
644  rrs_s = (double *) calloc(s->nfitbands, sizeof (double));
645 
646  //Intialize output solution array
647  for (iw=0;iw<s->nfree;iw++){
648  s->popt[iw] = BAD_FLT;
649  }
650 
651  //Compute subsurface remote-sensing reflectance
652  for (iw = 0; iw < s->nfitbands; iw++) {
653  rrs_s[iw] = s->Rrs_a[iw]/(0.52 + 1.7*s->Rrs_a[iw]);
654  }
655 
656  //model slope of bbp
657  s->bbp_s = 2.0 * (1.0 - 1.2 * expf(-0.9 * (rrs_s[s->i440]/rrs_s[s->i555])));
658 
659  //model slope of adg
660  s->adg_s = -(0.01447 + 0.00033*s->Rrs_a[s->i490]/s->Rrs_a[s->i555]);
661 
662  //Compute the spectral shapes adstar and bbpstar
663  calc_adgstar(s, 0);
664  calc_bbpstar(s, 0);
665 
666  /*Optimisation control parameters*/
667  opts[0] = 1E-3; //1E-3 LM_INIT_MU; /* scale factor for initial mu */
668  opts[1] = 1E-10; //1E-10; convergence in gradient
669  opts[2] = 1E-5; //1E-5 relative error desired in the approximate solution
670  opts[3] = 1E-7; //1E-7 relative error desired in the sum of squares
671 
672  /*Run L-M optimization with analytical jacobian*/
673  statusLM = dlevmar_der(sdp_qssa, sdp_jacqssa, p, rrs_s, s->nfree,
674  s->nfitbands, maxit, opts, info, NULL, NULL, (void *) s);
675 
676  free(rrs_s);
677 
678  //return solution if one exists
679  if (statusLM <= 0) {
680  return 1;
681  } else {
682  for (iw=0;iw<s->nfree;iw++){
683  s->popt[iw] = p[iw];
684  }
685  return 0;
686  }
687 }
688 
689 /* ------------------------------------------------------------------- */
690 /* check_l2_flags() - checks to see if l2 flags are set that indicate */
691 /* poor quality Rrs */
692 
693 /* Flags checked are: ATMFAIL(1), LAND (2), HIGLINT(4), PRODFAIL (31) */
694 /* STRAYLIGHT (9), CLOUD (10) */
695 
696 
697 int check_l2_flags(l2str *l2rec, int ip) {
698 
699  int i;
700  int nflag = 6;
701  int bit_position[6] = {1,2,4,9,10,31};
702 
703  int status = 0;
704 
705  //Check bit position in l2flags
706  for (i=0;i<nflag;i++) {
707  if (l2rec->l1rec->flags[ip] & (1<<bit_position[i])) {
708  status = 1;
709  break;
710  } else {
711  status = 0;
712  }
713  }
714 
715  return status;
716 }
717 
718 /* ------------------------------------------------------------------- */
719 /* set_sdp_flags() - checks range of derived pigments and sets quality */
720 /* flags */
721 
722 int set_sdp_flags(float pigvalue, int pigid, int ip, int16 *flag) {
723 
724  if (pigvalue < pigmin || pigvalue > pigmax) {
725 
726  if (pigid == 0) {
727  *flag |= BADTCHL;
728  }
729  if (pigid == 1) {
730  *flag |= BADZEA;
731  }
732  if (pigid == 2) {
733  *flag |= BADDVCHLA;
734  }
735  if (pigid == 3) {
736  *flag |= BADBUTOFUCO;
737  }
738  if (pigid == 4) {
739  *flag |= BADHEXOFUCO;
740  }
741  if (pigid == 5) {
742  *flag |= BADALLO;
743  }
744  if (pigid == 6) {
745  *flag |= BADMVCHLB;
746  }
747  if (pigid == 7) {
748  *flag |= BADNEO;
749  }
750  if (pigid == 8) {
751  *flag |= BADVIOLA;
752  }
753  if (pigid == 9) {
754  *flag |= BADFUCO;
755  }
756  if (pigid == 10) {
757  *flag |= BADCHLC12 ;
758  }
759  if (pigid == 11) {
760  *flag |= BASCHLC3;
761  }
762  if (pigid == 12) {
763  *flag |= BADPERID;
764  }
765 
766  pigvalue = BAD_FLT;
767 
768  return 1;
769 
770  } else {
771 
772  return 0;
773  }
774 }
775 
776 
777 /*----------------------------------------------------------------------*/
778 /* extract_band_3d() - utility function to extract selected band subset */
779 
780 /*----------------------------------------------------------------------*/
781 static void extract_band_3d(float *in_buf, float *out_buf, int numPixels, int numBands) {
782  float * out_ptr = out_buf;
783  for (int pix = 0; pix < numPixels; pix++) {
784  float *in_ptr = in_buf + pix * numBands;
785  for (int band_3d = 0; band_3d < input->nwavelengths_3d; band_3d++) {
786  int band = input->wavelength_3d_index[band_3d];
787  *out_ptr = in_ptr[band];
788  out_ptr++;
789  }
790  }
791 }
792 
793 
794 /* ------------------------------------------------------------------- */
795 /* run_sdp() - computes the pigment concentrations from 2nd derivativ */
796 /* of the model residuals */
797 
798 /* ------------------------------------------------------------------- */
799 
800 int run_sdp(l2str *l2rec) {
801 
802  static int firstCall = 1;
803 
804  static sdpstr sdpdata;
805  static sdpstr* s = &sdpdata;
806 
807  int iw, iwx, ig, ip, irep, ipb, iwxf;
808  int i400, i700, badrrs;
809 
810  float finewave = 1.0;
811  float band_temp;
812  double pigtemp_med;
813  double pigtemp_std;
814  double pigtemp;
815 
816  l1str *l1rec = l2rec->l1rec;
817  float *bands = l1rec->l1file->fwave;
818  int32_t npix = l1rec->npix;
819  int32_t nbands = l1rec->l1file->nbands;
820 
821  static float *vbands;
822  static double *fmrrs_s;
823  static float *fmRrs_a;
824  static float *RrsDiff;
825  static float *deriv2;
826  static float *fRrs_a;
827  static float *fsmthRrs_a;
828  static float *Rrs_a;
829  static double *pigtempArr;
830 
831  if (firstCall) {
832 
833  firstCall = 0;
834 
835  //Check sensor compatability. Must be hyperspectral sensor.
836  switch (l1rec->l1file->sensorID) {
837  case OCI:
838  break;
839  case OCIS:
840  break;
841  case HICO:
842  break;
843  case PRISM:
844  break;
845  default:
846  printf(
847  "-E- %s, %d: Spectral Derivative Pigment model only applicable to hypersectral data\n",
848  __FILE__, __LINE__);
849  exit(1);
850  break;
851  }
852 
853  s->nfree = 3;
854  i400 = windex(400, bands, nbands);
855  i700 = windex(700, bands, nbands);
856  s->nvbands = i700 - i400 + 1; //Number of sensor visible bands
857  s->nfitbands = nfitbands; //Number of fit bands
858  s->nfbands = (705 - 395 + 1)/finewave; //Number of fine resolution bands
859 
860  //Allocate output arrays
861  if ((pigArr =(float *) calloc(npix * npig, sizeof (float))) == NULL) {
862  printf("-E- %s line %d : error allocating memory for SDP.\n",
863  __FILE__, __LINE__);
864  exit(1);
865  }
866  if ((pigArr_unc =(float *) calloc(npix * npig, sizeof (float))) == NULL) {
867  printf("-E- %s line %d : error allocating memory for SDP.\n",
868  __FILE__, __LINE__);
869  exit(1);
870  }
871  //Allocate output arrays
872  if ((sdpflags =(int16 *) calloc(npix, sizeof (int16))) == NULL) {
873  printf("-E- %s line %d : error allocating memory for SDP.\n",
874  __FILE__, __LINE__);
875  exit(1);
876  }
877  if ((vbands = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
878  printf("-E- %s line %d : error allocating memory for SDP.\n",
879  __FILE__, __LINE__);
880  exit(1);
881  }
882  if ((s->g = (float *) calloc(2, sizeof (float))) == NULL) {
883  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
884  __FILE__, __LINE__);
885  exit(1);
886  }
887  if ((s->bindx = (int *) calloc(s->nvbands, sizeof (int))) == NULL) {
888  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
889  __FILE__, __LINE__);
890  exit(1);
891  }
892  if ((s->fbindx = (int *) calloc(nfitbands, sizeof (int))) == NULL) {
893  printf("-E- %s line %d : error allocating memory in init_iop_flag.\n",
894  __FILE__, __LINE__);
895  exit(1);
896  }
897  if ((s->bands = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
898  printf("-E- %s line %d : error allocating memory for SDP.\n",
899  __FILE__, __LINE__);
900  exit(1);
901  }
902  if ((s->fbands = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
903  printf("-E- %s line %d : error allocating memory for SDP.\n",
904  __FILE__, __LINE__);
905  exit(1);
906  }
907  if ((s->aw = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
908  printf("-E- %s line %d : error allocating memory for SDP.\n",
909  __FILE__, __LINE__);
910  exit(1);
911  }
912  if ((s->bbw = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
913  printf("-E- %s line %d : error allocating memory for SDP.\n",
914  __FILE__, __LINE__);
915  exit(1);
916  }
917  if ((s->aphA = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
918  printf("-E- %s line %d : error allocating memory for SDP.\n",
919  __FILE__, __LINE__);
920  exit(1);
921  }
922  if ((s->aphB = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
923  printf("-E- %s line %d : error allocating memory for SDP.\n",
924  __FILE__, __LINE__);
925  exit(1);
926  }
927  if ((s->faw = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
928  printf("-E- %s line %d : error allocating memory for SDP.\n",
929  __FILE__, __LINE__);
930  exit(1);
931  }
932  if ((s->fbbw = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
933  printf("-E- %s line %d : error allocating memory for SDP.\n",
934  __FILE__, __LINE__);
935  exit(1);
936  }
937  if ((s->faphA = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
938  printf("-E- %s line %d : error allocating memory for SDP.\n",
939  __FILE__, __LINE__);
940  exit(1);
941  }
942  if ((s->faphB = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
943  printf("-E- %s line %d : error allocating memory for SDP.\n",
944  __FILE__, __LINE__);
945  exit(1);
946  }
947  if ((s->adgstar = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
948  printf("-E- %s line %d : error allocating memory for SDP.\n",
949  __FILE__, __LINE__);
950  exit(1);
951  }
952  if ((s->bbpstar = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
953  printf("-E- %s line %d : error allocating memory for SDP.\n",
954  __FILE__, __LINE__);
955  exit(1);
956  }
957  if ((s->fadgstar = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
958  printf("-E- %s line %d : error allocating memory for SDP.\n",
959  __FILE__, __LINE__);
960  exit(1);
961  }
962  if ((s->fbbpstar = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
963  printf("-E- %s line %d : error allocating memory for SDP.\n",
964  __FILE__, __LINE__);
965  exit(1);
966  }
967  if ((s->Rrs_a = (float *) calloc(s->nvbands, sizeof (float))) == NULL) {
968  printf("-E- %s line %d : error allocating memory for SDP.\n",
969  __FILE__, __LINE__);
970  exit(1);
971  }
972  if ((s->popt = (double*) calloc(s->nfree, sizeof (double))) == NULL) {
973  printf("-E- %s line %d : error allocating memory for SDP.\n",
974  __FILE__, __LINE__);
975  exit(1);
976  }
977  if ((s->popt = (double*) calloc(s->nfree, sizeof (double))) == NULL) {
978  printf("-E- %s line %d : error allocating memory for SDP.\n",
979  __FILE__, __LINE__);
980  exit(1);
981  }
982  if ((pigtempArr = (double*) calloc(nreps, sizeof (double))) == NULL) {
983  printf("-E- %s line %d : error allocating memory for SDP.\n",
984  __FILE__, __LINE__);
985  exit(1);
986  }
987  if ((mRrs_out = calloc(npix * nbands, sizeof (float))) == NULL) {
988  printf("-E- %s line %d : error allocating memory for SDP.\n",
989  __FILE__, __LINE__);
990  exit(1);
991  }
992  if ((Rrs_diff_out = calloc(npix * nbands, sizeof (float))) == NULL) {
993  printf("-E- %s line %d : errorn allocating memory for SDP.\n",
994  __FILE__, __LINE__);
995  exit(1);
996  }
997  if ((d2Rrs_out = calloc(npix * nbands, sizeof (float))) == NULL) {
998  printf("-E- %s line %d : error allocating memory for SDP.\n",
999  __FILE__, __LINE__);
1000  exit(1);
1001  }
1002 
1003 
1004  //2D array of PCA coefficients (13 x 29000). Replicates are interleved in each row.
1005  s->Acoeff = (double **) allocate2d_double(npig, nreps*nAcoeff);
1006 
1007  //2D rray of PCA intercept coefficients (13 x 100)
1008  s->Ccoeff = (double **) allocate2d_double(npig, nreps);
1009 
1010 
1011  //Get the indices for the model fit bands
1012  for (iw = 0; iw < s->nfitbands; iw++) {
1013  iwx = windex(fitbands[iw], bands, nbands);
1014  s->bindx[iw] = iwx;
1015  }
1016 
1017 
1018  //Populate model spectral constants at model fit bands
1019  for (iw=0; iw<s->nfitbands; iw++) {
1020 
1021  //Subset of the model fit bands
1022  s->bands[iw] = bands[s->bindx[iw]];
1023 
1024  //Get phytoplankton absorption model coefficients
1025  aph_kramer_2022(s->bands[iw], &s->aphA[iw], &s->aphB[iw]);
1026 
1027  //get pure water absorption coefficients at fit bands
1028  get_sdp_aw(s->bands[iw], &s->aw[iw]) ;
1029 
1030  }
1031 
1032  //Compute reflectance model coefficients at fine resolution (400 - 700 nm, 1nm res)
1033  for (iw=0;iw<s->nfbands;iw++) {
1034 
1035  if (!iw){
1036  s->fbands[0] = 395;
1037  } else {
1038  s->fbands[iw] += s->fbands[iw -1] + 1;
1039  }
1040 
1041  //Get phytoplankton absorption model coefficients at 1nm
1042  aph_kramer_2022(s->fbands[iw], &s->faphA[iw], &s->faphB[iw]);
1043 
1044  //Get the water absorption coefficient at 1 nm
1045  get_sdp_aw(s->fbands[iw], &s->faw[iw]) ;
1046  }
1047 
1048  //Get bands at which model fits are calculated
1049  for (iw = 0; iw < s->nfitbands; iw++) {
1050  iwxf = windex(fitbands[iw], s->fbands, s->nfbands);
1051  s->fbindx[iw] = iwxf;
1052  }
1053 
1054  //Find array index for required wavelengths
1055  s->i440 = windex(440, s->bands, s->nvbands);
1056  s->i490 = windex(490, s->bands, s->nvbands);
1057  s->i555 = windex(550, s->bands, s->nvbands);
1058  s->i440f = windex(440, s->fbands, s->nfbands);
1059 
1060  //Get the SPD PCA model coefficients
1061  get_sdp_coeff(s->Acoeff, s->Ccoeff, pignames, npig, nAcoeff);
1062 
1063  //set Gordon rrs model coefficients
1064  s->g[0] = 0.0949;
1065  s->g[1] = 0.0794;
1066 
1067  }
1068 
1069 
1070  if ((fmrrs_s = (double *) calloc(s->nfbands, sizeof (double))) == NULL) {
1071  printf("-E- %s line %d : error allocating memory for SDP.\n",
1072  __FILE__, __LINE__);
1073  exit(1);
1074  }
1075  if ((Rrs_a = (float *) calloc(nbands, sizeof (float))) == NULL) {
1076  printf("-E- %s line %d : error allocating memory for SDP.\n",
1077  __FILE__, __LINE__);
1078  exit(1);
1079  }
1080  if ((fRrs_a = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
1081  printf("-E- %s line %d : error allocating memory for SDP.\n",
1082  __FILE__, __LINE__);
1083  exit(1);
1084  }
1085  if ((fsmthRrs_a = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
1086  printf("-E- %s line %d : error allocating memory for SDP.\n",
1087  __FILE__, __LINE__);
1088  exit(1);
1089  }
1090  if ((fmRrs_a = (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
1091  printf("-E- %s line %d : error allocating memory for SDP.\n",
1092  __FILE__, __LINE__);
1093  exit(1);
1094  }
1095  if ((RrsDiff= (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
1096  printf("-E- %s line %d : error allocating memory for SDP.\n",
1097  __FILE__, __LINE__);
1098  exit(1);
1099  }
1100  if ((deriv2= (float *) calloc(s->nfbands, sizeof (float))) == NULL) {
1101  printf("-E- %s line %d : error allocating memory for SDP.\n",
1102  __FILE__, __LINE__);
1103  exit(1);
1104  }
1105 
1106 
1107 
1108  //Initialize outputs for scanline
1109  for (ip=0; ip < npix; ip++) {
1110 
1111  sdpflags[ip] = 0;
1112 
1113  for (ig=0; ig < npig; ig++) {
1114  pigArr[ip*npig + ig] = BAD_FLT;
1115  pigArr_unc[ip*npig + ig] = BAD_FLT;
1116  }
1117 
1118  for (iw=0;iw<nbands;iw++){
1119 
1120  mRrs_out[ip*nbands + iw] = BAD_FLT;
1121  Rrs_diff_out[ip*nbands + iw] = BAD_FLT;
1122  d2Rrs_out[ip*nbands + iw] = BAD_FLT;
1123 
1124  }
1125 
1126  }
1127 
1128  //Process pixels across a scan line
1129  for (ip = 0; ip<l1rec->npix; ip++) {
1130 
1131  ipb = ip*l1rec->l1file->nbands;
1132 
1133  for (iw=0;iw<nbands;iw++){
1134  Rrs_a[iw] = l2rec->Rrs[ipb + iw];
1135  }
1136 
1137  //Interpolate the observed Rrs to 1 nm resolution
1138  sdp_interp(bands, Rrs_a, s->fbands, fRrs_a, nbands, s->nfbands);
1139 
1140  //Smooth the interpolated input Rrs data using 5nm median filter
1141  //sdp_filter(fRrs_a, fsmthRrs_a, s->nfbands, 5);
1142 
1143  sdp_mean_filter(fRrs_a, fsmthRrs_a, s->nfbands, 5);
1144 
1145  //Iniitalize bad rrs flag (helps us flag invalid Rrs_a values)
1146  badrrs=0;
1147 
1148  //Initialize the seawater absorption coefficient at 1nm resolution
1149  for (iw = 0; iw < s->nfitbands; iw++) {
1150 
1151  //Load seawater absorption and backscattering coefficients
1152  //at model fit bands
1153  //s->aw[iw] = l1rec->sw_a[ipb + s->bindx[iw]];
1154 
1155  get_sdp_aw(bands[s->bindx[iw]], &s->aw[iw]) ;
1156 
1157  s->bbw[iw] = l1rec->sw_bb[ipb + s->bindx[iw]];
1158 
1159  //Get the smoothed Rrs values at model fit bands
1160  s->Rrs_a[iw] = fsmthRrs_a[s->fbindx[iw]];
1161 
1162  //Check for bad Rrs
1163  if (s->Rrs_a[iw] == BAD_FLT){
1164  badrrs += 1;
1165  }
1166  }
1167 
1168  //Calculate bbw at fine 1nm resolution
1169  for (iw = 0; iw < s->nfbands; iw++) {
1170  s->fbbw[iw] = seawater_bb(s->fbands[iw], l1rec->sstref[ip],
1171  l1rec->sssref[ip], 0.039);
1172  }
1173 
1174 
1175  //Run the inverse semianalytical model
1176  if(run_sdp_saa(s)) {
1177  l2rec->l1rec->flags[ip] |= PRODFAIL;
1178  };
1179 
1180 
1181  //Compute the spectral shapes adsgtar and bbpstar at 1nm resolution
1182  calc_adgstar(s, 1);
1183  calc_bbpstar(s, 1);
1184 
1185  //Run forward model to generatd Rrs at fine resolution (e.g., 1nm)
1186  sdp_qssaf(s->popt, fmrrs_s, s->nfree, s->nfbands, s);
1187 
1188 
1189  //Compute the difference between observed and modelled above water Rrs
1190  //at 1 nm resolution
1191  for (iw=0;iw<s->nfbands;iw++) {
1192  fmRrs_a[iw] = (0.52 * fmrrs_s[iw]) / (1.0 - 1.7 * fmrrs_s[iw]);
1193  RrsDiff[iw] = fsmthRrs_a[iw] - fmRrs_a[iw];
1194  }
1195 
1196  //compute the 2nd derivative
1197  sdp_deriv2(s->fbands, RrsDiff, deriv2, s->nfbands);
1198 
1199  //Use all 100 x 13 A pigment coefficients. Use all 100 C coefficients.
1200  //Compute the pigment concentrations for each 100 coefficient replicates
1201  for (ig=0;ig<npig;ig++){
1202 
1203  pigtemp_med = 0.0;
1204  pigtemp_std = 0.0;
1205 
1206  for (irep=0;irep<nreps;irep++) {
1207 
1208  //reset temporary variables
1209  pigtemp=0.0;
1210  pigtempArr[irep] = 0.0;
1211 
1212  for (iw=0;iw<nAcoeff;iw++) {
1213  //Compute PCA sum 401 - 699 nm
1214  //Note: PCA model coefficients start at 401nm and derivative Rrs starts at 399nm.
1215  //So we need to start deriv2 multiplication at +6nm.
1216  pigtemp += deriv2[iw+6] * s->Acoeff[ig][irep*nAcoeff + iw];
1217  }
1218 
1219  //Add the intercept for the ith replicate pigment
1220  pigtempArr[irep] = pigtemp + s->Ccoeff[ig][irep];
1221 
1222  }
1223 
1224  //Compute median and standard deviation of the 100 replicate pigment computations
1225  pigtemp_med = gsl_stats_median(pigtempArr,1, nreps);
1226  pigtemp_std = sqrt(gsl_stats_variance(pigtempArr, 1, nreps));
1227 
1228  //check range of pigments and and flag s necessary
1229  set_sdp_flags(pigtemp_med, ig, ip, &sdpflags[ip]);
1230 
1231  pigArr[ip*npig + ig] = pigtemp_med;
1232  pigArr_unc[ip*npig + ig] = pigtemp_std;
1233  }
1234 
1235  //Extract the modelled Rrs, Rrs_diff, and 2nd derivative for output
1236  //at bands closest to OCI
1237  for (iw = 0; iw < nbands; iw++) {
1238 
1239  band_temp = bands[iw];
1240  iwx = windex(band_temp, s->fbands, s->nfbands);
1241 
1242  if (band_temp > 400 && band_temp < 700) {
1243  mRrs_out[ip*nbands + iw] = fmRrs_a[iwx];
1244  Rrs_diff_out[ip*nbands + iw] = RrsDiff[iwx];
1245  d2Rrs_out[ip*nbands + iw] = deriv2[iwx];
1246 
1247  } else {
1248  mRrs_out[ip*nbands + iw] = BAD_FLT;
1249  Rrs_diff_out[ip*nbands + iw] = BAD_FLT;
1250  d2Rrs_out[ip*nbands + iw] = BAD_FLT;
1251  }
1252 
1253  }
1254 
1255  }
1256 
1257  LastRecNum = l1rec->iscan;
1258 
1259  free(Rrs_a);
1260  free(fmrrs_s);
1261  free(fRrs_a);
1262  free(fmRrs_a);
1263  free(RrsDiff);
1264  free(deriv2);
1265  return 0;
1266 }
1267 
1268 
1269 /* ------------------------------------------------------------------- */
1270 /* get_sdp() - returns requested spectral derivative pigment products */
1271 
1272 /* ------------------------------------------------------------------- */
1273 void get_sdp(l2str *l2rec, l2prodstr *p, float prod[]) {
1274 
1275  l1str *l1rec = l2rec->l1rec;
1276  int ib = p->prod_ix;
1277 
1278  int32_t ip, ipg, ipb;
1279  int32_t npix = l2rec->l1rec->npix;
1280 
1281  //Check if calculations have been computed already
1282  if (!sdp_ran(l1rec->iscan)) {
1283  run_sdp(l2rec);
1284  }
1285 
1286  if(p->rank == 3) {
1287 
1288  switch (p->cat_ix) {
1289 
1290  case CAT_mRrs_sdp:
1291  extract_band_3d( mRrs_out, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
1292  break;
1293  case CAT_mRrs_diff_sdp:
1294  extract_band_3d( Rrs_diff_out, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
1295  break;
1296  case CAT_d2Rrs_diff_sdp:
1297  extract_band_3d( d2Rrs_out, prod, l2rec->l1rec->npix, l2rec->l1rec->l1file->nbands);
1298  break;
1299 
1300  }
1301 
1302  } else {
1303 
1304  for (ip = 0; ip < npix; ip++) {
1305 
1306  ipg = ip* npig;
1307  ipb = ip * l1rec->l1file->nbands + ib;
1308 
1309  switch (p->cat_ix) {
1310  case CAT_tchl_sdp:
1311  prod[ip] = (float) pigArr[ipg ];
1312  break;
1313  case CAT_zea_sdp:
1314  prod[ip] = (float) pigArr[ipg + 1];
1315  break;
1316  case CAT_dvchla_sdp:
1317  prod[ip] = (float) pigArr[ipg + 2];
1318  break;
1319  case CAT_butfuco_sdp:
1320  prod[ip] = (float) pigArr[ipg + 3];
1321  break;
1322  case CAT_hexfuco_sdp:
1323  prod[ip] = (float) pigArr[ipg + 4];
1324  break;
1325  case CAT_allo_sdp:
1326  prod[ip] = (float) pigArr[ipg+ 5];
1327  break;
1328  case CAT_mvchlb_sdp:
1329  prod[ip] = (float) pigArr[ipg + 6];
1330  break;
1331  case CAT_neo_sdp:
1332  prod[ip] = (float) pigArr[ipg + 7];
1333  break;
1334  case CAT_viola_sdp:
1335  prod[ip] = (float) pigArr[ipg+ 8];
1336  break;
1337  case CAT_fuco_sdp:
1338  prod[ip] = (float) pigArr[ipg+ 9];
1339  break;
1340  case CAT_chlc12_sdp:
1341  prod[ip] = (float) pigArr[ipg + 10];
1342  break;
1343  case CAT_chlc3_sdp:
1344  prod[ip] = (float) pigArr[ipg + 11];
1345  break;
1346  case CAT_perid_sdp:
1347  prod[ip] = (float) pigArr[ipg + 12];
1348  break;
1349  case CAT_flags_sdp:
1350  prod[ip] = sdpflags[ip];
1351  break;
1352  case CAT_mRrs_sdp:
1353  prod[ip] = (float) mRrs_out[ipb];
1354  break;
1355  case CAT_mRrs_diff_sdp:
1356  prod[ip] = (float) Rrs_diff_out[ipb];
1357  break;
1358  case CAT_d2Rrs_diff_sdp:
1359  prod[ip] = (float) d2Rrs_out[ipb];
1360  break;
1361  default:
1362  printf("Error: %s : Unknown product specifier: %d\n", __FILE__, p->cat_ix);
1363  exit(FATAL_ERROR);
1364  break;
1365  }
1366  }
1367  }
1368 
1369 }
sdp_mean_filter
int sdp_mean_filter(float *y, float *yfilt, int n, int window)
Definition: get_sdp.c:362
CAT_hexfuco_sdp
#define CAT_hexfuco_sdp
Definition: l2prod.h:567
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
table_read_r4
float * table_read_r4(char *input_filename, int m, int n)
Definition: table_io.cpp:2540
nonscience_parameters::success
integer success
Definition: nonscience_parameters.f90:8
j
int j
Definition: decode_rs.h:73
sdp_qssa
void sdp_qssa(double *p, double *mrrs_s, int m, int n, void *data)
Definition: get_sdp.c:466
status
int status
Definition: l1_czcs_hdf.c:32
table_column_count
int table_column_count(char *input_filename)
Definition: table_io.cpp:2524
CAT_neo_sdp
#define CAT_neo_sdp
Definition: l2prod.h:570
CAT_d2Rrs_diff_sdp
#define CAT_d2Rrs_diff_sdp
Definition: l2prod.h:579
OCI
#define OCI
Definition: sensorDefs.h:42
CAT_perid_sdp
#define CAT_perid_sdp
Definition: l2prod.h:575
maxit
int maxit
Definition: pml_iop_tables.c:42
extract_band_3d
float * extract_band_3d(float *in_buf, float *out_buf)
Definition: prodgen.c:67
table
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT table(e.g. m1) to be used for interpolation. The table values will be linearly interpolated using the tables corresponding to the node times bracketing the granule time. If the granule time falls before the time of the first node or after the time of the last node
allocateMemory
void * allocateMemory(size_t numBytes, const char *name)
Definition: allocateMemory.c:7
C
const double C
Definition: calc_par.c:102
NULL
#define NULL
Definition: decode_rs.h:63
calc_bbpstar
void calc_bbpstar(sdpstr *s, int opt)
Definition: get_sdp.c:558
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
band
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second band
Definition: HOWTO_Add_a_product.txt:42
sdp_interp
int sdp_interp(float *x, float *y, float *xinterp, float *yinterp, int m, int n)
Definition: get_sdp.c:327
CAT_mRrs_sdp
#define CAT_mRrs_sdp
Definition: l2prod.h:577
spline
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
sdp_ran
int sdp_ran(int recnum)
Definition: get_sdp.c:67
BADZEA
#define BADZEA
Definition: get_sdp.h:10
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
sdp_filter
int sdp_filter(float *y, float *yfilt, int n, int window)
Definition: get_sdp.c:396
table_row_count
int table_row_count(char *input_filename)
Definition: table_io.cpp:2528
input
instr * input
Definition: main_l2binmatch.cpp:27
BADNEO
#define BADNEO
Definition: get_sdp.h:16
l12_proto.h
CAT_chlc3_sdp
#define CAT_chlc3_sdp
Definition: l2prod.h:574
A
const float A
Definition: calc_par.c:100
BADHEXOFUCO
#define BADHEXOFUCO
Definition: get_sdp.h:13
get_sdp_coeff
int get_sdp_coeff(double **A, double **C, char **pigments, int n, int m)
Definition: get_sdp.c:245
CAT_tchl_sdp
#define CAT_tchl_sdp
Definition: l2prod.h:563
CAT_flags_sdp
#define CAT_flags_sdp
Definition: l2prod.h:576
BADFUCO
#define BADFUCO
Definition: get_sdp.h:18
recnum
read recnum
Definition: HOWTO_Add_a_sensor.txt:72
hico.value_locate.x
list x
Definition: value_locate.py:64
BADMVCHLB
#define BADMVCHLB
Definition: get_sdp.h:15
run_sdp_saa
int run_sdp_saa(sdpstr *s)
Definition: get_sdp.c:633
color_dtdb.y
y
Definition: color_dtdb.py:430
CAT_dvchla_sdp
#define CAT_dvchla_sdp
Definition: l2prod.h:565
BADVIOLA
#define BADVIOLA
Definition: get_sdp.h:17
OCIS
#define OCIS
Definition: sensorDefs.h:43
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
PRISM
#define PRISM
Definition: sensorDefs.h:31
filename
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
BADPERID
#define BADPERID
Definition: get_sdp.h:21
ih
int32_t ih
Definition: atrem_corl1.h:161
check_l2_flags
int check_l2_flags(l2str *l2rec, int ip)
Definition: get_sdp.c:697
run_sdp
int run_sdp(l2str *l2rec)
Definition: get_sdp.c:800
aph_kramer_2022
int aph_kramer_2022(float wave, float *A, float *B)
Definition: get_sdp.c:81
data
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
BADDVCHLA
#define BADDVCHLA
Definition: get_sdp.h:11
get_sdp_aw
int get_sdp_aw(float wave, float *A)
Definition: get_sdp.c:154
sdp_qssaf
void sdp_qssaf(double *p, double *mrrs_s, int m, int n, void *data)
Definition: get_sdp.c:497
CAT_mRrs_diff_sdp
#define CAT_mRrs_diff_sdp
Definition: l2prod.h:578
jac
int32_t jac
Definition: atrem_corl1v2.h:350
CAT_mvchlb_sdp
#define CAT_mvchlb_sdp
Definition: l2prod.h:569
BADBUTOFUCO
#define BADBUTOFUCO
Definition: get_sdp.h:12
B
const float B
Definition: calc_par.c:101
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
BADALLO
#define BADALLO
Definition: get_sdp.h:14
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
u
data_t u
Definition: decode_rs.h:74
sdp_deriv2
int sdp_deriv2(float *x, float *y, float *deriv2, int n)
Definition: get_sdp.c:432
get_sdp.h
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
splint
subroutine splint(xa, ya, y2a, n, x, y)
Definition: atrem_app_refl_f90_cubeio.f:3270
allocate2d_double
double ** allocate2d_double(size_t h, size_t w)
Allocate a two-dimensional array of type double of a given size.
Definition: allocate2d.c:145
strncpy
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second and prod_ix set to PARAM_TYPE_INT name_prefix is compared with the beginning of the product name If name_suffix is not empty the it must match the end of the product name The characters right after the prefix are read as an integer and prod_ix is set to that number strncpy(l2prod->name_prefix, "myprod", UNITLEN)
CAT_butfuco_sdp
#define CAT_butfuco_sdp
Definition: l2prod.h:566
BADTCHL
#define BADTCHL
Definition: get_sdp.h:9
s
data_t s[NROOTS]
Definition: decode_rs.h:75
BASCHLC3
#define BASCHLC3
Definition: get_sdp.h:20
create_images.bands
bands
Definition: create_images.py:91
HICO
#define HICO
Definition: sensorDefs.h:25
CAT_viola_sdp
#define CAT_viola_sdp
Definition: l2prod.h:571
CAT_chlc12_sdp
#define CAT_chlc12_sdp
Definition: l2prod.h:573
get_sdp
void get_sdp(l2str *l2rec, l2prodstr *p, float prod[])
Definition: get_sdp.c:1273
i
int i
Definition: decode_rs.h:71
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
CAT_allo_sdp
#define CAT_allo_sdp
Definition: l2prod.h:568
sdp_jacqssa
void sdp_jacqssa(double *p, double *jac, int n, int m, void *data)
Definition: get_sdp.c:594
CAT_zea_sdp
#define CAT_zea_sdp
Definition: l2prod.h:564
seawater_bb
float seawater_bb(float wave, float sst, float sss, double delta)
Definition: seawater.c:176
calc_adgstar
void calc_adgstar(sdpstr *s, int opt)
Definition: get_sdp.c:525
BADCHLC12
#define BADCHLC12
Definition: get_sdp.h:19
npix
int npix
Definition: get_cmp.c:28
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
CAT_fuco_sdp
#define CAT_fuco_sdp
Definition: l2prod.h:572
set_sdp_flags
int set_sdp_flags(float pigvalue, int pigid, int ip, int16 *flag)
Definition: get_sdp.c:722


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