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
l1_czcs.c
Go to the documentation of this file.
1 /*
2  * W. Robinson, SAIC, 10 Dec 2004 I/O routines for CZCS
3  */
4 #include "l1.h"
5 #include <hdf4utils.h>
6 #include "mfhdf.h"
7 #include "l1_czcs.h"
8 #include <math.h>
9 #include <libnav.h>
10 
11 #define NREC_IN_BUF 1
12 #define NBND_CZCS 5
13 #define POS_ERR_THRESH 2000. /* orbit position error tolerence */
14 
15 int cz_posll_2_satang(float *, int, float *, float *, float *, float *);
16 void matrix_mult(double[3], double[3][3], double[3]);
17 void cross_prod(double *, double *, double *);
18 
19 static int syear, sday; /* data start date */
20 static int32 smsec; /* data start time */
21 static int16 eyear, eday; /* data end date */
22 static int32 emsec; /* data end time */
23 static int32 nscan; /* number of scans */
24 static int32 npix; /* number pixels per scan */
25 
26 static char dtype[8];
27 
28 static int32 nsta;
29 static int32 ninc;
30 
31 static uint8 *counts, cz_band_present;
32 static int32 *msec;
33 static int16 *gain;
34 static float32 *tilt, *att_ang, *slope, *intercept;
35 static float32 *ctl_pt_lat, *ctl_pt_lon, *pos, *pos_err;
36 static int32 nctl_pt, *ctl_pt_cols;
37 static float *ctl_pt_vx, *ctl_pt_vy, *ctl_pt_vz, *y2_vx, *y2_vy, *y2_vz, *ctl_pt_x;
38 static float *lt750; /* internal 750 mn data source */
39 static char *ring_sat; /* set to 1 if 443, 520 or 550 are saturated for ringing
40  mask computation */
41 
42 int czcs_ring(int gain, float *lt750, char *ring_sat, l1str *l1rec)
43 /*******************************************************************
44 
45  czcs_ring
46 
47  purpose: use the Mueller ringing flagging algorithm to set the stray
48  light mask for CZCS
49 
50  Returns type: int - 0 if OK
51 
52  Parameters: (in calling order)
53  Type Name I/O Description
54  ---- ---- --- -----------
55  int gain I gain for the line of data
56  float * lt750 I 750 nm total rads
57  char * ring_sat I 1 if bands 1,2,or 3 are
58  saturated = use the 750
59  excess rad in this case
60  l1str * l1rec I/O L1 struc including tot rads,
61  flags
62 
63  Modification history:
64  Programmer Date Description of change
65  ---------- ---- ---------------------
66  W. Robinson 31 May 2005 Original development
67  W. Robinson 31 Jul 2007 repaired problem that occurs in SGIs
68  if bsum = 0, avoid log and set ring
69  distance to 0
70 
71  Based on the algorithm of Mueller, J. L., Nimbus-7 CZCS: electronic
72  overshoot due to cloud reflectance, Appl. Optics, Vol 27, No 3, 1 Feb,
73  1988, pp 438 - 440
74  Modified to only include the excess radiance in 750 for pixels where
75  bands 1, 2, or 3 have saturated. This is to reduce the masking from
76  coastal areas where the vis is not as bright as for cloud
77 
78  *******************************************************************/ {
79  static float gtrans[4] = {1.0, 1.25, 1.5, 2.1};
80  /* note that the values below have +- of 9.7 and 3.3 in paper */
81  /* standard values followed by low, high confidence values
82  static float mueller_int = 3.9, mueller_slp = 30.8;
83  static float mueller_int = -6.2, mueller_slp = 27.5;
84  static float mueller_int = 13.6, mueller_slp = 34.1;
85  */
86  static float mueller_int = 3.9, mueller_slp = 30.8;
87 
88  static float l0_750 = 2.45;
89 
90  float gfac, lthresh, gfac_ln, excess_brit[1968], bsub, bsum;
91  int i, ipx, btarg, last_btarg, igrp, iavg, pxloc, n_ring;
92 
93  /*
94  * set up gfac - gain factor for line, lthresh, brightness threshold
95  * and flags for current and last pixel above threshold radiance
96  */
97  gfac = gtrans[ gain - 1];
98  lthresh = l0_750 / gfac;
99  gfac_ln = log(gfac);
100  btarg = 0;
101  last_btarg = 0;
102 
103  /*
104  * go through pixels and set the mask for pixels downstream of bright target
105  */
106  for (i = 0; i < npix; i++) {
107  if (*(lt750 + i) > lthresh) {
108  btarg = 1;
109  excess_brit[i] =
110  (*(ring_sat + i) == 1) ? *(lt750 + i) - lthresh : 0.;
111  l1rec->stlight[i] = 1;
112  } else
113  excess_brit[i] = 0.;
114 
115  if ((btarg == 0) && (last_btarg == 1)) {
116  /*
117  * sum up excess brightness in 5 groups of 10, weighted by f(dist)
118  */
119  bsum = 0.;
120  for (igrp = 0; igrp < 5; igrp++) {
121  bsub = 0.;
122  for (iavg = 0; iavg < 10; iavg++) {
123  pxloc = i - 1 - (iavg + igrp * 10);
124  if (pxloc >= 0)
125  bsub += excess_brit[ pxloc ];
126  }
127  bsub = bsub / 10.;
128  bsum += bsub * exp(-0.32 * (igrp + 1));
129  }
130  /*
131  * compute # pixels to mask and mark the pixels in that range
132  */
133  bsum = (bsum > 450.) ? 450. : bsum;
134  if (bsum > 0.) {
135  n_ring = mueller_int + mueller_slp * (log(bsum) + gfac_ln);
136  n_ring = (n_ring < 0) ? 0 : n_ring;
137  } else
138  n_ring = 0;
139 
140  for (ipx = 0; ipx < n_ring; ipx++) {
141  pxloc = ipx + i;
142  if (pxloc < npix)
143  l1rec->stlight[pxloc] = 1;
144  }
145  }
146  /*
147  * lastly, reset to go to next pixel
148  */
149  last_btarg = btarg;
150  btarg = 0;
151  }
152  return 0;
153 }
154 
155 #define NBND 4
156 #define NGAIN 4
157 #define NEPOCH 5
158 
159 int get_czcscal(char *file, int orbit, int16 year, int16 day, int32 msec, short l1acnt[], float slope750, float intercept750, int16 igain, float32 l1brads[])
160 /******************************************************************
161 
162  get_czcscal
163  purpose: replicate the Evans & Gordon calibration
164 
165  Modification history:
166  Programmer Date Description of change
167  ---------- ---- ---------------------
168  Sean Bailey 08 Feb 2006 Original development
169 
170  Based on the original czcscaleg.f
171  Modified to use cal_czcs.hdf
172 
173  *******************************************************************/ {
174  static int firstCall = 1;
175  static float slope_coeff[NGAIN][NBND_CZCS];
176  static float offsets[NGAIN][NBND_CZCS];
177  static float rk[NGAIN][NBND_CZCS];
178  static float tdecay[NEPOCH][NBND_CZCS];
179  static int orbitepoch[NEPOCH];
180  static float timeconst[2];
181  static float time_dep[NBND_CZCS][NGAIN][3];
182 
183  int epochidx = 0;
184  int i, j, iorbit;
185  float rden, rnum;
186  float tcorr, slp, S;
187  static int orbepoch670 = 6750;
188  float g[NBND];
189  float64 jsec = yds2unix(year, day, ((double) (msec)) / 1000.0);
190  float64 refjsec = yds2unix(1978, 278, (double) 0.0);
191  float64 decayday = (jsec - refjsec) / 86400.;
192 
193  igain -= 1;
194 
195  if (firstCall) {
196  char name [H4_MAX_NC_NAME] = "";
197  char sdsname[H4_MAX_NC_NAME] = "";
198  int32 sd_id;
199  int32 sds_id;
200  int32 rank;
201  int32 nt;
202  int32 dims[H4_MAX_VAR_DIMS];
203  int32 nattrs;
204  int32 start[3] = {0, 0, 0};
205  int32 status;
206  /* get the current calibration */
207  if (file == NULL) {
208  fprintf(stderr,
209  "-E %s Line %d: No calibration file specified.\n",
210  __FILE__, __LINE__);
211  exit(1);
212  }
213 
214  if (want_verbose)
215  printf("Loading caltable: %s\n", file);
216  /* Open the file */
217  sd_id = SDstart(file, DFACC_RDONLY);
218  if (sd_id == FAIL) {
219  fprintf(stderr, "-E- %s line %d: SDstart(%s, %d) failed.\n",
220  __FILE__, __LINE__, file, DFACC_RDONLY);
221  exit(1);
222  }
223 
224  strcpy(sdsname, "slope_coeff");
225  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
226  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
227  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) slope_coeff);
228  if (status != 0) {
229  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
230  __FILE__, __LINE__, sdsname, file);
231  exit(1);
232  } else {
233  status = SDendaccess(sds_id);
234  }
235 
236  strcpy(sdsname, "offsets");
237  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
238  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
239  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) offsets);
240  if (status != 0) {
241  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
242  __FILE__, __LINE__, sdsname, file);
243  exit(1);
244  } else {
245  status = SDendaccess(sds_id);
246  }
247 
248  strcpy(sdsname, "rk");
249  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
250  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
251  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) rk);
252  if (status != 0) {
253  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
254  __FILE__, __LINE__, sdsname, file);
255  exit(1);
256  } else {
257  status = SDendaccess(sds_id);
258  }
259 
260  strcpy(sdsname, "dec");
261  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
262  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
263  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) tdecay);
264  if (status != 0) {
265  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
266  __FILE__, __LINE__, sdsname, file);
267  exit(1);
268  } else {
269  status = SDendaccess(sds_id);
270  }
271 
272  strcpy(sdsname, "iorbdec");
273  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
274  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
275  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) orbitepoch);
276  if (status != 0) {
277  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
278  __FILE__, __LINE__, sdsname, file);
279  exit(1);
280  } else {
281  status = SDendaccess(sds_id);
282  }
283 
284  strcpy(sdsname, "timeconst");
285  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
286  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
287  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) timeconst);
288  if (status != 0) {
289  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
290  __FILE__, __LINE__, sdsname, file);
291  exit(1);
292  } else {
293  status = SDendaccess(sds_id);
294  }
295 
296  strcpy(sdsname, "time_dep");
297  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
298  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
299  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) time_dep);
300  if (status != 0) {
301  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
302  __FILE__, __LINE__, sdsname, file);
303  exit(1);
304  } else {
305  status = SDendaccess(sds_id);
306  }
307  /* terminate access to the SD interface and close the file */
308  status = SDend(sd_id);
309 
310  firstCall = 0;
311  }
312 
313  iorbit = MIN(orbit, 39000);
314  for (j = 0; j < NEPOCH - 1; j++) {
315  if ((iorbit >= orbitepoch[j]) && (iorbit <= orbitepoch[j + 1]))
316  epochidx = j;
317  }
318 
319  rden = (float) orbitepoch[epochidx + 1] - (float) orbitepoch[epochidx];
320 
321  for (i = 0; i < NBND; i++) {
322  /* Calculate time dependent correction... */
323  /* double exp
324  tcorr = time_dep[i][igain][0] - time_dep[i][igain][1] * (1.0 - exp(timeconst[0]*decayday)) - time_dep[i][igain][2] *(1.0 - exp(timeconst[1]*decayday));
325  */
326  /* single exp */
327  tcorr = time_dep[i][igain][0] - time_dep[i][igain][1] * (1.0 - exp(-1 * time_dep[i][igain][2] * decayday));
328 
329  rnum = tdecay[epochidx + 1][i] - tdecay[epochidx][i];
330  slp = rnum / rden;
331 
332  /* Antione modification for 670nm */
333  if (i == 3 && iorbit > orbepoch670) {
334  rden = (float) orbitepoch[4] - (float) orbepoch670;
335  rnum = tdecay[4][i] - tdecay[3][i];
336  slp = rnum / rden;
337  S = (float) iorbit - (float) orbepoch670;
338  g[i] = 1.0 / (S * slp + tdecay[3][i]);
339 
340  } else {
341  S = (float) iorbit - (float) orbitepoch[epochidx];
342  g[i] = 1.0 / (S * slp + tdecay[epochidx][i]);
343  }
344  /*fprintf(stderr,"i: %d g: %8.6f tcorr: %8.6f\n",i,g[i],tcorr);*/
345 
346  l1brads[i] = (tcorr * g[i] * rk[igain][i] * slope_coeff[igain][i] * (float) l1acnt[i]) + offsets[igain][i];
347 
348  }
349 
350  l1brads[4] = slope750 * (float) l1acnt[4] + intercept750;
351 
352  return 0;
353 }
354 
355 /*
356  W. Robinson, SAIC, 6 Jan 2006 add code to read the position and
357  position error SDSes
358  */
359 
360 int openl1_czcs(filehandle *file) {
361 
362  /* */
363  /* get_l1a_open interface */
364  /* */
365  int32 fileID;
366  int32 status;
367  int16 sy, sd;
368  int i;
369 
370  /* open the file and get the file ID */
371  fileID = SDstart(file->name, DFACC_RDONLY);
372 
373  if (fileID < 0) {
374  fprintf(stderr,
375  "-E %s Line %d: Error opening %s for reading.\n",
376  __FILE__, __LINE__, file->name);
377  return (1);
378  }
379  status = getHDFattr(fileID, "Start Year", "", (VOIDP) & sy);
380  syear = (int32) sy;
381  status = getHDFattr(fileID, "Start Day", "", (VOIDP) & sd);
382  sday = (int32) sd;
383  status = getHDFattr(fileID, "Start Millisec", "", (VOIDP) & smsec);
384  status = getHDFattr(fileID, "End Year", "", (VOIDP) & eyear);
385  status = getHDFattr(fileID, "End Day", "", (VOIDP) & eday);
386  status = getHDFattr(fileID, "End Millisec", "", (VOIDP) & emsec);
387  status = getHDFattr(fileID, "Number of Scan Lines", "", (VOIDP) & nscan);
388  status = getHDFattr(fileID, "Data Type", "", (VOIDP) & dtype);
389 
390  status = getHDFattr(fileID, "Pixels per Scan Line", "", (VOIDP) & npix);
391  status = getHDFattr(fileID, "LAC Pixel Start Number", "", (VOIDP) & nsta);
392  status = getHDFattr(fileID, "LAC Pixel Subsampling", "", (VOIDP) & ninc);
393 
394  status = getHDFattr(fileID, "Orbit Number", "", (VOIDP) & file->orbit_number);
395  /* WDR not in czcs l1 file, so comment for now
396  status = getHDFattr(fileID, "Orbit Node Longitude", "",
397  (VOIDP)&file->orbit_node_lon);
398  status = getHDFattr(fileID, "Node Crossing Time", "",
399  (VOIDP)&file->node_crossing_time[0]);
400  */
401  status = getHDFattr(fileID, "Number of Pixel Control Points", "",
402  (VOIDP) & nctl_pt);
403  status = getHDFattr(fileID, "Parameter Presence Code", "",
404  (VOIDP) & cz_band_present);
405 
406  /* call cdata.f to initialize global FORTRAN common block data */
407  cdata_();
408 
409 
410  file->npix = npix;
411  file->nscan = nscan;
412  file->sensorID = CZCS;
413  file->sd_id = fileID;
414  strcpy(file->spatialResolution, "825 m");
415 
416  msec = (int32 *) calloc(nscan, sizeof (int32));
417  status = rdSDS(file->sd_id, "msec", 0, 0, 0, 0, (VOIDP) msec);
418 
419  att_ang = (float32 *) calloc(nscan * 3, sizeof ( float32));
420  status = rdSDS(file->sd_id, "att_ang", 0, 0, 0, 0, (VOIDP) att_ang);
421 
422  ctl_pt_cols = (int32 *) calloc(nctl_pt, sizeof ( int32));
423  status = rdSDS(file->sd_id, "cntl_pt_cols", 0, 0, 0, 0,
424  (VOIDP) ctl_pt_cols);
425  ctl_pt_x = (float *) calloc(nctl_pt, sizeof ( float));
426  for (i = 0; i < nctl_pt; i++)
427  ctl_pt_x[i] = (float) ctl_pt_cols[i] - 1.;
428 
429  tilt = (float32 *) calloc(nscan, sizeof (float32));
430  status = rdSDS(file->sd_id, "tilt", 0, 0, 0, 0, (VOIDP) tilt);
431 
432  slope = (float32 *) calloc(nscan * 6, sizeof (float32));
433  status = rdSDS(file->sd_id, "slope", 0, 0, 0, 0, (VOIDP) slope);
434  intercept = (float32 *) calloc(nscan * 6, sizeof (float32));
435  status = rdSDS(file->sd_id, "intercept", 0, 0, 0, 0, (VOIDP) intercept);
436  /* get nav orbit data */
437  pos = (float32 *) malloc(nscan * 3 * sizeof (float32));
438  status = rdSDS(file->sd_id, "orb_vec", 0, 0, 0, 0, (VOIDP) pos);
439  pos_err = (float32 *) malloc(nscan * sizeof (float32));
440  status = rdSDS(file->sd_id, "pos_err", 0, 0, 0, 0, (VOIDP) pos_err);
441 
442  gain = (int16 *) malloc(nscan * sizeof ( int16));
443  counts = (uint8 *) malloc(npix * NBND_CZCS * sizeof ( uint8));
444  ctl_pt_lat = (float32 *) malloc(nctl_pt * sizeof ( float32));
445  ctl_pt_lon = (float32 *) malloc(nctl_pt * sizeof ( float32));
446  ctl_pt_vx = (float32 *) malloc(nctl_pt * sizeof ( float32));
447  ctl_pt_vy = (float32 *) malloc(nctl_pt * sizeof ( float32));
448  ctl_pt_vz = (float32 *) malloc(nctl_pt * sizeof ( float32));
449  y2_vx = (float32 *) malloc(npix * sizeof ( float32));
450  y2_vy = (float32 *) malloc(npix * sizeof ( float32));
451  y2_vz = (float32 *) malloc(npix * sizeof ( float32));
452  lt750 = (float *) malloc(npix * sizeof ( float));
453  ring_sat = (char *) malloc(npix * sizeof ( char));
454 
455  return (status);
456 }
457 
458 /*
459  W. Robinson 31 May 2005 add ringing masking call
460  W. Robinson, SAIC, 6 Jan 2005 add code to compute sat angles with
461  cz_posll_2_satang if pos_err is acceptable
462  J. Gales 23 Sep 2005 add kludge for subsetted pixel range in
463  satang
464  (note that satang may not give proper values in pixel subsets with this
465  but it may be academic with...)
466  W. Robinson, SAIC, 10 Sep 2010 use vectors instead of lat, lon when
467  doing interpolation of ctl pt lat, lon to
468  all samples
469  */
470 
471 int readl1_czcs(filehandle *file, int32_t recnum, l1str *l1rec) {
472  /*void czcscal_( int *, float[], float[], short *, int *, float * );*/
473  void satang_(double *, double *, float *, float *, float *, float *,
474  float *, float *, float *, float *);
475  void sunangs_(int *, int *, float *, float *, float *, float *, float *);
476  int ll2vec(float *, float *);
477  int vec2ll(float *, float *);
478  short cnt_vec[NBND_CZCS];
479  float lt_lcl[NBND_CZCS], yout1, yout2, yout3, llvec[2], vec[3], gmt;
480  int ipx, ibnd, orbit, i, tpix;
481  uint8 cal_sum[5], cal_scan[6];
482  int32 status, navbad, ltsat;
483  double pi, radeg;
484  int32_t ib;
485 
486  int32_t nwave = l1rec->l1file->nbands;
487  int32_t *bindx = l1rec->l1file->bindx;
488  char *cal_path = l1_input->calfile;
489  /* load scan times */
490  /* int year = *l1rec->year;
491  int day = *l1rec->day;*/
492  /* int32 msec = l1rec->msec[recnum];*/
493 
494  float lonbuf[1968], latbuf[1968], senzbuf[1968], senabuf[1968];
495 
496  /*
497  * read in the line of counts, latitude, longitude
498  */
499  status = rdSDS(file->sd_id, "band1", recnum, 0, 1, npix, (VOIDP) counts);
500  status = rdSDS(file->sd_id, "band2", recnum, 0, 1, npix,
501  (VOIDP) (counts + npix));
502  status = rdSDS(file->sd_id, "band3", recnum, 0, 1, npix,
503  (VOIDP) (counts + 2 * npix));
504  status = rdSDS(file->sd_id, "band4", recnum, 0, 1, npix,
505  (VOIDP) (counts + 3 * npix));
506  status = rdSDS(file->sd_id, "band5", recnum, 0, 1, npix,
507  (VOIDP) (counts + 4 * npix));
508 
509  status = rdSDS(file->sd_id, "gain", recnum, 0, 1, 1, (VOIDP) gain);
510  status = rdSDS(file->sd_id, "latitude", recnum, 0, 1, nctl_pt,
511  (VOIDP) ctl_pt_lat);
512  status = rdSDS(file->sd_id, "longitude", recnum, 0, 1, nctl_pt,
513  (VOIDP) ctl_pt_lon);
514 
515  status = rdSDS(file->sd_id, "cal_sum", recnum, 0, 1, 5, (VOIDP) cal_sum);
516  status = rdSDS(file->sd_id, "cal_scan", recnum, 0, 1, 6, (VOIDP) cal_scan);
517  /*
518  * flag setting for entire line: set navfail if cal_sum shows problems
519  * (suspect it is * bad ephemeris or atitude) and set hilt if bands
520  * missing
521  */
522  ltsat = 0;
523  navbad = 0;
524  if ((cz_band_present & 0XF8) != 0XF8)
525  ltsat = 1;
526  else {
527  if ((cal_sum[3] != 0) || (cal_sum[4] != 0))
528  navbad = 1;
529  if ((cal_scan[0] != 0) || (cal_scan[1] != 0) ||
530  (cal_scan[2] != 0) || (cal_scan[3] != 0) ||
531  (cal_scan[4] != 0)) {
532  ltsat = 1;
533  navbad = 1;
534  }
535  }
536  /*
537  * calibrate the radiances and set flags for pixels
538  */
539  for (ipx = 0; ipx < npix; ipx++) {
540  for (ibnd = 0; ibnd < NBND_CZCS; ibnd++) {
541  cnt_vec[ibnd] = *(counts + ipx + ibnd * npix);
542  if ((cnt_vec[ibnd] == 255) || (ltsat == 1)) l1rec->hilt[ipx] = 1;
543  }
544  *(ring_sat + ipx) = ((cnt_vec[0] == 255) || (cnt_vec[1] == 255)
545  || (cnt_vec[2] == 255)) ? 1 : 0;
546  /*
547  * call the fortran calibration routine
548  */
549  orbit = (int) file->orbit_number;
550  /*czcscal_( &orbit, slope + recnum * 6,
551  intercept + recnum * 6, cnt_vec, &lgain, lt_lcl );*/
552  status = get_czcscal(cal_path, orbit, syear, sday, msec[recnum],
553  cnt_vec, slope[4], intercept[4], gain[0], lt_lcl);
554  /*
555  * assign the calibrated radiances
556  */
557  for (ibnd = 0; ibnd < nwave; ibnd++) {
558  ib = bindx[ ibnd ];
559  l1rec->Lt[ ipx * nwave + ib ] = lt_lcl[ibnd];
560  }
561  /*
562  * save the 750 nm data here
563  */
564  *(lt750 + ipx) = lt_lcl[4];
565 
566  if (navbad == 1) l1rec->navfail[ipx] = 1;
567  }
568  /*
569  * set the ringing mask
570  */
571  czcs_ring(gain[0], lt750, ring_sat, l1rec);
572  /*
573  * get navigation at all points from control point lat and lon values
574  * use spline interpolation to fill in, do in vector space to avoid date
575  * line discontinuity
576  */
577  for (i = 0; i < nctl_pt; i++) {
578  llvec[0] = ctl_pt_lat[i];
579  llvec[1] = ctl_pt_lon[i];
580  if (ll2vec(llvec, vec) == 0) {
581  ctl_pt_vx[i] = *vec;
582  ctl_pt_vy[i] = *(vec + 1);
583  ctl_pt_vz[i] = *(vec + 2);
584  } else {
585  fprintf(stderr, "-E %s Line %d: ll2vec failure.\n",
586  __FILE__, __LINE__);
587  exit(1);
588  }
589  }
590  spline(ctl_pt_x, ctl_pt_vx, nctl_pt, 1e30, 1e30, y2_vx);
591  spline(ctl_pt_x, ctl_pt_vy, nctl_pt, 1e30, 1e30, y2_vy);
592  spline(ctl_pt_x, ctl_pt_vz, nctl_pt, 1e30, 1e30, y2_vz);
593  for (i = 0; i < npix; i++) {
594  tpix = i * ninc /*+ spix*/;
595  splint(ctl_pt_x, ctl_pt_vx, y2_vx, nctl_pt, tpix, &yout1);
596  splint(ctl_pt_x, ctl_pt_vy, y2_vy, nctl_pt, tpix, &yout2);
597  splint(ctl_pt_x, ctl_pt_vz, y2_vz, nctl_pt, tpix, &yout3);
598 
599  *vec = yout1;
600  *(vec + 1) = yout2;
601  *(vec + 2) = yout3;
602  vec2ll(vec, llvec);
603 
604  l1rec->lon[i] = llvec[1];
605  l1rec->lat[i] = llvec[0];
606  }
607  /*
608  * calculate geometry. For sensor angles, use the orbit info if the
609  * position error is available and within error threshold
610  */
611  if ((*(pos_err + recnum) >= 0.) &&
612  (*(pos_err + recnum) < POS_ERR_THRESH)) {
613  cz_posll_2_satang((pos + 3 * recnum), npix, l1rec->lat, l1rec->lon,
614  l1rec->senz, l1rec->sena);
615  } else {
616  pi = PI;
617  radeg = RADEG;
618  for (i = 0; i < 1968; i++) {
619  lonbuf[i] = 0.0;
620  latbuf[i] = 0.0;
621  senzbuf[i] = 0.0;
622  senabuf[i] = 0.0;
623  }
624 
625  for (i = 0; i < npix; i++) {
626  lonbuf[i] = l1rec->lon[i];
627  latbuf[i] = l1rec->lat[i];
628  }
629 
630  satang_(&pi, &radeg, tilt + recnum, att_ang + 3 * recnum + 1,
631  att_ang + 3 * recnum + 2, att_ang + 3 * recnum, lonbuf,
632  latbuf, senzbuf, senabuf);
633 
634  for (i = 0; i < npix; i++) {
635  l1rec->senz[i] = senzbuf[i];
636  l1rec->sena[i] = senabuf[i];
637  }
638  }
639 
640 
641  /* for( i = spix; i < epix; i = i + ninc )*/
642  for (i = 0; i < npix; i++) {
643  gmt = (float) msec[recnum] / (1000. * 3600);
644  sunangs_(&syear, &sday, &gmt, (l1rec->lon) + i, (l1rec->lat) + i,
645  (l1rec->solz) + i, (l1rec->sola) + i);
646  }
647  /*
648  * set scan time
649  */
650  double secs = (double) (msec[recnum] / 1000.0);
651  int16_t year = syear;
652  int16_t day = sday;
653 
654  if (recnum > 0 && msec[recnum] < msec[recnum - 1]) { /* adjust for day rollover */
655  day += 1;
656  if (day > (365 + (year % 4 == 0))) {
657  year += 1;
658  day = 1;
659  }
660  }
661  l1rec->scantime = yds2unix(year, day, secs);
662 
663  return (status);
664 }
665 
666 /* W. Robinson, SAIC, 6 Jan 2005 include the pos, pos_err arrays in free */
667 
668 int closel1_czcs(filehandle *file) {
669  free(msec);
670  free(tilt);
671  free(att_ang);
672  free(counts);
673  free(ctl_pt_lat);
674  free(ctl_pt_lon);
675  free(ctl_pt_vx);
676  free(ctl_pt_vy);
677  free(ctl_pt_vz);
678  free(y2_vx);
679  free(y2_vy);
680  free(y2_vz);
681  free(ctl_pt_x);
682  free(ctl_pt_cols);
683  free(slope);
684  free(intercept);
685  free(lt750);
686  free(ring_sat);
687  free(pos);
688  free(pos_err);
689 
690  SDend(file->sd_id);
691 
692  return (0);
693 }
694 
695 #define re 6378.137
696 #define f 1. / 298.257
697 #define omf2 ( 1 - f ) * ( 1 - f )
698 
699 int cz_posll_2_satang(float *pos, int npix, float *lat, float *lon,
700  float *senz, float *sena)
701 /*******************************************************************
702 
703  cz_posll_2_satang
704 
705  purpose: convert satellite position, lat, lon into satellite zenith
706  and azimuth for a line of czcs data
707 
708  Returns type: int - 0 if no problems
709 
710  Parameters: (in calling order)
711  Type Name I/O Description
712  ---- ---- --- -----------
713  float * pos I position (x, y, z) of
714  satellite in meters
715  int npix I # pixels of lat, lon to find
716  sat angles for
717  float * lat I latitude array in degrees E
718  float * lon I longitude array in degrees E
719  float * senz O sensor zenith in degrees
720  float * sena O sensor azimuth in degrees
721 
722  Modification history:
723  Programmer Date Description of change
724  ---------- ---- ---------------------
725  W. Robinson, SAIC 5-Jan-2006 Original development
726 
727  *******************************************************************/ {
728  double up[3], ea[3], no[3], radeg, upxy, xlmat[3][3], xlatg, gv[3];
729  double r, rh[3], rl[3];
730  int i, j;
731 
732  radeg = 90. / asin(1.);
733 
734  for (i = 0; i < npix; i++) {
735  /*
736  * Compute the local vertical, East and North unit vectors
737  */
738  up[0] = cos(*(lat + i) / radeg) * cos(*(lon + i) / radeg);
739  up[1] = cos(*(lat + i) / radeg) * sin(*(lon + i) / radeg);
740  up[2] = sin(*(lat + i) / radeg);
741 
742  upxy = sqrt(up[0] * up[0] + up[1] * up[1]);
743 
744  ea[0] = -up[1] / upxy;
745  ea[1] = up[0] / upxy;
746  ea[2] = 0.;
747 
748  cross_prod(up, ea, no);
749 
750  /*
751  * Store vectors in 3x3 array
752  */
753  for (j = 0; j < 3; j++) {
754  xlmat[0][j] = ea[j];
755  xlmat[1][j] = no[j];
756  xlmat[2][j] = up[j];
757  }
758  /*
759  * Compute geocentric position vector
760  */
761  xlatg = radeg * atan(tan(*(lat + i) / radeg) * omf2);
762  gv[0] = cos(xlatg / radeg) * cos(*(lon + i) / radeg);
763  gv[1] = cos(xlatg / radeg) * sin(*(lon + i) / radeg);
764  gv[2] = sin(xlatg / radeg);
765 
766  r = re * (1. - f) /
767  sqrt(1. - (2. - f) * f * pow(cos(xlatg / radeg), 2));
768 
769  /* note that pos is in m but rest of constants are in km, hence the
770  1000 factor */
771  for (j = 0; j < 3; j++)
772  rh[j] = pos[j] / 1000. - r * gv[j];
773  /*
774  * Transform the pixel-to-spacecraft vector into the local frame
775  */
776  matrix_mult(rh, xlmat, rl);
777  /*
778  * Compute the sensor zenith and azimuth
779  * if zenith close to 0, set azimuth to 0 (and normalize azimuth)
780  */
781  *(senz + i) = radeg * atan2(sqrt(rl[0] * rl[0] + rl[1] * rl[1]),
782  rl[2]);
783  if (*(senz + i) > 0.05)
784  *(sena + i) = radeg * atan2(rl[0], rl[1]);
785  else
786  *(sena + i) = 0.;
787 
788  if (*(sena + i) < 0.) *(sena + i) += 360.;
789  }
790  /*
791  * and end
792  */
793  return 0;
794 }
795 
796 void matrix_mult(double vecin[3], double matrix[3][3], double vecout[3])
797 /*******************************************************************
798 
799  matrix_mult
800 
801  purpose: multiply a vector by a matrix and produce the output vector
802 
803  Returns type: void
804 
805  Parameters: (in calling order)
806  Type Name I/O Description
807  ---- ---- --- -----------
808  double[3] vecin I input vector
809  double[3][3] matrix I rotation matrix
810  double[3] vecout O rotated vector
811 
812  Modification history:
813  Programmer Date Description of change
814  ---------- ---- ---------------------
815  W. Robinson, SAIC 2-Nov-2005 Original development
816 
817  *******************************************************************/ {
818  int i, j;
819  for (i = 0; i < 3; i++) {
820  vecout[i] = 0.;
821  for (j = 0; j < 3; j++) {
822  vecout[i] += matrix[i][j] * vecin[j];
823  }
824  }
825 }
826 
827 void cross_prod(double *v1, double *v2, double *vout)
828 /*******************************************************************
829 
830  cross_prod
831 
832  purpose: produce the cross product of 2 vectors
833 
834  Returns type: void
835 
836  Parameters: (in calling order)
837  Type Name I/O Description
838  ---- ---- --- -----------
839  double * v1 I vector 1
840  double * v2 I vector 1
841  double * vout O v1 X v2
842 
843  Modification history:
844  Programmer Date Description of change
845  ---------- ---- ---------------------
846  W. Robinson, SAIC 4-Jan-2006 Original development
847 
848  *******************************************************************/ {
849  *vout = v1[1] * v2[2] - v1[2] * v2[1];
850  *(vout + 1) = v1[2] * v2[0] - v1[0] * v2[2];
851  *(vout + 2) = v1[0] * v2[1] - v1[1] * v2[0];
852 }
NGAIN
#define NGAIN
Definition: l1_czcs.c:156
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
bindx
const int bindx[3]
Definition: DbLutNetcdf.cpp:28
MIN
#define MIN(x, y)
Definition: rice.h:169
r
int r
Definition: decode_rs.h:73
ll2vec
int ll2vec(float *ll, float *vec)
Definition: ll2vec.c:6
g
e g
Definition: HOWTO_Add_a_product.txt:5
j
int j
Definition: decode_rs.h:73
rdSDS
int rdSDS(int32_t fileID, const char sdsname[], int32_t start1, int32_t start2, int32_t edges1, int32_t edges2, void *array_data)
day
int32_t day
Definition: atrem_corl1v2.h:308
status
int status
Definition: l1_czcs_hdf.c:32
yds2unix
double yds2unix(int16_t year, int16_t day, double secs)
Definition: yds2unix.c:7
orbit
*********************************************HISTORY FILE for MOD_PR02AQUA **Version ****Point of no algorithm change is made in this so the and when the blackbody temperature is above a threshold Since the LWIR FPA temperature oscillates on orbit
Definition: HISTORY.txt:106
fileID
int16_t fileID
Definition: main_l1det2det.c:44
FAIL
#define FAIL
Definition: ObpgReadGrid.h:18
lt750
float * lt750
Definition: l1_czcs_hdf.c:38
NULL
#define NULL
Definition: decode_rs.h:63
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
POS_ERR_THRESH
#define POS_ERR_THRESH
Definition: l1_czcs.c:13
f
#define f
Definition: l1_czcs.c:696
pi
const double pi
Definition: VcstCmnConsts.h:473
ctl_pt_lat
float32 * ctl_pt_lat
Definition: l1_czcs_hdf.c:35
ctl_pt_cols
int32 * ctl_pt_cols
Definition: l1_czcs_hdf.c:36
spline
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
ctl_pt_vy
float * ctl_pt_vy
Definition: l1_czcs_hdf.c:37
pos
float32 * pos
Definition: l1_czcs_hdf.c:35
matrix
float ** matrix(long nrl, long nrh, long ncl, long nch)
Definition: nrutil.c:60
cz_band_present
uint8 cz_band_present
Definition: l1_czcs_hdf.c:30
ctl_pt_vz
float * ctl_pt_vz
Definition: l1_czcs_hdf.c:37
get_czcscal
int get_czcscal(char *file, int orbit, int16 year, int16 day, int32 msec, short l1acnt[], float slope750, float intercept750, int16 igain, float32 l1brads[])
Definition: l1_czcs.c:159
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
file
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
NBND
#define NBND
Definition: l1_czcs.c:155
PI
#define PI
Definition: l3_get_org.c:6
run_local.lat
lat
Definition: run_local.py:69
sunangs_
void sunangs_(int *, int *, float *, float *, float *, float *, float *)
Definition: sunangs.c:25
y2_vz
float * y2_vz
Definition: l1_czcs_hdf.c:37
recnum
read recnum
Definition: HOWTO_Add_a_sensor.txt:72
y2_vx
float * y2_vx
Definition: l1_czcs_hdf.c:37
run_local.lon
lon
Definition: run_local.py:69
cz_posll_2_satang
int cz_posll_2_satang(float *, int, float *, float *, float *, float *)
Definition: l1_czcs.c:699
readl1_czcs
int readl1_czcs(filehandle *file, int32_t recnum, l1str *l1rec)
Definition: l1_czcs.c:471
hico.proc_hico.sena
sena
Definition: proc_hico.py:241
satang_
void satang_(double *, double *, float *, float *, float *, float *, float *, float *, float *, float *)
slope
float32 slope[]
Definition: l2lists.h:30
cdata_
void cdata_()
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
getHDFattr
int getHDFattr(int32_t fileID, const char attrname[], const char sdsname[], void *data)
want_verbose
int want_verbose
Definition: genutils_globals.c:3
ring_sat
char * ring_sat
Definition: l1_czcs_hdf.c:39
intercept
float32 intercept[]
Definition: l2lists.h:44
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
nctl_pt
int32 nctl_pt
Definition: l1_czcs_hdf.c:36
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
matrix_mult
void matrix_mult(double[3], double[3][3], double[3])
Definition: l1_czcs.c:796
vec2ll
int vec2ll(float *vec, float *ll)
Definition: ll2vec.c:53
ctl_pt_x
float * ctl_pt_x
Definition: l1_czcs_hdf.c:37
ctl_pt_lon
float32 * ctl_pt_lon
Definition: l1_czcs_hdf.c:35
re
#define re
Definition: l1_czcs.c:695
v2
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame the required RAM for each execution is MB on the DEC ALPHA and MB on the SGI Octane v2
Definition: HISTORY.txt:728
dtype
dtype
Definition: DDataset.hpp:31
y2_vy
float * y2_vy
Definition: l1_czcs_hdf.c:37
rank
Extra metadata that will be written to the HDF4 file l2prod rank
Definition: HOWTO_Add_a_product.txt:80
CZCS
#define CZCS
Definition: sensorDefs.h:17
cross_prod
void cross_prod(double *, double *, double *)
Definition: l1_czcs.c:827
splint
subroutine splint(xa, ya, y2a, n, x, y)
Definition: atrem_app_refl_f90_cubeio.f:3270
NBND_CZCS
#define NBND_CZCS
Definition: l1_czcs.c:12
openl1_czcs
int openl1_czcs(filehandle *file)
Definition: l1_czcs.c:360
pos_err
float32 * pos_err
Definition: l1_czcs_hdf.c:35
att_ang
float32 * att_ang
Definition: l1_czcs_hdf.c:34
omf2
#define omf2
Definition: l1_czcs.c:697
NEPOCH
#define NEPOCH
Definition: l1_czcs.c:157
l1_czcs.h
tilt
int16_t * tilt
Definition: l2bin.cpp:80
i
int i
Definition: decode_rs.h:71
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
czcs_ring
int czcs_ring(int gain, float *lt750, char *ring_sat, l1str *l1rec)
Definition: l1_czcs.c:42
npix
int npix
Definition: get_cmp.c:28
ctl_pt_vx
float * ctl_pt_vx
Definition: l1_czcs_hdf.c:37
closel1_czcs
int closel1_czcs(filehandle *file)
Definition: l1_czcs.c:668


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