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
swim.c
Go to the documentation of this file.
1 /* =================================================================== */
2 /* module swim.c - Shallow Water Inversion Model */
3 /* */
4 /* This module contains the functions to optimize and evaluate the */
5 /* SWIM Bio-optical model. The SWIM algorithm was designed to improve */
6 /* retrievals in optically shallow waters by considering both water */
7 /* column depth and benthic albedo */
8 /* */
9 /* References: */
10 /* McKinna et al (2015) A semianalytical ocean color */
11 /* inversion algorithm with explicit water-column depth and substrate */
12 /* reflectance parameterization, JGR Oceans, 120, */
13 /* doi: 10.1002/2014JC010224 */
14 /* */
15 /* Reichstetter et al (2014) Seafloor brightness map of the Great */
16 /* Barrier Reef, Australia, derived from biodiversity data, PANGEA, */
17 /* doi:10.1594/PANGAEA.835979 */
18 /* */
19 /* Implementation: */
20 /* L. McKinna and D. Shea, NASA/OBPG/SAIC, March 2015 */
21 /* */
22 /* Notes: */
23 /* In lieu of user-supplied data: */
24 /* - Uses ETOP01 as default bathymetry. */
25 /* - Uses single albedo spectrum as default. */
26 /* =================================================================== */
27 
28 #include <stdio.h>
29 #include <math.h>
30 #include <levmar.h>
31 #include <netcdf.h>
32 
33 #include "l12_proto.h"
34 
35 #define NPAR 3 //Number of free model parameters
36 #define MAXITR 1500 //Maximum number of LM optimisation iterations
37 static int32_t swimScanNum = -1; //last scan num that the model was calculated for
38 static int32_t LastRecNum = -1;
39 
40 static int32_t* lambda; //array holding wavelengths for this sensor
41 static int32_t nbandVis; //number of visible bands for a given sensor
42 
43 //output product pointers
44 static float *adg; // adg value for the scan line adg[pix, band]
45 static float *bbp; // bbp value for the scan line bbp[pix, band]
46 static float *aph; // aph value for the scan line aph[pix, band]
47 static float *atot; // atot value for the scan line atot[pix, band]
48 static float *bbtot; // bbtot value for the scan line bbtot[pix, band]
49 static int16 *iter;
50 
51 //Bio-optical model parameters
52 static float *aphstar; //spectral shape of phyto abs (aphstar[band])
53 static float *adgstar; //spectral shape of detritus + cdom abs (adgstar[band])
54 static float *bbpstar; //spectral shape of particle backscat (bbpstar[band])
55 static float *aw; //pure water abs spectra (aw[band])
56 static float *bbw; //pure water backscat spectra (bbw[band])
57 
58 //static float grd1 = 0.0949; //the "g0" coefficient in rrs function Lee et al 1998
59 //static float grd2 = 0.0794; //the "g1" coefficient in rrs function Lee et al 1998
60 
61 static const double grd1 = 0.08945; //the "g0" coefficient in rrs function Lee et al 2002
62 static const double grd2 = 0.1247; //the "g1" coefficient in rrs function Lee et al 2002
63 
64 static double invCosSolz; // 1.0 / cos(solar_zenith_angle)
65 static double invCosSenz; // 1.0/ cos(sensor_viewing_angle)
66 static double depth; // depth at this pixel
67 
68 // variables for reading the benthic netCDF file
69 static int ncfile; // netCDF file ID
70 static int benthicProportionId; // netCDF variable id for benthicProportion
71 static int benthicROutOfBounds; // binary flag, pixel not within supplied benthic reflectance lat/lon grid
72 static int benthicRFileExist; // binary flag, benthic reflectance file supplied to l2gen or not.
73 static size_t numLat; // number of latitudes in the benthicProportion array
74 static size_t numLon; // number of longitudes in the benthicProportion array
75 static size_t numBottomTypes; // number of bottom types in the benthicProportion array
76 static double upperLat; // latitude of benthicProportion[0][0][bottomType]
77 static double lowerLat; // latitude of benthicProportion[numLat-1][0][bottomType]
78 static double deltaLat; // latitude grid spacing
79 static double leftLon; // longitude of benthicProportion[0][0][bottomType]
80 static double rightLon; // longitude of benthicProportion[0][numLon-1][bottomType]
81 static double deltaLon; // longitude grid spacing
82 
83 static double* reflectance; // bottom reflectance for each sensor wavelength (reflectance[bottomType][band])
84 static double* benthicRefl; // bottom percent for this pixel (benthicRefl[band])
85 
86 //Hard code in a regional chlorophyll-specific phytoplankton spectra
87 //for tropical waters (collected in southern GBR).
88 static float taphw[10] = {412.0, 443.0, 469.0, 488.0, 531.0, 551.0, 555.0,
89  645.0, 667.0, 678.0};
90 static float taphs[10] = {0.840353901, 1.0, 0.821449702, 0.689387045,
91  0.217028284, 0.14646583, 0.134941382, 0.122233711, 0.291530253,
92  0.356272348};
93 
94 //Hard code in a default benthic albedo in the case there is no default global albedo file found,
95 //or data is out of lat/lon bounds.
96 static float reflw[10] = {412.0, 443.0, 469.0, 488.0, 531.0, 551.0, 555.0,
97  645.0, 667.0, 678.0};
98 static float refls[10] = {0.167663599, 0.197465145, 0.224745351, 0.24005827,
99  0.28106442, 0.297352629, 0.30027302, 0.337329977, 0.309064323, 0.310894269};
100 
101 /* ------------------------------------------------------------------------*/
102 /* swim_ran() - determine if swim has been run for scanline */
103 /* ------------------------------------------------------------------------*/
104 int swim_ran(int recnum)
105 {
106  if ( recnum == LastRecNum )
107  return 1;
108  else
109  return 0;
110 }
111 
112 /* ------------------------------------------------------------------------*/
113 /* getAttr() - get an attribute from the global ncfile. Use NC_GLOBAL for */
114 /* varid to get global attributes */
115 
116 /* ------------------------------------------------------------------------*/
117 double getAttr(int varid, char* attrName) {
118  double val;
119  if (nc_get_att_double(ncfile, varid, attrName, &val) != NC_NOERR) {
120  fprintf(stderr,
121  "-E- %s line %d: could get attribute %s from netCDF File.\n",
122  __FILE__, __LINE__, attrName);
123  exit(1);
124  }
125  return val;
126 }
127 
128 /* -----------------------------------------------------------------------*/
129 /* getVarId() - get variable ID from the global ncfile. */
130 
131 /* -----------------------------------------------------------------------*/
132 int getVarId(char* name) {
133  int varid;
134  if (nc_inq_varid(ncfile, name, &varid) != NC_NOERR) {
135  fprintf(stderr, "-E- %s line %d: could not find %s in netCDF File.\n",
136  __FILE__, __LINE__, name);
137  exit(1);
138  }
139  return varid;
140 }
141 
142 /* -----------------------------------------------------------------------*/
143 /* getDimensionId() - get the dimension IDs from a variable */
144 
145 /* -----------------------------------------------------------------------*/
146 void getDimensionIds(int varid, int* dimIds) {
147  if (nc_inq_vardimid(ncfile, varid, dimIds) != NC_NOERR) {
148  char name[NC_MAX_NAME];
149  nc_inq_varname(ncfile, varid, name);
150  fprintf(stderr,
151  "-E- %s line %d: could not find dimension Ids of %s in netCDF File.\n",
152  __FILE__, __LINE__, name);
153  exit(1);
154  }
155 }
156 
157 /* -----------------------------------------------------------------------*/
158 /* getDimensionLength() - get the length of the dimension ID */
159 
160 /* -----------------------------------------------------------------------*/
161 size_t getDimensionLength(int dimId) {
162  size_t length;
163  if (nc_inq_dimlen(ncfile, dimId, &length) != NC_NOERR) {
164  char name[NC_MAX_NAME];
165  nc_inq_dim(ncfile, dimId, name, &length);
166  fprintf(stderr,
167  "-E- %s line %d: could not get size of dimension \"%s\" in netCDF File.\n",
168  __FILE__, __LINE__, name);
169  exit(1);
170  }
171  return length;
172 }
173 
174 /* -----------------------------------------------------------------------*/
175 /* initBenthicFile() - init the global variables from the bottom */
176 /* reflectance netCDF file */
177 
178 /* -----------------------------------------------------------------------*/
179 void initBenthicFile(char* fileName) {
180  int varid;
181  int dimIds[3];
182  size_t bottom;
183  int band;
184  size_t numWavelengths; // number of wavelengths in reflectance array
185  double firstWavelength; // wavelength of reflectance[0]
186  double deltaWavelength; // delta between each wavelength of reflectance
187  double reflectanceScale; // reflectance scale factor
188  double reflectanceOffset; // reflectance offset
189 
190  //Was a benthic reflectance filename supplied?
191  if (strlen(fileName) != 0) {
192  printf("\n");
193  printf("Reading benthic reflectance from: %s \n", fileName);
194  printf("\n");
195  benthicRFileExist = 1;
196  } else {
197  printf("\n");
198  printf("-E- No benthic reflectance file supplied: use default benthic albedo spectra \n");
199  printf("\n");
200  benthicRFileExist = 0;
201  //Return out of function
202  return;
203  }
204 
205  //If file exists and cannot be read, exit and print error
206  if (benthicRFileExist) {
207  if (nc_open(fileName, NC_NOWRITE, &ncfile) != NC_NOERR) {
208  fprintf(stderr, "-E- %s line %d: could not open netCDF File \"%s\".\n",
209  __FILE__, __LINE__, fileName);
210  exit(1);
211  }
212  }
213 
214  lowerLat = getAttr(NC_GLOBAL, "lower_lat");
215  upperLat = getAttr(NC_GLOBAL, "upper_lat");
216  leftLon = getAttr(NC_GLOBAL, "left_lon");
217  rightLon = getAttr(NC_GLOBAL, "right_lon");
218 
219  // normalize the lons
220  while (leftLon > 180.0)
221  leftLon -= 360.0;
222  while (leftLon < -180.0)
223  leftLon += 360.0;
224  while (rightLon > 180.0)
225  rightLon -= 360.0;
226  while (rightLon < -180.0)
227  rightLon += 360.0;
228  while (rightLon <= leftLon)
229  rightLon += 360.0;
230 
231  firstWavelength = getAttr(NC_GLOBAL, "firstWavelength");
232  deltaWavelength = getAttr(NC_GLOBAL, "deltaWavelength");
233 
234  // make sure deltaWavelength is valid
235  if (deltaWavelength <= 0) {
236  fprintf(stderr,
237  "-E- %s line %d: deltaWavelength must be > 0 in netCDF File \"%s\".\n",
238  __FILE__, __LINE__, fileName);
239  exit(1);
240  }
241 
242  varid = getVarId("reflectance");
243 
244  // get reflectance scale and offset
245  reflectanceScale = 1.0 / getAttr(varid, "scale_factor");
246  reflectanceOffset = getAttr(varid, "add_offset");
247 
248  getDimensionIds(varid, dimIds);
249  numBottomTypes = getDimensionLength(dimIds[0]);
250  numWavelengths = getDimensionLength(dimIds[1]);
251 
252  // allocate global variable
253  reflectance = (double*) allocateMemory(
254  numBottomTypes * nbandVis * sizeof (double), "reflectance");
255 
256  size_t start[2];
257  size_t count[2];
258  count[0] = 1;
259 
260  // load up reflectance for each bottom type for each sensor wavelength
261  for (bottom = 0; bottom < numBottomTypes; bottom++) {
262  start[0] = bottom;
263  for (band = 0; band < nbandVis; band++) {
264 
265  // average an 11nm band width
266  if (lambda[band] < firstWavelength) {
267  reflectance[bottom * nbandVis + band] = 0.0;
268  continue;
269  }
270  int firstIndex = round(
271  (lambda[band] - 5 - firstWavelength) / deltaWavelength);
272  if (firstIndex < 0)
273  firstIndex = 0;
274  if (firstIndex >= numWavelengths) {
275  reflectance[bottom * nbandVis + band] = 0.0;
276  continue;
277  }
278  int numIndex = round(11.0 / deltaWavelength);
279  if (firstIndex + numIndex >= numWavelengths) {
280  numIndex = numWavelengths - firstIndex;
281  }
282  if (numIndex < 1)
283  numIndex = 1;
284 
285  start[1] = firstIndex;
286  count[1] = numIndex;
287 
288  ushort tmpShort[numIndex];
289  if (nc_get_vara_ushort(ncfile, varid, start, count,
290  tmpShort) != NC_NOERR) {
291  fprintf(stderr,
292  "-E- %s line %d: could not read values of reflectance in netCDF File \"%s\".\n",
293  __FILE__, __LINE__, fileName);
294  exit(1);
295  }
296 
297  // average the 11nm wide readings and scale
298  int i;
299  double sum = 0.0;
300  for (i = 0; i < numIndex; i++)
301  sum += tmpShort[i];
302  reflectance[bottom * nbandVis + band] = sum / numIndex
303  * reflectanceScale + reflectanceOffset;
304 
305  } // for bands
306  } // for bottoms
307 
308  // get dimensions of the benthicProportion variable
309  benthicProportionId = getVarId("benthicProportion");
310 
311  getDimensionIds(benthicProportionId, dimIds);
312  numLat = getDimensionLength(dimIds[0]);
313  numLon = getDimensionLength(dimIds[1]);
314  size_t tmpSize = getDimensionLength(dimIds[2]);
315 
316  if (tmpSize != numBottomTypes) {
317  fprintf(stderr,
318  "-E- %s line %d: numBottomTypes in benthicProportion is not the same as in reflectance in netCDF File \"%s\".\n",
319  __FILE__, __LINE__, fileName);
320  exit(1);
321  }
322 
323  deltaLat = (upperLat - lowerLat) / numLat;
324  deltaLon = (rightLon - leftLon) / numLon;
325 
326 }
327 
328 /* -----------------------------------------------------------------------*/
329 /* getDefaultBenthicR() - get default benthic reflectance spectra and */
330 /* interpolate to sensor wavelengths */
331 
332 /* -----------------------------------------------------------------------*/
333 void getDefaultBenthicR() {
334 
335  static int32_t refln = -1;
336  int band;
337 
338  /* read default benthic reflectance table size */
339  refln = 0;
340  for (refln = 0; refln < nbandVis; refln++) {
341  if (reflw[refln] < 0) {
342  break;
343  }
344  }
345 
346  /*Interploate the hard-coded benthic albedo to sensor wavelenghth*/
347  for (band = 0; band < nbandVis; band++) {
348  benthicRefl[band] = linterp(reflw, refls, refln, lambda[band]);
349  }
350 
351 }
352 
353 /* -----------------------------------------------------------------------*/
354 /* getBottomReflectance() - get the benthic reflectance for this pixel */
355 /* into benthicRefl global array */
356 
357 /* -----------------------------------------------------------------------*/
358 void getBottomReflectance(float lat, float lon) {
359  int latIndex;
360  int lonIndex;
361  int latFlag = 0;
362  int lonFlag = 0;
363 
364  int band;
365  int bottom;
366  double sum;
367 
368  // see if coordinate is covered by the file
369  if (lat < lowerLat || lat > upperLat) {
370  latFlag = 1;
371  }
372 
373  // see if coordinate is covered by the file
374  if (lon < leftLon || lon > rightLon) {
375  lonFlag = 1;
376  }
377 
378  //If the pixel is outside the lat/lon range of the user input benthic albedo file,
379  // a default hard-coded sand benthic albedo is used.
380  if (latFlag || lonFlag) {
381  //Set benthic out of bounds flag to 1;
382  benthicROutOfBounds = 1;
383 
384  //Call function to interpolate default benthic reflectance
385  getDefaultBenthicR();
386 
387  //Break out of getBottomReflectance function
388  return;
389 
390  } else {
391  /*Else the pixel is within the user-supplied benthic albedo file*/
392  //Set benthic out of bounds flag to 0;
393  benthicROutOfBounds = 0;
394  }
395 
396  //printf("benthic out of bounds = %d \n",benthicROutOfBounds);
397  latIndex = (lat - lowerLat) / deltaLat;
398  lonIndex = (lon - leftLon) / deltaLon;
399 
400  size_t start[3];
401  start[0] = latIndex;
402  start[1] = lonIndex;
403  start[2] = 0;
404  size_t count[3];
405  count[0] = 1;
406  count[1] = 1;
407  count[2] = numBottomTypes;
408 
409  uint8_t benthicPercent[numBottomTypes];
410  if (nc_get_vara_uchar(ncfile, benthicProportionId, start, count,
411  benthicPercent) != NC_NOERR) {
412  fprintf(stderr,
413  "-E- %s line %d: could not read values from benthicProportion in netCDF File.\n",
414  __FILE__, __LINE__);
415  exit(1);
416  }
417 
418  for (band = 0; band < nbandVis; band++) {
419  sum = 0.0;
420  for (bottom = 0; bottom < numBottomTypes; bottom++) {
421  sum += reflectance[bottom * nbandVis + band] * benthicPercent[bottom]
422  / 100.0;
423  }
424 
425  benthicRefl[band] = sum;
426  }
427 }
428 
429 /* ------------------------------------------------------------------- */
430 /* swim_func() - optically shallow semianlytical reflectance model */
431 
432 /* ------------------------------------------------------------------- */
433 void swim_func(double *initialParams, double *rrsTotal, int numParams,
434  int numBands, void* dataPtr) {
435 
436  int iw; //iterator
437  double rrsDeep; //deep water rrs term
438  double rrsColumn; //water column rrs term
439  double rrsBenthos; //bottom contribution rrs term
440  double bb; //total backscattering coefficient
441  double ac; //total absorption coefficient
442  double kappa; //kappa = ac + bb
443  double u; //u = bb_total/(kappa)
444  double duC; // Upward attenuation term
445  double duB; // Downward attenuation term
446 
447  //-------------------------------------------//
448  /* Shallow water sub-surface reflectance model*/
449  //-------------------------------------------//
450  // Iterate over wavelengths
451  for (iw = 0; iw < numBands; iw++) {
452 
453  //Calculate the bulk IOPs
454  ac = aw[iw] + initialParams[0] * aphstar[iw]
455  + initialParams[1] * adgstar[iw];
456  bb = bbw[iw] + initialParams[2] * bbpstar[iw];
457 
458 
459  //Attenuation coefficients
460  kappa = bb + ac;
461  u = bb / kappa;
462 
463  //Pathlength elongation factors
464  duC = 1.03 * pow((1 + 2.48 * u), 0.5);
465  duB = 1.04 * pow((1.54 + 5.4 * u), 0.5);
466 
467  //Deep water reflectance
468  rrsDeep = grd1 * u + grd2 * pow(u, 2);
469 
470  //Water column term
471  rrsColumn = rrsDeep
472  * (1.0
473  + exp(
474  -1.0 * kappa * depth
475  * (invCosSolz + (duC * invCosSenz))));
476 
477  //Bottom contribution term
478  rrsBenthos = (benthicRefl[iw] / M_PI)
479  * exp(-1.0 * kappa * depth * (invCosSolz + (duB * invCosSenz)));
480 
481  rrsTotal[iw] = rrsColumn + rrsBenthos;
482 
483  }
484 }
485 
486 /* ---------------------------------------------------------------------- */
487 /* run_swim() - calculates the SWIM model for full scanline */
488 
489 /* ---------------------------------------------------------------------- */
490 void run_swim(l2str *l2rec) {
491 
492  //Declare variables
493  static int32_t taphn = -1;
494  //static float *taphw = NULL;
495  //static float *taphs = NULL;
496  static float adg_s = -1;
497  static float bbp_s = -1;
498 
499  //Iterators indexes
500  int i, iw;
501  int validr;
502  int statusLM;
503  l1str *l1rec = l2rec->l1rec;
504  filehandle *l1file = l1rec->l1file;
505  double Rrs; //above-water observed remote sensing reflectance
506  double rrs_s[l1rec->l1file->nbands]; //Array holding sub-surface observed remote sensing reflectance
507  // double rrs_a[l1rec->l1file->nbands]; //Array holding above-water observed remote sensing reflectance
508  double fitParams[NPAR]; //Array for holding initial guesses for aph440, adg440, bbp550
509  double opts[LM_OPTS_SZ]; //Atray for holding Levmar optionsfloat get_benthic_r(float wave, float chl)
510  double info[LM_INFO_SZ]; //Array for holding Levmar info
511 
512  float aph443; //the value of phyto abs at 443 nm
513  float bbp443; //the value of detritus + cdom abs at 443 nm
514  float adg443; //the value of particle backscatter 443 nm
515 
516  int defaultRefl; //binary flag, was default reflectance used?
517 
518  // first check to see if this is the first time we have run
519  // so we can allocate memory and init stuff
520  if (swimScanNum == -1) {
521 
522  /* limit to visible wavelengths */
523  nbandVis = rdsensorinfo(l1file->sensorID, l1_input->evalmask,
524  "NbandsVIS", NULL);
525 
526  /*Get the wavelengths - needed for IOP spectral models*/
527  rdsensorinfo(l1file->sensorID, l1_input->evalmask, "Lambda",
528  (void**) &lambda);
529 
530  /*Allocate memory*/
531  adg = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
532  "adg");
533  bbp = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
534  "bbp");
535  aph = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
536  "aph");
537  atot = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
538  "atot");
539  bbtot = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
540  "bbtot");
541  iter = (int16*) allocateMemory(l1rec->npix * nbandVis * sizeof (int16),
542  "iter");
543  aphstar = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
544  "aphstar");
545  adgstar = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
546  "adgstar");
547  bbpstar = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
548  "bbpstar");
549  aw = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
550  "aw");
551  bbw = (float*) allocateMemory(l1rec->npix * nbandVis * sizeof (float),
552  "bbw");
553  benthicRefl = (double*) allocateMemory(
554  l1rec->npix * nbandVis * sizeof (double), "benthicRefl");
555 
556  /*--------------------------------------------------------*/
557  /*IOP models ar defined once here and allocated to memory */
558  /*--------------------------------------------------------*/
559 
560  /*Need to load up the model coefficient shapes (aph*, adg*, bbp*)*/
561  /* set adg slope and bbp power law static coefficients */
562  adg_s = 0.017;
563  bbp_s = 1.0;
564 
565  /* read aphstar table size */
566  taphn = 0;
567  for (taphn = 0; taphn < nbandVis; taphn++)
568  if (taphw[taphn] < 0)
569  break;
570 
571  /*Loop over wavelengths */
572  for (iw = 0; iw < nbandVis; iw++) {
573 
574  /*Phytoplankton absorption shape*/
575  aphstar[iw] = linterp(taphw, taphs, taphn, lambda[iw]);
576 
577  //printf("aphi = %f @ lambda = %d \n", aphstar[iw], lambda[iw] );
578 
579  /*Coloured dissolved and detrital matter absorption shape*/
580  adgstar[iw] = exp(-adg_s * (lambda[iw] - 443.0));
581 
582  /*Particulate matter backscatter shape*/
583  bbpstar[iw] = pow((443.0 / lambda[iw]), bbp_s);
584 
585  /*Pure water absorption coefficient*/
586  aw[iw] = aw_spectra(lambda[iw], BANDW);
587 
588  /*Pure water backscatter coefficient*/
589  bbw[iw] = bbw_spectra(lambda[iw], BANDW);
590  }
591 
592  /*Begin reading benthic reflectance data*/
593  initBenthicFile(input->breflectfile);
594  }
595 
596  /*record that we have run the model on this scan line*/
597  swimScanNum = l1rec->iscan;
598 
599  /*----------------------------------------*/
600  /* Iterate through pixels */
601  /*----------------------------------------*/
602  for (i = 0; i < l1rec->npix; i++) {
603 
604  // init the spectral outputs
605  for (iw = 0; iw < nbandVis; iw++) {
606  adg[i * nbandVis + iw] = BAD_FLT;
607  bbp[i * nbandVis + iw] = BAD_FLT;
608  aph[i * nbandVis + iw] = BAD_FLT;
609  atot[i * nbandVis + iw] = BAD_FLT;
610  bbtot[i * nbandVis + iw] = BAD_FLT;
611  iter[i * nbandVis + iw] = 0;
612  }
613 
614  //Set status to default value of 0.
615  validr = 0;
616 
617  //Use the quality control from GSM to flag bad pixelsparamFit
618  if (l1rec->mask[i] == 0) {
619 
620  /* Solar and Sensor Viewing Geometry */
621  double solzRad = (l1rec->solz[i] * M_PI) / 180.0; /*Convert to radians*/
622  double senzRad = (l1rec->senz[i] * M_PI) / 180.0; /*convert to radians*/
623  invCosSolz = 1.0 / cos(solzRad); /*Check rad or deg */
624  invCosSenz = 1.0 / cos(senzRad); /*Check rad or deg */
625 
626  depth = 0 - l1rec->dem[i];
627 
628  for (iw = 0; iw < nbandVis; iw++) {
629  Rrs = l2rec->Rrs[i * l1file->nbands + iw];
630  // rrs_a[iw] = Rrs;
631  //Sub-surface remote sensing reflectance
632  rrs_s[iw] = Rrs / (0.52 + 1.7 * Rrs);
633  if (rrs_s[iw] >= 0.0)
634  validr = validr + 1;
635  }
636 
637  /* if less than 3 valid radiances, fail pixel */
638  if (validr <= 3) {
639  l1rec->flags[i] |= PRODFAIL;
640  continue;
641  }
642 
643  /*---------------------------------------------------*/
644  /* fill up the benthicRefl global array */
645  defaultRefl = -1;
646  if (!benthicRFileExist) {
647 
648  if (defaultRefl < 0) {
649  //call function to interpolate default reflectance
650  getDefaultBenthicR();
651 
652  //Assign this flag so the default reflectance does not
653  //need to be calculated multiple times.
654  defaultRefl = 1;
655  }
656 
657  } else {
658  getBottomReflectance(l1rec->lat[i], l1rec->lon[i]);
659  }
660 
661 
662  /*If the hard-coded reflectance is used, set l2flag to PRODWARN
663  if ( benthicROutOfBounds ) {
664  l2rec->flags[i] |= PRODWARN;
665  };*/
666 
667  /*---------------------------------------------------*/
668  /* fit model for this pixel using levmar optimisation*/
669  /*---------------------------------------------------*/
670 
671  /*Initial guess (LM starting points)*/
672  fitParams[0] = 0.2;
673  fitParams[1] = 0.01;
674  fitParams[2] = 0.001;
675 
676  /*Optimisation control parameters*/
677  opts[0] = 1E-3; //1E-3 LM_INIT_MU; /* scale factor for initial mu */
678  opts[1] = 1E-10; //1E-10 LM_INIT_MU; convergence in gradient
679  opts[2] = 1E-5; //1E-5 relative error desired in the approximate solution
680  opts[3] = 1E-7; //1E-7 relative error desired in the sum of squares
681  opts[4] = 1E-6; //1E-6 LM_DIFF_DELTA; /*step used in difference approximation to the Jacobian*/
682 
683  /*Set box constraints*/
684  /*L/U Bounds inferred from Werdell et al. 2013*/
685  double lowerBounds[NPAR] = {-0.0003755, -0.0003755, -0.0001195625}; //Lower bounds for three parameters
686  double upperBounds[NPAR] = {5.0, 5.0, 5.0}; //Upper bounds for three parameters
687  int m = NPAR;
688  int n = nbandVis;
689  int maxit = MAXITR;
690 
691  /*Run optimization*/
692  /*levmar routine (dlevmar), box contstraints (bc), numerical partials (diff)*/
693  //For further details, see source code in: lmbc_core.c //
694  statusLM = dlevmar_bc_dif(swim_func, fitParams, rrs_s, m, n,
695  lowerBounds, upperBounds, NULL, maxit, opts, info, NULL,
696  NULL, NULL);
697 
698  /*---------------------------------------------------*/
699  /* Write products to arrays */
700  /*---------------------------------------------------*/
701 
702  /* save results and flag failure */
703  //printf("num iteration %f \n",info[5]);
704  //if (statusLM <= 0 || info[5] == maxit) {
705  if (statusLM <= 0) {
706  l1rec->flags[i] |= PRODFAIL;
707 
708  } else {
709  /*Output results and calculate IOPs*/
710  aph443 = fitParams[0]; //absorption of phytoplankton at 443 nm band
711  adg443 = fitParams[1]; //absorption of detritus + cdom at 443 nm band
712  bbp443 = fitParams[2]; //particle backscattering at 443 nm band
713 
714  /*calculate spectral IOPs for output*/
715  for (iw = 0; iw < nbandVis; iw++) {
716 
717  //absorption of phytoplankton
718  //Note: have normalised aphi440=1.0 - fix later with regional aphi
719  //Free parameter in model thus aph443 rather than Chl
720  aph[i * nbandVis + iw] = aph443 * aphstar[iw];
721  adg[i * nbandVis + iw] = adg443 * adgstar[iw]; //absorption detritus + cdom
722  bbp[i * nbandVis + iw] = bbp443 * bbpstar[iw]; //backscattering of particles
723  atot[i * nbandVis + iw] = aw[iw] + aph443 * aphstar[iw]
724  + adg443 * adgstar[iw]; //total absorption
725  bbtot[i * nbandVis + iw] = bbw[iw] + bbp443 * bbpstar[iw]; //total backscattering
726 
727  }
728  }
729 
730  }
731  }
732 
733  LastRecNum = l1rec->iscan;
734 
735 }
736 
737 /* -------------------------------------------------------------------- */
738 /* get_swim() - returns requested swim product for full scanline */
739 
740 /* ---------------------------------------------------------------------*/
741 void get_swim(l2str *l2rec, l2prodstr *p, float prod[]) {
742 
743  int32_t ip;
744  l1str *l1rec = l2rec->l1rec;
745 
746  if (!swim_ran(l1rec->iscan)) {
747  run_swim(l2rec);
748  }
749 
750  int bandNum = p->prod_ix;
751  if (bandNum < 0) {
752  printf("-E- %s line %d : prod_ix must be positive, prod_ix=%d\n",
753  __FILE__, __LINE__, p->prod_ix);
754  exit(1);
755  }
756  if (bandNum >= nbandVis) {
757  printf(
758  "-E- %s line %d : prod_ix must be less than numBands=%d, prod_ix=%d\n",
759  __FILE__, __LINE__, nbandVis, p->prod_ix);
760  exit(1);
761  }
762 
763  for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
764 
765  switch (p->cat_ix) {
766 
767  case CAT_a_swim:
768  prod[ip] = atot[ip * nbandVis + bandNum];
769  break;
770 
771  case CAT_bb_swim:
772  prod[ip] = bbtot[ip * nbandVis + bandNum];
773  break;
774 
775  case CAT_adg_swim:
776  prod[ip] = adg[ip * nbandVis + bandNum];
777  break;
778 
779  case CAT_aph_swim:
780  prod[ip] = aph[ip * nbandVis + bandNum];
781  break;
782 
783  case CAT_bbp_swim:
784  prod[ip] = bbp[ip * nbandVis + bandNum];
785  break;
786 
787  default:
788  printf("-E- %s line %d : erroneous product ID %d passed to swim.\n",
789  __FILE__, __LINE__, p->cat_ix);
790  exit(1);
791  }
792  }
793 }
794 
795 /* ------------------------------------------------------------------- */
796 /* interface to convl12() to return SWIM bulk iops */
797 
798 /* ------------------------------------------------------------------- */
799 void iops_swim(l2str *l2rec) {
800 
801  int ib, ip;
802  l1str *l1rec = l2rec->l1rec;
803 
804  if (!swim_ran(l1rec->iscan)) {
805  run_swim(l2rec);
806  }
807 
808  for (ip = 0; ip < l1rec->npix; ip++) {
809  for (ib = 0; ib < nbandVis; ib++) {
810  l2rec->a[ip * l1rec->l1file->nbands + ib] = atot[ip * nbandVis + ib];
811  l2rec->bb[ip * l1rec->l1file->nbands + ib] = bbtot[ip * nbandVis + ib];
812 
813  }
814  }
815 
816  return;
817 }
818 
819 /* ------------------------------------------------------------------- */
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
color_dtdb.bottom
bottom
Definition: color_dtdb.py:376
MAXITR
#define MAXITR
Definition: swim.c:36
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
NPAR
#define NPAR
Definition: swim.c:35
getBottomReflectance
void getBottomReflectance(float lat, float lon)
Definition: swim.c:358
maxit
int maxit
Definition: pml_iop_tables.c:42
getDimensionLength
size_t getDimensionLength(int dimId)
Definition: swim.c:161
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
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
initBenthicFile
void initBenthicFile(char *fileName)
Definition: swim.c:179
CAT_bbp_swim
#define CAT_bbp_swim
Definition: l2prod.h:280
PRODFAIL
#define PRODFAIL
Definition: l2_flags.h:41
lambda
data_t lambda[NROOTS+1]
Definition: decode_rs.h:75
get_swim
void get_swim(l2str *l2rec, l2prodstr *p, float prod[])
Definition: swim.c:741
CAT_a_swim
#define CAT_a_swim
Definition: l2prod.h:276
input
instr * input
Definition: main_l2binmatch.cpp:27
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
l12_proto.h
run_local.lat
lat
Definition: run_local.py:69
getAttr
double getAttr(int varid, char *attrName)
Definition: swim.c:117
CAT_aph_swim
#define CAT_aph_swim
Definition: l2prod.h:279
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
CAT_adg_swim
#define CAT_adg_swim
Definition: l2prod.h:278
ushort
unsigned short ushort
Definition: l1_viirs_l1b.c:16
getDefaultBenthicR
void getDefaultBenthicR()
Definition: swim.c:333
getDimensionIds
void getDimensionIds(int varid, int *dimIds)
Definition: swim.c:146
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
create_images.Rrs
def Rrs
Definition: create_images.py:92
ac
float * ac[MAX_BANDS]
Definition: pml_iop_tables.c:26
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
swim_ran
int swim_ran(int recnum)
Definition: swim.c:104
getVarId
int getVarId(char *name)
Definition: swim.c:132
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
u
data_t u
Definition: decode_rs.h:74
swim_func
void swim_func(double *initialParams, double *rrsTotal, int numParams, int numBands, void *dataPtr)
Definition: swim.c:433
rdsensorinfo
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
BANDW
#define BANDW
Definition: l1.h:53
CAT_bb_swim
#define CAT_bb_swim
Definition: l2prod.h:277
i
int i
Definition: decode_rs.h:71
iops_swim
void iops_swim(l2str *l2rec)
Definition: swim.c:799
run_swim
void run_swim(l2str *l2rec)
Definition: swim.c:490
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
count
int count
Definition: decode_rs.h:79


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