Due to the lapse in federal government funding, NASA is not updating this website. We sincerely regret this inconvenience.
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
gsm.c
Go to the documentation of this file.
1 /* =================================================================== */
2 /* module gsm.c - Garver-Seigel-Maritorena 2001 Bio-optical model */
3 /* */
4 /* This module contains the functions to optimize and evaluate the */
5 /* Garver-Seigel-Maritorena 2001 Bio-optical model. The optimization */
6 /* is performed for each set of Rrs values within a scanline, using */
7 /* the Amoeba simplex minimization technique (Numerical Recipes). */
8 /* */
9 /* Reference: */
10 /* Maritorena S., D.A. Siegel & A. Peterson, Optimization of a Semi- */
11 /* Analytical Ocean Color Model for Global Scale Applications, Applied */
12 /* Optics, 41(15): 2705-2714, 2002. */
13 /* */
14 /* Implementation: */
15 /* B. Franz, NASA/OCRPG/SAIC, September 2004 (rewrite of previous) */
16 /* */
17 /* =================================================================== */
18 
19 #include <stdio.h>
20 #include <math.h>
21 #include "l12_proto.h"
22 #include "amoeba.h"
23 #include <gsl/gsl_errno.h>
24 #include <gsl/gsl_vector.h>
25 #include <gsl/gsl_blas.h>
26 #include <gsl/gsl_multifit_nlin.h>
27 
28 #define NPAR 3
29 #define MAXITR 500
30 
31 #define GSMDEFAULT 0
32 #define CHESAPEAKE 1
33 
34 #define AMOEBA 0
35 #define LEVMARQ 1
36 
37 static float badval = BAD_FLT;
38 
39 static int32_t nbandsvis;
40 static float grd1 = 0.0949;
41 static float grd2 = 0.0794;
42 static float *aphstar;
43 static float *adgstar;
44 static float *bbpstar;
45 static float *lambda;
46 static float *aw;
47 static float *bbw;
48 
49 static int32_t coefset = GSMDEFAULT;
50 static int32_t pixday;
51 static float32 pixlon;
52 static float32 pixlat;
53 
54 static double *dFdx;
55 static double *dxda;
56 static double *dxdb;
57 
58 static int32_t GSMRecNum = -1;
59 static double *chl;
60 static double *adg;
61 static double *bbp;
62 static int16 *iter;
63 
64 static float merit;
65 
66 struct datastruct {
67  size_t n;
68  double * y;
69  double * sigma;
70 } data;
71 
72 int gsm_ran(int recnum) {
73  if (recnum == GSMRecNum)
74  return 1;
75  else
76  return 0;
77 }
78 
79 
80 /* ------------------------------------------------------------------- */
81 /* get_cb_model - retrieves model coefficients for Chesapeake region */
82 
83 /* ------------------------------------------------------------------- */
84 void set_cb_model(float wave[], int32_t nwave, float lon, float lat, int32_t day) {
85  static int firstCall = 1;
86  static int *iwtab;
87  static int *interp;
88  static int iseason = -1;
89 
90  static int32_t ntwave = 6;
91  static float twave[6] = {412.0, 443.0, 490.0, 510.0, 555.0, 670.0};
92  static char season[4][7] = {"Spring", "Summer", "Fall", "Winter"};
93 
94  static float tcoefs[5][4][6 + 2] = {
95  { /* Upper Bay */
96  {0.02119, 0.02509, 0.01282, 0.00910, 0.00427, 0.02087, 0.01218, 0.0}, /* Spring */
97  {0.02653, 0.02979, 0.01655, 0.01208, 0.00470, 0.02122, 0.01218, 0.0}, /* Summer */
98  {0.02259, 0.02573, 0.01372, 0.01019, 0.00376, 0.02066, 0.01218, 0.0}, /* Fall */
99  {0.02259, 0.02573, 0.01372, 0.01019, 0.00376, 0.02066, 0.01218, 0.0}
100  }, /* Winter */
101  { /* Mid-Bay */
102  {0.02001, 0.02212, 0.01279, 0.00974, 0.00449, 0.01588, 0.01385, 0.0}, /* Spring */
103  {0.03345, 0.03900, 0.02318, 0.01664, 0.00609, 0.02285, 0.01385, 0.0}, /* Summer */
104  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01385, 0.0}, /* Fall */
105  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01385, 0.0}
106  }, /* Winter */
107  { /* Lower Bay */
108  {0.02001, 0.02212, 0.01279, 0.00974, 0.00449, 0.01588, 0.01330, 0.0}, /* Spring */
109  {0.03345, 0.03900, 0.02318, 0.01664, 0.00609, 0.02285, 0.01330, 0.0}, /* Summer */
110  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01330, 0.0}, /* Fall */
111  {0.02758, 0.03080, 0.01826, 0.01438, 0.00691, 0.02112, 0.01330, 0.0}
112  }, /* Winter */
113  { /* Inshore MAB */
114  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}, /* Spring */
115  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}, /* Summer */
116  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}, /* Fall */
117  {0.07123, 0.08843, 0.06024, 0.04072, 0.01693, 0.03815, 0.01236, 0.0}
118  }, /* Winter */
119  { /* Offshore MAB */
120  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}, /* Spring */
121  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}, /* Summer */
122  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}, /* Fall */
123  {0.11331, 0.14678, 0.09832, 0.06048, 0.01920, 0.04349, 0.01646, 1.0}
124  } /* Winter */
125 
126  /* aph*412,aph*443,aph*490,aph*510,aph*555,aph*670, Sdg, Sbbp */
127  };
128 
129  float *taphstar, adg_s, bbp_s;
130  int iw, iregion;
131  if (firstCall == 1) {
132  firstCall = 0;
133  if ((iwtab = (int *) calloc(nbandsvis, sizeof (int))) == NULL) {
134  printf("-E- %s line %d : error allocating memory for iwtab in gsm:set_cb_model.\n",
135  __FILE__, __LINE__);
136  exit(1);
137  }
138  if ((interp = (int *) calloc(nbandsvis, sizeof (int))) == NULL) {
139  printf("-E- %s line %d : error allocating memory for interp in gsm:set_cb_model.\n",
140  __FILE__, __LINE__);
141  exit(1);
142  }
143  }
144 
145  if (iseason < 0) {
146 
147  /* first call, determine the season */
148  if (day >= 79 && day < 172) iseason = 0; /* Spring */
149  else if (day >= 172 && day < 265) iseason = 1; /* Summer */
150  else if (day >= 265 && day < 355) iseason = 2; /* Fall */
151  else iseason = 3; /* Winter */
152 
153  printf("\nInitializing Chesapeake GSM model for %s\n", season[iseason]);
154 
155  /* determine if we need to interpolate */
156  for (iw = 0; iw < nwave; iw++) {
157  iwtab[iw] = windex(wave[iw], twave, ntwave);
158  if (fabsf(wave[iw] - twave[iwtab[iw]]) > 0.5) {
159  printf("\nGSM coefficients will be interpolated for %5.1f nm band.\n", wave[iw]);
160  interp[iw] = 1;
161  } else
162  interp[iw] = 0;
163  }
164  }
165 
166  /* determine region */
167  iregion = -1;
168  if (lat < 37.0 && lon > (-(lat + 152) / 2.5))
169  iregion = 4; /* Offshore MAB */
170  else if (lat >= 37.0 && lon > ((lat - 153.31) / 1.5385))
171  iregion = 4; /* Offshore MAB */
172  else if (lon > -75.7)
173  iregion = 3; /* Inshore MAB */
174  else if (lat < 37.6)
175  iregion = 2; /* Lower Bay */
176  else if (lat < 38.6)
177  iregion = 1; /* Mid Bay */
178  else
179  iregion = 0; /* Upper Bay */
180 
181  /* compute coefficients */
182  taphstar = &tcoefs[iregion][iseason][0];
183  adg_s = tcoefs[iregion][iseason][ntwave + 0];
184  bbp_s = tcoefs[iregion][iseason][ntwave + 1];
185  for (iw = 0; iw < nwave; iw++) {
186  if (interp[iw])
187  aphstar[iw] = linterp(twave, taphstar, ntwave, wave[iw]);
188  else
189  aphstar[iw] = taphstar[iwtab[iw]];
190  adgstar[iw] = exp(-adg_s * (wave[iw] - 443.0));
191  bbpstar[iw] = pow((443.0 / wave[iw]), bbp_s);
192  }
193 
194  return;
195 }
196 
197 
198 /* ------------------------------------------------------------------- */
199 /* gsm_amb() GSM model, set-up to work with amoeba optimization */
200 /* */
201 /* Input: */
202 /* auxdata - amoeba interface structure */
203 /* par[] - model parameters */
204 /* */
205 /* Output: */
206 /* function returns chi-squared, which is the quantity minimized. */
207 /* */
208 
209 /* ------------------------------------------------------------------- */
210 double gsm_amb(FITSTRUCT *auxdata, double par[]) {
211 
212  int iw;
213  float x, F, bb, ac;
214  double merit = 0.0;
215  int nwave = auxdata->npnts;
216 
217  if (coefset == CHESAPEAKE) {
218  set_cb_model(lambda, nwave, pixlon, pixlat, pixday);
219  }
220 
221  for (iw = 0; iw < nwave; iw++) {
222 
223  ac = aw [iw] + aphstar[iw] * par[0] + par[1] * adgstar[iw];
224  bb = bbw[iw] + par[2] * bbpstar[iw];
225  x = bb / (ac + bb);
226  F = grd1 * x + grd2 * pow(x, 2);
227 
228  auxdata->yfit[iw] = F;
229 
230  merit += pow((auxdata->y[iw] - F), 2) * auxdata->wgt[iw];
231  }
232 
233  auxdata->merit = merit;
234 
235  return (merit);
236 }
237 
238 
239 /* ------------------------------------------------------------------- */
240 /* fit_gsm_amb() - run amoeba optimization of GSM01 for one spectrum */
241 
242 /* ------------------------------------------------------------------- */
243 int fit_gsm_amb(double Rrs[], double wts[], int32_t npts, double fitparms[],
244  double Rrs_fit[], int16 *itercnt) {
245  static float tol = 1.e-6; /* fractional change in chisqr */
246  static FITSTRUCT auxdata; /* amoeba interface structure */
247  static double init[NPAR * (NPAR + 1)]; /* initial simplex */
248  static double p0[] = {.1, .02, .0001}; /* starting parameters (C,adg,bb)*/
249  static double d0[] = {5., 1.0, 0.01}; /* characteristic length */
250 
251  int i, j;
252  short isml;
253  int status = 1;
254 
255  auxdata.niter = MAXITR; /* max number of iterations allowed */
256  auxdata.nfunc = NPAR; /* number of model parameters */
257  auxdata.npnts = npts; /* number of wavelengths (Rrs values) */
258  auxdata.y = Rrs; /* Input Rrs values */
259  auxdata.wgt = wts; /* Input weights on Rrs values */
260  auxdata.yfit = Rrs_fit; /* Output model predicted Rrs values */
261 
262  /* initialize simplex with first guess model parameters */
263  for (j = 0; j < NPAR + 1; j++) for (i = 0; i < NPAR; i++) init[j * NPAR + i] = p0[i];
264  for (i = 0; i < NPAR; i++) {
265  init[NPAR + i * (NPAR + 1)] += d0[i];
266  fitparms[i] = 0.0;
267  }
268 
269  /* run optimization */
270  isml = amoeba(init, &auxdata, gsm_amb, tol);
271 
272  /* check convergence and record parameter results */
273  if (auxdata.niter > MAXITR)
274  /* failed to converge */
275  status = 0;
276 
277  for (i = 0; i < NPAR; i++)
278  fitparms[i] = init[NPAR * isml + i];
279 
280  *itercnt = auxdata.niter;
281 
282  return (status);
283 }
284 
285 
286 /* ------------------------------------------------------------------- */
287 /* run_gsm_amb() - returns requested GSM product for full scanline */
288 
289 /* ------------------------------------------------------------------- */
290 void run_gsm_amb(l2str *l2rec) {
291 
292  int32_t ip, iw;
293  double model [NPAR];
294  double *Rrs;
295  double *wts;
296  double *Rrs_fit;
297  int16 itercnt;
298  int16 bndcnt;
299  int status;
300 
301  l1str *l1rec = l2rec->l1rec;
302 
303  if ((Rrs = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
304  printf("-E- %s line %d : error allocating memory for Rrs in gsm:run_gsm_amb.\n",
305  __FILE__, __LINE__);
306  exit(1);
307  }
308  if ((wts = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
309  printf("-E- %s line %d : error allocating memory for wts in gsm:run_gsm_amb.\n",
310  __FILE__, __LINE__);
311  exit(1);
312  }
313  if ((Rrs_fit = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
314  printf("-E- %s line %d : error allocating memory for Rrs_fit in gsm:run_gsm_amb.\n",
315  __FILE__, __LINE__);
316  exit(1);
317  }
318  /* if this is a new scanline, fit model for every pixel of the scanline */
319 
320  for (ip = 0; ip < l1rec->npix; ip++) {
321 
322  status = 1;
323  bndcnt = 0;
324  chl [ip] = badval;
325  adg [ip] = badval;
326  bbp [ip] = badval;
327  iter[ip] = 0;
328 
329  if (l1rec->mask[ip] == 0) {
330 
331  /* for CB model selection */
332  pixlon = l1rec->lon[ip];
333  pixlat = l1rec->lat[ip];
334  int16_t year, day;
335  double sec;
336  unix2yds(l1rec->scantime, &year, &day, &sec);
337  pixday = (int32_t) day;
338 
339  for (iw = 0; iw < nbandsvis; iw++) {
340  wts[iw] = 1.0;
341  Rrs[iw] = l2rec->Rrs[ip * l1rec->l1file->nbands + iw];
342  Rrs[iw] = Rrs[iw] / (0.52 + 1.7 * Rrs[iw]);
343  if (Rrs[iw] <= 0.0) status = 0;
344  else bndcnt++;
345  }
346 
347  /* if less than 3 valid radiances, fail pixel */
348  if (bndcnt < 3) {
349  l1rec->flags[ip] |= PRODFAIL;
350  continue;
351  }
352 
353  /* if any negative reflectance, set warning flag */
354  /* but allow processing to continue */
355  if (status == 0) {
356  l1rec->flags[ip] |= PRODWARN;
357  }
358 
359  /* fit model for this pixel */
360  status = fit_gsm_amb(Rrs, wts, nbandsvis, model, Rrs_fit, &itercnt);
361 
362  /* save results and flag failure */
363  if (status != 1 || model[0] <= 0.0) {
364  l1rec->flags[ip] |= PRODFAIL;
365  } else {
366  chl[ip] = model[0];
367  adg[ip] = model[1];
368  bbp[ip] = model[2];
369  }
370  iter[ip] = itercnt;
371  } else {
372  l1rec->flags[ip] |= PRODFAIL;
373  }
374  }
375 
376  GSMRecNum = l1rec->iscan;
377 
378  free(Rrs);
379  free(wts);
380  free(Rrs_fit);
381  return;
382 }
383 
384 
385 /* ------------------------------------------------------------------- */
386 /* gsm_lm_f() GSM01 model, set-up to work with L.M. optimization */
387 /* */
388 /* Input: */
389 /* par[] - model parameters */
390 /* data - data to fit */
391 /* f - gsl_multifit_function */
392 /* */
393 
394 /* ------------------------------------------------------------------- */
395 int gsm_lm_f(const gsl_vector *par, void *data, gsl_vector *f) {
396  static int first = 1;
397 
398  int iw;
399  float x, F, bb, ac;
400  size_t nwave = ((struct datastruct *) data)->n;
401  double *y = ((struct datastruct *) data)->y;
402  double *sigma = ((struct datastruct *) data)->sigma;
403 
404  if (first == 1) {
405 
406  if ((dFdx = (double *) calloc(nwave, sizeof (double))) == NULL) {
407  printf("-E- %s line %d : error allocating memory for dFdx in gsm:gsm_lm_f.\n",
408  __FILE__, __LINE__);
409  exit(1);
410  }
411  if ((dxda = (double *) calloc(nwave, sizeof (double))) == NULL) {
412  printf("-E- %s line %d : error allocating memory for dxda in gsm:gsm_lm_f.\n",
413  __FILE__, __LINE__);
414  exit(1);
415  }
416  if ((dxdb = (double *) calloc(nwave, sizeof (double))) == NULL) {
417  printf("-E- %s line %d : error allocating memory for dxdb in gsm:gsm_lm_f.\n",
418  __FILE__, __LINE__);
419  exit(1);
420  }
421  first = 0;
422  }
423  if (coefset == CHESAPEAKE) {
424  set_cb_model(lambda, nwave, pixlon, pixlat, pixday);
425  }
426 
427  merit = 0;
428  for (iw = 0; iw < nwave; iw++) {
429  ac = aw [iw] +
430  aphstar[iw] * gsl_vector_get(par, 0) + gsl_vector_get(par, 1) * adgstar[iw];
431  bb = bbw[iw] + gsl_vector_get(par, 2) * bbpstar[iw];
432  x = bb / (ac + bb);
433  F = grd1 * x + grd2 * pow(x, 2);
434 
435  dFdx[iw] = grd1 + 2 * grd2*x;
436  dxda[iw] = -x * x / bb;
437  dxdb[iw] = x / bb + dxda[iw];
438 
439  gsl_vector_set(f, iw, (F - y[iw]) / sigma[iw]);
440 
441  merit += pow((F - y[iw]) / sigma[iw], 2);
442  }
443 
444  return GSL_SUCCESS;
445 }
446 
447 
448 /* ------------------------------------------------------------------- */
449 /* gsm_lm_df() GSM01 model, set-up to work with L.M. optimization */
450 /* */
451 /* Input: */
452 /* par[] - model parameters */
453 /* data - data to fit */
454 /* J - Jacobian */
455 /* */
456 
457 /* ------------------------------------------------------------------- */
458 int gsm_lm_df(const gsl_vector *par, void *data, gsl_matrix *J) {
459  int iw;
460  size_t nwave = ((struct datastruct *) data)->n;
461  double *sigma = ((struct datastruct *) data)->sigma;
462 
463  if (coefset == CHESAPEAKE) {
464  set_cb_model(lambda, nwave, pixlon, pixlat, pixday);
465  }
466 
467  for (iw = 0; iw < nwave; iw++) {
468 
469  /* Jacobian matrix J(iw,j) = dfi / dxj, */
470  gsl_matrix_set(J, iw, 0, dFdx[iw] * dxda[iw] * aphstar[iw] / sigma[iw]);
471  gsl_matrix_set(J, iw, 1, dFdx[iw] * dxda[iw] * adgstar[iw] / sigma[iw]);
472  gsl_matrix_set(J, iw, 2, dFdx[iw] * dxdb[iw] * bbpstar[iw] / sigma[iw]);
473  }
474 
475  return GSL_SUCCESS;
476 }
477 
478 int gsm_lm_fdf(const gsl_vector *par, void *data,
479  gsl_vector *f, gsl_matrix *J) {
480  gsm_lm_f(par, data, f);
481  gsm_lm_df(par, data, J);
482 
483  return GSL_SUCCESS;
484 }
485 
486 
487 /* ---------------------------------------------------------------------- */
488 /* run_gsm_lm() - returns requested GSM product for full scanline */
489 
490 /* ---------------------------------------------------------------------- */
491 void run_gsm_lm(l2str *l2rec) {
492  static int firstCall = TRUE;
493 
494  int32_t ip, iw;
495  double *Rrs;
496  double *wts;
497  double *sigma;
498  int16 itercnt;
499  int16 bndcnt;
500  int status;
501  l1str *l1rec = l2rec->l1rec;
502 
503  const gsl_multifit_fdfsolver_type *T;
504  gsl_multifit_fdfsolver *s;
505 
506  static gsl_multifit_function_fdf f;
507  double x_init[3] = {0.0, 0.0, 0.0};
508  gsl_vector_view x = gsl_vector_view_array(x_init, NPAR);
509  if ((Rrs = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
510  printf("-E- %s line %d : error allocating memory for Rrs in gsm:run_gsm_amb.\n",
511  __FILE__, __LINE__);
512  exit(1);
513  }
514  if ((wts = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
515  printf("-E- %s line %d : error allocating memory for wts in gsm:run_gsm_amb.\n",
516  __FILE__, __LINE__);
517  exit(1);
518  }
519  if ((sigma = (double *) calloc(l1rec->l1file->nbands, sizeof (double))) == NULL) {
520  printf("-E- %s line %d : error allocating memory for sigma in gsm:run_gsm_amb.\n",
521  __FILE__, __LINE__);
522  exit(1);
523  }
524 
525  if (firstCall == TRUE) {
526 
527  firstCall = FALSE;
528 
529  /* Set up data structure */
530  data.n = nbandsvis;
531  data.y = Rrs;
532  data.sigma = sigma;
533 
534  /* Set up multifit function structure */
535  f.f = &gsm_lm_f;
536  f.df = &gsm_lm_df;
537  f.fdf = &gsm_lm_fdf;
538  f.n = nbandsvis;
539  f.p = NPAR;
540  f.params = &data;
541  }
542 
543  T = gsl_multifit_fdfsolver_lmsder;
544  s = gsl_multifit_fdfsolver_alloc(T, nbandsvis, NPAR);
545 
546  /* if this is a new scanline, fit model for every pixel of the scanline */
547 
548  for (ip = 0; ip < l1rec->npix; ip++) {
549 
550  status = 1;
551  bndcnt = 0;
552  chl [ip] = badval;
553  adg [ip] = badval;
554  bbp [ip] = badval;
555  iter[ip] = 0;
556 
557  if (l1rec->mask[ip] == 0) {
558 
559  /* for CB model selection */
560  pixlon = l1rec->lon[ip];
561  pixlat = l1rec->lat[ip];
562  int16_t year, day;
563  double sec;
564  unix2yds(l1rec->scantime, &year, &day, &sec);
565  pixday = (int32_t) day;
566 
567 
568  for (iw = 0; iw < nbandsvis; iw++) {
569  wts[iw] = 1.0;
570  sigma[iw] = 1 / sqrt(wts[iw]);
571  Rrs[iw] = l2rec->Rrs[ip * l1rec->l1file->nbands + iw];
572  Rrs[iw] = Rrs[iw] / (0.52 + 1.7 * Rrs[iw]);
573  if (Rrs[iw] <= 0.0) status = 0;
574  else bndcnt++;
575  }
576 
577  /* if less than 3 valid radiances, fail pixel */
578  if (bndcnt < 3) {
579  l1rec->flags[ip] |= PRODFAIL;
580  continue;
581  }
582 
583  /* if any negative reflectance, set warning flag */
584  /* but allow processing to continue */
585  if (status == 0) {
586  l1rec->flags[ip] |= PRODWARN;
587  }
588 
589  /* fit model for this pixel */
590  gsl_multifit_fdfsolver_set(s, &f, &x.vector);
591  itercnt = 0;
592  do {
593  itercnt++;
594  status = gsl_multifit_fdfsolver_iterate(s);
595 
596  //printf ("status = %s\n", gsl_strerror (status));
597 
598  if (status)
599  break;
600 
601  status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
602 
603  } while (status == GSL_CONTINUE && itercnt < MAXITR);
604 
605 
606  /* save results and flag failure */
607  if (itercnt == MAXITR || gsl_vector_get(s->x, 0) <= 0.0) {
608  l1rec->flags[ip] |= PRODFAIL;
609  } else {
610  chl[ip] = gsl_vector_get(s->x, 0);
611  adg[ip] = gsl_vector_get(s->x, 1);
612  bbp[ip] = gsl_vector_get(s->x, 2);
613  }
614  iter[ip] = itercnt;
615  } else {
616  l1rec->flags[ip] |= PRODFAIL;
617  }
618  }
619 
620  GSMRecNum = l1rec->iscan;
621  gsl_multifit_fdfsolver_free(s);
622 
623  free(Rrs);
624  free(wts);
625  free(sigma);
626  return;
627 }
628 
629 
630 /* ---------------------------------------------------------------------- */
631 /* run_gsm - runs optimization using requested method */
632 
633 /* ---------------------------------------------------------------------- */
634 void run_gsm(l2str *l2rec) {
635  static int firstCall = 1;
636  static int last_npix;
637 
638  l1str *l1rec = l2rec->l1rec;
639  filehandle *l1file = l1rec->l1file;
640  int32_t nbands = l1file->nbands;
641 
642  /* return requested product, return -1 for abs/scatter in NIR */
643 
644  if (firstCall) {
645  firstCall = 0;
646 
647  int32_t iw;
648  int32_t taphn = -1;
649  float *taphw = NULL;
650  float *taphs = NULL;
651  float adg_s = -1;
652  float bbp_s = -1;
653 
654  int16_t year, day;
655  double sec;
656  unix2yds(l2rec->l1rec->scantime, &year, &day, &sec);
657  pixday = day;
658 
659  if ((aphstar = (float *) calloc(nbands, sizeof (float))) == NULL) {
660  printf("-E- %s line %d : error allocating memory for gsm.\n",
661  __FILE__, __LINE__);
662  exit(1);
663  }
664  if ((adgstar = (float *) calloc(nbands, sizeof (double))) == NULL) {
665  printf("-E- %s line %d : error allocating memory for gsm.\n",
666  __FILE__, __LINE__);
667  exit(1);
668  }
669  if ((bbpstar = (float *) calloc(nbands, sizeof (double))) == NULL) {
670  printf("-E- %s line %d : error allocating memory for gsm.\n",
671  __FILE__, __LINE__);
672  exit(1);
673  }
674  if ((lambda = (float *) calloc(nbands, sizeof (float))) == NULL) {
675  printf("-E- %s line %d : error allocating memory for gsm.\n",
676  __FILE__, __LINE__);
677  exit(1);
678  }
679  if ((aw = (float *) calloc(nbands, sizeof (float))) == NULL) {
680  printf("-E- %s line %d : error allocating memory for gsm.\n",
681  __FILE__, __LINE__);
682  exit(1);
683  }
684  if ((bbw = (float *) calloc(nbands, sizeof (float))) == NULL) {
685  printf("-E- %s line %d : error allocating memory for gsm.\n",
686  __FILE__, __LINE__);
687  exit(1);
688  }
689  if ((chl = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
690  printf("-E- %s line %d : error allocating memory for gsm.\n",
691  __FILE__, __LINE__);
692  exit(1);
693  }
694  if ((adg = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
695  printf("-E- %s line %d : error allocating memory for gsm.\n",
696  __FILE__, __LINE__);
697  exit(1);
698  }
699  if ((bbp = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
700  printf("-E- %s line %d : error allocating memory for gsm.\n",
701  __FILE__, __LINE__);
702  exit(1);
703  }
704  if ((iter = (int16 *) calloc(l1rec->npix, sizeof (int16))) == NULL) {
705  printf("-E- %s line %d : error allocating memory for gsm.\n",
706  __FILE__, __LINE__);
707  exit(1);
708  }
709  last_npix = l1rec->npix;
710 
711  /* set model coefficient set */
712  coefset = input->gsm_opt;
713 
714  /* build array of sensor wavelengths */
715  for (iw = 0; iw < l1rec->l1file->nbands; iw++) {
716  lambda[iw] = l1rec->l1file->fwave[iw];
717  }
718 
719  /* limit to visible wavelengths */
720  nbandsvis = rdsensorinfo(l1rec->l1file->sensorID, l1_input->evalmask, "NbandsVIS", NULL);
721 
722  /* set static coefficients */
723  if (coefset == GSMDEFAULT) {
724  adg_s = input->gsm_adg_s;
725  bbp_s = input->gsm_bbp_s;
726  taphw = input->gsm_aphw; /* aphstar table wavelengths */
727  taphs = input->gsm_aphs; /* aphstar table coefficients */
728  taphn = 0; /* aphstar table size */
729  for (taphn = 0; taphn < nbandsvis; taphn++)
730  if (taphw[taphn] < 0) break;
731  for (iw = 0; iw < nbandsvis; iw++) {
732  aphstar[iw] = linterp(taphw, taphs, taphn, lambda[iw]);
733  adgstar[iw] = exp(-adg_s * (lambda[iw] - 443.0));
734  bbpstar[iw] = pow((443.0 / lambda[iw]), bbp_s);
735  }
736  }
737 
738  /* save associated irradiances and water absorption/backscatter */
739  /* use band-averaged vales if no nLw out-of-band correction applied */
740 
741  if (l1_input->outband_opt >= 2) {
742  for (iw = 0; iw < nbandsvis; iw++) {
743  aw [iw] = aw_spectra(lambda[iw], BANDW);
744  bbw[iw] = bbw_spectra(lambda[iw], BANDW);
745  }
746  } else {
747  for (iw = 0; iw < nbandsvis; iw++) {
748  aw [iw] = l1rec->l1file->aw[iw];
749  bbw[iw] = l1rec->l1file->bbw[iw];
750  }
751  }
752  }
753 
754  // see if the number of pixels got bigger
755  if (l1rec->npix > last_npix) {
756  free(chl);
757  if ((chl = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
758  printf("-E- %s line %d : error allocating memory for gsm.\n",
759  __FILE__, __LINE__);
760  exit(1);
761  }
762  free(adg);
763  if ((adg = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
764  printf("-E- %s line %d : error allocating memory for gsm.\n",
765  __FILE__, __LINE__);
766  exit(1);
767  }
768  free(bbp);
769  if ((bbp = (double *) calloc(l1rec->npix, sizeof (double))) == NULL) {
770  printf("-E- %s line %d : error allocating memory for gsm.\n",
771  __FILE__, __LINE__);
772  exit(1);
773  }
774  free(iter);
775  if ((iter = (int16 *) calloc(l1rec->npix, sizeof (int16))) == NULL) {
776  printf("-E- %s line %d : error allocating memory for gsm.\n",
777  __FILE__, __LINE__);
778  exit(1);
779  }
780  last_npix = l1rec->npix;
781  }
782 
783  switch (input->gsm_fit) {
784  case AMOEBA:
785  run_gsm_amb(l2rec);
786  break;
787  case LEVMARQ:
788  run_gsm_lm(l2rec);
789  break;
790  default:
791  printf("%s Line %d: Unknown optimization method for GSM %d\n",
792  __FILE__, __LINE__, input->gsm_fit);
793  exit(1);
794  break;
795  }
796 }
797 
798 
799 /* ------------------------------------------------------------------- */
800 /* get_gsm() - returns requested GSM product for full scanline */
801 
802 /* ------------------------------------------------------------------- */
803 void get_gsm(l2str *l2rec, l2prodstr *p, float prod[]) {
804  int prodID = p->cat_ix;
805  int band = p->prod_ix;
806  int32_t ip, iw;
807 
808  l1str *l1rec = l2rec->l1rec;
809  filehandle *l1file = l1rec->l1file;
810 
811  if (band >= 0) {
812  if (l1file->fwave[band] > 700.0) {
813  for (ip = 0; ip < l1rec->npix; ip++)
814  prod[ip] = -1.0;
815  return;
816  }
817  }
818 
819  iw = p->prod_ix;
820 
821  if (!gsm_ran(l1rec->iscan))
822  run_gsm(l2rec);
823 
824  for (ip = 0; ip < l1rec->npix; ip++) {
825 
826  switch (prodID) {
827 
828  case CAT_chl_gsm:
829  prod[ip] = (float) chl[ip];
830  break;
831 
832  case CAT_adg_gsm:
833  if (band < 0)
834  prod[ip] = (float) adg[ip];
835  else
836  prod[ip] = (float) adg[ip] * adgstar[iw];
837  break;
838 
839  case CAT_bbp_gsm:
840  if (band < 0)
841  prod[ip] = (float) bbp[ip];
842  else
843  prod[ip] = (float) bbp[ip] * bbpstar[iw];
844  break;
845 
846  case CAT_aph_gsm:
847  prod[ip] = (float) aphstar[iw] * chl[ip];
848  break;
849 
850  case CAT_a_gsm:
851  prod[ip] = (float) aw[iw] + aphstar[iw] * chl[ip] + adg[ip] * adgstar[iw];
852  break;
853 
854  case CAT_bb_gsm:
855  prod[ip] = (float) bbw[iw] + bbp[ip] * bbpstar[iw];
856  break;
857 
858  default:
859  printf("-E- %s line %d : erroneous product ID %d passed to gsm.\n",
860  __FILE__, __LINE__, prodID);
861  exit(1);
862  }
863  }
864 
865  return;
866 }
867 
868 
869 /* ------------------------------------------------------------------- */
870 /* get_gsm_iter() - returns iteration count */
871 
872 /* ------------------------------------------------------------------- */
873 int16 *get_iter_gsm(l2str *l2rec) {
874  coefset = input->gsm_opt;
875  if (!gsm_ran(l2rec->l1rec->iscan))
876  run_gsm(l2rec);
877 
878  return iter;
879 }
880 
881 /* Interface to convl12() to return GSM iops */
882 
883 void iops_gsm(l2str *l2rec) {
884  int32_t ib, ip, ipb;
885  int32_t nbands = l2rec->l1rec->l1file->nbands;
886 
887  coefset = input->gsm_opt;
888  if (!gsm_ran(l2rec->l1rec->iscan))
889  run_gsm(l2rec);
890 
891  for (ip = 0; ip < l2rec->l1rec->npix; ip++)
892  for (ib = 0; ib < nbands; ib++) {
893  ipb = ip * nbands + ib;
894  l2rec->a [ipb] = (float) aw[ib] + aphstar[ib] * chl[ip] + adg[ip] * adgstar[ib];
895  l2rec->bb[ipb] = (float) bbw[ib] + bbp[ip] * bbpstar[ib];
896  }
897 
898  return;
899 }
900 
CAT_bb_gsm
#define CAT_bb_gsm
Definition: l2prod.h:94
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
j
int j
Definition: decode_rs.h:73
CAT_adg_gsm
#define CAT_adg_gsm
Definition: l2prod.h:64
datastruct::y
double * y
Definition: gsm.c:68
day
int32_t day
Definition: atrem_corl1v2.h:308
status
int status
Definition: l1_czcs_hdf.c:32
MAXITR
#define MAXITR
Definition: gsm.c:29
CAT_aph_gsm
#define CAT_aph_gsm
Definition: l2prod.h:95
FITSTRUCT::niter
int niter
Definition: amoeba.h:9
FALSE
#define FALSE
Definition: rice.h:164
NULL
#define NULL
Definition: decode_rs.h:63
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
PRODWARN
#define PRODWARN
Definition: l2_flags.h:13
flags::PRODWARN
@ PRODWARN
TRUE
#define TRUE
Definition: rice.h:165
run_gsm_lm
void run_gsm_lm(l2str *l2rec)
Definition: gsm.c:491
get_gsm
void get_gsm(l2str *l2rec, l2prodstr *p, float prod[])
Definition: gsm.c:803
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
lambda
data_t lambda[NROOTS+1]
Definition: decode_rs.h:75
amoeba
short amoeba(double *pnts, FITSTRUCT *auxdata, double(*func)(FITSTRUCT *, double[]), float tol)
Definition: amoeba.c:14
FITSTRUCT::merit
double merit
Definition: amoeba.h:4
input
instr * input
Definition: main_l2binmatch.cpp:27
F
const double F
l12_proto.h
data
struct datastruct data
get_iter_gsm
int16 * get_iter_gsm(l2str *l2rec)
Definition: gsm.c:873
run_gsm_amb
void run_gsm_amb(l2str *l2rec)
Definition: gsm.c:290
set_cb_model
void set_cb_model(float wave[], int32_t nwave, float lon, float lat, int32_t day)
Definition: gsm.c:84
init
int init(int32_t ipr, int32_t jpr, char *efile, char *pfile)
Definition: proj_report.c:51
J
int J
Definition: Usds.c:62
run_local.lat
lat
Definition: run_local.py:69
fit_gsm_amb
int fit_gsm_amb(double Rrs[], double wts[], int32_t npts, double fitparms[], double Rrs_fit[], int16 *itercnt)
Definition: gsm.c:243
f
double precision function f(R1)
Definition: tmd.lp.f:1454
interp::interp
Interp< T, U, K > interp(T *lut, std::vector< U * > &grid, Args &&...dims)
Definition: Interpolation.hpp:278
NPAR
#define NPAR
Definition: gsm.c:28
recnum
read recnum
Definition: HOWTO_Add_a_sensor.txt:72
run_local.lon
lon
Definition: run_local.py:69
linterp
float linterp(float xin[], float yin[], int32_t nin, float xout)
Definition: lspline.c:59
hico.value_locate.x
list x
Definition: value_locate.py:64
FITSTRUCT::wgt
double * wgt
Definition: amoeba.h:6
color_dtdb.y
y
Definition: color_dtdb.py:430
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
bbw_spectra
float bbw_spectra(int32_t wl, int32_t width)
Definition: water_spectra.c:142
tol
float tol
Definition: pml_iop_tables.c:45
create_images.Rrs
def Rrs
Definition: create_images.py:92
GSMDEFAULT
#define GSMDEFAULT
Definition: gsm.c:31
unix2yds
void unix2yds(double usec, short *year, short *day, double *secs)
ac
float * ac[MAX_BANDS]
Definition: pml_iop_tables.c:26
FITSTRUCT::y
double * y
Definition: amoeba.h:5
run_gsm
void run_gsm(l2str *l2rec)
Definition: gsm.c:634
gsm_ran
int gsm_ran(int recnum)
Definition: gsm.c:72
datastruct::sigma
double * sigma
Definition: gsm.c:69
AMOEBA
#define AMOEBA
Definition: gsm.c:34
FITSTRUCT::nfunc
short int nfunc
Definition: amoeba.h:2
iops_gsm
void iops_gsm(l2str *l2rec)
Definition: gsm.c:883
FITSTRUCT::yfit
double * yfit
Definition: amoeba.h:8
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
CAT_a_gsm
#define CAT_a_gsm
Definition: l2prod.h:93
datastruct::n
size_t n
Definition: gsm.c:67
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
datastruct
Definition: gsm.c:66
windex
int windex(float wave, float twave[], int ntwave)
Definition: windex.c:73
CAT_bbp_gsm
#define CAT_bbp_gsm
Definition: l2prod.h:65
FITSTRUCT::npnts
short int npnts
Definition: amoeba.h:3
s
data_t s[NROOTS]
Definition: decode_rs.h:75
gsm_lm_f
int gsm_lm_f(const gsl_vector *par, void *data, gsl_vector *f)
Definition: gsm.c:395
rdsensorinfo
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
model
int32_t model
Definition: atrem_corl1.h:132
BANDW
#define BANDW
Definition: l1.h:53
gsm_lm_fdf
int gsm_lm_fdf(const gsl_vector *par, void *data, gsl_vector *f, gsl_matrix *J)
Definition: gsm.c:478
gsm_amb
double gsm_amb(FITSTRUCT *auxdata, double par[])
Definition: gsm.c:210
flags::PRODFAIL
@ PRODFAIL
CAT_chl_gsm
#define CAT_chl_gsm
Definition: l2prod.h:63
i
int i
Definition: decode_rs.h:71
gsm_lm_df
int gsm_lm_df(const gsl_vector *par, void *data, gsl_matrix *J)
Definition: gsm.c:458
CHESAPEAKE
#define CHESAPEAKE
Definition: gsm.c:32
LEVMARQ
#define LEVMARQ
Definition: gsm.c:35
p
float p[MODELMAX]
Definition: atrem_corl1.h:131


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