OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
l1a_osmi.c
Go to the documentation of this file.
1 #include <hdf4utils.h>
2 #include "l1a.h"
3 #include "eng_qual.h"
4 #include "l1.h"
5 #include "l1a_seawifs.h"
6 #include "call1a_proto_osmi.h"
7 #include "getcal_proto_osmi.h"
8 #include "orbit.h"
9 #include "lunsol.h"
10 
11 #include <hdf.h>
12 #include <mfhdf.h>
13 
14 #define LAC_PIXEL_NUM 1285
15 #define NREC_IN_BUF 10
16 #define NBANDS_OSMI 8
17 
18 static int16_t syear, sday; /* data start date */
19 static int32_t smsec; /* data start time */
20 static int16_t eyear, eday; /* data end date */
21 static int32_t emsec; /* data end time */
22 static int32_t nscan; /* number of scans */
23 static int32_t npix; /* number pixels per scan */
24 
25 static int16_t cal_year;
26 static int16_t cal_day;
27 static int32_t cal_msec;
28 
29 static float eoffset;
30 static float egain[8], temps[8], mirror[8];
31 static float slope[BANDS_DIMS_1A * 96];
32 static float dc[BANDS_DIMS_1A * 96];
33 static float sm[1044];
34 
35 static float t_const[BANDS_DIMS_1A];
36 static float t_linear[BANDS_DIMS_1A];
37 static float t_quadratic[BANDS_DIMS_1A];
38 
39 static char cal_path_tab[128];
40 
41 static int32_t nsta;
42 static int32_t ninc;
43 
44 static int16_t *l1a_data = NULL; /* raw radiances band-int-by-pix */
45 
46 static float *l1b_buffer = NULL; /* l1b radiances band-int-by-pix */
47 static int32_t *msec;
48 static float *sc_ana;
49 static double *sc_ttag;
50 static double *sc_ttag_inc;
51 
52 static char dtype[8];
53 static char cal_path_tab[128];
54 
55 static short l2_flags_buffer[LAC_PIXEL_NUM];
56 
57 static float ylat[LAC_PIXEL_NUM];
58 static float xlon[LAC_PIXEL_NUM];
59 static float solz[LAC_PIXEL_NUM];
60 static float sola[LAC_PIXEL_NUM];
61 static float senz[LAC_PIXEL_NUM];
62 static float sena[LAC_PIXEL_NUM];
63 
64 int32_t get_l1a_rec_osmi(int32_t, int32_t, cal_mod_struc *, float **, int16_t **);
65 
66 void intpos_(double *epoch1, double pos1[3], double vel1[3],
67  double *epoch2, double pos2[3], double vel2[3],
68  double *epoch, double pos [3], double vel [3]);
69 
70 int openl1a_osmi(filehandle *file) {
71 
72  /* */
73  /* input files */
74  /* */
75  char *cal_path = l1_input->calfile;
76 
77  /* */
78  /* get_l1a_open interface */
79  /* */
80  int32_t fileID;
81  int32_t status;
82 
83  /* open the file and get the file ID */
84  fileID = SDstart(file->name, DFACC_RDONLY);
85 
86  if (fileID < 0) {
87  fprintf(stderr,
88  "-E %s Line %d: Error opening %s for reading.\n",
89  __FILE__, __LINE__, file->name);
90  return (1);
91  }
92 
93  status = getHDFattr(fileID, "Start Year", "", (VOIDP) & syear);
94  status = getHDFattr(fileID, "Start Day", "", (VOIDP) & sday);
95  status = getHDFattr(fileID, "Start Millisec", "", (VOIDP) & smsec);
96  status = getHDFattr(fileID, "End Year", "", (VOIDP) & eyear);
97  status = getHDFattr(fileID, "End Day", "", (VOIDP) & eday);
98  status = getHDFattr(fileID, "End Millisec", "", (VOIDP) & emsec);
99  status = getHDFattr(fileID, "Number of Scan Lines", "", (VOIDP) & nscan);
100  status = getHDFattr(fileID, "Data Type", "", (VOIDP) & dtype);
101 
102  status = getHDFattr(fileID, "Pixels per Scan Line", "", (VOIDP) & npix);
103  status = getHDFattr(fileID, "LAC Pixel Start Number", "", (VOIDP) & nsta);
104  status = getHDFattr(fileID, "LAC Pixel Subsampling", "", (VOIDP) & ninc);
105 
106 
107  /* get the current calibration model solution */
108  if (cal_path == NULL) {
109  fprintf(stderr,
110  "-E %s Line %d: No calibration file specified.\n",
111  __FILE__, __LINE__);
112  exit(1);
113  }
114  status = get_cal_osmi(cal_path, syear, sday, eday, smsec,
115  &cal_year, &cal_day, &cal_msec,
116  &eoffset, &egain[0], &temps[0],
117  &mirror[0], &t_const[0], &t_linear[0],
118  &t_quadratic[0], &slope[0], &dc[0],
119  &sm[0]);
120  if (status < 0) {
121  fprintf(stderr,
122  "-E- %s line %d: Error applying calibration table \"%s\".\n",
123  __FILE__, __LINE__, cal_path);
124  exit(status);
125  }
126 
127 
128  file->npix = npix;
129  file->nscan = nscan;
130  file->sensorID = OSMI;
131  file->sd_id = fileID;
132  strcpy(file->spatialResolution, "850 m");
133 
134  strcpy(cal_path_tab, cal_path);
135 
136  l1a_data = (int16_t *) calloc(npix * NBANDS_OSMI*NREC_IN_BUF, sizeof (int16_t));
137  l1b_buffer = (float *) calloc(npix*NBANDS_OSMI, sizeof (float));
138 
139  msec = (int32_t *) calloc(nscan, sizeof (int32_t));
140  status = rdSDS(file->sd_id, "msec", 0, 0, 0, 0, (VOIDP) msec);
141 
142  sc_ana = (float *) calloc(nscan * 40, sizeof (float));
143  status = rdSDS(file->sd_id, "sc_ana", 0, 0, 0, 0, (VOIDP) sc_ana);
144 
145  sc_ttag = (double *) calloc(nscan, sizeof (double));
146  status = rdSDS(file->sd_id, "sc_ttag", 0, 0, 0, 0, (VOIDP) sc_ttag);
147 
148  sc_ttag_inc = (double *) calloc(nscan, sizeof (double));
149  status = rdSDS(file->sd_id, "sc_ttag_inc", 0, 0, 0, 0, (VOIDP) sc_ttag_inc);
150 
151  return (status);
152 }
153 
154 int readl1a_osmi(filehandle *file, int32_t recnum, l1str *l1rec) {
155  /*recnum is 0-based */
156 
157  /* */
158  /* get_l1a_rec interface */
159  /* */
160  static cal_mod_struc cal_mod; /* cal modification structure */
161  static int16_t *l2_flags; /* radiance quality flags for L2 */
162  static float *l1b_data;
163  static int chindx[6] = {0, 1, 2, 4, 6, 7}; /* active bands in L1B data block*/
164 
165  int32_t status = 0;
166 
167  int32_t nwave = l1rec->l1file->nbands;
168  int32_t *bindx = l1rec->l1file->bindx;
169 
170  /* */
171  /* local vars */
172  /* */
173  int32_t ipix; /* pixel number */
174  int32_t iw, ib;
175 
176 
177 
178  /* Get l1a data */
179  /* ------------ */
180  status = get_l1a_rec_osmi(file->sd_id, recnum + 1, &cal_mod,
181  &l1b_data, &l2_flags);
182 
183 
184  /* */
185  /* Copy scan geolocation and view geometry */
186  /* */
187  memcpy(l1rec->lat, ylat, npix * sizeof (float));
188  memcpy(l1rec->lon, xlon, npix * sizeof (float));
189  memcpy(l1rec->solz, solz, npix * sizeof (float));
190  memcpy(l1rec->sola, sola, npix * sizeof (float));
191  memcpy(l1rec->senz, senz, npix * sizeof (float));
192  memcpy(l1rec->sena, sena, npix * sizeof (float));
193 
194  /* */
195  /* Copy L1B radiances, pixel interlaced by band. Add per-band */
196  /* view angles. */
197  /* */
198  for (ipix = 0; ipix < file->npix; ipix++) {
199 
200  if (l1rec->sena[ipix] < -180.0) l1rec->sena[ipix] += 360.0;
201 
202  for (iw = 0; iw < nwave; iw++) {
203  ib = bindx[iw];
204  l1rec->Lt [ipix * nwave + ib] = l1b_data[chindx[iw] * npix + ipix];
205  }
206  l1rec->stlight[ipix] = ((l2_flags[ipix] & STRAYLIGHT) > 0);
207  l1rec->hilt [ipix] = ((l2_flags[ipix] & HILT) > 0);
208  if (recnum < 96) l1rec->hilt[ipix] = 1; /* bad position interp */
209  if (ipix < 15) l1rec->hilt[ipix] = 1; /* instrument shading? */
210  if (ipix > 1028) l1rec->hilt[ipix] = 1; /* instrument shading? */
211  }
212 
213  /* */
214  /* Set scan time in L1B output record */
215  /* */
216  int16_t year = (int16_t) syear;
217  int16_t day = (int16_t) sday;
218 
219  if ((msec[recnum] < msec[recnum + 1]) && (recnum < nscan - 1)) {
220  /* adjust for day rollover */
221  day += 1;
222  if (day > (365 + (year % 4 == 0))) {
223  year += 1;
224  day = 1;
225  }
226  }
227 
228  l1rec->scantime = yds2unix(year, day, (double) (msec[recnum] / 1.e3));
229 
230  return (status);
231 }
232 
233 int closel1a_osmi(filehandle *file) {
234 
235  free(l1a_data);
236  free(l1b_buffer);
237  free(msec);
238  free(sc_ana);
239  free(sc_ttag);
240  free(sc_ttag_inc);
241 
242  SDend(file->sd_id);
243 
244  return (0);
245 }
246 
247 int32_t get_l1a_rec_osmi(int32_t sd_id, int32_t recno, cal_mod_struc *cal_mod,
248  float **l1b_data, int16_t **l2_flags) {
249  /* recno is 1-based */
250 
251  /* */
252  /* local vars */
253  /* */
254  int32_t ipix; /* pixel number */
255  int32_t idet; /* detector number */
256  int32_t irec;
257  int32_t start[3] = {0, 0, 0};
258  int32_t edges[3];
259 
260  int16_t band_buf[LAC_PIXEL_NUM * NBANDS_OSMI];
261 
262  static int32_t crec;
263  static int16_t n_read = 0;
264  static int32_t max_rec_in_rdbuf;
265  static int32_t offset;
266  static int32_t i, j;
267 
268  static byte first = 1;
269 
270  int16_t i16dum;
271 
272  int32_t recno1, recno2;
273  double uepoch, epoch, epoch1, epoch2;
274  double pos1[3], vel1[3];
275  double pos2[3], vel2[3];
276  double pos [3], vel [3];
277 
278  struct GEODETIC obs;
279  struct TOPOCENTRIC top;
280 
281  /* */
282  /* Read L1A data scan by scan, build L1B record, and write. */
283 
284  if (first) {
285  max_rec_in_rdbuf = recno - 1;
286  }
287 
288  crec = recno;
289 
290  if (crec > nscan) return 0;
291 
292 
293  if (crec > max_rec_in_rdbuf || crec < offset) {
294 
295  n_read = (nscan - recno >= NREC_IN_BUF) ?
296  NREC_IN_BUF : nscan - recno + 1;
297 
298 
299  edges[0] = 1;
300  start[1] = 0;
301  edges[1] = npix;
302  edges[2] = 8;
303  for (irec = 0; irec < n_read; irec++) {
304  start[0] = recno - 1 + irec;
305  SDreaddata(SDselect(sd_id, SDnametoindex(sd_id, "l1a_data")),
306  start, NULL, edges, (VOIDP) band_buf);
307  for (ipix = 0; ipix < npix; ipix++)
308  for (idet = 0; idet < 8; idet++)
309  l1a_data[irec * (NBANDS_OSMI * npix) + idet * npix + ipix] =
310  band_buf[ipix * NBANDS_OSMI + idet];
311  }
312 
313  offset = recno;
314  max_rec_in_rdbuf = offset + n_read - 1;
315 
316  }
317 
318  if ((crec - offset) >= 0) {
319  j = crec - offset;
320 
322  msec[recno - 1], dtype, nsta, npix,
323  (recno - 1) % 96, &i16dum, i16dum, 0,
325  cal_mod);
326 
327 
328  rdSDS(sd_id, "lat", recno - 1, 0, 1, npix, (VOIDP) ylat);
329  rdSDS(sd_id, "lon", recno - 1, 0, 1, npix, (VOIDP) xlon);
330 
331  /* Calculate sunz/a and senz/a */
332 
333  uepoch = itojul(700101, 0); /* unix time epoch */
334 
335  if (recno <= 96) {
336  recno1 = recno - 1;
337  recno2 = recno - 1;
338  } else {
339  recno1 = recno - 1; /* current scan */
340  recno2 = recno - 96 - 1; /* next scan in time */
341  }
342 
343  /* get time and ecef s/c pos/vel for bounding records */
344  epoch1 = uepoch + sc_ttag[recno1] / 86400.0;
345  pos1[0] = sc_ana[40 * (recno1) + 9] * 1000.0;
346  pos1[1] = sc_ana[40 * (recno1) + 10] * 1000.0;
347  pos1[2] = sc_ana[40 * (recno1) + 11] * 1000.0;
348  vel1[0] = sc_ana[40 * (recno1) + 12] * 1000.0;
349  vel1[1] = sc_ana[40 * (recno1) + 13] * 1000.0;
350  vel1[2] = sc_ana[40 * (recno1) + 14] * 1000.0;
351 
352  epoch2 = uepoch + sc_ttag[recno2] / 86400.0;
353  pos2[0] = sc_ana[40 * (recno2) + 9] * 1000.0;
354  pos2[1] = sc_ana[40 * (recno2) + 10] * 1000.0;
355  pos2[2] = sc_ana[40 * (recno2) + 11] * 1000.0;
356  vel2[0] = sc_ana[40 * (recno2) + 12] * 1000.0;
357  vel2[1] = sc_ana[40 * (recno2) + 13] * 1000.0;
358  vel2[2] = sc_ana[40 * (recno2) + 14] * 1000.0;
359 
360  for (i = 0; i < npix; i++) {
361 
362  epoch = uepoch + (sc_ttag[recno - 1] +
363  /*
364  ((npix - 1 - i) * sc_ttag_inc[recno-1]))/86400.0;
365  */
366  ((i) * sc_ttag_inc[recno - 1])) / 86400.0;
367 
368  obs.lat = DTOR(ylat[i]);
369  obs.lon = DTOR(xlon[i]);
370  obs.height = 0;
371 
372  sunpos(tconv(epoch - AU / (VLIGHT * 86400.0), "utc:tdt", NULL), pos, NULL);
373  ctotc(pos, &obs, gmha(epoch), &top, NULL);
374  solz[i] = (float) (90.0 - RTOD(top.el));
375  sola[i] = (float) RTOD(top.az);
376  if (sola[i] < 0.0) sola[i] = sola[i] + 360.0;
377 
378  intpos_(&epoch1, pos1, vel1, &epoch2, pos2, vel2, &epoch, pos, vel);
379 
380  /*
381  printf("%d %lf %lf %lf %lf\n",i,epoch1,pos1[0],pos1[1],pos1[2]);
382  printf("%d %lf %lf %lf %lf\n",i,epoch ,pos [0],pos [1],pos [2]);
383  printf("%d %lf %lf %lf %lf\n",i,epoch2,pos2[0],pos2[1],pos2[2]);
384  */
385 
386  ctotc(pos, &obs, gmha(epoch), &top, NULL);
387  senz[i] = (float) (90.0 - RTOD(top.el));
388  sena[i] = (float) RTOD(top.az);
389  if (sena[i] < 0.0) sena[i] = sena[i] + 360.0;
390 
391  }
392 
393  }
394 
395  *l1b_data = l1b_buffer;
396  *l2_flags = l2_flags_buffer;
397  first = 0;
398 
399  return (0);
400 }
401 
int32_t get_l1a_rec_osmi(int32_t, int32_t, cal_mod_struc *, float **, int16_t **)
Definition: l1a_osmi.c:247
const int bindx[3]
Definition: DbLutNetcdf.cpp:28
int16 eday
Definition: l1_czcs_hdf.c:17
int j
Definition: decode_rs.h:73
int rdSDS(int32_t fileID, const char sdsname[], int32_t start1, int32_t start2, int32_t edges1, int32_t edges2, void *array_data)
int32_t day
double t_const[BANDS_DIMS_1A]
Definition: l1a_seawifs.c:49
int status
Definition: l1_czcs_hdf.c:32
double yds2unix(int16_t year, int16_t day, double secs)
Definition: yds2unix.c:7
double height
Definition: orbit.h:119
int16_t fileID
char cal_path_tab[128]
Definition: l1a_seawifs.c:70
#define OSMI
Definition: sensorDefs.h:16
#define NULL
Definition: decode_rs.h:63
void ctotc(r, struct GEODETIC *obs, double gha, struct TOPOCENTRIC *top, daer)
Definition: ctotc.c:43
void intpos_(double *epoch1, double pos1[3], double vel1[3], double *epoch2, double pos2[3], double vel2[3], double *epoch, double pos[3], double vel[3])
read l1rec
#define STRAYLIGHT
Definition: l2_flags.h:19
int readl1a_osmi(filehandle *file, int32_t recnum, l1str *l1rec)
Definition: l1a_osmi.c:154
double itojul(int32_t date, int32_t time)
Definition: time-utils.c:100
int16_t * l1a_data
Definition: l1a_seawifs.c:84
const int chindx[8]
Definition: DbLutNetcdf.cpp:27
#define NBANDS_OSMI
Definition: l1a_osmi.c:16
float ylat[LAC_PIXEL_NUM]
Definition: l1a_seawifs.c:92
float32 * pos
Definition: l1_czcs_hdf.c:35
int32 * msec
Definition: l1_czcs_hdf.c:31
int16 eyear
Definition: l1_czcs_hdf.c:17
int syear
Definition: l1_czcs_hdf.c:15
short l2_flags_buffer[LAC_PIXEL_NUM]
Definition: l1a_seawifs.c:79
int32 nscan
Definition: l1_czcs_hdf.c:19
int32 smsec
Definition: l1_czcs_hdf.c:16
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:15
int32 nsta
Definition: l1_czcs_hdf.c:27
int32_t get_cal_osmi(char *cal_path, int16_t syear, int16_t sday, int16_t eday, int32_t msec, int16_t *cal_year, int16_t *cal_day, int32_t *cal_msec, float *eoffset, float *egain, float *temp, float *mirror, float *t_const, float *t_linear, float *t_quadratic, float *slope, float *dc, float *sm)
Definition: get_cal_osmi.c:99
int openl1a_osmi(filehandle *file)
Definition: l1a_osmi.c:70
#define DTOR(X)
Definition: orbit.h:51
read recnum
#define BANDS_DIMS_1A
Definition: level_1a_index.h:8
float32 slope[]
Definition: l2lists.h:30
int32 ninc
Definition: l1_czcs_hdf.c:28
double gmha(double tjd)
Definition: gmha.c:61
#define VLIGHT
Definition: lunsol.h:23
l1_input_t * l1_input
Definition: l1_options.c:9
int getHDFattr(int32_t fileID, const char attrname[], const char sdsname[], void *data)
double lon
Definition: orbit.h:118
int closel1a_osmi(filehandle *file)
Definition: l1a_osmi.c:233
#define HILT
Definition: l2_flags.h:15
float xlon[LAC_PIXEL_NUM]
Definition: l1a_seawifs.c:93
double tconv(double tjd, char conv[], char **errstr)
int32 emsec
Definition: l1_czcs_hdf.c:18
float * l1b_buffer
Definition: l1a_seawifs.c:86
int32_t calibrate_l1a_osmi(char *cal_path, int16_t syear, int16_t sday, int32_t smsec, int16_t eday, int32_t msec, char *dtype, int32_t st_samp, int32_t nsamp, int32_t fpixel, int16_t gain[4], int16_t offset, int16_t scan_temp, int16_t *l1a_data, float *l1b_data, cal_mod_struc *cal_mod)
#define LAC_PIXEL_NUM
Definition: l1a_osmi.c:14
dtype
Definition: DDataset.hpp:31
#define AU
Definition: lunsol.h:19
l2prod offset
#define NREC_IN_BUF
Definition: l1a_osmi.c:15
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")
double lat
Definition: orbit.h:117
int sunpos(double tjd, double r[3], char **errstr)
int npix
Definition: get_cmp.c:27
#define RTOD(X)
Definition: orbit.h:54