ocssw V2020
water_spectra.c
Go to the documentation of this file.
1 /* ============================================================================================== */
2 /* module water.c - functions to read and return absorption and backscatter of pure sea water */
3 /* */
4 /* B. Franz, NASA/GSFC Ocean Color Discipline Processing Group, Sep. 2004 */
5 /* */
6 /* ============================================================================================== */
7 
8 #include "l12_proto.h"
9 
10 #define MINAWTAB 200
11 #define MAXAWTAB 2449
12 #define INTAWTAB 1
13 #define NAWTAB ((MAXAWTAB-MINAWTAB)/INTAWTAB + 1)
14 
15 static double awtab [NAWTAB];
16 static double bbwtab[NAWTAB];
17 static int32_t ntab = NAWTAB;
18 static int min_wl = MINAWTAB;
19 static int del_wl = INTAWTAB;
20 
21 
22 /* ---------------------------------------------------------------------------------------------- */
23 /* read_water_spectra() - called once to load look-up table static arrays */
24 /* ---------------------------------------------------------------------------------------------- */
25 void read_water_spectra(void) {
26  // TODO: Dynamically allocate awtab and bbwtab arrays
27  static int firstCall = 1;
28 
29  FILE *fp;
30  char *line;
31  int32_t i, status, netcdf_input;
32  float wl = 0, aw = 0, bw = 0;
33  int ncid, varid;
34 
35  if (!firstCall) return;
36 
37  /* Does the file exist? */
38  if (access(input->water_spectra_file, F_OK) || access(input->water_spectra_file, R_OK)) {
39  printf("-E- %s: water_spectra_file '%s' does not exist or cannot open.\n",
40  __FILE__, input->water_spectra_file);
41  exit(EXIT_FAILURE);
42  }
43 
44  /* test for NetCDF input file */
45  status = nc_open(input->water_spectra_file, NC_NOWRITE, &ncid);
46  netcdf_input = (status == NC_NOERR);
47  if (netcdf_input) {
48  /* Get the aw data. */
49  if ((nc_inq_varid(ncid, "aw", &varid)) == NC_NOERR) {
50  nc_get_var_double(ncid, varid, &awtab[0]);
51  } else {
52  printf("-E- %s: water_spectra_file '%s' does not have aw.\n",
53  __FILE__, input->water_spectra_file);
54  exit(EXIT_FAILURE);
55  }
56  /* Get the bw data. */
57  if ((nc_inq_varid(ncid, "bw", &varid)) == NC_NOERR) {
58  nc_get_var_double(ncid, varid, &bbwtab[0]);
59  for (i = 0; i < ntab; i++) {
60  bbwtab[i] *= 0.5;
61  }
62  } else {
63  printf("-E- %s: water_spectra_file '%s' does not have bw.\n",
64  __FILE__, input->water_spectra_file);
65  exit(EXIT_FAILURE);
66  }
67 
68  nc_close(ncid);
69  } else {
70  /* This reader is for the "old" SeaBASS formatted ASCII input
71  * It still suffers from the "can't handle blank lines or missing
72  * comment character" problem
73  * Deprecate as soon as netCDF version is official
74  */
75  if ((fp = fopen(input->water_spectra_file, "r")) == NULL) {
76  fprintf(stderr,
77  "-E- %s line %d: unable to open %s for reading\n", __FILE__, __LINE__, input->water_spectra_file);
78  exit(1);
79  }
80 
81  line = (char *) calloc(80, sizeof (char));
82  i = 0;
83  while (i < ntab) {
84  if (fgets(line, 80, fp) == NULL) {
85  fprintf(stderr, "-E- %s line %d: error reading %s at line %d\n", __FILE__, __LINE__, input->water_spectra_file, i);
86  exit(1);
87  }
88  if (line[0] != '/' && line[0] != '!') {
89  sscanf(line, "%f %f %f", &wl, &aw, &bw);
90  awtab [i] = aw;
91  bbwtab[i] = bw / 2.0;
92  i++;
93  }
94  }
95  free(line);
96  }
97  firstCall = 0;
98 }
99 
100 
101 /* ---------------------------------------------------------------------------------------------- */
102 /* aw_spectra() - returns water absorption at wavelength, wl, averaged over bandwidth, width */
103 
104 /* ---------------------------------------------------------------------------------------------- */
105 float aw_spectra(int32_t wl, int32_t width) {
106  static int firstCall = 1;
107 
108  int32_t itab = (wl - min_wl) / del_wl;
109  int32_t imin = MAX(itab - width / 2 / del_wl, 0);
110  int32_t imax = MIN(itab + width / 2 / del_wl, ntab);
111  float aw = 0;
112  int32_t i;
113 
114  if (firstCall) {
116  firstCall = 0;
117  }
118 
119  if (itab < 0) {
120  //fprintf(stderr,
121  // "-W- %s line %d: wavelength of %d outside aw table range.\n",__FILE__,__LINE__,wl);
122  itab = 0;
123  } else if (itab > ntab - 1) {
124  //fprintf(stderr,
125  // "-W- %s line %d: wavelength of %d outside aw table range.\n",__FILE__,__LINE__,wl);
126  itab = ntab - 1;
127  }
128 
129  for (i = imin; i <= imax; i++)
130  aw += (float) awtab[i];
131 
132  aw /= (imax - imin + 1);
133 
134  return (aw);
135 }
136 
137 
138 /* ---------------------------------------------------------------------------------------------- */
139 /* bbw_spectra() - returns water backscatter at wavelength, wl, averaged over bandwidth, width */
140 
141 /* ---------------------------------------------------------------------------------------------- */
142 float bbw_spectra(int32_t wl, int32_t width) {
143  static int firstCall = 1;
144 
145  int32_t itab = (wl - min_wl) / del_wl;
146  int32_t imin = MAX(itab - width / 2 / del_wl, 0);
147  int32_t imax = MIN(itab + width / 2 / del_wl, ntab);
148  float bbw = 0;
149  int32_t i;
150 
151  if (firstCall) {
153  firstCall = 0;
154  }
155 
156  if (itab < 0) {
157  //fprintf(stderr,
158  // "-W- %s line %d: wavelength of %d outside bbw table range.\n",__FILE__,__LINE__,wl);
159  itab = 0;
160  } else if (itab > ntab - 1) {
161  //fprintf(stderr,
162  // "-W- %s line %d: wavelength of %d outside bbw table range.\n",__FILE__,__LINE__,wl);
163  itab = ntab - 1;
164  }
165 
166  for (i = imin; i <= imax; i++)
167  bbw += (float) bbwtab[i];
168 
169  bbw /= (imax - imin + 1);
170 
171  return (bbw);
172 }
173 
174 
175 /* ---------------------------------------------------------------------------------------------- */
176 /* returns aw and bbw at specified "sensor" wavelengths, appropriate to the derived nLw */
177 
178 /* ---------------------------------------------------------------------------------------------- */
179 void get_aw_bbw(l2str *l2rec, float wave[], int nwave, float *aw, float *bbw) {
180  int ib, iw;
181 
182  if (input->outband_opt >= 2) {
183  for (ib = 0; ib < nwave; ib++) {
184  aw [ib] = l2rec->l1rec->sw_a[ib];
185  bbw[ib] = l2rec->l1rec->sw_bb[ib];
186  }
187  } else {
188  float *senaw = l2rec->l1rec->l1file->aw;
189  float *senbbw = l2rec->l1rec->l1file->bbw;
190  for (ib = 0; ib < nwave; ib++) {
191  iw = bindex_get(wave[ib]);
192  aw [ib] = senaw [iw];
193  bbw[ib] = senbbw[iw];
194  }
195  }
196 }
#define MAX(A, B)
Definition: swl0_utils.h:26
#define MIN(x, y)
Definition: rice.h:169
int status
Definition: l1_czcs_hdf.c:31
#define MINAWTAB
Definition: water_spectra.c:10
#define NULL
Definition: decode_rs.h:63
float bbw_spectra(int32_t wl, int32_t width)
#define INTAWTAB
Definition: water_spectra.c:12
void get_aw_bbw(l2str *l2rec, float wave[], int nwave, float *aw, float *bbw)
int bindex_get(int32_t wave)
Definition: windex.c:43
instr * input
void read_water_spectra(void)
Definition: water_spectra.c:25
Definition: common.h:9
#define NAWTAB
Definition: water_spectra.c:13
float aw_spectra(int32_t wl, int32_t width)
int i
Definition: decode_rs.h:71