OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
bin_climatology.c
Go to the documentation of this file.
1 #include "l12_proto.h"
2 
3 static int32 *basebin;
4 static int32 *numbin;
5 static int32 totbins;
6 static int32 nrows;
7 
8 #define MAXNVDATA 50
9 
10 #define MAXCPROD 3
11 #define CHLORA 0
12 #define NLW412 1
13 #define ANGSTROM 2
14 
15 static int32 fid;
16 static int32 sdfid;
17 static int32 vgid;
18 static int32 vdata_id[MAXNVDATA];
19 
20 static int32_t *bins;
21 static float *data [MAXCPROD];
22 static int32 n_records;
23 
24 static char pnames[MAXCPROD][20] = {"chlor_a", "nLw_412", "angstrom"};
25 static int loaded[MAXCPROD] = {0, 0, 0};
26 
27 /* ---------------------------------------------------------------------------------------- */
28 
29 /* ---------------------------------------------------------------------------------------- */
30 int32 open_l3b(char *l3file, int32_t day) {
31  static int nmonths = 12;
32  static int firstCall = 1;
33  static int clim_day [12] = {1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335};
34  static char clim_files[12][19] ={"S001031_L3b_MC.hdf", "S032059_L3b_MC.hdf", "S060090_L3b_MC.hdf",
35  "S091120_L3b_MC.hdf", "S121151_L3b_MC.hdf", "S152181_L3b_MC.hdf",
36  "S182212_L3b_MC.hdf", "S213243_L3b_MC.hdf", "S244273_L3b_MC.hdf",
37  "S274304_L3b_MC.hdf", "S305334_L3b_MC.hdf", "S335365_L3b_MC.hdf"};
38 
39  char *tmp_str;
40 
41  intn i, j;
42  int32 vg_ref;
43  int32 tag;
44  int32 ref;
45  int32 vdid;
46  char nam_buf[80];
47  char cls_buf[80];
48  char file[FILENAME_MAX];
49 
50 
51  if (l3file != NULL) {
52 
53  strcpy(file, l3file);
54  if (firstCall) printf("\nOpening L3 data file %s\n", file);
55  } else {
56 
57  /* Determine filename from day number */
58  for (i = 1; i < nmonths; i++)
59  if (clim_day[i] > day)
60  break;
61  i--;
62 
63  if ((tmp_str = getenv("OCDATAROOT")) == NULL) {
64  printf("OCDATAROOT environment variable is not defined.\n");
65  exit(1);
66  }
67  strcpy(file, tmp_str);
68  strcat(file, "/common/");
69  strcat(file, clim_files[i]);
70 
71  if (firstCall) printf("\nOpening climatology file %s\n", file);
72  }
73 
74  firstCall = 0;
75 
76 
77  /* Open the file and start the SD interface */
78  if ((fid = Hopen(file, DFACC_RDONLY, 0)) < 0) {
79  printf("-E- %s line %d: Cannot open input HDF file on Hopen\n", __FILE__, __LINE__);
80  printf("-E- Level-3 data file - %s\n", file);
81  return (0);
82  }
83  Vstart(fid);
84  if ((sdfid = SDstart(file, DFACC_RDONLY)) < 0) {
85  printf("-E- %s line %d: Cannot open input HDF file on SDstart- %s\n", __FILE__, __LINE__, file);
86  exit(1);
87  }
88 
89  /* Locate the Binned Data vgroup */
90 
91  vgid = -1;
92  for (i = 0; i < MAXNVDATA; i++) {
93  vdata_id[i] = -1;
94  }
95 
96  vg_ref = Vfind(fid, "Level-3 Binned Data");
97  if (vg_ref < 0) {
98  printf("-E- %s line %d: Cannot locate Level-3 Binned Data in %s\n", __FILE__, __LINE__, file);
99  exit(1);
100  }
101 
102  vgid = Vattach(fid, vg_ref, "r");
103 
104 
105  /* Locate all the required vdata in the vgroup */
106 
107  j = 0;
108  for (i = 0; i < Vntagrefs(vgid); i++) {
109 
110  Vgettagref(vgid, i, &tag, &ref);
111  vdid = VSattach(fid, ref, "r");
112 
113  if (vdid == -1) {
114  printf("-E- %s line %d: Problem opening Vdata (reference # %d)\n", __FILE__, __LINE__, ref);
115  exit(1);
116  }
117 
118  VSgetname(vdid, nam_buf);
119  VSgetclass(vdid, cls_buf);
120 
121  if (strcmp(cls_buf, "DataMain") == 0) {
122  vdata_id[0] = vdid;
123  }
124 
125  if (strcmp(cls_buf, "Index") == 0) {
126  vdata_id[1] = vdid;
127  nrows = VSelts(vdata_id[1]);
128  }
129 
130  if (strcmp(cls_buf, "DataSubordinate") == 0) {
131  vdata_id[2 + j++] = vdid;
132  }
133  }
134 
135  return 1;
136 }
137 
138 
139 /* ---------------------------------------------------------------------------------------- */
140 
141 /* ---------------------------------------------------------------------------------------- */
142 int32 close_l3b() {
143  intn i;
144 
145  if (vgid != -1) {
146 
147  Vdetach(vgid);
148 
149  for (i = 0; i < MAXNVDATA; i++) {
150  if (vdata_id[i] != -1) VSdetach(vdata_id[i]);
151  }
152  }
153 
154  SDend(sdfid);
155  Vend(fid);
156  Hclose(fid);
157 
158  return SUCCEED;
159 }
160 
161 
162 /* ---------------------------------------------------------------------------------------- */
163 /* binary search for nearest bin in bin list */
164 
165 /* ---------------------------------------------------------------------------------------- */
166 int32 nearest_bin(int32 bin_num) {
167  int32 jl = -1, ju = n_records, jm = 0;
168  int ascnd;
169 
170  ascnd = (bins[n_records - 1] >= bins[0]);
171  while (ju - jl > 1) {
172  jm = (ju + jl) / 2;
173  if (ascnd == (bin_num >= bins[jm]))
174  jl = jm;
175  else
176  ju = jm;
177  }
178  if (bin_num == bins[jl]) return (jl);
179  if (jl + 1 < n_records && bin_num == bins[jl + 1]) return (jl + 1);
180  if (jl > 0 && bin_num == bins[jl - 1]) return (jl - 1);
181  if (bin_num == bins[0]) return (0);
182  if (bin_num == bins[n_records - 1]) return (n_records - 1);
183  return (-1);
184 }
185 
186 /* ---------------------------------------------------------------------------------------- */
187 
188 /* ---------------------------------------------------------------------------------------- */
189 void initbin(void) {
190  int row;
191  float32 latbin;
192 
193  /* ----------------- */
194  /* Set up bin arrays */
195  /* ----------------- */
196  if ((numbin = (int32 *) calloc(nrows, sizeof (int32))) == NULL) {
197  printf("-E- %s:%d Error allocating memory to numbin\n", __FILE__, __LINE__);
198  exit(FATAL_ERROR);
199  }
200  if ((basebin = (int32 *) calloc(nrows + 1, sizeof (int32))) == NULL) {
201  printf("-E- %s:%d Error allocating memory to basebin\n", __FILE__, __LINE__);
202  exit(FATAL_ERROR);
203  }
204  basebin[0] = 1;
205  for (row = 0; row < nrows; row++) {
206  latbin = ((row + 0.5)*180.0 / nrows) - 90.0;
207  numbin[row] = (int32) (2 * nrows * cos(latbin * PI / 180.0) + 0.5);
208  if (row > 0) {
209  basebin[row] = basebin[row - 1] + numbin[row - 1];
210  }
211  }
212  totbins = basebin[nrows - 1] + numbin[nrows - 1] - 1;
213 }
214 
215 /* ---------------------------------------------------------------------------------------- */
216 
217 /* ---------------------------------------------------------------------------------------- */
218 void closebin(void) {
219  free(numbin);
220  free(basebin);
221 }
222 
223 /* ---------------------------------------------------------------------------------------- */
224 
225 /* ---------------------------------------------------------------------------------------- */
226 double constrain_lat(double lat) {
227  if (lat > 90) lat = 90;
228  if (lat < -90) lat = -90;
229  return (lat);
230 }
231 
232 /* ---------------------------------------------------------------------------------------- */
233 
234 /* ---------------------------------------------------------------------------------------- */
235 double constrain_lon(double lon) {
236  while (lon < -180) lon += 360;
237  while (lon > 180) lon -= 360;
238  return (lon);
239 }
240 
241 /* ---------------------------------------------------------------------------------------- */
242 
243 /* ---------------------------------------------------------------------------------------- */
244 int16 lat2row(double lat) {
245  int16 row;
246 
247  row = (int16) ((90 + lat) * nrows / 180.0);
248  if (row >= nrows) row = nrows - 1;
249  return (row);
250 }
251 
252 /* ---------------------------------------------------------------------------------------- */
253 
254 /* ---------------------------------------------------------------------------------------- */
255 int32 rowlon2bin(int16 row, double lon) {
256  int16 col;
257  int32 bin;
258 
259  lon = constrain_lon(lon);
260  col = (int16) ((lon + 180.0) * numbin[row] / 360.0);
261  if (col >= numbin[row]) col = numbin[row] - 1;
262  bin = basebin[row] + col;
263  return (bin);
264 }
265 
266 /* ---------------------------------------------------------------------------------------- */
267 
268 /* ---------------------------------------------------------------------------------------- */
269 int32 latlon2bin(double lat, double lon) {
270  int16 row, col;
271  int32 bin;
272 
273  /* Constrain latitudes to [-90,90] and longitudes to [-180,180]. */
274  lat = constrain_lat(lat);
275  lon = constrain_lon(lon);
276 
277  row = lat2row(lat);
278  col = (int16) ((lon + 180.0) * numbin[row] / 360.0);
279  if (col >= numbin[row]) col = numbin[row] - 1;
280  bin = basebin[row] + col;
281  return (bin);
282 }
283 
284 /* ---------------------------------------------------------------------------------------- */
285 
286 /* ---------------------------------------------------------------------------------------- */
287 float bin_climatology(char *l3file, int32_t day, float lon, float lat, char *prodname) {
288  static int firstCall = 1;
289  static char lastprod[20] = "";
290  static int nofile = 0;
291  static int32 prodID = -1;
292  static int32 interlace0;
293 
294  int32 interlace;
295  int32 vdata_size;
296  int32 prodinx;
297  char fields[60];
298  char vdata_name[H4_MAX_NC_NAME];
299  char buf[200];
300  int32 vsize;
301  float32 wgt;
302 
303  int status;
304  int32 nread;
305  int32 nprod;
306  int32_t i;
307  int32 bin_num;
308  int32 bin_idx;
309 
310  float32 *sum_buf = NULL;
311  uint8 *bin_buf = NULL;
312 
313 
314  if (strncmp(prodname, lastprod, 8) != 0 || prodID < 0) {
315  prodID = -1;
316  for (i = 0; i < MAXCPROD; i++)
317  if (strncmp(prodname, pnames[i], 8) == 0) {
318  prodID = i;
319  break;
320  }
321  if (prodID == -1) {
322  printf("-E- %s line %d: Unknown climatology product\n", __FILE__, __LINE__);
323  exit(1);
324  }
325  }
326 
327  if (firstCall) {
328 
329  /* get bin information */
330 
331  if (!open_l3b(l3file, day)) {
332  printf("-E- Level-3 %s needed for the correction\n", pnames[prodID]);
333  nofile = 1;
334  firstCall = 0;
335  exit(1);
336  }
337 
338  status = VSsetfields(vdata_id[0], "bin_num,nobs,weights,flags_set");
339  status = VSsetfields(vdata_id[1], "start_num,begin,extent,max");
340 
341  /* allocate space for bin numbers */
342  VSinquire(vdata_id[0], &n_records, &interlace0, fields, &vdata_size, vdata_name);
343  if ((bins = (int32_t *) malloc(n_records * sizeof (int32_t))) == NULL) {
344  printf("-E- %s line %d: Error allocating memory to L3 file bins\n", __FILE__, __LINE__);
345  exit(1);
346  }
347 
348  /* allocate temp space for bin_num (4), nobs (4), weights (4), flags_set (2) */
349  vsize = 2 * sizeof (int32) + sizeof (float32) + sizeof (int16);
350  if ((bin_buf = (uint8 *) calloc(n_records, vsize)) == NULL) {
351  printf("-E- %s line %d: Error allocating memory for L3 file bin info\n", __FILE__, __LINE__);
352  exit(1);
353  }
354 
355  /* get bin info numbers for each filled bin */
356  nread = VSread(vdata_id[0], (uint8 *) bin_buf, n_records, interlace0);
357  if (nread == -1) {
358  printf("-E- %s line %d: Unable to read bin numbers.\n", __FILE__, __LINE__);
359  exit(FATAL_ERROR);
360  }
361 
362  /* copy bin numbers */
363  for (i = 0; i < n_records; i++) {
364  memcpy((void *) &bins[i], (const void *) &bin_buf[i * vsize], 4);
365  }
366 
367  close_l3b();
368  free(bin_buf);
369 
370  firstCall = 0;
371  }
372 
373  if (nofile) return (-999.0);
374 
375  if (!loaded[prodID]) {
376 
377  /* load product into static array */
378 
379  if (l3file != NULL)
380  printf("Loading %s coverage.\n", pnames[prodID]);
381  else
382  printf("Loading %s climatology.\n", pnames[prodID]);
383 
384 
385  if (!open_l3b(l3file, day)) {
386  printf("-E- %s line %d: Error opening climatology file\n", __FILE__, __LINE__);
387  exit(1);
388  }
389 
390  initbin();
391 
392  /* Calculate the number of L3 products */
393  /* ----------------------------------- */
394  nprod = MAXNVDATA - 2;
395  for (i = 0; i < MAXNVDATA; i++) {
396  if (vdata_id[i] == -1) nprod--;
397  }
398 
399  if (nprod < 1) {
400  printf("-E- %s: There are too few products in the L3 file\n", __FILE__);
401  exit(1);
402  }
403 
404  status = VSsetfields(vdata_id[0], "bin_num,nobs,weights,flags_set");
405  status = VSsetfields(vdata_id[1], "start_num,begin,extent,max");
406 
407  /* locate product vdata & set field names */
408  for (i = 0; i < nprod; i++) {
409  VSgetname(vdata_id[2 + i], fields);
410  if (strncmp(fields, pnames[prodID], 8) == 0) {
411  prodinx = i;
412  strcpy(buf, fields);
413  strcat(buf, "_sum,");
414  strcat(buf, fields);
415  strcat(buf, "_sum_sq");
416  status = VSsetfields(vdata_id[2 + i], buf);
417  if (status != 0) {
418  printf("-E- %s: Failed to set field for %s.\n", __FILE__, pnames[prodID]);
419  exit(1);
420  }
421  break;
422  }
423  }
424  if (i == nprod) {
425  printf("-E- %s: Failed to locate %s.\n", __FILE__, pnames[prodID]);
426  exit(1);
427  }
428 
429  /* allocate space for product mean */
430  VSinquire(vdata_id[2 + prodinx], &n_records, &interlace, fields, &vdata_size, vdata_name);
431  if ((data[prodID] = (float32 *) malloc(n_records * sizeof (float32))) == NULL) {
432  printf("-E- %s: Error allocating memory to L3 product %s\n", __FILE__, prodname);
433  exit(1);
434  }
435 
436  /* allocate temp space for sum and sum of squares */
437  if ((sum_buf = (float32 *) calloc(2 * n_records, sizeof (float32))) == NULL) {
438  printf("-E- %s: Error allocating memory to L3 file products\n", __FILE__);
439  exit(1);
440  }
441 
442  /* allocate temp space for bin numbers, nobs, wts, etc. */
443  vsize = 2 * sizeof (int32) + sizeof (float32) + sizeof (int16);
444  if ((bin_buf = (uint8 *) calloc(n_records, vsize)) == NULL) {
445  printf("-E- %s: Error allocating memory for L3 file bin info\n", __FILE__);
446  exit(1);
447  }
448 
449  /* get bin info for each filled bin */
450  nread = VSread(vdata_id[0], (uint8 *) bin_buf, n_records, interlace0);
451  if (nread == -1) {
452  printf("-E- %s line %d: Unable to read bin numbers.\n", __FILE__, __LINE__);
453  exit(1);
454  }
455 
456 
457  /* get product sum and sum_sq for each filled bin */
458  /*
459  VSseek(vdata_id[2+prodinx], 0);
460  */
461  nread = VSread(vdata_id[2 + prodinx], (uint8 *) sum_buf, n_records, interlace);
462  if (nread == -1) {
463  printf("-E- %s line %d: Unable to read sum/sum_sq\n", __FILE__, __LINE__);
464  exit(1);
465  }
466 
467  for (i = 0; i < n_records; i++) {
468  memcpy((void *) &wgt, (const void *) &bin_buf[i * vsize + 6], 4);
469  data[prodID][i] = (float32) sum_buf[2 * i] / wgt;
470  /*
471  printf("%d %f %f %f\n",bins[i],sum_buf[2*i],wgt,data[prodID][i]);
472  */
473  }
474 
475  close_l3b();
476  free(bin_buf);
477  free(sum_buf);
478  closebin();
479 
480  loaded[prodID] = 1;
481  }
482 
483 
484  /* now, find the bin of interest for this product */
485 
486  bin_num = latlon2bin(lat, lon);
487  bin_idx = nearest_bin(bin_num);
488 
489  if (bin_idx >= 0) return (data[prodID][bin_idx]);
490  else return (-999.0);
491 }
492 
493 
494 
integer, parameter int16
Definition: cubeio.f90:3
int j
Definition: decode_rs.h:73
int32_t day
int status
Definition: l1_czcs_hdf.c:32
int32 latlon2bin(double lat, double lon)
float bin_climatology(char *l3file, int32_t day, float lon, float lat, char *prodname)
#define NULL
Definition: decode_rs.h:63
float * lat
int32 open_l3b(char *l3file, int32_t day)
int32 close_l3b()
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed file
Definition: HISTORY.txt:413
#define PI
Definition: l3_get_org.c:6
double constrain_lon(double lon)
int16 lat2row(double lat)
void initbin(void)
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE fields
Definition: HISTORY.txt:400
#define FATAL_ERROR
Definition: swl0_parms.h:5
int32_t l3file(int32_t sdfid, int32_t c_sdfid, int32_t *nbins, int32_t *c_nbins, char *ptype)
Definition: l3stat_chk.c:347
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
void closebin(void)
int32 nearest_bin(int32 bin_num)
int32 rowlon2bin(int16 row, double lon)
#define MAXNVDATA
float * lon
int i
Definition: decode_rs.h:71
#define MAXCPROD
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
double constrain_lat(double lat)