OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
calib_calibrate_l1a.c
Go to the documentation of this file.
1 /*-----------------------------------------------------------------------------
2  File : calibrate_l1a.c
3 
4  Contents:
5  calibrate_l1a - takes an array of level-1A raw counts and returns
6  a corresponding arry of radiance values after
7  applying the sensor calibration
8 
9  Other relevant files:
10  cal_l1a.h - various #defined constants, other include files
11  (get_cal.h, getcal_proto.h, call1a_proto.h) and
12  also includes hdf.h
13  get_cal.c - a higher layer of calibration input functions
14  get_cal_misc.c - a lower layer of calibration input functions
15  tcal_l1a.c - a test routine to test calibrate_l1a
16 
17  Notes:
18  HDF library used - HDF3.3r3p1
19 
20  Modification history:
21  Programmer Organization Date Description of change
22  -------------- ------------ -------- ---------------------
23  Lakshmi Kumar Hughes STX 06/07/94 Original development
24  Lakshmi Kumar Hughes STX 10/30/95 Filtering scan_mod corr.
25  to data of dtype = SOL &
26  TDI.
27  Lakshmi Kumar Hughes STX 11/20/95 Filtering scan_mod corr.
28  to data of dtype = IGC
29  Lakshmi Kumar Hughes STX 01/25/96 Modified local variable
30  "rads" to global var
31  Gene Eplee GSC 05/24/96 Implemented quadratic
32  temporal calibration
33  correction (gain & offset)
34  Lakshmi Kumar Hughes STX 10/17/96 Removed cal_year, cal_day
35  output arguments.
36  Ref. V5.0 I/O specs.
37  Gene Eplee GSC 01/14/97 Reworked populating of
38  g_f at the calibration
39  knees.
40  Gene Eplee GSC 04/25/97 Modified system gain and
41  offset corrections to
42  allow calling routine to
43  override cal table values.
44  Gene Eplee GSC 09/09/97 Modified dark restore
45  subtraction to use
46  mean values for GAC, LAC,
47  and HRPT data.
48 -----------------------------------------------------------------------------*/
49 
50 #include "calib_cal_l1a.h"
51 #include "calib_call1a_proto.h"
52 
53 /*-----------------------------------------------------------------------------
54  Function: calibrate_l1a
55 
56  Returns: int32 (status)
57  Returns a status code of 0 when successful. Otherwise returns
58  -1 - to indicate calibration file open/close error
59  -2 - to indicate read error
60  -3 - to indicate time error (if the given time cannot be found)
61  -4 - to indicate insufficient memory error
62 
63  Description:
64  The function calibrate_l1a takes an array of level-1A raw counts
65  and returns a corresponding array of radiance values (level-1B data)
66  after applying the sensor calibration.
67 
68  Arguments: (in calling order)
69  Type Name I/O Description
70  ---- ---- --- -----------
71  char * cal_path I calibration file path
72  int16 syear I year of data start time
73  int16 sday I day-of-year for data start time
74  int32 smsec I milliseconds of the day for data
75  start tm as returned by get_l1a_open
76  int16 eday I day-of-year for data end time
77  int32 msec I millisecs of the day for data start
78  time as returned by get_l1a_record
79  char * dtype I data type flag
80  int32 st_samp I start pixel/sample number to process
81  int32 nsamp I samples per scan line
82  int16 * dark_rest I dark restore values for all 8 bands
83  float32 * dark_mean I mean dark restore values for all 8
84  bands
85  int16 * gain I gains for all 8 bands; as returned
86  by get_l1a_record
87  int16 * tdi I input TDI for all 8 bands
88  int16 * scan_temp I digitized scan temperatures
89  int16 * side I mirror side of scan line
90  int16 * l1a_data I Level-1A data; Band interleave by pixel
91  float32 * l1b_data O Sensor calibrated radiance values
92  corresponding to l1a_data
93  struct cal_mod_struc I Structure for override control of
94  cal_mod system gain and offset
95 
96  Notes:
97 
98  Modification history:
99  Programmer Organization Date Description of change
100  -------------- ------------ -------- ---------------------
101  Lakshmi Kumar Hughes STX 06/07/94 Original development
102  Lakshmi Kumar Hughes STX 02/07/95 Bug fix--Applied time
103  dependent correction.
104  Added cal_year & cal_day
105  arguments.
106  (ref v4.2 I/O specs.)
107  Lakshmi Kumar Hughes STX 10/30/95 Filtering scan_mod corr.
108  to data of dtype = SOL &
109  TDI.
110  Lakshmi Kumar Hughes STX 11/20/95 Filtering scan_mod corr.
111  to data of dtype = IGC
112  Lakshmi Kumar Hughes STX 01/25/96 Modified loc var "rads"
113  to global var "cal_rads"
114  to make it accessible to
115  l1a_read rtn.
116  Gene Eplee GSC 05/24/96 Implemented quadratic
117  temporal calibration
118  correction (gain & offset)
119  Lakshmi Kumar Hughes STX 10/17/96 Removed cal_year, cal_day
120  output arguments.
121  Ref. V5.0 I/O specs.
122  Gene Eplee GSC 01/14/97 Reworked populating of
123  g_f at the calibration
124  knees.
125  Gene Eplee GSC 04/25/97 Modified system gain and
126  offset corrections to
127  allow calling routine to
128  override cal table values.
129  Lakshmi Kumar GSC 05/02/97 Removed non prototype
130  defns.
131  Gene Eplee GSC 09/09/97 Modified dark restore
132  subtraction to use
133  mean values for GAC, LAC,
134  and HRPT data.
135 -----------------------------------------------------------------------------*/
136 float32 cal_counts[NBANDS][4][5]; /*digital cnts (zero-offs corrected */
137 float32 cal_rads[NBANDS][4][5]; /*radiances corresponding to knees */
138 
139 int32 calibrate_l1a(char *cal_path, int16 syear, int16 sday, int32 smsec,
140  int16 eday, int32 msec, char *dtype, int32 st_samp,
141  int32 nsamp, int16 *dark_rest, float32 *dark_mean, int16 *gain,
142  int16 *tdi, int16 *scan_temp, int16 side, int16 *l1a_data,
143  float32 *l1b_data, cal_mod_struc *cal_mod) {
144 
145  static int32 called_get_call = 0, first_call = 1;
146  static int16 pr_syear = 0, pr_sday = 0;
147  static int32 pr_smsec = 0;
148  static int16 pr_gain[NBANDS] = {-1, -1, -1, -1, -1, -1, -1, -1};
149  static int16 pr_tdi[NBANDS] = {-1, -1, -1, -1, -1, -1, -1, -1};
150  static float32 pr_cal_gain[NBANDS] = {-1, -1, -1, -1, -1, -1, -1, -1};
151 
152 
153  int16 cal_year, cal_day; /* calibrate entry year & day */
154  int16 status, n;
155  int16 ref_year, ref_day, ref_min; /* calibration model reference time */
156  int16 band, knee, count, pixel; /* local indeces */
157  int16 count1, count2, tdi_flag, gain_flag; /* temporary storage fields */
158  int16 l1_data; /*temp space for dark-restore corrected cnt*/
159  int32 ref_jday, data_jday; /* julian reference day and data day */
160  float32 slope; /*loc_slope for converting l1a cnt and rad */
161  static float32 cal_gain[NBANDS]; /* calibration model system gain */
162  float64 delta_min; /* calibration model timescales */
163  static float64 delta_t;
164  static float32 g_f[NBANDS][1024]; /*gain factors look up table */
165 
166  /*get_cal output parameters */
167  static float32 temps[256][NBANDS]; /*temp correction co_efficients */
168  static float32 scan_mod[2][1285]; /*scan modulation corr factors */
169  static float32 mirror[2][NBANDS]; /*mirror side-0 and 1 corr factors */
170  static float64 tfactor_const[NBANDS]; /* time dependent constant term */
171  static float64 tfactor_linear[NBANDS]; /* time dependent linear term */
172  static float64 tfactor_quadratic[NBANDS]; /* time dependent quadratic term */
173  static float32 cal_offset[NBANDS]; /* calibration model system offset */
174 
175  for (band = 0; band < NBANDS; band++) {
176  if (tdi[band] < 0)
177  tdi[band] = 0;
178  if (tdi[band] > 255)
179  tdi[band] = 255;
180  if (scan_temp[band] < 0)
181  scan_temp[band] = 0;
182  if (scan_temp[band] > 255)
183  scan_temp[band] = 255;
184  }
185 
186  for (band = 0; band < NBANDS && tdi[band] == pr_tdi[band]; band++);
187 
188  if (band < NBANDS)
189  tdi_flag = 1;
190  else
191  tdi_flag = 0;
192 
193  if (first_call || syear != pr_syear || sday != pr_sday || smsec != pr_smsec
194  || tdi_flag) {
195 
196  first_call = 0;
197  pr_syear = syear;
198  pr_sday = sday;
199  pr_smsec = smsec;
200  for (band = 0; band < NBANDS; band++)
201  pr_tdi[band] = tdi[band];
202 
203  if ((status = get_cal(cal_path, syear, sday, eday, msec, dtype, tdi,
204  &cal_year, &cal_day, &ref_year, &ref_day, &ref_min, temps,
205  scan_mod, mirror, tfactor_const, tfactor_linear,
206  tfactor_quadratic, cal_offset, cal_counts, cal_rads)) < 0)
207  return status;
208  called_get_call = 1;
209 #ifdef DEBUG
210  printf("\ncalled get_cal\n");
211 #endif
212 
213  /* Define the timescale for the calibration system gain. */
214  ref_jday = jul2jday(ref_year, ref_day);
215  data_jday = jul2jday(syear, sday);
216  delta_min = (float64) ((data_jday - ref_jday)*1440 - ref_min);
217  delta_t = delta_min + (float64) msec / 60000.0;
218 
219  /* Compute the system gain */
220  for (band = 0; band < NBANDS; band++) {
221  cal_gain[band] = (float32) (tfactor_const[band]
222  + delta_t * tfactor_linear[band] + delta_t * delta_t * tfactor_quadratic[band]);
223  if (cal_gain[band] != pr_cal_gain[band]) {
224  printf("band %d: dt = %dmin cal_offset = %g cal_gain = %7.5f + %.3e * dt + %.3e * dt * dt = %g\n",
225  band, (int32_t) delta_t, cal_offset[band], tfactor_const[band], tfactor_linear[band], tfactor_quadratic[band], cal_gain[band]);
226  pr_cal_gain[band] = cal_gain[band];
227  }
228  }
229 
230  }
231 
232  for (band = 0; band < NBANDS && gain[band] == pr_gain[band]; band++);
233 
234  if (band < NBANDS)
235  gain_flag = 1;
236  else
237  gain_flag = 0;
238 
239  if (called_get_call || gain_flag) {
240  called_get_call = 0;
241  for (band = 0; band < NBANDS; band++)
242  pr_gain[band] = gain[band];
243  for (band = 0; band < NBANDS; band++) {
244 #ifdef DEBUG
245  printf("TDI band gain knee cnt1 cnt2 slope\n");
246  printf("-----------------------------------------\n\n");
247 #endif
248  for (knee = 1; knee <= 4; knee++) {
249  n = 1;
250  while (((int16) cal_counts[band][gain[band]][knee] ==
251  (int16) cal_counts[band][gain[band]][knee - n]) && n <= knee)
252  n++;
253  count1 = (int16) cal_counts[band][gain[band]][knee - n] + 1;
254  count2 = (int16) cal_counts[band][gain[band]][knee];
255  if (knee == 1)
256  count1 = 0;
257  if (knee == 4)
258  count2 = 1023;
259  slope = (cal_rads[band][gain[band]][knee] -
260  cal_rads[band][gain[band]][knee - n]) /
261  (cal_counts[band][gain[band]][knee] -
262  cal_counts[band][gain[band]][knee - n]);
263 #ifdef DEBUG
264  printf("%3d %3d %3d %4d %5d %5d %10.8f\n",
265  tdi[band], band, gain[band], knee, count1, count2, slope);
266 #endif
267  for (count = count1; count <= count2; count++)
268  g_f[band][count] =
269  slope * (count - cal_counts[band][gain[band]][knee - n]) +
270  cal_rads[band][gain[band]][knee - n];
271  }
272  }
273  }
274 
275 #ifdef DEBUG
276  printf("\n\n------------- g_f table values --------------\n");
277  printf("\nindex b1\tb2\tb3\tb4\tb5\tb6\tb7\tb8\n");
278  for (i = 0; i < 1024; i++) {
279  printf("\n%4d", i);
280  for (j = 0; j < NBANDS; j++)
281  printf("%9.5f", g_f[j][i]);
282  }
283 
284 #endif
285 
286 
287  if (st_samp < 1)
288  st_samp = 1;
289  for (band = 0; band < BANDS; band++) {
290 
291  /* Check for override on system gain and offset */
292  /* and pass back the gain, offset used */
293  if (cal_mod) {
294  if (cal_mod->flag == 1) {
295  cal_gain[band] = (float32) cal_mod->gain[band];
296  cal_mod->offset[band] = cal_offset[band];
297  } else if (cal_mod->flag == 2) {
298  cal_offset[band] = (float32) cal_mod->offset[band];
299  cal_mod->gain[band] = cal_gain[band];
300  } else if (cal_mod->flag == 3) {
301  cal_gain[band] = (float32) cal_mod->gain[band];
302  cal_offset[band] = (float32) cal_mod->offset[band];
303  } else {
304  cal_mod->gain[band] = cal_gain[band];
305  cal_mod->offset[band] = cal_offset[band];
306  }
307  }
308 
309  for (pixel = st_samp - 1; pixel <= st_samp + nsamp - 2; pixel++) {
310 
311  /* Use dark restore values from each scan line for calibration data.
312  Use mean dark restore values for imaging data. */
313  if (!strcmp(dtype, "SOL") ||
314  !strcmp(dtype, "TDI") ||
315  !strcmp(dtype, "IGC") ||
316  !strcmp(dtype, "LUN"))
317  l1_data = l1a_data[pixel * BANDS + band] - dark_rest[band];
318  else
319  l1_data = l1a_data[pixel * BANDS + band] - (int16) (dark_mean[band] + 0.5);
320 
321  if (l1_data < 0) l1_data = 0;
322  if (l1_data > 1023) l1_data = 1023;
323  l1b_data[pixel * BANDS + band] = g_f[band][l1_data] * mirror[side][band] * temps[scan_temp[band]][band];
324  if (strcmp(dtype, "SOL") && strcmp(dtype, "TDI") && strcmp(dtype, "IGC"))
325  l1b_data[pixel * BANDS + band] *= scan_mod[band % 2][pixel];
326  l1b_data[pixel * BANDS + band] =
327  l1b_data[pixel * BANDS + band] * cal_gain[band] + cal_offset[band];
328  if (l1b_data[pixel * BANDS + band] > 100)
329  printf("l1b_data %g l1a_data %d cal_gain %g g_f %g mirror %g temps %g scan_mod %g\n", l1b_data[pixel * BANDS + band], l1a_data[pixel * BANDS + band], cal_gain[band], g_f[band][l1_data], mirror[side][band], temps[scan_temp[band]][band], scan_mod[band % 2][pixel]);
330 
331  }
332  }
333  return SUCCEED;
334 }
335 
336 /*---------------------------------------------------------------------------*/
337 
338 /*------------------------------------------------------------------------------
339 
340  Function: jul2jday
341 
342  Returns type: int
343 
344  Description:
345 
346  Convert year and day-of-year pair to Julian Day. This routine uses
347  the same cal2jday routine by just specifying the month to be 1 and
348  the day-of-month to be day-of-year.
349 
350  Parameters: (in calling order)
351  Type Name I/O Description
352  ---- ---- --- -----------
353  int year I year
354  int yday I day of year [1,366]
355 
356  Modification history:
357  Programmer Date Description of change
358  ---------- ---- ---------------------
359  Frank Chen 30-Jul-1993 Original development
360 
361 ------------------------------------------------------------------------------*/
362 
363 int jul2jday(int year, int yday) {
364  int jday;
365 
366  jday = cal2jday(year, 1, yday);
367  return (jday);
368 }
369 
370 /*------------------------------------------------------------------------------
371 
372  Function: cal2jday
373 
374  Returns type: int
375 
376  Description:
377 
378  This function converts a calendar date to the corresponding Julian
379  day starting at noon on the calendar date. The algorithm used is
380  from Van Flandern and Pulkkinen, Ap. J. Supplement Series 41,
381  November 1979, p.400
382  This will also work when month is 1 and day-of-month is specified
383  as day-of-year.
384 
385  Parameters: (in calling order)
386  Type Name I/O Description
387  ---- ---- --- -----------
388  int year I year
389  int month I month [1,12]
390  int mday I day of month [1,31]
391 
392  Modification history:
393  Programmer Date Description of change
394  ---------- ---- ---------------------
395  Fred Patt 04-Nov-1992 Original written in FORTRAN(jd)
396  Frank Chen 30-Jul-1993 Rewrite in C
397 
398 ------------------------------------------------------------------------------*/
399 
400 int cal2jday(int year, int month, int mday) {
401  int jday;
402 
403  jday = 367 * year
404  - 7 * (year + (month + 9) / 12) / 4
405  + 275 * month / 9
406  + mday
407  + 1721014;
408  /* additional calculation is needed only for dates outside of */
409  /* the period March 1, 1900 to Feburary 28, 2100 */
410  /*
411  jday = jday + 15 - 3 * ((year + (month - 9) / 7) / 100 + 1) / 4;
412  */
413  return (jday);
414 }
415 
cubeio::int16
integer, parameter int16
Definition: cubeio.f90:3
eday
int16 eday
Definition: l1_czcs_hdf.c:17
j
int j
Definition: decode_rs.h:73
dark_mean
float dark_mean[8]
Definition: l1a_seawifs.c:34
status
int status
Definition: l1_czcs_hdf.c:32
gain
int16 * gain
Definition: l1_czcs_hdf.c:33
calib_call1a_proto.h
band
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second band
Definition: HOWTO_Add_a_product.txt:42
calibrate_l1a
int32 calibrate_l1a(char *cal_path, int16 syear, int16 sday, int32 smsec, int16 eday, int32 msec, char *dtype, int32 st_samp, int32 nsamp, int16 *dark_rest, float32 *dark_mean, int16 *gain, int16 *tdi, int16 *scan_temp, int16 side, int16 *l1a_data, float32 *l1b_data, cal_mod_struc *cal_mod)
Definition: calib_calibrate_l1a.c:139
l1a_data
int16_t * l1a_data
Definition: l1a_seawifs.c:84
dark_rest
int16_t * dark_rest
Definition: l1a_seawifs.c:89
msec
int32 * msec
Definition: l1_czcs_hdf.c:31
syear
int syear
Definition: l1_czcs_hdf.c:15
cal_counts
float32 cal_counts[NBANDS][4][5]
Definition: calib_calibrate_l1a.c:136
jday
int32_t jday(int16_t i, int16_t j, int16_t k)
Definition: jday.c:4
smsec
int32 smsec
Definition: l1_czcs_hdf.c:16
spacetrack_tles.year
year
Definition: spacetrack_tles.py:170
st_samp
int st_samp
Definition: get_cmp.c:33
sday
int sday
Definition: l1_czcs_hdf.c:15
get_cal
int32 get_cal(char *cal_path, int16 syear, int16 sday, int16 eday, int32 msec, char *dtype, int16 *tdi, int16 *cal_year, int16 *cal_day, int16 *ref_year, int16 *ref_day, int16 *ref_min, float32 temps[256][BANDS], float32 scan_mod[2][1285], float32 mirror[2][BANDS], float64 *t_const, float64 *t_linear, float64 *t_quadratic, float32 *cal_offs, float32 counts[BANDS][4][5], float32 rads[BANDS][4][5])
Definition: calib_get_cal.c:98
jul2jday
int jul2jday(int year, int yday)
Definition: calib_calibrate_l1a.c:363
calib_cal_l1a.h
slope
float32 slope[]
Definition: l2lists.h:30
tdi
int16_t tdi[BANDS_DIMS_1A]
Definition: l1a_seawifs.c:37
pixel
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution pixel
Definition: HISTORY.txt:192
ref_day
int16_t ref_day
Definition: l1a_seawifs.c:42
cal2jday
int cal2jday(int year, int month, int mday)
Definition: calib_calibrate_l1a.c:400
dtype
dtype
Definition: DDataset.hpp:31
side
int16_t * side
Definition: l1a_seawifs.c:88
ref_year
int16_t ref_year
Definition: l1a_seawifs.c:41
scan_mod
float scan_mod[2][1285]
Definition: l1a_seawifs.c:45
i
int i
Definition: decode_rs.h:71
NBANDS
@ NBANDS
Definition: make_L3_v1.1.c:53
cal_rads
float32 cal_rads[NBANDS][4][5]
Definition: calib_calibrate_l1a.c:137
BANDS
#define BANDS
count
int count
Definition: decode_rs.h:79


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

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 Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster