NASA Logo
Ocean Color Science Software

ocssw V2022
l1_viirs_l1b.c
Go to the documentation of this file.
1 
7 #include "l1_viirs_l1b.h"
8 #include "l1.h"
9 
10 #include "nc4utils.h"
11 #include <libnav.h>
12 #include <float.h>
13 #include "calibrate_viirs.h"
14 #include <math.h>
15 
16 typedef unsigned short ushort;
17 
18 /***********************************************************************
19  Calibration
20  ***********************************************************************/
21 static double ***f_cal_corr = NULL; /* f table correction [band][det][ms] */
22 
23 /***********************************************************************
24  Initialization
25  ***********************************************************************/
26 static float *Fobar; // reflectance to radiance conversion factors
27 static int extract_pixel_start = 0;
28 static int extract_pixel_stop = 0;
29 
30 /* Variables to load for each input file and group */
31 typedef struct {
32  size_t index;
33  char *name;
34 } varlist;
35 
36 /*----- GEO file -----*/
37 
38 /* Scan Line attributes, to be read on initialization.
39  All values are dimensioned as [number_of_scans]. */
40 static const char* SCN_GRP = "scan_line_attributes";
41 
42 enum scn_var {
43  SCN_STIME, /*< double scan_start_time */
44  //SCN_ETIME, /*< double scan_end_time */
45  //SCN_MTIME, /*< double ev_mid_time */
46  SCN_MSIDE, /*< ubyte HAM_side */
47  SCN_MODE, /*< ubyte sensor_mode */
48  SCN_QUAL, /*< short scan_quality */
50 };
51 static const varlist VARLIST_SCN[] = {
52  { SCN_STIME, "scan_start_time"},
53  //{ SCN_ETIME, "scan_end_time" },
54  //{ SCN_MTIME, "ev_mid_time" },
55  { SCN_MSIDE, "HAM_side"},
56  { SCN_MODE, "sensor_mode"},
57  { SCN_QUAL, "scan_quality"}
58 };
59 static var_str_nc *scn[NVARS_SCN];
60 
61 /* Navigation data, to be read on initialization.
62  All values are float; dimensioned [number_of_scans,3] except as noted. */
63 static const char* NAV_GRP = "navigation_data";
64 
65 enum nav_var {
66  //NAV_QUAT, /*< att_quat_ev [number_of_scans, quaternion_elements] */
67  NAV_ANG, /*< att_ang */
68  NAV_POS, /*< orb_pos_ev */
69  NAV_VEL, /*< orb_vel_ev */
70  //NAV_SOLJ, /*< solar_j2000 */
71  //NAV_SOLV, /*< solar_inst */
72  //NAV_SOLD, /*< earth_sun_distance [number_of_scans] */
73  //NAV_LUNJ, /*< lunar_j2000 */
74  //NAV_LUNV, /*< lunar_inst */
75  //NAV_LUND, /*< earth_moon_distance [number_of_scans] */
77 };
78 static const varlist VARLIST_NAV[] = {
79  //{ NAV_QUAT, "att_quat_ev" },
80  { NAV_ANG, "att_ang"}, // (roll, pitch, yaw: degrees)
81  { NAV_POS, "orb_pos_ev"}, // (ECR: meters)
82  { NAV_VEL, "orb_vel_ev"}, // (ECR: meters/second)
83  //{ NAV_SOLJ, "solar_j2000"},
84  //{ NAV_SOLV, "solar_inst"},
85  //{ NAV_SOLD, "earth_sun_distance"}, // AU
86  //{ NAV_LUNJ, "lunar_j2000" },
87  //{ NAV_LUNV, "lunar_inst" },
88  //{ NAV_LUND, "earth_moon_distance" }
89 };
90 static var_str_nc *nav[NVARS_NAV]; /* read on initialization */
91 
92 /* Geolocation data, to be read line-by-line.
93  All values are dimensioned as [number_of_lines, number_of_pixels]. */
94 static const char* GEO_GRP = "geolocation_data";
95 
96 enum geo_var {
97  GEO_LAT, /*< float latitude */
98  GEO_LON, /*< float longitude */
99  GEO_HGT, /*< short height (-> float) */
100  //GEO_RNG, /*< short range (-> float) */
101  GEO_SENA, /*< short sensor_azimuth (-> float) */
102  GEO_SENZ, /*< short sensor_zenith (-> float) */
103  GEO_SOLA, /*< short solar_azimuth (-> float) */
104  GEO_SOLZ, /*< short solar_zenith (-> float) */
105  //GEO_MASK, /*< ubyte land_water_mask */
106  GEO_QUAL, /*< ubyte quality_flag */
108 };
109 static const varlist VARLIST_GEO[] = {
110  { GEO_LAT, "latitude"},
111  { GEO_LON, "longitude"},
112  { GEO_HGT, "height"},
113  //{ GEO_RNG, "range" },
114  { GEO_SENA, "sensor_azimuth"},
115  { GEO_SENZ, "sensor_zenith"},
116  { GEO_SOLA, "solar_azimuth"},
117  { GEO_SOLZ, "solar_zenith"},
118  //{ GEO_MASK, "land_water_mask" },
119  { GEO_QUAL, "quality_flag"}
120 };
121 static var_str_nc *geo[NVARS_GEO];
122 
123 /*----- L1B file -----*/
124 
125 /* Scan Line attributes, to be read on initialization.
126  All values are dimensioned as [number_of_scans]. */
128  L1BSCN_STIME, /*< double scan_start_time */
129  //L1BSCN_ETIME, /*< double scan_end_time */
130  //L1BSCN_MTIME, /*< double ev_mid_time */
131  L1BSCN_MODE, /*< ubyte scan_state_flags */
132  L1BSCN_QUAL, /*< ubyte scan_quality_flags */
134 };
135 static const varlist VARLIST_L1BSCN[] = {
136  { L1BSCN_STIME, "scan_start_time"},
137  //{ L1BSCN_ETIME, "scan_end_time" },
138  //{ L1BSCN_MTIME, "ev_mid_time" },
139  { L1BSCN_MODE, "scan_state_flags"},
140  { L1BSCN_QUAL, "scan_quality_flags"}
141 };
142 static var_str_nc *l1bscn[NVARS_L1BSCN];
143 
144 /* Radiometric data, to be read line-by-line.
145  All values are dimensioned as [number_of_lines, number_of_pixels]. */
146 static const char* L1B_GRP = "observation_data";
147 
148 enum bandtypes {
150 };
151 
152 static const varlist VARLIST_L1B[] = {
153  /*** Reflective Solar Bands (aka VIS/SWIR) ***/
154  { RSB, "M01"}, // 410 nm
155  { RSB, "M02"}, // 443 nm
156  { RSB, "M03"}, // 486 nm (blue)
157  { RSB, "M04"}, // 551 nm (green)
158  { RSB, "M05"}, // 671 nm (red)
159  { RSB, "M06"}, // 745 nm
160  { RSB, "M07"}, // 862 nm
161  { RSB, "M08"}, // 1238 nm
162  { CIR, "M09"}, // 1378 nm (cirrus)
163  { RSB, "M10"}, // 1601 nm
164  { RSB, "M11"}, // 2257 nm
165  /*** Thermal Emissive Bands ***/
166  { TEB, "M12"}, // 3700 nm
167  { TEB, "M13"}, // 4050 nm
168  { TEB, "M14"}, // 8550 nm
169  { TEB, "M15"}, // 10763 nm
170  { TEB, "M16"} // 12013 nm
171 };
172 #define MAXBANDS 16
173 static var_str_nc *l1b[MAXBANDS]; /* ushort -> float */
174 static var_str_nc *l1bq[MAXBANDS]; /* ubyte; name = band+"_quality_flags" */
175 
176 static float latGeoFillValue = -999.9;
177 static float lonGeoFillValue = -999.9;
178 static short senzGeoFillValue = -32768;
179 static short senaGeoFillValue = -32768;
180 static short solzGeoFillValue = -32768;
181 static short solaGeoFillValue = -32768;
182 
183 static float latL2FillValue = -999.0;
184 static float lonL2FillValue = -999.0;
185 static float senzL2FillValue = -32767;
186 static float senaL2FillValue = -32767;
187 static float solzL2FillValue = -32767;
188 static float solaL2FillValue = -32767;
189 
190 /*---------------------------------------------------------------------*/
191 
192 /* Info for any VIIRS input file */
193 
194 typedef struct {
195  int32_t id; /*< NetCDF4 file ID */
196  char file[FILENAME_MAX]; /*< file path */
197  char title[FILENAME_MAX];
198  size_t nscans;
199  size_t nlines;
200  size_t npixls;
202  char start_time[25];
203  char end_time[25];
205 } viirs_file;
206 
207 void print_viirs_file(const viirs_file info) {
208  printf("file:\t%s\n", info.file);
209  printf("title:\t%s\n", info.title);
210  printf("time range:\t%s %s\n", info.start_time, info.end_time);
211  printf("orbit: \t%d\n", info.orbit_number);
212  printf("nscans =\t%d\n", (int) info.nscans);
213  printf("nlines =\t%d\n", (int) info.nlines);
214  printf("npixls =\t%d\n", (int) info.npixls);
215 }
216 
217 int init_viirs_file(const char filename[FILENAME_MAX],
218  viirs_file *info) {
219  int32_t grpid, varid, dimid;
220  size_t attlen; // text attribute length
221  int status;
222 
223  /* open file */
224  bzero(info, sizeof(*info));
225  strcpy(info->file, filename);
226  status = nc_open(info->file, NC_NOWRITE, &info->id);
227  if(status != NC_NOERR) {
228  printf("-E- init_viirs_file - Could not open netCDF file \"%s\"\n", info->file);
229  exit(EXIT_FAILURE);
230  }
231 
232 
233  /* load dimensions */
234  nc_inq_dimid(info->id, "number_of_scans", &dimid);
235  nc_inq_dimlen(info->id, dimid, &info->nscans);
236  nc_inq_dimid(info->id, "number_of_lines", &dimid);
237  nc_inq_dimlen(info->id, dimid, &info->nlines);
238  nc_inq_dimid(info->id, "number_of_pixels", &dimid);
239  nc_inq_dimlen(info->id, dimid, &info->npixls);
240 
241  /* load global attributes */
242  nc_get_att_text(info->id, NC_GLOBAL,
243  "title", (char*) &info->title);
244  nc_inq_attlen(info->id, NC_GLOBAL, "title", &attlen);
245  info->title[attlen] = '\0'; // null terminate
246 
247  nc_get_att_text(info->id, NC_GLOBAL,
248  "time_coverage_start", (char*) &info->start_time);
249  nc_inq_attlen(info->id, NC_GLOBAL, "time_coverage_start", &attlen);
250  info->start_time[attlen] = '\0';
251 
252  nc_get_att_text(info->id, NC_GLOBAL,
253  "time_coverage_end", (char*) &info->end_time);
254  nc_inq_attlen(info->id, NC_GLOBAL, "time_coverage_end", &attlen);
255  info->end_time[attlen] = '\0';
256 
257  // yeah, borrowed this from the example code on UCAR's website...
258  nc_type vr_type; /* attribute type */
259  size_t vr_len; /* attribute length */
260  if ((nc_inq_att(info->id, NC_GLOBAL, "orbit_number", &vr_type, &vr_len) == NC_NOERR)) {
261  nc_get_att_int(info->id, NC_GLOBAL,
262  "orbit_number", &info->orbit_number);
263  } else {
264  nc_get_att_int(info->id, NC_GLOBAL,
265  "OrbitNumber", &info->orbit_number);
266  }
267 
268  if ((nc_inq_att(info->id, NC_GLOBAL, "extract_pixel_start", &vr_type, &vr_len) == NC_NOERR)) {
269  nc_get_att_int(info->id, NC_GLOBAL,
270  "extract_pixel_start", &extract_pixel_start);
271  extract_pixel_start--; // Attribute is one-based
272  nc_get_att_int(info->id, NC_GLOBAL,
273  "extract_pixel_stop", &extract_pixel_stop);
274  extract_pixel_stop--; // Attribute is one-based
275  }
276 
277  /* load start time for each scan */
278  nc_inq_grp_ncid(info->id, "scan_line_attributes", &grpid);
279  nc_inq_varid(grpid, "scan_start_time", &varid);
280  TRYMEM(__FILE__, __LINE__,
281  (info->scan_start_time = malloc(info->nscans *
282  sizeof(*info->scan_start_time))));
283  nc_get_var_double(grpid, varid, info->scan_start_time);
284 
285  return SUCCESS;
286 }
287 
288 /*---------------------------------------------------------------------*/
289 
290 /* Info for VIIRS L1B and GEO files */
291 
293  size_t i;
294  const char* title =
295  "VIIRS M-band Reflected Solar Band and Thermal Emissive Band Data";
296  grp_str_nc file = { 0 };
297  char qaname[FILENAME_MAX + 1];
298 
299  /* Verify file type from global attributes */
300  if (strcmp(l1binfo.title, title) != 0) {
301  fprintf(stderr,
302  "Input file %s has unexpected product title!\n"
303  "Expected:\t\"%s\"\n"
304  "Actual: \t\"%s\"\n",
305  l1binfo.file, title, l1binfo.title);
306  return -1;
307  }
308 
309  /* Determine file contents */
310  file.id = l1binfo.id;
311  load_grp_nc(&file);
312 
313  /* Populate array holding info and data for L1BSCN variables */
314  for (i = 0; i < NVARS_L1BSCN; i++) {
315  l1bscn[i] = find_var_byname_nc(file, VARLIST_L1BSCN[i].name, SCN_GRP);
316  readall_var(l1bscn[i]);
317  }
318 
319  /* For each L1B variable, */
320  for (i = 0; i < MAXBANDS; i++) {
321 
322  /* Populate array holding info for L1B variables */
323  l1b[i] = find_var_byname_nc(file, VARLIST_L1B[i].name, L1B_GRP);
324  // returns NULL if var not found.
325 
326  /* also load QA flag info (unscaled) */
327  sprintf(qaname, "%s%s", VARLIST_L1B[i].name, "_quality_flags");
328  l1bq[i] = find_var_byname_nc(file, qaname, L1B_GRP);
329 
330  }
331 
332  /* TO DO: free unused memory allocated to grp_str_nc file. */
333  return SUCCESS;
334 }
335 
337  size_t i;
338  const char* title = "VIIRS M-band Geolocation Data";
339  grp_str_nc file = { 0 };
340 
341  /* Verify file type from global attributes */
342  if (strcmp(geoinfo.title, title) != 0) {
343  fprintf(stderr,
344  "Input file %s has unexpected product title!\n"
345  "Expected:\t\"%s\"\n"
346  "Actual: \t\"%s\"\n",
347  geoinfo.file, title, geoinfo.title);
348  return -1;
349  }
350 
351  /* Determine file contents */
352  file.id = geoinfo.id;
353  load_grp_nc(&file);
354 
355  /* Populate array holding info for GEO variables */
356  for (i = 0; i < NVARS_GEO; i++) {
357  geo[i] = find_var_byname_nc(file, VARLIST_GEO[i].name, GEO_GRP);
358  }
359 
360  TRY_NC(__FILE__, __LINE__,
361  nc_get_att_float(geo[GEO_LAT]->grpid, geo[GEO_LAT]->id, "_FillValue", &latGeoFillValue));
362  TRY_NC(__FILE__, __LINE__,
363  nc_get_att_float(geo[GEO_LON]->grpid, geo[GEO_LON]->id, "_FillValue", &lonGeoFillValue));
364  TRY_NC(__FILE__, __LINE__,
365  nc_get_att_short(geo[GEO_SENZ]->grpid, geo[GEO_SENZ]->id, "_FillValue", &senzGeoFillValue));
366  TRY_NC(__FILE__, __LINE__,
367  nc_get_att_short(geo[GEO_SENA]->grpid, geo[GEO_SENA]->id, "_FillValue", &senaGeoFillValue));
368  TRY_NC(__FILE__, __LINE__,
369  nc_get_att_short(geo[GEO_SOLZ]->grpid, geo[GEO_SOLZ]->id, "_FillValue", &solzGeoFillValue));
370  TRY_NC(__FILE__, __LINE__,
371  nc_get_att_short(geo[GEO_SOLA]->grpid, geo[GEO_SOLA]->id, "_FillValue", &solaGeoFillValue));
372 
373  // grab the fill values for L2 products
374  int status;
375  productInfo_t* info = allocateProductInfo();
376  status = findProductInfo("lat", VIIRSN, info);
377  if (status)
378  latL2FillValue = info->fillValue;
379  status = findProductInfo("lon", VIIRSN, info);
380  if (status)
381  lonL2FillValue = info->fillValue;
382  status = findProductInfo("sena", VIIRSN, info);
383  if (status)
384  senaL2FillValue = info->fillValue;
385  status = findProductInfo("senz", VIIRSN, info);
386  if (status)
387  senzL2FillValue = info->fillValue;
388  status = findProductInfo("sola", VIIRSN, info);
389  if (status)
390  solaL2FillValue = info->fillValue;
391  status = findProductInfo("solz", VIIRSN, info);
392  if (status)
393  solzL2FillValue = info->fillValue;
394  freeProductInfo(info);
395 
396  /* Populate array holding info and data for SCN variables */
397  for (i = 0; i < NVARS_SCN; i++) {
398  scn[i] = find_var_byname_nc(file, VARLIST_SCN[i].name, SCN_GRP);
399  readall_var(scn[i]);
400  }
401 
402  /* Populate array holding info and data for NAV variables */
403  for (i = 0; i < NVARS_NAV; i++) {
404  nav[i] = find_var_byname_nc(file, VARLIST_NAV[i].name, NAV_GRP);
405  if (!nav[i]) {
406  char *tmpStr = malloc(strlen(VARLIST_NAV[i].name) + 5);
407  sprintf(tmpStr, "%s_mid", VARLIST_NAV[i].name);
408  nav[i] = find_var_byname_nc(file, tmpStr, NAV_GRP);
409  if (!nav[i]) {
410  printf("-E- init_viirs_geofile - could not find variable %s in the geo file\n", VARLIST_NAV[i].name);
411  exit(EXIT_FAILURE);
412  }
413  free(tmpStr);
414  }
415  readall_var(nav[i]);
416  }
417 
418  /* TO DO: free unused memory allocated to grp_str_nc file. */
419  return SUCCESS;
420 }
421 
422 /*---------------------------------------------------------------------*/
423 
424 /* Open and validate L1B and GEO files */
425 
426 int openl1_viirs_l1b(filehandle *l1file) {
427  viirs_file l1binfo, geoinfo;
428  int iscan;
429 
430  // make sure a GEO file was passed in
431  if(!l1file->geofile || !l1file->geofile[0]) {
432  printf("-E- openl1b_viirs_nc - VIIRS L1B files require a GEO file.\n");
433  exit(EXIT_FAILURE);
434  }
435 
436  /*----- Initialize L1B and GEO info -----*/
437  init_viirs_file(l1file->name, &l1binfo);
438  init_viirs_file(l1file->geofile, &geoinfo);
439 
440  /*----- Populate filehandle structure -----*/
441  l1file->ndets = l1binfo.nlines / l1binfo.nscans;
442  l1file->nscan = l1binfo.nlines;
443  l1file->npix = l1binfo.npixls;
444  l1file->terrain_corrected = 1; // presumed.
445  l1file->orbit_number = l1binfo.orbit_number;
446  strcpy(l1file->spatialResolution, "750 m");
447 
448  /*----- Check that input files are compatible -----*/
449 
450  /* dimensions */
451  if (
452  (l1binfo.nscans != geoinfo.nscans) ||
453  (l1binfo.nlines != geoinfo.nlines) ||
454  (l1binfo.npixls != geoinfo.npixls)) {
455  fprintf(stderr, "Geometry mismatch!\n");
456  fprintf(stderr, "L1B: nscans = %3d, nlines = %4d, npixls = %4d\n",
457  (int) l1binfo.nscans, (int) l1binfo.nlines, (int) l1binfo.npixls);
458  fprintf(stderr, "GEO: nscans = %3d, nlines = %4d, npixls = %4d\n",
459  (int) geoinfo.nscans, (int) geoinfo.nlines, (int) geoinfo.npixls);
460  return -1;
461  }
462 
463  /* time and orbit */
464  if (
465  (l1binfo.orbit_number != geoinfo.orbit_number) ||
466  (strcmp(l1binfo.start_time, geoinfo.start_time) != 0) ||
467  (strcmp(l1binfo.end_time, geoinfo.end_time) != 0)) {
468  fprintf(stderr, "Time coverage mismatch!\n");
469  fprintf(stderr, "L1B: Orbit %6d, %s - %s\n",
470  l1binfo.orbit_number, l1binfo.start_time, l1binfo.end_time);
471  fprintf(stderr, "GEO: Orbit %6d, %s - %s\n",
472  geoinfo.orbit_number, geoinfo.start_time, geoinfo.end_time);
473  return -1;
474  }
475 
476  /* scan start times */
477  for (iscan = 0; iscan < l1binfo.nscans; iscan++)
478  if (fabs(l1binfo.scan_start_time[iscan] -
479  geoinfo.scan_start_time[iscan]) > DBL_EPSILON) {
480  fprintf(stderr, "Scan time mismatch!\n");
481  fprintf(stderr, "Scan %d\tL1B=%f\tGEO=%f\tdiff=%f\n",
482  iscan,
483  l1binfo.scan_start_time[iscan],
484  geoinfo.scan_start_time[iscan],
485  l1binfo.scan_start_time[iscan] -
486  geoinfo.scan_start_time[iscan]);
487  return -1;
488  }
489 
490  /*----- Load L1B and GEO info into global variables -----*/
491  init_viirs_l1bfile(l1binfo); // l1b, l1bq
492  init_viirs_geofile(geoinfo); // scn, nav, geo
493 
494  /*----- Ancillary Data -----*/
495  // TODO: remove this?
496  rdsensorinfo(l1file->sensorID, l1_input->evalmask,
497  "Fobar", (void **) &Fobar);
498 
499  /*----- Calibration LUT -----*/
500  if (l1_input->calfile[0]) {
501  double grantime = l1binfo.scan_start_time[0]; // granule start time
502  grantime -= leapseconds_since_1993(grantime); // convert TAI93 to "UTC93"
503  grantime += 1104537600.0; // convert "UTC93" to UTC (1958) (verify)
504  grantime *= 1000000.0; // convert to IET
505  load_fcal_lut(l1_input->calfile, (int64_t) grantime, &f_cal_corr);
506  }
507 
508  free(l1binfo.scan_start_time); // free memory allocated
509  free(geoinfo.scan_start_time); // by init_viirs_file()
510  return SUCCESS;
511 }
512 
513 /***********************************************************************
514  Read Lines
515  ***********************************************************************/
516 
517 int read_var_1line(var_str_nc *var, size_t iline) {
518  size_t typesize;
519  size_t npixl = 0;
520  size_t start[2] = {0, 0}; // assume dims [nline,npixl]
521  size_t count[2] = {1, 1};
522 
523  if (var->ndims == 2) {
524 
525  /* determine # bytes to read */
526  npixl = var->dim[1].len;
527  TRY_NC(__FILE__, __LINE__,
528  nc_inq_type(var->id, var->type, NULL, &typesize));
529 
530  /* allocate one line as needed (first call) */
531  if (var->data == NULL) {
532  TRYMEM(__FILE__, __LINE__,
533  (var->data = calloc(npixl, typesize)));
534  var->dim[0].len = 1; // to represent data stored
535  }
536  /* read data */
537  start[0] = iline; // read one line
538  count[1] = npixl; // read all pixels
539  TRY_NC(__FILE__, __LINE__,
540  nc_get_vara(var->grpid, var->id, start, count, var->data));
541  }
542  return npixl;
543 }
544 
545 int scale_short(var_str_nc *var, float* dest) {
546  /* cheating here: assume scaling is always short to float. */
547 
548  size_t ipix;
549  size_t npix = var->dim[1].len;
550  short tmpval, fillval;
551  short minval, maxval;
552  float scale, offset;
553 
554  /* load scaling factors */
555  TRY_NC(__FILE__, __LINE__,
556  nc_get_att_short(var->grpid, var->id, "_FillValue", &fillval));
557  TRY_NC(__FILE__, __LINE__,
558  nc_get_att_short(var->grpid, var->id, "valid_min", &minval));
559  TRY_NC(__FILE__, __LINE__,
560  nc_get_att_short(var->grpid, var->id, "valid_max", &maxval));
561  TRY_NC(__FILE__, __LINE__,
562  nc_get_att_float(var->grpid, var->id, "scale_factor", &scale));
563  TRY_NC(__FILE__, __LINE__,
564  nc_get_att_float(var->grpid, var->id, "add_offset", &offset));
565 
566  for (ipix = 0; ipix < npix; ipix++) {
567  tmpval = ((short *) var->data)[ipix];
568  if ((tmpval == fillval) ||
569  (tmpval < minval) ||
570  (tmpval > maxval))
571  dest[ipix] = (float) tmpval;
572  else
573  dest[ipix] = scale * (float) tmpval + offset;
574  }
575 
576  return SUCCESS;
577 }
578 
579 int scale_l1bvals(l1str *l1rec) {
580 
581  uint32_t tmpval, fillval;
582  uint32_t minval, maxval;
583  float scale, offset;
584  double f_corr;
585 
586  size_t ipb, ipix, npix;
587  size_t iband;
588  size_t irsb = 0;
589  size_t iteb = 0;
590  enum bandtypes bandtype;
591  //char flag; // ubyte
592 
593  /* Initializations */
594  npix = l1rec->npix;
595 
596  /* Step through all bands */
597  for (iband = 0; iband < MAXBANDS; iband++) {
598  bandtype = VARLIST_L1B[iband].index;
599 
600  if (l1b[iband] == NULL) {
601  if (bandtype == TEB) iteb++;
602  if (bandtype == RSB) irsb++;
603  continue; // skip rest of processing for missing band
604  }
605 
606  TRY_NC(__FILE__, __LINE__,
607  nc_get_att_uint(l1b[iband]->grpid, l1b[iband]->id,
608  "_FillValue", &fillval));
609  TRY_NC(__FILE__, __LINE__,
610  nc_get_att_uint(l1b[iband]->grpid, l1b[iband]->id,
611  "valid_min", &minval));
612  TRY_NC(__FILE__, __LINE__,
613  nc_get_att_uint(l1b[iband]->grpid, l1b[iband]->id,
614  "valid_max", &maxval));
615  TRY_NC(__FILE__, __LINE__,
616  nc_get_att_float(l1b[iband]->grpid, l1b[iband]->id,
617  "scale_factor", &scale));
618  TRY_NC(__FILE__, __LINE__,
619  nc_get_att_float(l1b[iband]->grpid, l1b[iband]->id,
620  "add_offset", &offset));
621 
622  /* get specific f table cal correction */
623  f_corr = (f_cal_corr == NULL) ? 1.0
624  : f_cal_corr[iband][l1rec->detnum][l1rec->mside];
625 
626  /*** Thermal Emissive Bands ***/
627  if (bandtype == TEB) {
628  for (ipix = 0; ipix < npix; ipix++) {
629  //ipb = ipix * l1rec->l1file->nbandsir + iteb; // band varies fastest
630  ipb = ipix * NBANDSIR + iteb; // band varies fastest
631  l1rec->Ltir[ipb] = 0; // init to fill value
632 
633  /* skip out-of-bounds values */
634  if (strcmp(VARLIST_L1B[iband].name, "M13") == 0) {
635  tmpval = ((uint32_t *) l1b[iband]->data)[ipix];
636  } else {
637  tmpval = ((ushort *) l1b[iband]->data)[ipix];
638  }
639  if ((tmpval == fillval) ||
640  (tmpval < minval) ||
641  (tmpval > maxval))
642  continue;
643 
644  /* apply radiance scaling */
645  l1rec->Ltir[ipb] = scale * (float) tmpval + offset;
646 
647  /* Convert radiance from W/m^2/um/sr to mW/cm^2/um/sr */
648  l1rec->Ltir[ipb] /= 10.0;
649 
650  /* Apply F-factor */
651  l1rec->Ltir[ipb] *= f_corr;
652 
653  } // ipix
654  iteb++;
655  } // TEB
656 
657  /*** Cirrus Band ***/
658  else if (bandtype == CIR) {
659  for (ipix = 0; ipix < npix; ipix++) {
660  l1rec->rho_cirrus[ipix] = BAD_FLT; // init to fill value
661 
662  /* skip out-of-bounds values */
663  tmpval = ((ushort *) l1b[iband]->data)[ipix];
664  if ((tmpval == fillval) ||
665  (tmpval < minval) ||
666  (tmpval > maxval))
667  continue;
668 
669  /* apply reflectance scaling */
670  l1rec->rho_cirrus[ipix] = scale * (float) tmpval + offset;
671 
672  /* Normalize reflectance by solar zenith angle */
673  l1rec->rho_cirrus[ipix] /= cos(l1rec->solz[ipix] / RADEG);
674 
675  /* Apply F-factor */
676  l1rec->rho_cirrus[ipix] *= f_corr;
677 
678  } // ipix
679  } // CIR
680 
681  /*** Reflective Solar Bands ***/
682  else { // if (bandtype == RSB)
683 
684  l1rec->Fo[irsb] = Fobar[irsb] * l1rec->fsol;
685 
686  for (ipix = 0; ipix < npix; ipix++) {
687  ipb = ipix * l1rec->l1file->nbands + irsb; // band varies fastest
688  l1rec->Lt[ipb] = BAD_FLT; // init to fill value
689 
690  /* skip night data for visible bands */
691  if (l1rec->solz[ipix] > SOLZNIGHT)
692  continue;
693 
694  /* skip out-of-bounds values */
695  tmpval = ((ushort *) l1b[iband]->data)[ipix];
696  if ((tmpval == fillval) ||
697  (tmpval < minval) ||
698  (tmpval > maxval))
699  continue;
700 
701  /* apply reflectance scaling */
702  l1rec->Lt[ipb] = scale * (float) tmpval + offset;
703 
704  /* convert from reflectance to radiance */
705  l1rec->Lt[ipb] *= l1rec->Fo[irsb] / PI;
706 
707  /* Apply F-factor */
708  l1rec->Lt[ipb] *= f_corr;
709 
710  /* TO DO: flag any suspect values */
711  /* TO DO: flag hilt */
712  //flag = ((char*)l1bq[iband]->data)[ipix];
713  /* if ((flag & 4) && (irsb < 13)) ; */
714 
715  } // ipix
716  irsb++;
717  } // RSB
718 
719  } // iband
720 
721  return SUCCESS;
722 }
723 
724 int readl1_viirs_l1b(filehandle *l1file,
725  const int32_t iline,
726  l1str *l1rec,
727  int lonlat) {
728  float nbytes = l1file->npix * sizeof(float);
729  size_t ivar;
730  size_t ipix;
731  size_t iscan;
732  size_t i;
733 
734  float ang[3]; // degrees
735  float pos[3]; // km
736  float vel[3]; // km/sec
737  float sen_mat[3][3], coeff[10]; // for ocorient & mnorm
738  double mnorm[3];
739 
740  /*----- Basic info -----*/
741  // l1rec->sensorID = l1file->sensorID;
742  l1rec->npix = l1file->npix;
743  l1rec->detnum = iline % l1file->ndets;
744 
745  /*----- Scan-line values -----*/
746  iscan = iline / l1file->ndets;
747 
748  /* Scan Start Time */
749  double TAI93sec = ((double*) scn[SCN_STIME]->data)[iscan];
750  if (TAI93sec < 0.0) {
751  l1rec->scantime = BAD_FLT;
752  } else {
753  l1rec->scantime = tai93_to_unix(TAI93sec);
754  }
755 
756  /* Mirror Side */
757  l1rec->mside = (int32_t) ((char*) scn[SCN_MSIDE]->data)[iscan];
758 
759  /* Earth-sun distance correction for this scan */
760  int16_t year, day;
761  double dsec;
762  unix2yds(l1rec->scantime, &year, &day, &dsec);
763  int32_t yr = (int32_t) year;
764  int32_t dy = (int32_t) day;
765  int32_t msec = (int32_t) (dsec * 1000.0);
766  double esdist = esdist_(&yr, &dy, &msec);
767  l1rec->fsol = pow(1.0 / esdist, 2);
768 
769  /* TO DO: check SCN_MODE, SCN_QUAL for both GEO and L1B */
770  int8_t scanqual = ((int8_t*) l1bscn[L1BSCN_QUAL]->data)[iscan];
771 
772  if (scanqual & 1) {
773  l1file->sv_with_moon = 1;
774  }
775 
776  // if lonlat mode, read only lon, lat and solz
777  if (lonlat) {
778  read_var_1line(geo[GEO_LAT], iline);
779  read_var_1line(geo[GEO_LON], iline);
780  read_var_1line(geo[GEO_SOLZ], iline);
781  }
782  // else, read all the variables
783  else {
784  /*----- Geolocation swath values -----*/
785  for (ivar = 0; ivar < NVARS_GEO; ivar++)
786  read_var_1line(geo[ivar], iline);
787  }
788 
789  // copy lon, lat and solz into l1rec
790  memcpy(l1rec->lat, geo[GEO_LAT]->data, nbytes);
791  memcpy(l1rec->lon, geo[GEO_LON]->data, nbytes);
792  scale_short(geo[GEO_SOLZ], l1rec->solz);
793 
794  // stop reading if lonlat mode, dont need anything else
795  if (lonlat)
796  return LIFE_IS_GOOD;
797 
798  scale_short(geo[GEO_HGT], l1rec->height);
799  scale_short(geo[GEO_SOLA], l1rec->sola);
800  scale_short(geo[GEO_SENZ], l1rec->senz);
801  scale_short(geo[GEO_SENA], l1rec->sena);
802 
803  /* Load Angles */
804  for (i = 0; i < 3; i++) {
805  ang[i] = ((float*) nav[NAV_ANG]->data)[iscan * 3 + i]; // degrees
806  pos[i] = ((float*) nav[NAV_POS]->data)[iscan * 3 + i] / 1000.; // m -> km
807  vel[i] = ((float*) nav[NAV_VEL]->data)[iscan * 3 + i] / 1000.; // m/s -> km/s
808  }
809 
810  /* Check for non-nominal roll, pitch, or yaw */
811  /* badatt = */
812  /* (fabs(ang[0]) > MAX_ATTERR) || */
813  /* (fabs(ang[1]) > MAX_ATTERR) || */
814  /* (fabs(ang[2]) > MAX_ATTERR); */
815 
816  /* Compute polarization rotation angles */
817  ocorient_(pos, vel, ang, sen_mat, coeff);
818  for (i = 0; i < 3; i++)
819  mnorm[i] = sen_mat[i][0];
820  compute_alpha(l1rec->lon, l1rec->lat,
821  l1rec->senz, l1rec->sena,
822  mnorm, l1rec->npix, l1rec->alpha);
823 
824  //check for moon in spaceview port
825  //l1file->sv_with_moon = 1; // as in l1_viirs_h5.c
826 
827  /*----- Radiometric swath values -----*/
828  for (ivar = 0; ivar < MAXBANDS; ivar++) {
829  if (l1b[ivar] == NULL) continue; // skip missing band
830  read_var_1line(l1b[ivar], iline);
831  read_var_1line(l1bq[ivar], iline);
832  }
833  scale_l1bvals(l1rec); // get reflectance & radiance
834  radiance2bt(l1rec, -1); // calculate brightness temperature
835 
836  /*----- Check pixel values -----*/
837  for (ipix = 0; ipix < l1rec->npix; ipix++) {
838  l1rec->pixnum[ipix] = ipix + extract_pixel_start;
839  if (l1rec->scantime < 0)
840  l1rec->flags[ipix] |= NAVFAIL;
841  if (scanqual & 4)
842  l1rec->flags[ipix] |= NAVFAIL;
843 
844  // 1 Input_invalid
845  // 2 Pointing_bad
846  // 4 Terrain_bad
847  // check for Input_invalid or Pointing_bad (==3)
848  if (((unsigned char*) (geo[GEO_QUAL]->data))[ipix] & 3)
849  l1rec->flags[ipix] |= NAVFAIL;
850 
851  // check for Terrain_bad (==4)
852  if (((unsigned char*) (geo[GEO_QUAL]->data))[ipix] & 4)
853  l1rec->flags[ipix] |= NAVWARN;
854 
855  flag_bowtie_deleted(l1rec, ipix, extract_pixel_start);
856  //l1rec->navwarn[ipix] |= (l1file->sv_with_moon == 1);
857 
858  // convert the fill values to the proper value
859  if (l1rec->lat[ipix] == latGeoFillValue)
860  l1rec->lat[ipix] = latL2FillValue;
861  if (l1rec->lon[ipix] == lonGeoFillValue)
862  l1rec->lon[ipix] = lonL2FillValue;
863  if (l1rec->sena[ipix] == senaGeoFillValue)
864  l1rec->sena[ipix] = senaL2FillValue;
865  if (l1rec->senz[ipix] == senzGeoFillValue)
866  l1rec->senz[ipix] = senzL2FillValue;
867  if (l1rec->sola[ipix] == solaGeoFillValue)
868  l1rec->sola[ipix] = solaL2FillValue;
869  if (l1rec->solz[ipix] == solzGeoFillValue)
870  l1rec->solz[ipix] = solzL2FillValue;
871  }
872 
873  return SUCCESS;
874 }
875 
876 void flag_bowtie_deleted(l1str *l1rec, size_t ipix, int extract_offset) {
877  int pix = ipix + extract_offset;
878  if (pix < AGZONE1 || pix >= AGZONE5) {
879  if (l1rec->detnum < 2 || l1rec->detnum > 13) {
880  l1rec->flags[ipix] |= BOWTIEDEL;
881  }
882  } else if ((pix >= AGZONE1 && pix < AGZONE2)
883  || (pix >= AGZONE4 && pix < AGZONE5)) {
884  if (l1rec->detnum == 0 || l1rec->detnum == 15) {
885  l1rec->flags[ipix] |= BOWTIEDEL;
886  }
887  }
888 }
889 
890 /***********************************************************************
891  Cleanup
892  ***********************************************************************/
893 
895  /* int32_t i; */
896 
897  /* Free memory allocated for global variables */
898  /* for (i = 0; i < NVARS_SCN; i++) */
899  /* free_vars_nc(scn[i],1) ; */
900  /* for (i = 0; i < NVARS_NAV; i++) */
901  /* free_vars_nc(nav[i],1) ; */
902  /* for (i = 0; i < NVARS_GEO; i++) */
903  /* free_vars_nc(geo[i],1) ; */
904  /* for (i = 0; i < NVARS_L1BSCN; i++) */
905  /* free_vars_nc(l1bscn[i],1) ; */
906  /* for (i = 0; i < MAXBANDS; i++) { */
907  /* free_vars_nc(l1b[i], 1) ; */
908  /* free_vars_nc(l1bq[i],1) ; */
909  /* } */
910 
911  /* Free any memory allocated for input files */
912  /* End NetCDF access */
913 
914  /* TO DO: make grp_str_nc l1bfile & geofile global?
915  Then use free_grp_nc & nc_close. Won't need free_vars_nc calls above,
916  since they point to same var_str_nc memory. */
917 
918  return SUCCESS;
919 }
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
size_t len
Definition: nc4utils.h:44
#define SUCCESS
Definition: ObpgReadGrid.h:15
bandtypes
Definition: l1_viirs_l1b.c:148
#define TRY_NC(file, line, ncstat)
Definition: nc4utils.h:31
void freeProductInfo(productInfo_t *info)
int32_t day
int status
Definition: l1_czcs_hdf.c:32
void print_viirs_file(const viirs_file info)
Definition: l1_viirs_l1b.c:207
double * scan_start_time
Definition: l1_viirs_l1b.c:204
double mnorm[3]
#define NBANDSIR
Definition: filehandle.h:23
@ NVARS_GEO
Definition: l1_viirs_l1b.c:107
@ SCN_STIME
Definition: l1_viirs_l1b.c:43
void flag_bowtie_deleted(l1str *l1rec, size_t ipix, int extract_offset)
Definition: l1_viirs_l1b.c:876
size_t index
Definition: l1_viirs_l1b.c:32
int grpid
Definition: nc4utils.h:55
nav_var
Definition: l1_viirs_l1b.c:65
#define NULL
Definition: decode_rs.h:63
@ SCN_MODE
Definition: l1_viirs_l1b.c:47
void ocorient_(float *pos, float *vel, float *att, float(*)[3], float *coef)
int closel1_viirs_l1b()
Definition: l1_viirs_l1b.c:894
@ SCN_QUAL
Definition: l1_viirs_l1b.c:48
read l1rec
@ NAV_ANG
Definition: l1_viirs_l1b.c:67
#define VIIRSN
Definition: sensorDefs.h:23
l1bscn_var
Definition: l1_viirs_l1b.c:127
@ NAV_POS
Definition: l1_viirs_l1b.c:68
#define NAV_GRP
Definition: goci_slot.c:12
float32 * pos
Definition: l1_czcs_hdf.c:35
int32 * msec
Definition: l1_czcs_hdf.c:31
int init_viirs_l1bfile(viirs_file l1binfo)
Definition: l1_viirs_l1b.c:292
@ GEO_SENZ
Definition: l1_viirs_l1b.c:102
double esdist_(int32_t *year, int32_t *day, int32_t *msec)
int readl1_viirs_l1b(filehandle *l1file, const int32_t iline, l1str *l1rec, int lonlat)
Definition: l1_viirs_l1b.c:724
@ SCN_MSIDE
Definition: l1_viirs_l1b.c:46
char end_time[25]
Definition: l1_viirs_l1b.c:203
#define AGZONE5
Definition: l1.h:69
void bzero()
int orbit_number
Definition: l1_viirs_l1b.c:201
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
char start_time[25]
Definition: l1_viirs_l1b.c:202
character(len=1000) if
Definition: names.f90:13
#define LIFE_IS_GOOD
Definition: passthebuck.h:4
#define PI
Definition: l3_get_org.c:6
nc_type type
Definition: nc4utils.h:58
int openl1_viirs_l1b(filehandle *l1file)
Definition: l1_viirs_l1b.c:426
@ L1BSCN_STIME
Definition: l1_viirs_l1b.c:128
void compute_alpha(float lon[], float lat[], float senz[], float sena[], double mnorm[3], int npix, float alpha[])
Definition: compute_alpha.c:8
@ GEO_HGT
Definition: l1_viirs_l1b.c:99
int load_grp_nc(grp_str_nc *grp)
Definition: nc4utils.c:423
int ndims
Definition: nc4utils.h:59
@ GEO_LAT
Definition: l1_viirs_l1b.c:97
subroutine lonlat(alon, alat, xlon, ylat)
Definition: lonlat.f:2
size_t npixls
Definition: l1_viirs_l1b.c:200
int leapseconds_since_1993(double tai93time)
Definition: leapsecond.c:58
@ L1BSCN_MODE
Definition: l1_viirs_l1b.c:131
unsigned short ushort
Definition: l1_viirs_l1b.c:16
productInfo_t * allocateProductInfo()
double tai93_to_unix(double tai93)
Definition: yds2tai.c:16
int id
Definition: nc4utils.h:56
@ NVARS_NAV
Definition: l1_viirs_l1b.c:76
var_str_nc * find_var_byname_nc(grp_str_nc nc, const char *varname, const char *grpname)
Definition: nc4utils.c:314
int scale_short(var_str_nc *var, float *dest)
Definition: l1_viirs_l1b.c:545
char file[FILENAME_MAX]
Definition: l1_viirs_l1b.c:196
int read_var_1line(var_str_nc *var, size_t iline)
Definition: l1_viirs_l1b.c:517
void * data
Definition: nc4utils.h:64
#define AGZONE4
Definition: l1.h:68
@ GEO_SENA
Definition: l1_viirs_l1b.c:101
l1_input_t * l1_input
Definition: l1_options.c:9
@ GEO_LON
Definition: l1_viirs_l1b.c:98
@ NAV_VEL
Definition: l1_viirs_l1b.c:69
char * name
Definition: l1_viirs_l1b.c:33
void unix2yds(double usec, short *year, short *day, double *secs)
size_t nlines
Definition: l1_viirs_l1b.c:199
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
@ NVARS_SCN
Definition: l1_viirs_l1b.c:49
int32_t id
Definition: l1_viirs_l1b.c:195
#define RADEG
Definition: czcs_ctl_pt.c:5
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
#define SOLZNIGHT
Definition: l1.h:58
int findProductInfo(const char *productName, int sensorId, productInfo_t *info)
int scale_l1bvals(l1str *l1rec)
Definition: l1_viirs_l1b.c:579
dim_str_nc * dim
Definition: nc4utils.h:60
#define MAXBANDS
Definition: l1_viirs_l1b.c:172
#define BAD_FLT
Definition: jplaeriallib.h:19
@ NVARS_L1BSCN
Definition: l1_viirs_l1b.c:133
int load_fcal_lut(char *calfile, int64_t UTC58usec, double ****ftable)
#define AGZONE2
Definition: l1.h:66
@ CIR
Definition: l1_viirs_l1b.c:149
void radiance2bt(l1str *l1rec, int resolution)
Definition: brightness.c:170
geo_var
Definition: l1_viirs_l1b.c:96
#define BOWTIEDEL
Definition: l2_flags.h:39
@ GEO_SOLZ
Definition: l1_viirs_l1b.c:104
#define fabs(a)
Definition: misc.h:93
int32_t iscan
size_t nscans
Definition: l1_viirs_l1b.c:198
int readall_var(var_str_nc *var)
Definition: nc4utils.c:334
@ L1BSCN_QUAL
Definition: l1_viirs_l1b.c:132
@ RSB
Definition: l1_viirs_l1b.c:149
int init_viirs_file(const char filename[FILENAME_MAX], viirs_file *info)
Definition: l1_viirs_l1b.c:217
real *8 function esdist(iyr, iday, msec)
Definition: esdist.f:3
#define AGZONE1
Definition: l1.h:65
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
@ GEO_SOLA
Definition: l1_viirs_l1b.c:103
l2prod offset
#define TRYMEM(file, line, memstat)
Definition: l1_hmodis_hdf.c:23
@ GEO_QUAL
Definition: l1_viirs_l1b.c:106
#define NAVWARN
Definition: l2_flags.h:27
int i
Definition: decode_rs.h:71
scn_var
Definition: l1_viirs_l1b.c:42
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int npix
Definition: get_cmp.c:28
int init_viirs_geofile(viirs_file geoinfo)
Definition: l1_viirs_l1b.c:336
@ TEB
Definition: l1_viirs_l1b.c:149
char title[FILENAME_MAX]
Definition: l1_viirs_l1b.c:197
#define NAVFAIL
Definition: l2_flags.h:36
int count
Definition: decode_rs.h:79