ocssw V2020
l1_czcs_hdf.c
Go to the documentation of this file.
1 /*
2  * W. Robinson, SAIC, 10 Dec 2004 I/O routines for CZCS
3  */
4 #include "hdf.h"
5 #include "mfhdf.h"
6 #include "l12_proto.h"
7 #include "l1_czcs_hdf.h"
8 #include <math.h>
9 
10 #define NREC_IN_BUF 1
11 #define NBND_CZCS 5
12 #define POS_ERR_THRESH 2000. /* orbit position error tolerence */
13 
14 int syear, sday; /* data start date */
15 int32 smsec; /* data start time */
16 int16 eyear, eday; /* data end date */
17 int32 emsec; /* data end time */
18 int32 nscan; /* number of scans */
19 int32 npix; /* number pixels per scan */
20 int32 spix; /* start pixel number (from 0) */
21 int32 dpix; /* LAC pixel increment */
22 int32 epix;
23 
24 char dtype[8];
25 
26 int32 nsta;
27 int32 ninc;
28 
30 int32 *msec;
33 float32 *tilt, *att_ang, *slope, *intercept;
37 float *lt750; /* internal 750 mn data source */
38 char *ring_sat; /* set to 1 if 443, 520 or 550 are saturated for ringing
39  mask computation */
40 
41 int czcs_ring(int gain, float *lt750, char *ring_sat, l1str *l1rec)
42 /*******************************************************************
43 
44  czcs_ring
45 
46  purpose: use the Mueller ringing flagging algorithm to set the stray
47  light mask for CZCS
48 
49  Returns type: int - 0 if OK
50 
51  Parameters: (in calling order)
52  Type Name I/O Description
53  ---- ---- --- -----------
54  int gain I gain for the line of data
55  float * lt750 I 750 nm total rads
56  char * ring_sat I 1 if bands 1,2,or 3 are
57  saturated = use the 750
58  excess rad in this case
59  l1str * l1rec I/O L1 struc including tot rads,
60  flags
61 
62  Modification history:
63  Programmer Date Description of change
64  ---------- ---- ---------------------
65  W. Robinson 31 May 2005 Original development
66  W. Robinson 31 Jul 2007 repaired problem that occurs in SGIs
67  if bsum = 0, avoid log and set ring
68  distance to 0
69 
70  Based on the algorithm of Mueller, J. L., Nimbus-7 CZCS: electronic
71  overshoot due to cloud reflectance, Appl. Optics, Vol 27, No 3, 1 Feb,
72  1988, pp 438 - 440
73  Modified to only include the excess radiance in 750 for pixels where
74  bands 1, 2, or 3 have saturated. This is to reduce the masking from
75  coastal areas where the vis is not as bright as for cloud
76 
77  *******************************************************************/ {
78  static float gtrans[4] = {1.0, 1.25, 1.5, 2.1};
79  /* note that the values below have +- of 9.7 and 3.3 in paper */
80  /* standard values followed by low, high confidence values
81  static float mueller_int = 3.9, mueller_slp = 30.8;
82  static float mueller_int = -6.2, mueller_slp = 27.5;
83  static float mueller_int = 13.6, mueller_slp = 34.1;
84  */
85  static float mueller_int = 3.9, mueller_slp = 30.8;
86 
87  static float l0_750 = 2.45;
88 
89  float gfac, lthresh, gfac_ln, excess_brit[1968], bsub, bsum;
90  int i, ipx, btarg, last_btarg, igrp, iavg, pxloc, n_ring;
91  extern int32 npix;
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  spix = nsta - 1;
456  dpix = ninc;
457  epix = spix + npix * dpix;
458 
459 
460  return (status);
461 }
462 
463 /*
464  W. Robinson 31 May 2005 add ringing masking call
465  W. Robinson, SAIC, 6 Jan 2005 add code to compute sat angles with
466  cz_posll_2_satang if pos_err is acceptable
467  J. Gales 23 Sep 2005 add kludge for subsetted pixel range in
468  satang
469  (note that satang may not give proper values in pixel subsets with this
470  but it may be academic with...)
471  W. Robinson, SAIC, 10 Sep 2010 use vectors instead of lat, lon when
472  doing interpolation of ctl pt lat, lon to
473  all samples
474  */
475 
476 int readl1_czcs(filehandle *file, int32_t recnum, l1str *l1rec) {
477  /*void czcscal_( int *, float[], float[], short *, int *, float * );*/
478  void satang_(double *, double *, float *, float *, float *, float *,
479  float *, float *, float *, float *);
480  void sunangs_(int *, int *, float *, float *, float *, float *, float *);
481  int ll2vec(float *, float *);
482  int vec2ll(float *, float *);
483  short cnt_vec[NBND_CZCS];
484  float lt_lcl[NBND_CZCS], yout1, yout2, yout3, llvec[2], vec[3], gmt;
485  int ipx, ibnd, orbit, i, tpix;
486  uint8 cal_sum[5], cal_scan[6];
487  int32 status, navbad, ltsat;
488  double pi, radeg;
489  int32_t ib;
490 
491  int32_t nwave = l1rec->l1file->nbands;
492  int32_t *bindx = l1rec->l1file->bindx;
493  char *cal_path = file->calfile;
494  /* load scan times */
495  /* int year = *l1rec->year;
496  int day = *l1rec->day;*/
497  /* int32 msec = l1rec->msec[recnum];*/
498 
499  float lonbuf[1968], latbuf[1968], senzbuf[1968], senabuf[1968];
500 
501  /*
502  * read in the line of counts, latitude, longitude
503  */
504  status = rdSDS(file->sd_id, "band1", recnum, 0, 1, npix, (VOIDP) counts);
505  status = rdSDS(file->sd_id, "band2", recnum, 0, 1, npix,
506  (VOIDP) (counts + npix));
507  status = rdSDS(file->sd_id, "band3", recnum, 0, 1, npix,
508  (VOIDP) (counts + 2 * npix));
509  status = rdSDS(file->sd_id, "band4", recnum, 0, 1, npix,
510  (VOIDP) (counts + 3 * npix));
511  status = rdSDS(file->sd_id, "band5", recnum, 0, 1, npix,
512  (VOIDP) (counts + 4 * npix));
513 
514  status = rdSDS(file->sd_id, "gain", recnum, 0, 1, 1, (VOIDP) gain);
515  status = rdSDS(file->sd_id, "latitude", recnum, 0, 1, nctl_pt,
516  (VOIDP) ctl_pt_lat);
517  status = rdSDS(file->sd_id, "longitude", recnum, 0, 1, nctl_pt,
518  (VOIDP) ctl_pt_lon);
519 
520  status = rdSDS(file->sd_id, "cal_sum", recnum, 0, 1, 5, (VOIDP) cal_sum);
521  status = rdSDS(file->sd_id, "cal_scan", recnum, 0, 1, 6, (VOIDP) cal_scan);
522  /*
523  * flag setting for entire line: set navfail if cal_sum shows problems
524  * (suspect it is * bad ephemeris or atitude) and set hilt if bands
525  * missing
526  */
527  ltsat = 0;
528  navbad = 0;
529  if ((cz_band_present & 0XF8) != 0XF8)
530  ltsat = 1;
531  else {
532  if ((cal_sum[3] != 0) || (cal_sum[4] != 0))
533  navbad = 1;
534  if ((cal_scan[0] != 0) || (cal_scan[1] != 0) ||
535  (cal_scan[2] != 0) || (cal_scan[3] != 0) ||
536  (cal_scan[4] != 0)) {
537  ltsat = 1;
538  navbad = 1;
539  }
540  }
541  /*
542  * calibrate the radiances and set flags for pixels
543  */
544  for (ipx = 0; ipx < npix; ipx++) {
545  for (ibnd = 0; ibnd < NBND_CZCS; ibnd++) {
546  cnt_vec[ibnd] = *(counts + ipx + ibnd * npix);
547  if ((cnt_vec[ibnd] == 255) || (ltsat == 1)) l1rec->hilt[ipx] = 1;
548  }
549  *(ring_sat + ipx) = ((cnt_vec[0] == 255) || (cnt_vec[1] == 255)
550  || (cnt_vec[2] == 255)) ? 1 : 0;
551  /*
552  * call the fortran calibration routine
553  */
554  orbit = (int) file->orbit_number;
555  /*czcscal_( &orbit, slope + recnum * 6,
556  intercept + recnum * 6, cnt_vec, &lgain, lt_lcl );*/
557  status = get_czcscal(cal_path, orbit, syear, sday, msec[recnum],
558  cnt_vec, slope[4], intercept[4], gain[0], lt_lcl);
559  /*
560  * assign the calibrated radiances
561  */
562  for (ibnd = 0; ibnd < nwave; ibnd++) {
563  ib = bindx[ ibnd ];
564  l1rec->Lt[ ipx * nwave + ib ] = lt_lcl[ibnd];
565  }
566  /*
567  * save the 750 nm data here
568  */
569  *(lt750 + ipx) = lt_lcl[4];
570 
571  if (navbad == 1) l1rec->navfail[ipx] = 1;
572  }
573  /*
574  * set the ringing mask
575  */
577  /*
578  * get navigation at all points from control point lat and lon values
579  * use spline interpolation to fill in, do in vector space to avoid date
580  * line discontinuity
581  */
582  for (i = 0; i < nctl_pt; i++) {
583  llvec[0] = ctl_pt_lat[i];
584  llvec[1] = ctl_pt_lon[i];
585  if (ll2vec(llvec, vec) == 0) {
586  ctl_pt_vx[i] = *vec;
587  ctl_pt_vy[i] = *(vec + 1);
588  ctl_pt_vz[i] = *(vec + 2);
589  } else {
590  fprintf(stderr, "-E %s Line %d: ll2vec failure.\n",
591  __FILE__, __LINE__);
592  exit(1);
593  }
594  }
595  spline(ctl_pt_x, ctl_pt_vx, nctl_pt, 1e30, 1e30, y2_vx);
596  spline(ctl_pt_x, ctl_pt_vy, nctl_pt, 1e30, 1e30, y2_vy);
597  spline(ctl_pt_x, ctl_pt_vz, nctl_pt, 1e30, 1e30, y2_vz);
598  for (i = 0; i < npix; i++) {
599  tpix = i * ninc /*+ spix*/;
600  splint(ctl_pt_x, ctl_pt_vx, y2_vx, nctl_pt, tpix, &yout1);
601  splint(ctl_pt_x, ctl_pt_vy, y2_vy, nctl_pt, tpix, &yout2);
602  splint(ctl_pt_x, ctl_pt_vz, y2_vz, nctl_pt, tpix, &yout3);
603 
604  *vec = yout1;
605  *(vec + 1) = yout2;
606  *(vec + 2) = yout3;
607  vec2ll(vec, llvec);
608 
609  l1rec->lon[i] = llvec[1];
610  l1rec->lat[i] = llvec[0];
611  }
612  /*
613  * calculate geometry. For sensor angles, use the orbit info if the
614  * position error is available and within error threshold
615  */
616  if ((*(pos_err + recnum) >= 0.) &&
617  (*(pos_err + recnum) < POS_ERR_THRESH)) {
618  cz_posll_2_satang((pos + 3 * recnum), npix, l1rec->lat, l1rec->lon,
619  l1rec->senz, l1rec->sena);
620  } else {
621  pi = PI;
622  radeg = RADEG;
623  for (i = 0; i < 1968; i++) {
624  lonbuf[i] = 0.0;
625  latbuf[i] = 0.0;
626  senzbuf[i] = 0.0;
627  senabuf[i] = 0.0;
628  }
629 
630  for (i = 0; i < npix; i++) {
631  lonbuf[i] = l1rec->lon[i];
632  latbuf[i] = l1rec->lat[i];
633  }
634 
635  satang_(&pi, &radeg, tilt + recnum, att_ang + 3 * recnum + 1,
636  att_ang + 3 * recnum + 2, att_ang + 3 * recnum, lonbuf,
637  latbuf, senzbuf, senabuf);
638 
639  for (i = 0; i < npix; i++) {
640  l1rec->senz[i] = senzbuf[i];
641  l1rec->sena[i] = senabuf[i];
642  }
643  }
644 
645 
646  /* for( i = spix; i < epix; i = i + ninc )*/
647  for (i = 0; i < npix; i++) {
648  gmt = (float) msec[recnum] / (1000. * 3600);
649  sunangs_(&syear, &sday, &gmt, (l1rec->lon) + i, (l1rec->lat) + i,
650  (l1rec->solz) + i, (l1rec->sola) + i);
651  }
652  /*
653  * set scan time
654  */
655  double secs = (double) (msec[recnum] / 1000.0);
656  int16_t year = syear;
657  int16_t day = sday;
658 
659  if (msec[recnum] < msec[recnum - 1]) { /* adjust for day rollover */
660  day += 1;
661  if (day > (365 + (year % 4 == 0))) {
662  year += 1;
663  day = 1;
664  }
665  }
666  l1rec->scantime = yds2unix(year, day, secs);
667 
668  return (status);
669 }
670 
671 /* W. Robinson, SAIC, 6 Jan 2005 include the pos, pos_err arrays in free */
672 
673 int closel1_czcs(filehandle *file) {
674  free(msec);
675  free(tilt);
676  free(att_ang);
677  free(counts);
678  free(ctl_pt_lat);
679  free(ctl_pt_lon);
680  free(ctl_pt_vx);
681  free(ctl_pt_vy);
682  free(ctl_pt_vz);
683  free(y2_vx);
684  free(y2_vy);
685  free(y2_vz);
686  free(ctl_pt_x);
687  free(ctl_pt_cols);
688  free(slope);
689  free(intercept);
690  free(lt750);
691  free(ring_sat);
692  free(pos);
693  free(pos_err);
694 
695  SDend(file->sd_id);
696 
697  return (0);
698 }
699 
700 #define re 6378.137
701 #define f 1. / 298.257
702 #define omf2 ( 1 - f ) * ( 1 - f )
703 
704 int cz_posll_2_satang(float *pos, int npix, float *lat, float *lon,
705  float *senz, float *sena)
706 /*******************************************************************
707 
708  cz_posll_2_satang
709 
710  purpose: convert satellite position, lat, lon into satellite zenith
711  and azimuth for a line of czcs data
712 
713  Returns type: int - 0 if no problems
714 
715  Parameters: (in calling order)
716  Type Name I/O Description
717  ---- ---- --- -----------
718  float * pos I position (x, y, z) of
719  satellite in meters
720  int npix I # pixels of lat, lon to find
721  sat angles for
722  float * lat I latitude array in degrees E
723  float * lon I longitude array in degrees E
724  float * senz O sensor zenith in degrees
725  float * sena O sensor azimuth in degrees
726 
727  Modification history:
728  Programmer Date Description of change
729  ---------- ---- ---------------------
730  W. Robinson, SAIC 5-Jan-2006 Original development
731 
732  *******************************************************************/ {
733  double up[3], ea[3], no[3], radeg, upxy, xlmat[3][3], xlatg, gv[3];
734  double r, rh[3], rl[3];
735  int i, j;
736 
737  radeg = 90. / asin(1.);
738 
739  for (i = 0; i < npix; i++) {
740  /*
741  * Compute the local vertical, East and North unit vectors
742  */
743  up[0] = cos(*(lat + i) / radeg) * cos(*(lon + i) / radeg);
744  up[1] = cos(*(lat + i) / radeg) * sin(*(lon + i) / radeg);
745  up[2] = sin(*(lat + i) / radeg);
746 
747  upxy = sqrt(up[0] * up[0] + up[1] * up[1]);
748 
749  ea[0] = -up[1] / upxy;
750  ea[1] = up[0] / upxy;
751  ea[2] = 0.;
752 
753  cross_prod(up, ea, no);
754 
755  /*
756  * Store vectors in 3x3 array
757  */
758  for (j = 0; j < 3; j++) {
759  xlmat[0][j] = ea[j];
760  xlmat[1][j] = no[j];
761  xlmat[2][j] = up[j];
762  }
763  /*
764  * Compute geocentric position vector
765  */
766  xlatg = radeg * atan(tan(*(lat + i) / radeg) * omf2);
767  gv[0] = cos(xlatg / radeg) * cos(*(lon + i) / radeg);
768  gv[1] = cos(xlatg / radeg) * sin(*(lon + i) / radeg);
769  gv[2] = sin(xlatg / radeg);
770 
771  r = re * (1. - f) /
772  sqrt(1. - (2. - f) * f * pow(cos(xlatg / radeg), 2));
773 
774  /* note that pos is in m but rest of constants are in km, hence the
775  1000 factor */
776  for (j = 0; j < 3; j++)
777  rh[j] = pos[j] / 1000. - r * gv[j];
778  /*
779  * Transform the pixel-to-spacecraft vector into the local frame
780  */
781  matrix_mult(rh, xlmat, rl);
782  /*
783  * Compute the sensor zenith and azimuth
784  * if zenith close to 0, set azimuth to 0 (and normalize azimuth)
785  */
786  *(senz + i) = radeg * atan2(sqrt(rl[0] * rl[0] + rl[1] * rl[1]),
787  rl[2]);
788  if (*(senz + i) > 0.05)
789  *(sena + i) = radeg * atan2(rl[0], rl[1]);
790  else
791  *(sena + i) = 0.;
792 
793  if (*(sena + i) < 0.) *(sena + i) += 360.;
794  }
795  /*
796  * and end
797  */
798  return 0;
799 }
800 
801 void matrix_mult(double vecin[3], double matrix[3][3], double vecout[3])
802 /*******************************************************************
803 
804  matrix_mult
805 
806  purpose: multiply a vector by a matrix and produce the output vector
807 
808  Returns type: void
809 
810  Parameters: (in calling order)
811  Type Name I/O Description
812  ---- ---- --- -----------
813  double[3] vecin I input vector
814  double[3][3] matrix I rotation matrix
815  double[3] vecout O rotated vector
816 
817  Modification history:
818  Programmer Date Description of change
819  ---------- ---- ---------------------
820  W. Robinson, SAIC 2-Nov-2005 Original development
821 
822  *******************************************************************/ {
823  int i, j;
824  for (i = 0; i < 3; i++) {
825  vecout[i] = 0.;
826  for (j = 0; j < 3; j++) {
827  vecout[i] += matrix[i][j] * vecin[j];
828  }
829  }
830 }
831 
832 void cross_prod(double *v1, double *v2, double *vout)
833 /*******************************************************************
834 
835  cross_prod
836 
837  purpose: produce the cross product of 2 vectors
838 
839  Returns type: void
840 
841  Parameters: (in calling order)
842  Type Name I/O Description
843  ---- ---- --- -----------
844  double * v1 I vector 1
845  double * v2 I vector 1
846  double * vout O v1 X v2
847 
848  Modification history:
849  Programmer Date Description of change
850  ---------- ---- ---------------------
851  W. Robinson, SAIC 4-Jan-2006 Original development
852 
853  *******************************************************************/ {
854  *vout = v1[1] * v2[2] - v1[2] * v2[1];
855  *(vout + 1) = v1[2] * v2[0] - v1[0] * v2[2];
856  *(vout + 2) = v1[0] * v2[1] - v1[1] * v2[0];
857 }
integer, parameter int16
Definition: cubeio.f90:3
const int bindx[3]
Definition: DbLutNetcdf.cpp:31
#define MIN(x, y)
Definition: rice.h:169
int r
Definition: decode_rs.h:73
int16 eday
Definition: l1_czcs_hdf.c:16
int czcs_ring(int gain, float *lt750, char *ring_sat, l1str *l1rec)
Definition: l1_czcs_hdf.c:41
int ll2vec(float *ll, float *vec)
Definition: ll2vec.c:5
int j
Definition: decode_rs.h:73
int32_t day
int status
Definition: l1_czcs_hdf.c:31
double yds2unix(int16_t year, int16_t day, double secs)
Definition: yds2unix.c:7
float32 * tilt
Definition: l1_czcs_hdf.c:33
constexpr double RADEG
int16 * gain
Definition: l1_czcs_hdf.c:32
*********************************************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
int16_t fileID
#define FAIL
Definition: ObpgReadGrid.h:18
float * lt750
Definition: l1_czcs_hdf.c:37
#define NULL
Definition: decode_rs.h:63
read l1rec
#define NBND_CZCS
Definition: l1_czcs_hdf.c:11
#define NBND
Definition: l1_czcs_hdf.c:155
const double pi
float32 * ctl_pt_lat
Definition: l1_czcs_hdf.c:34
int32 * ctl_pt_cols
Definition: l1_czcs_hdf.c:35
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
float * ctl_pt_vy
Definition: l1_czcs_hdf.c:36
float32 * intercept
Definition: l1_czcs_hdf.c:33
#define VOIDP
Definition: hdf5utils.h:11
float32 * pos
Definition: l1_czcs_hdf.c:34
float * lat
int32 * msec
Definition: l1_czcs_hdf.c:30
uint8 * counts
Definition: l1_czcs_hdf.c:29
void cross_prod(double *v1, double *v2, double *vout)
Definition: l1_czcs_hdf.c:832
int16 eyear
Definition: l1_czcs_hdf.c:16
float ** matrix(long nrl, long nrh, long ncl, long nch)
Definition: nrutil.c:60
uint8 cz_band_present
Definition: l1_czcs_hdf.c:29
int syear
Definition: l1_czcs_hdf.c:14
float * ctl_pt_vz
Definition: l1_czcs_hdf.c:36
int32 nscan
Definition: l1_czcs_hdf.c:18
int32 smsec
Definition: l1_czcs_hdf.c:15
int readl1_czcs(filehandle *file, int32_t recnum, l1str *l1rec)
Definition: l1_czcs_hdf.c:476
#define POS_ERR_THRESH
Definition: l1_czcs_hdf.c:12
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
int sday
Definition: l1_czcs_hdf.c:14
#define PI
Definition: l3_get_org.c:6
int32 nsta
Definition: l1_czcs_hdf.c:26
dtype
Definition: dataset.h:17
int cz_posll_2_satang(float *pos, int npix, float *lat, float *lon, float *senz, float *sena)
Definition: l1_czcs_hdf.c:704
void sunangs_(int *, int *, float *, float *, float *, float *, float *)
float * y2_vz
Definition: l1_czcs_hdf.c:36
read recnum
float * y2_vx
Definition: l1_czcs_hdf.c:36
int closel1_czcs(filehandle *file)
Definition: l1_czcs_hdf.c:673
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_hdf.c:159
void satang_(double *, double *, float *, float *, float *, float *, float *, float *, float *, float *)
#define re
Definition: l1_czcs_hdf.c:700
int rdSDS(int32 fileID, const char sdsname[], int32 start1, int32 start2, int32 edges1, int32 edges2, VOIDP array_data)
int32 ninc
Definition: l1_czcs_hdf.c:27
#define NGAIN
Definition: l1_czcs_hdf.c:156
void cdata_()
int openl1_czcs(filehandle *file)
Definition: l1_czcs_hdf.c:360
int want_verbose
char * ring_sat
Definition: l1_czcs_hdf.c:38
int32 nctl_pt
Definition: l1_czcs_hdf.c:35
#define NEPOCH
Definition: l1_czcs_hdf.c:157
int vec2ll(float *vec, float *ll)
Definition: ll2vec.c:52
int32 emsec
Definition: l1_czcs_hdf.c:17
#define omf2
Definition: l1_czcs_hdf.c:702
float * ctl_pt_x
Definition: l1_czcs_hdf.c:36
float32 * ctl_pt_lon
Definition: l1_czcs_hdf.c:34
void matrix_mult(double vecin[3], double matrix[3][3], double vecout[3])
Definition: l1_czcs_hdf.c:801
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
int32 dpix
Definition: l1_czcs_hdf.c:21
int32 spix
Definition: l1_czcs_hdf.c:20
char * name
Definition: Granule.c:1234
float * y2_vy
Definition: l1_czcs_hdf.c:36
Extra metadata that will be written to the HDF4 file l2prod rank
#define CZCS
Definition: sensorDefs.h:17
subroutine splint(xa, ya, y2a, n, x, y)
float32 * slope
Definition: l1_czcs_hdf.c:33
float * lon
int32 epix
Definition: l1_czcs_hdf.c:22
int lgain
Definition: l1_czcs_hdf.c:31
float32 * pos_err
Definition: l1_czcs_hdf.c:34
#define f
Definition: l1_czcs_hdf.c:701
float32 * att_ang
Definition: l1_czcs_hdf.c:33
int getHDFattr(int32 fileID, const char attrname[], const char sdsname[], VOIDP data)
int i
Definition: decode_rs.h:71
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int32 npix
Definition: l1_czcs_hdf.c:19
float * ctl_pt_vx
Definition: l1_czcs_hdf.c:36