OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
l1_sgli.c
Go to the documentation of this file.
1 /* ============================================================================ */
2 /* module l1_sgli.c - functions to read SGLI to MSL12
3  * W. Robinson, SAIC 28 Sep 2016
4  * ============================================================================ */
5 
6 #include <math.h>
7 
8 #include "l1.h"
9 #include <h5io.h>
10 #include "sgli.h"
11 #include "l1_sgli.h"
12 
13 #define SGLI_NVNIR 9
14 #define SGLI_NVNIR_LG 2
15 #define SGLI_NVNIR_GEOM 2
16 #define SGLI_NSWIR 4
17 
18 /* swirir_exist indicates (1) if the swir/ir file is available and
19  resqkm indicates (1) if data is at quarter km resolution, else 1 km */
20 static int32_t *scan_times, swirir_exist, resqkm, resample, resamp_1km;
21 static int32_t n_vnir = SGLI_NVNIR, n_vnir_lg = SGLI_NVNIR_LG;
22 static int32_t n_vnir_geom = SGLI_NVNIR_GEOM, n_swir = SGLI_NSWIR;
23 static float *flt_buf, vnir_scale[SGLI_NVNIR], vnir_off[SGLI_NVNIR];
24 static float swir_scale[SGLI_NSWIR], swir_off[SGLI_NVNIR];
25 static float vnir_sat[SGLI_NVNIR], swir_sat[SGLI_NSWIR];
26 static float vnir_geom_scale[SGLI_NVNIR_GEOM], vnir_geom_off[SGLI_NVNIR_GEOM];
27 static float vnir_lg_scale[SGLI_NVNIR_LG], vnir_lg_off[SGLI_NVNIR_LG];
28 static float vnir_lg_sat[SGLI_NVNIR_LG];
29 static uint16_t vnir_mx_dn[SGLI_NVNIR], vnir_lg_mx_dn[SGLI_NVNIR_LG];
30 static uint16_t swir_mx_dn[SGLI_NSWIR];
31 static uint16_t *ui16_buf_1km;
32 static int32_t npix, npix_1km, nscan, npix_tie, nscan_tie;
33 static int32_t npix_tie_1km, nscan_tie_1km;
34 static h5io_str vnir_fid, swir_ir_fid, vnir_hg_dsid[SGLI_NVNIR];
35 static h5io_str vnir_lg_dsid[SGLI_NVNIR_LG], gen_dsid;
36 static h5io_str vnir_geom_dsid[SGLI_NVNIR_GEOM], swir_dsid[SGLI_NSWIR];
37 static double st_day_unix;
38 static char *geo_q, *swir_q, *sen_q, *sol_q;
39 static gsl_interp_accel *accel_x, *accel_y;
40 static gsl_spline2d *geo_int_id[3];
41 static double *xa, *ya, *xa_1km, *ya_1km;
42 static double *geo_x, *geo_y, *geo_z;
43 static float *tie_el, *tie_az;
44 static grid_res_str grid_res[2];
45 static band_geom_str band_geom_vnir[SGLI_NVNIR]; /* for sensor geom_per_band */
46 static band_geom_str band_geom_vnir_lg[SGLI_NVNIR_LG];
47 static band_geom_str band_geom_swir[SGLI_NSWIR];
48 static band_geom_str band_geom_nominal[2]; /* for nominal sen, sol angles */
49 
50 sgli_t* createPrivateData_sgli() {
51 
52  sgli_t* data = (sgli_t*) calloc(1, sizeof (sgli_t));
53  if (data == NULL) {
54  fprintf(stderr, "-E- %s line %d: unable to allocate private data for sgli\n",
55  __FILE__, __LINE__);
56  exit(1);
57  }
58 
59  return data;
60 }
61 
62 /*
63  *
64  */
65 int sgli_file_ver(h5io_str *fid) {
66  /*
67  sgli_file_ver
68 
69  purpose: do a basic validity check on a SGLI file
70 
71  Returns 0 if all checks are OK
72  Parameters: (in calling order)
73  Type Name I/O Description
74  ---- ---- --- -----------
75  h5io_str * fid I file id of file to check
76 
77  Modification history:
78  Programmer Date Description of change
79  ---------- ---- ---------------------
80  Wayne Robinson 14 Oct 2016 Original development
81 
82  -----------------------------------------------------------------------*/
83  int32_t stage_pass = 0;
84  h5io_str g_id;
85  char str_store[255];
86  /*
87  This will check (in group Global_attributes/attributes):
88  Satellite: Global Change Observation Mission - Climate (GCOM-C)
89  Sensor: Second-generation Global Imager (SGLI)
90  Product_level: Level-1B
91  Product_name: Top of atmosphere radiance (reflectance)
92  */
93 
94  if (h5io_set_grp(fid, "Global_attributes", &g_id) == 0) {
95  if (h5io_attr_exist(&g_id, "Satellite") == 0) {
96  if (h5io_rd_attr(&g_id, "Satellite", (void *) str_store) == 0) {
97  if (strncmp(str_store, "Global Change Observation Mission - Climate (GCOM-C)", 52) == 0) {
98  stage_pass = 1;
99  }
100  }
101  }
102  if ((stage_pass == 1) && (h5io_attr_exist(&g_id, "Sensor") == 0)) {
103  stage_pass = 0;
104  if (h5io_rd_attr(&g_id, "Sensor", (void *) str_store) == 0) {
105  if (strncmp(str_store, "Second-generation Global Imager (SGLI)", 38) == 0) {
106  stage_pass = 1;
107  }
108  }
109  }
110  if ((stage_pass == 1) && (h5io_attr_exist(&g_id, "Product_level") == 0)) {
111  stage_pass = 0;
112  if (h5io_rd_attr(&g_id, "Product_level", (void *) str_store) == 0) {
113  if (strncmp(str_store, "Level-1B", 8) == 0) {
114  stage_pass = 1;
115  }
116  }
117  }
118  if ((stage_pass == 1) && (h5io_attr_exist(&g_id, "Product_name") == 0)) {
119  stage_pass = 0;
120  if (h5io_rd_attr(&g_id, "Product_name", (void *) str_store) == 0) {
121  if (strncmp(str_store, "Top of atmosphere radiance (reflectance)", 40) == 0) {
122  stage_pass = 1;
123  }
124  }
125  }
126  h5io_close(&g_id);
127  }
128  return ( stage_pass == 0) ? 1 : 0;
129 }
130 
131 int sgli_rad_info(h5io_str *fid, char *ds_nam, h5io_str *dsid,
132  int *dim_siz, uint16_t *mx_dn, float *scale, float *off, float *sat)
133 /*-----------------------------------------------------------------------
134  sgli_rad_info
135 
136  purpose: get information about the radiance arrays in SGLI files
137 
138  Returns 0 if all checks are OK
139  Parameters: (in calling order)
140  Type Name I/O Description
141  ---- ---- --- -----------
142  h5io_str * fid I file id of file
143  char * ds_nam I dataset name in file
144  h5io_str * dsid O returned id of opened dataset
145  int * dim_siz O size of the dimensions
146  uint16_t * mx_dn O maximum DN for band counts
147  float * scale O scaling value for the data
148  float * off O offset for scaling
149  float * sat O saturation value of rad
150 
151  Modification history:
152  Programmer Date Description of change
153  ---------- ---- ---------------------
154  W. Robinson 20 Oct 2016 Original development
155 
156  -----------------------------------------------------------------------*/ {
157  H5T_class_t h5_class;
158  hid_t h5_native_typ;
159  int ndim, sto_len;
160 
161  if (h5io_set_ds(fid, ds_nam, dsid) != 0) {
162  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
163  __LINE__, ds_nam);
164  return 1;
165  }
166  /* get the sizes: npix, nscan and the resolution */
167  if (h5io_info(dsid, NULL, &h5_class, &h5_native_typ, &ndim, dim_siz,
168  &sto_len) != 0) {
169  fprintf(stderr, "%s, %d, E - h5io_info failure, dataset: %s\n",
170  __FILE__, __LINE__, ds_nam);
171  return 1;
172  }
173  /* get the max dn, scale, offset for these bands */
174  if (h5io_rd_attr(dsid, "Maximum_valid_DN", (uint16_t *) mx_dn) != 0) {
175  fprintf(stderr, "%s, %d, E - Cannot read Maximum_valid_DN for %s\n",
176  __FILE__, __LINE__, ds_nam);
177  return 1;
178  }
179  if (h5io_rd_attr(dsid, "Slope", (float *) scale) != 0) {
180  fprintf(stderr, "%s, %d, E - Cannot read Slope for %s\n",
181  __FILE__, __LINE__, ds_nam);
182  return 1;
183  }
184  if (h5io_rd_attr(dsid, "Offset", (float *) off) != 0) {
185  fprintf(stderr, "%s, %d, E - Cannot read Offset for %s\n",
186  __FILE__, __LINE__, ds_nam);
187  return 1;
188  }
189  if (h5io_rd_attr(dsid, "Saturation_radiance", (float *) sat) != 0) {
190  fprintf(stderr, "%s, %d, E - Cannot read Saturation_radiance for %s\n",
191  __FILE__, __LINE__, ds_nam);
192  return 1;
193  }
194  return 0;
195 }
196 
197 int view_ang_tie_init2(band_geom_str *band_geom, int32_t nbands )
198 /*-----------------------------------------------------------------------
199  view_ang_tie_init2
200 
201  purpose: initialize tie point view angle data for interpolation
202 
203  Returns 0
204 
205  Parameters: (in calling order)
206  Type Name I/O Description
207  ---- ---- --- -----------
208  band_geom_str * band_geom I/O array of structures for tie
209  point data
210  int32_t nbands I # bands (above array size
211 
212  Modification history:
213  Programmer Date Description of change
214  ---------- ---- ---------------------
215  Wayne Robinson 9 Dec 2016 Original development
216  Wayne Robinson 6 Mar 2019 Update for file format change
217 
218  -----------------------------------------------------------------------*/ {
219  int32_t ib, ip, il, ipl, npix_tie, nlin_tie;
220  int16_t *tie_el_16, *tie_az_16;
221  int start[2], count[2];
222  double *geo[3], el_cvt, az_cvt, xy_mod, *xa, *ya;
223  double deg2rad(double);
224  /*
225  * initial set up
226  */
227  tie_el_16 = (int16_t *) tie_el;
228  tie_az_16 = (int16_t *) tie_az;
229  /*
230  * process each band, getting the interpolation object for a unit vector
231  * version of the sensor angle
232  */
233  for (ib = 0; ib < nbands; ib++) {
234  npix_tie = band_geom[ib].grd_desc->npix_tie;
235  nlin_tie = band_geom[ib].grd_desc->nscn_sub_tie;
236  xa = band_geom[ib].grd_desc->xa;
237  ya = band_geom[ib].grd_desc->ya;
238  if ((band_geom[ib].qual =
239  (char *) calloc(npix_tie * nlin_tie, sizeof ( char))) == NULL) {
240  fprintf(stderr,
241  "%s, %d, E - Failed to allocate vnir sensor, solar quality storage\n",
242  __FILE__, __LINE__);
243  return 1;
244  }
245  /*
246  * read each zenith, azimuth
247  */
248  start[0] = band_geom[ib].grd_desc->tie_st_lin;
249  start[1] = 0;
250  count[0] = nlin_tie;
251  count[1] = npix_tie;
252 
253  if ((h5io_rd_ds_slice(&band_geom[ib].dsid[1], start, count,
254  (void *) tie_el_16) != 0) ||
255  (h5io_rd_ds_slice(&band_geom[ib].dsid[0], start, count,
256  (void *) tie_az_16) != 0)) {
257  fprintf(stderr, "%s, %d, E Error reading vnir view angle tie points\n",
258  __FILE__, __LINE__);
259  return 1;
260  }
261 
262  for (il = 0; il < nlin_tie; il++) {
263  for (ip = 0; ip < npix_tie; ip++) {
264  ipl = ip + il * npix_tie;
265  /* scale to float */
266  el_cvt = band_geom[ib].offset[1] +
267  band_geom[ib].scale[1] * tie_el_16[ipl];
268  az_cvt = band_geom[ib].offset[0] +
269  band_geom[ib].scale[0] * tie_az_16[ipl];
270 
271  if ((el_cvt < 0.) || (el_cvt > 180.) ||
272  (az_cvt < -180.) || (az_cvt > 180.)) {
273  el_cvt = 0;
274  az_cvt = 0;
275  band_geom[ib].qual[ipl] = 1;
276  }
277  xy_mod = sin(deg2rad(el_cvt));
278  geo_x[ipl] = xy_mod * cos(deg2rad(az_cvt));
279  geo_y[ipl] = xy_mod * sin(deg2rad(az_cvt));
280  geo_z[ipl] = cos(deg2rad(el_cvt));
281  }
282  }
283  geo[0] = geo_x;
284  geo[1] = geo_y;
285  geo[2] = geo_z;
286  /* set up for the 2D spline interpolation and accelerator */
287  for (ip = 0; ip < 3; ip++) {
288  band_geom[ib].int_id_sen[ip] =
289  gsl_spline2d_alloc(gsl_interp2d_bicubic, npix_tie, nlin_tie);
290  if ((gsl_spline2d_init(band_geom[ib].int_id_sen[ip], xa,
291  ya, geo[ip], npix_tie, nlin_tie)) != 0) {
292  fprintf(stderr, "%s, %d, E vnir gsl_spline2d_init error\n",
293  __FILE__, __LINE__);
294  return 1;
295  }
296  }
297  }
298  /*
299  * interpolation coefficients are made
300  */
301  return 0;
302 }
303 
304 int view_ang_tie_eval(band_geom_str *band_geom, int32_t ip, int32_t scan,
305  float *azi_v, float *zen_v, char *qual)
306 /*-----------------------------------------------------------------------
307  view_ang_tie_eval
308 
309  purpose: interpolate view angle to a pixel
310 
311  Returns 0
312 
313  Parameters: (in calling order)
314  Type Name I/O Description
315  ---- ---- --- -----------
316  band_geom_str * band_geom I structure for tie
317  point data for that band
318  int32_t ip I specific pixel to evaluate
319  int32_t scan I specific scan to evaluate
320  float * azi_v O returned azimuth
321  float * zen_v O returned zenith
322  char * qual O quality 1 = bad
323 
324  Modification history:
325  Programmer Date Description of change
326  ---------- ---- ---------------------
327  Wayne Robinson 12 Dec 2016 Original development
328 
329  -----------------------------------------------------------------------*/ {
330  double xvec, yvec, zvec;
331  double rad2deg(double);
332  int32_t ipb, scan_tie, tie_en_line, npix_tie, nscan_tie;
333  size_t ip_tie, il_tie;
334 
335  npix_tie = band_geom->grd_desc->npix_tie;
336  nscan_tie = band_geom->grd_desc->nscn_sub_tie;
337  tie_en_line = band_geom->grd_desc->tie_st_lin + nscan_tie - 1;
338  scan_tie = ( ( scan >= band_geom->grd_desc->resamp * tie_en_line)) ?
339  band_geom->grd_desc->resamp * tie_en_line - 1 : scan;
340 
341  if ((gsl_spline2d_eval_e(band_geom->int_id_sen[0], ip,
342  scan, accel_x, accel_y, &xvec) != 0) ||
343  (gsl_spline2d_eval_e(band_geom->int_id_sen[1], ip,
344  scan, accel_x, accel_y, &yvec) != 0) ||
345  (gsl_spline2d_eval_e(band_geom->int_id_sen[2], ip,
346  scan, accel_x, accel_y, &zvec) != 0)) {
347  fprintf(stderr, "%s, %d, E Error in vnir_lg gsl_spline2d_eval\n",
348  __FILE__, __LINE__);
349  return 1;
350  }
351  *zen_v = (float) rad2deg(acos(zvec));
352  *azi_v = (float) rad2deg(atan2(yvec, xvec));
353  /* it is possible that tie point view geom is bad, so */
354  ip_tie = gsl_interp_accel_find(accel_x, band_geom->grd_desc->xa, npix_tie, ip);
355  il_tie = gsl_interp_accel_find(accel_y, band_geom->grd_desc->ya, nscan_tie, scan_tie);
356  ipb = ip_tie + il_tie * npix_tie;
357  *qual = 0;
358  if ((band_geom->qual[ipb] != 0) ||
359  (band_geom->qual[ ipb + 1 ] != 0) ||
360  (band_geom->qual[ ipb + npix_tie ] != 0) ||
361  (band_geom->qual[ ipb + 1 + npix_tie ] != 0))
362  *qual = 1;
363  return 0;
364 }
365 
366 int openl1_sgli(filehandle * l1file) {
367  char *cptr, *vis_nir_file, *swir_ir_file, ds_nam[FILENAME_MAX];
368  char str_store[FILENAME_MAX];
369  int32_t ibnd, ip;
370  H5T_class_t h5_class;
371  hid_t h5_native_typ;
372  h5io_str g_id;
373  int ndim, sto_len, dim_siz[4];
374  float resolution;
375  int st_year, st_mon, st_day;
376  char *vnir_hg_bnam[] ={"01", "02", "03", "04", "05", "06", "07", "09", "10"};
377  char *vnir_lg_bnam[] = {"08", "11"};
378  char *swir_bnam[] = {"01", "02", "03", "04"};
379  char *vnir_geom_nam[] = {"Latitude", "Longitude"};
380  char *view_nam[] = {"Sensor", "Solar"};
381  char *ang_nam[] = {"azimuth", "zenith"};
382  char *sen_ang[2] = {"Sensor_azimuth", "Sensor_zenith"};
383 
384  //printf("%s, %d - I - SGLI vis/nir file name: %s\n", __FILE__, __LINE__,
385  // l1file->name);
386 
387  sgli_t *data = l1file->private_data = createPrivateData_sgli();
388 
389  /*
390  * Derive the SWIR/IR file name from the l1 file if not set
391  */
392  vis_nir_file = l1file->name;
393  swir_ir_file = data->swir_ir_file;
394 
395  if (swir_ir_file[0] == 0) {
396  strcpy(swir_ir_file, vis_nir_file);
397  if ((cptr = strstr(swir_ir_file, "VNRD")) == NULL) {
398  printf("%s, %d, I - VIS / NIR file has non-standard name - cannot create SWIR/IR file name\n", __FILE__, __LINE__);
399  printf("VIS/NIR file name: %s\n", vis_nir_file);
400  printf("Processing will only be done with the VIS/NIR data\n");
401  } else {
402  memcpy(cptr, "IRS", 3);
403  }
404  }
405  /*
406  * open the vis/nir file and all the band datasets
407  */
408  if (h5io_openr(vis_nir_file, 0, &vnir_fid) != 0) {
409  fprintf(stderr, "%s, %d, E - Failure to open %s\n", __FILE__,
410  __LINE__, vis_nir_file);
411  return 1;
412  }
413  /* set bands for reading, get info and check consistency */
414  for (ibnd = 0; ibnd < n_vnir; ibnd++) {
415  sprintf(ds_nam, "Image_data/Lt_VN%s", vnir_hg_bnam[ibnd]);
416  if (sgli_rad_info(&vnir_fid, ds_nam, (vnir_hg_dsid + ibnd),
417  dim_siz, (vnir_mx_dn + ibnd), (vnir_scale + ibnd),
418  (vnir_off + ibnd), (vnir_sat + ibnd)) != 0)
419  return 1;
420 
421  if (ibnd == 0) {
422  l1file->npix = npix = dim_siz[1];
423  l1file->nscan = nscan = dim_siz[0];
424 
425  /* officially set the resolution */
426  if (h5io_rd_attr((vnir_hg_dsid + ibnd), "Spatial_resolution",
427  (void *) &resolution) != 0) {
428  fprintf(stderr, "%s, %d, E - Cannot read band resolution\n", __FILE__,
429  __LINE__);
430  return 1;
431  }
432  resqkm = (resolution == 1000.) ? 0 : 1;
433 
434  /* set the read buffer for rad info */
435  if ((flt_buf = (float *) malloc(npix * sizeof (float))) == NULL) {
436  fprintf(stderr, "%s, %d, E - Cannot allocate the SGLI read buffer\n",
437  __FILE__, __LINE__);
438  return 1;
439  }
440  } else {
441  if ((dim_siz[1] != npix) || (dim_siz[0] != nscan)) {
442  fprintf(stderr, "%s, %d, E - Vis/NIR band array size mismatch\n",
443  __FILE__, __LINE__);
444  fprintf(stderr, "for band %s\n", ds_nam);
445  return 1;
446  }
447  }
448  }
449  /* open the low gain 673 and 865 nm bands and check */
450  for (ibnd = 0; ibnd < n_vnir_lg; ibnd++) {
451  sprintf(ds_nam, "Image_data/Lt_VN%s", vnir_lg_bnam[ibnd]);
452 
453  if (sgli_rad_info(&vnir_fid, ds_nam, (vnir_lg_dsid + ibnd),
454  dim_siz, (vnir_lg_mx_dn + ibnd), (vnir_lg_scale + ibnd),
455  (vnir_lg_off + ibnd), (vnir_lg_sat + ibnd)) != 0)
456  return 1;
457 
458  if ((dim_siz[1] != npix) || (dim_siz[0] != nscan)) {
459  fprintf(stderr, "%s, %d, E - Vis/NIR band array size mismatch\n",
460  __FILE__, __LINE__);
461  fprintf(stderr, "for band %s\n", ds_nam);
462  return 1;
463  }
464  }
465  /* get the start time */
466  if (h5io_set_grp(&vnir_fid, "Global_attributes", &g_id) != 0) {
467  fprintf(stderr, "%s, %d, E - Cannot set Global_attributes group\n",
468  __FILE__, __LINE__);
469  return 1;
470  }
471  if (h5io_rd_attr(&g_id, "Scene_start_time", (void *) str_store) != 0) {
472  fprintf(stderr, "%s, %d, E - Cannot read Scene_start_time\n",
473  __FILE__, __LINE__);
474  return 1;
475  }
476  h5io_close(&g_id);
477  sscanf(str_store, "%4d%2d%2d", &st_year, &st_mon, &st_day);
478  st_day_unix = ymds2unix((int16_t) st_year, (int16_t) st_mon, (int16_t) st_day,
479  0);
480 
481  /* get the time for each line - read it all here */
482  strcpy(ds_nam, "Image_data/Line_msec");
483  if (h5io_set_ds(&vnir_fid, ds_nam, &gen_dsid) != 0) {
484  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
485  __LINE__, ds_nam);
486  return 1;
487  }
488  scan_times = (int32_t *) malloc(nscan * sizeof (int32_t));
489  if (h5io_rd_ds(&gen_dsid, (void *) scan_times) != 0) {
490  fprintf(stderr, "%s, %d, E - Error reading dataset: %s\n", __FILE__,
491  __LINE__, ds_nam);
492  return 1;
493  }
494  if (h5io_close(&gen_dsid) != 0) {
495  fprintf(stderr, "%s, %d, E - Error closing dataset: %s\n", __FILE__,
496  __LINE__, ds_nam);
497  return 1;
498  }
499 
500  /* open the view geometry and geolocation */
501  /* and get sizes */
502  for (ibnd = 0; ibnd < n_vnir_geom; ibnd++) {
503  sprintf(ds_nam, "Geometry_data/%s", vnir_geom_nam[ibnd]);
504  if (h5io_set_ds(&vnir_fid, ds_nam, (vnir_geom_dsid + ibnd)) != 0) {
505  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
506  __LINE__, ds_nam);
507  return 1;
508  }
509  /* for first band, get the sizes: npix, nscan and the resampling interval */
510  if (ibnd == 0) {
511  if (h5io_info((vnir_geom_dsid + ibnd), NULL, &h5_class, &h5_native_typ,
512  &ndim, dim_siz, &sto_len) != 0) {
513  fprintf(stderr, "%s, %d, E - h5io_info failure\n", __FILE__,
514  __LINE__);
515  return 1;
516  }
517  npix_tie = dim_siz[1];
518  nscan_tie = dim_siz[0];
519  if (h5io_rd_attr((vnir_geom_dsid + ibnd), "Resampling_interval",
520  (void *) &resample) != 0) {
521  fprintf(stderr,
522  "%s, %d, E - Failed to read resampling for geo data\n",
523  __FILE__, __LINE__);
524  return 1;
525  }
526  }
527  if ((h5io_rd_attr((vnir_geom_dsid + ibnd), "Slope",
528  (void *) (vnir_geom_scale + ibnd)) != 0) ||
529  (h5io_rd_attr((vnir_geom_dsid + ibnd), "Offset",
530  (void *) (vnir_geom_off + ibnd)) != 0)) {
531  fprintf(stderr,
532  "%s, %d, E - Failed to read scale or offset for geo data\n",
533  __FILE__, __LINE__);
534  return 1;
535  }
536 
537  }
538  /*
539  * For the nominal sensor, solar angles
540  */
541  for (ibnd = 0; ibnd < 2; ibnd++) {
542  for (ip = 0; ip < 2; ip++) {
543  sprintf(ds_nam, "Geometry_data/%s_%s",
544  view_nam[ibnd], ang_nam[ip]);
545  if (h5io_set_ds(&vnir_fid, ds_nam,
546  &band_geom_nominal[ibnd].dsid[ip]) != 0) {
547  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
548  __LINE__, ds_nam);
549  return 1;
550  }
551  if ((h5io_rd_attr(&band_geom_nominal[ibnd].dsid[ip], "Slope",
552  (void *) (&band_geom_nominal[ibnd].scale[ip])) != 0) ||
553  (h5io_rd_attr(&band_geom_nominal[ibnd].dsid[ip], "Offset",
554  (void *) (&band_geom_nominal[ibnd].offset[ip])) != 0)) {
555  fprintf(stderr,
556  "%s, %d, E - Failed to read scale or offset for nominal geo data\n",
557  __FILE__, __LINE__);
558  return 1;
559  }
560  /* as it is unsure that all geom_per_band senX sizes are same as
561  the tie point size, do a check to be sure */
562  if (h5io_info(&band_geom_nominal[ibnd].dsid[ip], NULL, &h5_class,
563  &h5_native_typ, &ndim, dim_siz, &sto_len) != 0) {
564  fprintf(stderr, "%s, %d, E - h5io_info failure for band_geom_nominal\n",
565  __FILE__, __LINE__);
566  return 1;
567  }
568  if ((dim_siz[1] != npix_tie) || (dim_siz[0] != nscan_tie)) {
569  fprintf(stderr,
570  "%s, %d, E - Unexpected band_geom_nominal sizes found\n",
571  __FILE__, __LINE__);
572  return 1;
573  }
574  }
575  band_geom_nominal[ibnd].grd_desc = grid_res;
576  }
577  /*
578  * if the geom per band is used, open all the sena, senz for vnir and
579  * low gain bands
580  */
581  if (l1_input->geom_per_band == 1) {
582  for (ibnd = 0; ibnd < n_vnir; ibnd++) {
583  for (ip = 0; ip < 2; ip++) {
584  sprintf(ds_nam, "Geometry_data/%s_VN%s",
585  sen_ang[ip], vnir_hg_bnam[ibnd]);
586  if (h5io_set_ds(&vnir_fid, ds_nam,
587  &band_geom_vnir[ibnd].dsid[ip]) != 0) {
588  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
589  __LINE__, ds_nam);
590  return 1;
591  }
592  if ((h5io_rd_attr(&band_geom_vnir[ibnd].dsid[ip], "Slope",
593  (void *) (&band_geom_vnir[ibnd].scale[ip])) != 0) ||
594  (h5io_rd_attr(&band_geom_vnir[ibnd].dsid[ip], "Offset",
595  (void *) (&band_geom_vnir[ibnd].offset[ip])) != 0)) {
596  fprintf(stderr,
597  "%s, %d, E - Failed to read scale or offset for geo data\n",
598  __FILE__, __LINE__);
599  return 1;
600  }
601  /* as it is unsure that all geom_per_band senX sizes are same as
602  the tie point size, do a check to e sure */
603  if (h5io_info(&band_geom_vnir[ibnd].dsid[ip], NULL, &h5_class,
604  &h5_native_typ, &ndim, dim_siz, &sto_len) != 0) {
605  fprintf(stderr, "%s, %d, E - h5io_info failure for geom_per_band\n",
606  __FILE__, __LINE__);
607  return 1;
608  }
609  if ((dim_siz[1] != npix_tie) || (dim_siz[0] != nscan_tie)) {
610  fprintf(stderr, "%s, %d, E - Unexpected geom_per_band sizes found\n", __FILE__, __LINE__);
611  return 1;
612  }
613  }
614  band_geom_vnir[ibnd].grd_desc = grid_res;
615  }
616  /* low gain band-dependent geometry */
617  for (ibnd = 0; ibnd < n_vnir_lg; ibnd++) {
618  for (ip = 0; ip < 2; ip++) {
619  sprintf(ds_nam, "Geometry_data/%s_VN%s",
620  sen_ang[ip], vnir_lg_bnam[ibnd]);
621  if (h5io_set_ds(&vnir_fid, ds_nam,
622  &band_geom_vnir_lg[ibnd].dsid[ip]) != 0) {
623  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
624  __LINE__, ds_nam);
625  return 1;
626  }
627  if ((h5io_rd_attr(&band_geom_vnir_lg[ibnd].dsid[ip], "Slope",
628  (void *) (&band_geom_vnir_lg[ibnd].scale[ip])) != 0) ||
629  (h5io_rd_attr(&band_geom_vnir_lg[ibnd].dsid[ip], "Offset",
630  (void *) (&band_geom_vnir_lg[ibnd].offset[ip])) != 0)) {
631  fprintf(stderr,
632  "%s, %d, E - Failed to read scale or offset for geo data\n",
633  __FILE__, __LINE__);
634  return 1;
635  }
636  /* as it is unsure that all geom_per_band senX sizes are same as
637  the tie point size, do a check to e sure */
638  if (h5io_info(&band_geom_vnir_lg[ibnd].dsid[ip], NULL, &h5_class,
639  &h5_native_typ, &ndim, dim_siz, &sto_len) != 0) {
640  fprintf(stderr,
641  "%s, %d, E - h5io_info failure for vis lg geom_per_band\n",
642  __FILE__, __LINE__);
643  return 1;
644  }
645  if ((dim_siz[1] != npix_tie) || (dim_siz[0] != nscan_tie)) {
646  fprintf(stderr,
647  "%s, %d, E - Unexpected vis_lg geom_per_band sizes found\n",
648  __FILE__, __LINE__);
649  return 1;
650  }
651  }
652  band_geom_vnir_lg[ibnd].grd_desc = grid_res;
653  }
654  }
655  /* now for the SWIR/IR file */
656  swirir_exist = 0;
657  if (swir_ir_file[0] != 0) {
658  /* open the swir file and make sure it is right type */
659  if (h5io_openr(swir_ir_file, 0, &swir_ir_fid) != 0) {
660  fprintf(stderr, "%s, %d, I - The SWIR/IR file: %s\n",
661  __FILE__, __LINE__, swir_ir_file);
662  fprintf(stderr, "does not exist. Using VIS/NIR file only\n");
663  } else {
664  swirir_exist = 1;
665  fprintf(stderr, "%s, %d, I - SGLI SWIR/IR file will be used\n",
666  __FILE__, __LINE__);
667  fprintf(stderr, " NAME: %s\n", swir_ir_file );
668  if (sgli_file_ver(&swir_ir_fid) != 0) {
669  fprintf(stderr,
670  "%s, %d, E - SWIR/IR file: %s is not a SGLI L1b file\n", __FILE__,
671  __LINE__, swir_ir_file);
672  return 1;
673  }
674  /* open the swir bands and check that the array sizes match */
675  for (ibnd = 0; ibnd < n_swir; ibnd++) {
676  sprintf(ds_nam, "Image_data/Lt_SW%s", swir_bnam[ibnd]);
677  if (sgli_rad_info(&swir_ir_fid, ds_nam, (swir_dsid + ibnd),
678  dim_siz, (swir_mx_dn + ibnd), (swir_scale + ibnd),
679  (swir_off + ibnd), (swir_sat + ibnd)) != 0)
680  return 1;
681  if ((resqkm == 1) && (ibnd != 2)) {
682  /* the SW01, 02, 04 are at 1 km res for the qkm dataset
683  get # pixels for this */
684  if (ibnd == 0) {
685  npix_1km = dim_siz[1];
686  } else {
687  if (dim_siz[1] != npix_1km) {
688  fprintf(stderr, "%s, %d, E - image data size mismatch between\n",
689  __FILE__, __LINE__);
690  fprintf(stderr, "VIS/NIR file: %s\nand SWIR/IR file: %s\n",
691  vis_nir_file, swir_ir_file);
692  return 1;
693  }
694  }
695  } else {
696  if ((dim_siz[1] != npix) || (dim_siz[0] != nscan)) {
697  fprintf(stderr,
698  "%s, %d, E - image data size mismatch between\n",
699  __FILE__, __LINE__);
700  fprintf(stderr, "VIS/NIR file: %s\nand SWIR/IR file: %s\n",
701  vis_nir_file, swir_ir_file);
702  return 1;
703  }
704  }
705  } /* end swir/ir radiance work */
706  /* also check the geoloc tie datasets to match */
707  strcpy(ds_nam, "Geometry_data/Latitude");
708  if (h5io_set_ds(&swir_ir_fid, ds_nam, &gen_dsid) != 0) {
709  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
710  __LINE__, ds_nam);
711  return 1;
712  }
713  if (h5io_info(&gen_dsid, NULL, &h5_class, &h5_native_typ,
714  &ndim, dim_siz, &sto_len) != 0) {
715  fprintf(stderr, "%s, %d, E - h5io_info failure\n", __FILE__,
716  __LINE__);
717  return 1;
718  }
719  if ((npix_tie != dim_siz[1]) || (nscan_tie != dim_siz[0])) {
720  fprintf(stderr,
721  "%s, %d, E - tie point data size mismatch between VIS/NIR file:\n",
722  __FILE__, __LINE__);
723  fprintf(stderr, "VIS/NIR file: %s\nand SWIR/IR file: %s\n",
724  vis_nir_file, swir_ir_file);
725  return 1;
726  }
727  h5io_close(&gen_dsid);
728 
729  /* if the geom per band is used, open all the sena, senz for swir bands */
730  if (l1_input->geom_per_band == 1) {
731  for (ibnd = 0; ibnd < n_swir; ibnd++) {
732  for (ip = 0; ip < 2; ip++) {
733  sprintf(ds_nam, "Geometry_data/%s_SW%s",
734  sen_ang[ip], swir_bnam[ibnd]);
735  if (h5io_set_ds(&swir_ir_fid, ds_nam,
736  &band_geom_swir[ibnd].dsid[ip]) != 0) {
737  fprintf(stderr, "%s, %d, E - Cannot open dataset: %s\n", __FILE__,
738  __LINE__, ds_nam);
739  return 1;
740  }
741  if ((h5io_rd_attr(&band_geom_swir[ibnd].dsid[ip], "Slope",
742  (void *) (&band_geom_swir[ibnd].scale[ip])) != 0) ||
743  (h5io_rd_attr(&band_geom_swir[ibnd].dsid[ip], "Offset",
744  (void *) (&band_geom_swir[ibnd].offset[ip])) != 0)) {
745  fprintf(stderr,
746  "%s, %d, E - Failed to read scale or offset for geo data\n",
747  __FILE__, __LINE__);
748  return 1;
749  }
750  /* check size of the tie point data */
751  if (h5io_info(&band_geom_swir[ibnd].dsid[ip], NULL, &h5_class,
752  &h5_native_typ, &ndim, dim_siz, &sto_len) != 0) {
753  fprintf(stderr,
754  "%s, %d, E - h5io_info failure for swir geom_per_band\n",
755  __FILE__, __LINE__);
756  return 1;
757  }
758  if ((resqkm == 0) || (ibnd == 2) ) {
759  if ((dim_siz[1] != npix_tie) || (dim_siz[0] != nscan_tie)) {
760  fprintf(stderr,
761  "%s, %d, E - Unexpected swir geom_per_band sizes found\n",
762  __FILE__, __LINE__);
763  return 1;
764  }
765  band_geom_swir[ibnd].grd_desc = grid_res;
766  }
767  else {
768  if( ( ibnd == 0 ) && ( ip == 0 ) ) {
769  npix_tie_1km = dim_siz[1];
770  nscan_tie_1km = dim_siz[0];
771  resamp_1km = resample * 4;
772  } else {
773  if ((dim_siz[1] != npix_tie_1km) || (dim_siz[0] != nscan_tie_1km)) {
774  fprintf(stderr,
775  "%s, %d, E - Unexpected swir geom_per_band sizes found\n",
776  __FILE__, __LINE__);
777  return 1;
778  }
779  }
780  band_geom_swir[ibnd].grd_desc = ( grid_res + 1 );
781  }
782  }
783  }
784  }
785  }
786  }
787  return (0);
788 }
789 
790 /*
791  readl1_sgli - read a line of sgli data in
792  */
793 int readl1_sgli(filehandle *file, int32_t scan, l1str *l1rec) {
794  int32_t ip_tie, il_tie, sscan, escan, tie_st_line, tie_st_lin_1km;
795  int32_t il, ipl, ip_1km, ip_off, scan_1km, isub, scan_tie, nscn_sub_tie;
796  uint16_t *ui16_buf = (uint16_t *) flt_buf;
797  int32_t start[2], count[2], ib, ip, ipb, ib_lg, ilr, iv;
798  int32_t lt_lg_state, lt_hg_state;
799  double lat_cvt, lon_cvt;
800  double xvec, yvec, zvec, xy_mod, *geo[3];
801  double rad2deg(double), deg2rad(double);
802  char nav_bad;
803  float zen_v, azi_v;
804  static int32_t firstcall = 1, cur_scan_1km[3] = {-1, -1, -1};
805  static int32_t tie_en_line, tie_en_lin_1km;
806 
807  enum lt_state {
808  LT_BAD, LT_SAT, LT_NORM
809  };
810  int32_t nbands = l1rec->l1file->nbands;
811  uint16_t lt_strip = pow(2, 14) - 1;
812  float w_m2_to_mw_cm2 = 0.1; /* to store initial W m^-2 as mW cm^-2 */
813  float lt_tmp;
814 
815  /*
816  * do some setup in the first call
817  */
818  /*
819  * if required for that record, set up the geom_per_band storage
820  */
821  if ((l1_input->geom_per_band == 1) && (l1rec->geom_per_band == NULL)) {
822  init_geom_per_band(l1rec);
823  }
824 
825  if (firstcall == 1) {
826  /* Set up for interpolating the geo and view information that is stored
827  in tie point format at resample rate: resample
828  Only set up for the lines to be processed and always have bounding
829  tie point lines so as not to extrapolate.
830  This is especially important so pad out the tie lines */
831  sscan = l1_input->sline - 1;
832  escan = l1_input->eline - 1;
833  tie_st_line = sscan / resample - 1; /* pad out start tie line */
834  tie_en_line = escan / resample + 2; /* pad out end tie line */
835 
836  tie_st_line = (tie_st_line < 0) ? 0 : tie_st_line;
837  tie_en_line = (tie_en_line >= nscan_tie) ? nscan_tie - 1 : tie_en_line;
838  /*
839  * there needs to be 4 lines of tie point data for the spline
840  * to work - this will assure that there are 4
841  */
842  if (tie_en_line - tie_st_line < 3) {
843  if (tie_en_line >= 3)
844  tie_st_line = tie_en_line - 3;
845  else
846  tie_en_line = tie_st_line + 3;
847  }
848  if ((tie_st_line < 0) || (tie_en_line > nscan_tie - 1)) {
849  fprintf(stderr, "%s, %d, E - L1 file has fewer than 4 tie point lines\n",
850  __FILE__, __LINE__);
851  fprintf(stderr, " Spline interpolation is not possible\n");
852  /* in the unlikely event of NEEDING to process such short files,
853  * a linear interpolation will need to be added as a fallback option
854  */
855  return 1;
856  }
857  /* set up the 1st grid description */
858  nscn_sub_tie = tie_en_line - tie_st_line + 1;
859  xa = (double *) malloc(npix_tie * sizeof ( double));
860  ya = (double *) malloc(nscn_sub_tie * sizeof ( double));
861  grid_res[0].nscn_sub_tie = nscn_sub_tie;
862  grid_res[0].xa = xa;
863  grid_res[0].ya = ya;
864  grid_res[0].npix_tie = npix_tie;
865  grid_res[0].tie_st_lin = tie_st_line;
866  grid_res[0].resamp = resample;
867 
868  /* First, latitude, longitude - read them in */
869  start[0] = tie_st_line;
870  start[1] = 0;
871  count[0] = nscn_sub_tie;
872  count[1] = npix_tie;
873 
874  if (((tie_el = (float *) malloc(npix_tie * nscn_sub_tie * sizeof ( float)))
875  == NULL) ||
876  ((tie_az = (float *) malloc(npix_tie * nscn_sub_tie * sizeof ( float)))
877  == NULL) ||
878  ((geo_x = (double *) malloc(npix_tie * nscn_sub_tie * sizeof ( double)))
879  == NULL) ||
880  ((geo_y = (double *) malloc(npix_tie * nscn_sub_tie * sizeof ( double)))
881  == NULL) ||
882  ((geo_z = (double *) malloc(npix_tie * nscn_sub_tie * sizeof ( double)))
883  == NULL) ||
884  ((geo_q = (char *) calloc(npix_tie * nscn_sub_tie, sizeof ( char)))
885  == NULL)) {
886  fprintf(stderr, "%s, %d, E - Failed to allocate lat, lon storage\n",
887  __FILE__, __LINE__);
888  return 1;
889  }
890 
891  if (h5io_rd_ds_slice(vnir_geom_dsid, start, count, (void *) tie_el)
892  != 0) {
893  fprintf(stderr, "%s, %d, E Error reading latitude tie points\n",
894  __FILE__, __LINE__);
895  return 1;
896  }
897  if (h5io_rd_ds_slice((vnir_geom_dsid + 1), start, count, (void *) tie_az)
898  != 0) {
899  fprintf(stderr, "%s, %d, E Error reading longitude tie points\n",
900  __FILE__, __LINE__);
901  return 1;
902  }
903 
904  /* convert to x,y,z unit vectors */
905  for (il = 0; il < nscn_sub_tie; il++) {
906  ya[il] = (il + start[0]) * resample;
907  for (ip = 0; ip < npix_tie; ip++) {
908  if (il == 0) xa[ip] = ip * resample;
909  ipl = ip + il * npix_tie;
910  /* lat, lon not scaled
911  lat_cvt = vnir_geom_off[0] + vnir_geom_scale[0] * tie_el[ipl];
912  lon_cvt = vnir_geom_off[1] + vnir_geom_scale[1] * tie_az[ipl];
913  */
914  lat_cvt = tie_el[ipl];
915  lon_cvt = tie_az[ipl];
916  if ((lat_cvt < -90.) || (lat_cvt > 90.) ||
917  (lon_cvt < -180.) || (lon_cvt > 180.)) {
918  lat_cvt = 0;
919  lon_cvt = 0;
920  geo_q[ipl] = 1;
921  }
922  xy_mod = cos(deg2rad(lat_cvt));
923  geo_x[ipl] = xy_mod * cos(deg2rad(lon_cvt));
924  geo_y[ipl] = xy_mod * sin(deg2rad(lon_cvt));
925  geo_z[ipl] = sin(deg2rad(lat_cvt));
926  }
927  }
928  geo[0] = geo_x;
929  geo[1] = geo_y;
930  geo[2] = geo_z;
931  /* set up for the 2D spline interpolation and accelerator */
932  for (ip = 0; ip < 3; ip++) {
933  geo_int_id[ip] =
934  gsl_spline2d_alloc(gsl_interp2d_bicubic, npix_tie, nscn_sub_tie);
935  if ((gsl_spline2d_init(geo_int_id[ip], xa, ya, geo[ip], npix_tie,
936  nscn_sub_tie)) != 0) {
937  fprintf(stderr, "%s, %d, E gsl_spline2d_init error\n", __FILE__,
938  __LINE__);
939  return 1;
940  }
941  }
942  accel_x = gsl_interp_accel_alloc();
943  accel_y = gsl_interp_accel_alloc();
944  /*
945  * set up sensor and solar zenith and azimuth
946  */
947  if (view_ang_tie_init2(band_geom_nominal, 2 ) != 0)
948  return 1;
949  /*
950  * include the band-dependent geometry, if requested
951  */
952  if (l1_input->geom_per_band == 1) {
953  /* set up the spline coefficients for band dependent vnir */
954  if (view_ang_tie_init2(band_geom_vnir, n_vnir) != 0)
955  return 1;
956  /*
957  * set up the spline coefficients for band dependent vnir low gain
958  */
959  if (view_ang_tie_init2(band_geom_vnir_lg, n_vnir_lg ) != 0)
960  return 1;
961  }
962 
963  /*
964  * for the SWIR, set up various things
965  * for replication, we only need 1 line, but need to remember it
966  * for 3/4 swir bands
967  */
968  if (swirir_exist == 1) {
969  /*
970  * prepare for handling the 1km SWIR in a qkm file -
971  * for replication, we only need 1 line, but need to remember it
972  * for 3/4 swir bands
973  */
974  if (resqkm == 1) {
975  if (((ui16_buf_1km = (uint16_t *) malloc(npix_1km * 3 *
976  sizeof ( uint16_t))) == NULL) ||
977  ((swir_q = (char *) calloc(npix_1km * 3, sizeof ( char)))
978  == NULL)) {
979  fprintf(stderr, "%s, %d, E - Failed to allocate lat, lon storage\n",
980  __FILE__, __LINE__);
981  return 1;
982  }
983  }
984  /*
985  * if band-dependent sensor angles needed, set up the interpolation
986  * for this
987  */
988  if (l1_input->geom_per_band == 1) {
989  /* WDR move to after the if( resqkm == 1 ...
990  and set-up of the grid_res[1] values???
991  if (view_ang_tie_init2(band_geom_swir, n_swir ) != 0)
992  return 1;
993 */
994  /* set up the grid information for the tie point grids that
995  are used for the 1 km data */
996  if( resqkm == 1 ) {
997  tie_st_lin_1km = sscan / resamp_1km - 1; /* pad out start tie line */
998  tie_en_lin_1km = escan / resamp_1km + 2; /* pad out end tie line */
999 
1000  tie_st_lin_1km = (tie_st_lin_1km < 0) ? 0 : tie_st_lin_1km;
1001  tie_en_lin_1km = (tie_en_lin_1km >= nscan_tie_1km) ? nscan_tie_1km - 1 : tie_en_lin_1km;
1002  /*
1003  * there needs to be 4 lines of tie point data for the spline
1004  * to work - this will assure that there are 4
1005  */
1006  if (tie_en_lin_1km - tie_st_lin_1km < 3) {
1007  if (tie_en_lin_1km >= 3)
1008  tie_st_lin_1km = tie_en_lin_1km - 3;
1009  else
1010  tie_en_lin_1km = tie_st_lin_1km + 3;
1011  }
1012  if ((tie_st_lin_1km < 0) || (tie_en_lin_1km > nscan_tie_1km - 1)) {
1013  fprintf(stderr, "%s, %d, E - L1 file has fewer than 4 tie point lines\n",
1014  __FILE__, __LINE__);
1015  fprintf(stderr, " Spline interpolation is not possible\n");
1016  /* in the unlikely event of NEEDING to process such short files,
1017  * a linear interpolation will need to be added as a fallback option
1018  */
1019  return 1;
1020  }
1021  /* set up the 2nd grid description */
1022  grid_res[1].nscn_sub_tie = tie_en_lin_1km - tie_st_lin_1km + 1;
1023  xa_1km = (double *) malloc(npix_tie_1km * sizeof(double) );
1024  for (ip = 0; ip < npix_tie_1km; ip++)
1025  xa_1km[ip] = ip * resamp_1km;
1026  ya_1km = (double *) malloc(grid_res[1].nscn_sub_tie
1027  * sizeof(double) );
1028  for (il = 0; il < grid_res[1].nscn_sub_tie; il++)
1029  ya_1km[il] = (il + tie_st_lin_1km) * resamp_1km;
1030  grid_res[1].xa = xa_1km;
1031  grid_res[1].ya = ya_1km;
1032  grid_res[1].npix_tie = npix_tie_1km;
1033  grid_res[1].tie_st_lin = tie_st_lin_1km;
1034  grid_res[1].resamp = resamp_1km;
1035  }
1036  if (view_ang_tie_init2(band_geom_swir, n_swir ) != 0)
1037  return 1;
1038  }
1039  }
1040  firstcall = 0;
1041  }
1042  /*
1043  * the use of the last scan line in gsl_interp_accel_find causes problems,
1044  * so, make this adjustment. scan_tie is modified scan so its not at
1045  * the end line of the tie point grid
1046  */
1047  scan_tie = (scan >= (resample * tie_en_line)) ?
1048  resample * tie_en_line - 1 : scan;
1049 
1050  /* get the VIS/NIR read and scaled */
1051  for (ib = 0; ib < n_vnir; ib++) {
1052  /* read in the scan line of Lt data */
1053  start[0] = scan;
1054  start[1] = 0;
1055  count[0] = 1;
1056  count[1] = npix;
1057  if (h5io_rd_ds_slice((vnir_hg_dsid + ib), start, count,
1058  (void *) ui16_buf) != 0) {
1059  fprintf(stderr, "%s, %d, E Error reading VIS/NIR hg, index %d\n",
1060  __FILE__, __LINE__, ib);
1061  return 1;
1062  }
1063  /* for data that is within the valid DN, scale it to mw m^-2...
1064  and if above the saturation radiance, set the hilt */
1065  for (ip = 0; ip < npix; ip++) {
1066  ipb = ip * nbands + ib;
1067  if (*(ui16_buf + ip) < *(vnir_mx_dn + ib)) {
1068  /* They say the 1st 2 bits are for stray light and sign of corr
1069  but examples have none of this. For now, strip off top 2 bits
1070  when applying the scale, offset */
1071  lt_tmp = *(vnir_off + ib) + *(vnir_scale + ib) *
1072  (*(ui16_buf + ip) & lt_strip);
1073  if ((lt_tmp >= *(vnir_sat + ib)) &&
1074  (ib != 6) && (ib != 8))
1075  l1rec->hilt[ip] = 1;
1076  l1rec->Lt[ipb] = lt_tmp * w_m2_to_mw_cm2;
1077  }
1078  }
1079  }
1080  /*
1081  * nominal view geometry
1082  */
1083  for (iv = 0; iv < 2; iv++) {
1084  for (ip = 0; ip < npix; ip++) {
1085  if (view_ang_tie_eval(&band_geom_nominal[iv], ip, scan,
1086  &azi_v, &zen_v, &nav_bad) != 0) return 1;
1087  if (iv == 0) {
1088  l1rec->senz[ip] = zen_v;
1089  l1rec->sena[ip] = azi_v;
1090  } else {
1091  l1rec->solz[ip] = zen_v;
1092  l1rec->sola[ip] = azi_v;
1093  }
1094  if (nav_bad != 0) l1rec->navfail[ip] = 1;
1095  }
1096  }
1097  /* high gain band-dependent sena */
1098  /* and solar band-dependent values (just mirror the nominal) */
1099  if (l1_input->geom_per_band == 1) {
1100  for (ib = 0; ib < n_vnir; ib++) {
1101  /*
1102  * get interpolated view vectors and azimuth, zenith
1103  */
1104  for (ip = 0; ip < npix; ip++) {
1105  ipl = ip * nbands + ib;
1106  //printf( "ib = %d, ip = %d\n", ib, ip );
1107  if (view_ang_tie_eval(&band_geom_vnir[ib], ip, scan,
1108  &azi_v, &zen_v, &nav_bad) != 0) return 1;
1109  l1rec->geom_per_band->senz[ipl] = zen_v;
1110  l1rec->geom_per_band->sena[ipl] = azi_v;
1111  if (nav_bad != 0) l1rec->navfail[ip] = 1;
1112 
1113  l1rec->geom_per_band->solz[ipl] = l1rec->solz[ip];
1114  l1rec->geom_per_band->sola[ipl] = l1rec->sola[ip];
1115  }
1116  }
1117  }
1118  /* for the low gain bands, fill any saturated values with low gain part */
1119  for (ib = 0; ib < n_vnir_lg; ib++) {
1120  ib_lg = (ib == 0) ? 6 : 8;
1121  if (h5io_rd_ds_slice((vnir_lg_dsid + ib), start, count,
1122  (void *) ui16_buf) != 0) {
1123  fprintf(stderr, "%s, %d, E Error reading VIS/NIR hg, index %d\n",
1124  __FILE__, __LINE__, ib);
1125  return 1;
1126  }
1127  /* for data that is within the valid DN, scale it to mw m^-2...
1128  and if above the saturation radiance, set the hilt */
1129 
1130  /* logic to get the correct data between the low-and high-gain bands */
1131  /* NOTE that due to observed non-linearities near saturation in the
1132  high gain, the switch-over to low agin data use is set to .75 of
1133  the high gain saturation - ...* 0.75... */
1134  for (ip = 0; ip < npix; ip++) {
1135  /* set status of the Lt for the high gain data */
1136  ipb = ip * nbands + ib_lg;
1137  if (l1rec->Lt[ipb] == BAD_FLT)
1138  lt_hg_state = LT_BAD;
1139  else {
1140  lt_hg_state = LT_NORM;
1141  if (l1rec->Lt[ipb] >= *(vnir_sat + ib_lg) * 0.75 *
1142  w_m2_to_mw_cm2)
1143  lt_hg_state = LT_SAT;
1144  }
1145  /* if the hg is normal, nothing left, otherwise check the low gain */
1146  if (lt_hg_state != LT_NORM) {
1147  /* get lg rad and lg rad state */
1148  if (*(ui16_buf + ip) < *(vnir_lg_mx_dn + ib)) {
1149  lt_lg_state = LT_NORM;
1150  lt_tmp = *(vnir_lg_off + ib) + *(vnir_lg_scale + ib) *
1151  (*(ui16_buf + ip) & lt_strip);
1152  if (lt_tmp >= *(vnir_lg_sat + ib))
1153  lt_lg_state = LT_SAT;
1154  lt_tmp *= w_m2_to_mw_cm2;
1155  /* get the low gain sensor zenith, azimuth - may be needed */
1156  if (l1_input->geom_per_band == 1) {
1157  if (view_ang_tie_eval(&band_geom_vnir_lg[ib], ip, scan,
1158  &azi_v, &zen_v, &nav_bad) != 0) return 1;
1159  }
1160  } else
1161  lt_lg_state = LT_BAD;
1162  /* go through hg, lg states to set Lt, detector used, and hilt */
1163  /* setting. if we have band-dependent view angles, make sure
1164  the lg angles are used */
1165  if (lt_lg_state == LT_NORM) {
1166  l1rec->Lt[ipb] = lt_tmp;
1167  if (l1_input->geom_per_band == 1) {
1168  l1rec->geom_per_band->sena[ipb] = azi_v;
1169  l1rec->geom_per_band->senz[ipb] = zen_v;
1170  l1rec->navfail[ip] = nav_bad;
1171  }
1172  } else if (lt_lg_state == LT_SAT) {
1173  l1rec->Lt[ipb] = lt_tmp;
1174  l1rec->hilt[ip] = 1;
1175  if (l1_input->geom_per_band == 1) {
1176  l1rec->geom_per_band->sena[ipb] = azi_v;
1177  l1rec->geom_per_band->senz[ipb] = zen_v;
1178  l1rec->navfail[ip] = nav_bad;
1179  }
1180  } else if (lt_hg_state == LT_SAT) {
1181  /* lg is bad, so stick with hg value */
1182  l1rec->hilt[ip] = 1;
1183  /* use hg detector position - already there as default */
1184  } else {
1185  /* both Lt bad, use bad (already there) and
1186  hg detector position - already there as default */
1187  }
1188  }
1189  }
1190  /* end logic for hg and lg bands */
1191  }
1192  /* set the time for the line */
1193  l1rec->scantime = st_day_unix + (double) *(scan_times + scan) / 1.e3;
1194  /* set the geolocation and view geometry */
1195  for (ip = 0; ip < npix; ip++) {
1196  if ((gsl_spline2d_eval_e(geo_int_id[0], ip, scan, accel_x,
1197  accel_y, &xvec) != 0) ||
1198  (gsl_spline2d_eval_e(geo_int_id[1], ip, scan, accel_x,
1199  accel_y, &yvec) != 0) ||
1200  (gsl_spline2d_eval_e(geo_int_id[2], ip, scan, accel_x,
1201  accel_y, &zvec) != 0)) {
1202  fprintf(stderr, "%s, %d, E Error in gsl_spline2d_eval\n",
1203  __FILE__, __LINE__);
1204  return 1;
1205  }
1206  l1rec->lat[ip] = (float) rad2deg(asin(zvec));
1207  l1rec->lon[ip] = (float) rad2deg(atan2(yvec, xvec));
1208  /* it is possible that tie point nav is bad, so */
1209  ip_tie = gsl_interp_accel_find(accel_x, xa, npix_tie, ip);
1210  il_tie = gsl_interp_accel_find(accel_y, ya, nscan_tie, scan_tie);
1211  ipb = ip_tie + il_tie * npix_tie;
1212  if ((geo_q[ipb] != 0) || (geo_q[ ipb + 1 ] != 0) ||
1213  (geo_q[ ipb + npix_tie ] != 0) ||
1214  (geo_q[ ipb + 1 + npix_tie ] != 0))
1215  l1rec->navfail[ip] = 1;
1216  }
1217  /*
1218  * read the SWIR bands, if available
1219  * this will do replication for the 1 km bands in a qkm file
1220  */
1221  if (swirir_exist == 1) {
1222  for (ib = 0; ib < n_swir; ib++) {
1223  if ((resqkm == 1) && (ib != 2)) {
1224  /* replicate the 1km data to the qkm size */
1225  ilr = ib;
1226  if (ib == 3) ilr = 2;
1227  scan_1km = scan / 4;
1228  if (cur_scan_1km[ilr] != scan_1km) {
1229  start[0] = scan_1km;
1230  start[1] = 0;
1231  count[0] = 1;
1232  count[1] = npix_1km;
1233  if ((h5io_rd_ds_slice((swir_dsid + ib), start, count,
1234  (void *) (ui16_buf_1km + ilr * npix_1km))) != 0) {
1235  fprintf(stderr, "%s, %d, E Error reading SWIR, band index %d\n",
1236  __FILE__, __LINE__, ib);
1237  return 1;
1238  }
1239  cur_scan_1km[ilr] = scan_1km;
1240  }
1241  for (ip_1km = 0; ip_1km < npix_1km; ip_1km++) {
1242  ip = ip_1km * 4;
1243  ipb = ip * nbands + (ib + n_vnir);
1244  ip_off = ip_1km + npix_1km * ilr;
1245  if (*(ui16_buf_1km + ip_off) < *(swir_mx_dn + ib)) {
1246  lt_tmp = *(swir_off + ib) + *(swir_scale + ib) *
1247  (*(ui16_buf_1km + ip_off) & lt_strip);
1248  if (lt_tmp >= *(swir_sat + ib))
1249  for (isub = 0; isub < 4; isub++)
1250  l1rec->hilt[ ip + isub ] = 1;
1251  lt_tmp *= w_m2_to_mw_cm2;
1252  for (isub = 0; isub < 4; isub++)
1253  l1rec->Lt[ipb + isub * nbands] = lt_tmp;
1254  }
1255  }
1256  } else {
1257  /* deal with bands at all the same resolution */
1258  start[0] = scan;
1259  start[1] = 0;
1260  count[0] = 1;
1261  count[1] = npix;
1262  if (h5io_rd_ds_slice((swir_dsid + ib), start, count,
1263  (void *) ui16_buf) != 0) {
1264  fprintf(stderr, "%s, %d, E Error reading VIS/NIR hg, index %d\n",
1265  __FILE__, __LINE__, ib);
1266  return 1;
1267  }
1268  /* check for valid data in band, note initialized to BAD_FLT */
1269  for (ip = 0; ip < npix; ip++) {
1270  ipb = ip * nbands + (ib + n_vnir);
1271  if (*(ui16_buf + ip) <= *(swir_mx_dn + ib)) {
1272  lt_tmp = *(swir_off + ib) + *(swir_scale + ib) *
1273  (*(ui16_buf + ip) & lt_strip);
1274  if (lt_tmp >= *(swir_sat + ib))
1275  l1rec->hilt[ip] = 1;
1276  l1rec->Lt[ipb] = lt_tmp * w_m2_to_mw_cm2;
1277  }
1278  }
1279  }
1280  }
1281  /*
1282  * get the geom_per_band for the swir, note that all sena/z arrays
1283  * are relative to the 1/4 km resolution data and so will not need
1284  * the same treatment as the Lt did above for 3/4 bands
1285  */
1286  if (l1_input->geom_per_band == 1) {
1287  for (ib = 0; ib < n_swir; ib++) {
1288  /*
1289  * get interpolated view vectors and azimuth, zenith
1290  * As the resolution of the grids changes for qkm SWIR
1291  * in bands 0, 2, and 3 and then back after, do these
1292  * calls to gsl_interp_accel_reset
1293  */
1294  if ( ( resqkm == 1 ) && ( ib != 1 ) ) {
1295  gsl_interp_accel_reset( accel_x );
1296  gsl_interp_accel_reset( accel_y );
1297  }
1298  for (ip = 0; ip < npix; ip++) {
1299  ipb = ip * nbands + (ib + n_vnir);
1300  if (view_ang_tie_eval(&band_geom_swir[ib], ip, scan,
1301  &azi_v, &zen_v, &nav_bad) != 0) return 1;
1302  l1rec->geom_per_band->senz[ipb] = zen_v;
1303  l1rec->geom_per_band->sena[ipb] = azi_v;
1304  if (nav_bad != 0) l1rec->navfail[ip] = 1;
1305  /* solar values, though constant are also per-band */
1306  l1rec->geom_per_band->solz[ipb] = l1rec->solz[ip];
1307  l1rec->geom_per_band->sola[ipb] = l1rec->sola[ip];
1308  }
1309  }
1310  if ( resqkm == 1 ) {
1311  gsl_interp_accel_reset( accel_x );
1312  gsl_interp_accel_reset( accel_y );
1313  }
1314  }
1315  } else {
1316  /* if geom per band is used but no swir is available, fill
1317  view geom for swir bands with the nominal values */
1318  if (l1_input->geom_per_band == 1) {
1319  for (ib = 0; ib < n_swir; ib++) {
1320  for (ip = 0; ip < npix; ip++) {
1321  ipb = ip * nbands + (ib + n_vnir);
1322  l1rec->geom_per_band->senz[ipb] = l1rec->senz[ip];
1323  l1rec->geom_per_band->sena[ipb] = l1rec->sena[ip];
1324  l1rec->geom_per_band->solz[ipb] = l1rec->solz[ip];
1325  l1rec->geom_per_band->sola[ipb] = l1rec->sola[ip];
1326  }
1327  }
1328  }
1329  }
1330 
1331  return (0);
1332 }
1333 
1334 int closel1_sgli(filehandle *file) {
1335  int32_t ibnd, iv;
1336  sgli_t *data = file->private_data;
1337 
1338  /*
1339  * free allocated space, close all IDs and files
1340  */
1341  free(data);
1342  free(flt_buf);
1343  free(tie_el);
1344  free(tie_az);
1345  free(geo_x);
1346  free(geo_y);
1347  free(geo_z);
1348  free(geo_q);
1349  free(sen_q);
1350  free(sol_q);
1351  free(xa);
1352  free(ya);
1353 
1354  for (ibnd = 0; ibnd < n_vnir; ibnd++)
1355  h5io_close(vnir_hg_dsid + ibnd);
1356  for (ibnd = 0; ibnd < n_vnir_lg; ibnd++)
1357  h5io_close(vnir_lg_dsid + ibnd);
1358  for (ibnd = 0; ibnd < n_swir; ibnd++)
1359  h5io_close(swir_dsid + ibnd);
1360  for (ibnd = 0; ibnd < n_vnir_geom; ibnd++)
1361  h5io_close(vnir_geom_dsid + ibnd);
1362 
1363  h5io_close(&vnir_fid);
1364 
1365  for (ibnd = 0; ibnd < 3; ibnd++)
1366  gsl_spline2d_free(geo_int_id[ibnd]);
1367  /*
1368  * free splines for nominal sen, sol angles
1369  */
1370  for (ibnd = 0; ibnd < 2; ibnd++)
1371  for (iv = 0; iv < 3; iv++)
1372  gsl_spline2d_free(band_geom_nominal[ibnd].int_id_sen[iv]);
1373  /*
1374  * free splines for band-dependent sensor angles
1375  */
1376  if (l1_input->geom_per_band == 1) {
1377  for (ibnd = 0; ibnd < n_vnir; ibnd++)
1378  for (iv = 0; iv < 3; iv++)
1379  gsl_spline2d_free(band_geom_vnir[ibnd].int_id_sen[iv]);
1380  for (ibnd = 0; ibnd < n_vnir_lg; ibnd++)
1381  for (iv = 0; iv < 3; iv++)
1382  gsl_spline2d_free(band_geom_vnir_lg[ibnd].int_id_sen[iv]);
1383  }
1384  /*
1385  * for the SWIR/IR file, free allocated space, close all IDs and files
1386  */
1387  if (swirir_exist) {
1388  if (resqkm == 1) {
1389  free(xa_1km);
1390  free(ya_1km);
1391  free(ui16_buf_1km);
1392  }
1393 
1394  if (l1_input->geom_per_band == 1)
1395  for (ibnd = 0; ibnd < n_swir; ibnd++)
1396  for (iv = 0; iv < 3; iv++)
1397  gsl_spline2d_free(band_geom_swir[ibnd].int_id_sen[iv]);
1398 
1399  for (ibnd = 0; ibnd < n_swir; ibnd++)
1400  h5io_close(swir_dsid + ibnd);
1401 
1402  h5io_close(&swir_ir_fid);
1403  }
1404  return (0);
1405 }
1406 
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
init_geom_per_band
int init_geom_per_band(l1str *l1rec)
Definition: geom_per_band.c:7
closel1_sgli
int closel1_sgli(filehandle *file)
Definition: l1_sgli.c:1334
h5io_openr
int h5io_openr(char *file, int opt, h5io_str *id)
Definition: h5io.c:4
h5io_attr_exist
int h5io_attr_exist(h5io_str *id, char *attr_name)
Definition: h5io.c:1451
view_ang_tie_eval
int view_ang_tie_eval(band_geom_str *band_geom, int32_t ip, int32_t scan, float *azi_v, float *zen_v, char *qual)
Definition: l1_sgli.c:304
SGLI_NVNIR_GEOM
#define SGLI_NVNIR_GEOM
Definition: l1_sgli.c:15
qual
int16_t * qual
Definition: l2bin.cpp:86
NULL
#define NULL
Definition: decode_rs.h:63
readl1_sgli
int readl1_sgli(filehandle *file, int32_t scan, l1str *l1rec)
Definition: l1_sgli.c:793
scan
MOD_PR01 Production producing one five minute granule of output data in each run It can be configured to produce as many as three five minute granules per run Each execution with one construction record and one date file for each dataset In normal these are created by which splits them out of the hour datasets For LANCE they are created by which merges all session MODIS L0 datasets overlapping the requested time and extracts from the merged data those packets which fall within that time period Each scan of data is stored in the L1A granule that covers the start time of that scan
Definition: MOD_PR01_pr.txt:19
resolution
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that resolution
Definition: HISTORY.txt:188
SGLI_NSWIR
#define SGLI_NSWIR
Definition: l1_sgli.c:16
deg2rad
const double deg2rad
Definition: VcstCmnConsts.h:478
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
h5io_info
int h5io_info(h5io_str *id, char *attr_name, H5T_class_t *class, hid_t *native_typ, int *ndim, int *dim_siz, int *sto_len)
Definition: h5io.c:173
createPrivateData_sgli
sgli_t * createPrivateData_sgli()
Definition: l1_sgli.c:50
l1_sgli.h
h5io_close
int h5io_close(h5io_str *id)
Definition: h5io.c:115
nscan
int32 nscan
Definition: l1_czcs_hdf.c:19
SGLI_NVNIR
#define SGLI_NVNIR
Definition: l1_sgli.c:13
file
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed file
Definition: HISTORY.txt:413
h5io_set_grp
int h5io_set_grp(h5io_str *id, char *path_name, h5io_str *grp_id)
Definition: h5io.c:369
h5io_set_ds
int h5io_set_ds(h5io_str *id, char *path_name, h5io_str *ds_id)
Definition: h5io.c:324
openl1_sgli
int openl1_sgli(filehandle *l1file)
Definition: l1_sgli.c:366
h5io_rd_ds_slice
int h5io_rd_ds_slice(h5io_str *ds_id, int *start, int *count, void *data)
Definition: h5io.c:602
h5io_rd_attr
int h5io_rd_attr(h5io_str *id, char *attr_name, void *data)
Definition: h5io.c:412
sgli_rad_info
int sgli_rad_info(h5io_str *fid, char *ds_nam, h5io_str *dsid, int *dim_siz, uint16_t *mx_dn, float *scale, float *off, float *sat)
Definition: l1_sgli.c:131
color_dtdb.scale
int scale
Definition: color_dtdb.py:459
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
data
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
view_ang_tie_init2
int view_ang_tie_init2(band_geom_str *band_geom, int32_t nbands)
Definition: l1_sgli.c:197
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
sgli_file_ver
int sgli_file_ver(h5io_str *fid)
Definition: l1_sgli.c:65
ymds2unix
double ymds2unix(short year, short month, short day, double secs)
for
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
offset
l2prod offset
Definition: HOWTO_Add_a_product.txt:70
SGLI_NVNIR_LG
#define SGLI_NVNIR_LG
Definition: l1_sgli.c:14
rad2deg
const double rad2deg
Definition: VcstCmnConsts.h:479
h5io_rd_ds
int h5io_rd_ds(h5io_str *ds_id, void *data)
Definition: h5io.c:505
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
npix
int npix
Definition: get_cmp.c:27
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