ocssw V2020
l1_hmodis_hdf.c
Go to the documentation of this file.
1 /* ============================================================================ */
2 /* module l1_hmodis_hdf.c - functions to read MODIS HIRES L1B for MSL12 */
3 /* */
4 /* Written By: B. Franz, NASA/SIMBIOS, January 2003. */
5 /* Modified By: J. Gales, NASA/OBPG, August 2005. */
6 /* Conversion for 16-band MODIS: B. Franz, J. Gales, 2006. */
7 /* Completely restructured by G. Fireman in 2014. */
8 /* */
9 /* ============================================================================ */
10 
11 #define _XOPEN_SOURCE_EXTENDED /* for strdup() */
12 #include "libnav.h"
13 #include "hdf4utils.h"
14 #include "l12_proto.h"
15 #undef _XOPEN_SOURCE_EXTENDED
16 
17 #define SCAN_TIME_INTERVAL 1.4771
18 
19 /***********************************************************************/
20 #define TRYMEM(file,line,memstat) { \
21  if (memstat == NULL) { \
22  fprintf(stderr, \
23  "-E- %s line %d: Memory allocation error.\n", \
24  file, line); \
25  exit(1); } \
26  }
27 /***********************************************************************/
28 
29 /* SDSs of interest in each type of input file */
30 typedef struct {
31  int32_t index;
32  int32_t scandim;
33  char *name;
34 } sdslist;
35 
36 /* Positioning info to preload from MODIS Geolocation file */
37 enum REF_SDS {
43 };
44 static const sdslist SDSLIST_REF[] = {
45  { GEO_TAISEC, 0, "EV start time"},
46  { GEO_MSIDE, 0, "Mirror side"},
47  { GEO_ANGLES, 0, "attitude_angles"},
48  { GEO_MNORM, 0, "T_inst2ECR"},
49  { 0, 0, NULL}
50 };
51 
52 /* MODIS Geolocation swath SDSs */
53 enum GEO_SDS {
62 };
63 static const sdslist SDSLIST_GEO[] = {
64  { GEO_LON, 0, "Longitude"}, /* float32 */
65  { GEO_LAT, 0, "Latitude"}, /* float32 */
66  { GEO_HGT, 0, "Height"}, /* int16 */
67  { GEO_SOLZ, 0, "SolarZenith"}, /* int16 */
68  { GEO_SOLA, 0, "SolarAzimuth"}, /* int16 */
69  { GEO_SENZ, 0, "SensorZenith"}, /* int16 */
70  { GEO_SENA, 0, "SensorAzimuth"}, /* int16 */
71  { 0, 0, NULL}
72 };
73 
74 enum GEO_COEFFS {
77 };
78 static const char* geo_coeff[] = {/* index with enum GEO_COEFFS */
79  "scale_factor"
80 };
81 
82 /* MODIS L1B Band groupings */
83 enum L1B_SDS {
90 };
91 static sdslist SDSLIST_1KM[] = {// 1KM file contents
92  { RSB_250, 1, "EV_250_Aggr1km_RefSB"},
93  { RSB_500, 1, "EV_500_Aggr1km_RefSB"},
94  { RSB_1KM, 1, "EV_1KM_RefSB"},
95  { CIR_1KM, 0, "EV_Band26"}, // 1st dim=scan, not band
96  { TEB_1KM, 1, "EV_1KM_Emissive"},
97  { 0, 0, NULL}
98 };
99 static sdslist SDSLIST_HKM[] = {// HKM file contents
100  { RSB_250, 1, "EV_250_Aggr500_RefSB"},
101  { RSB_500, 1, "EV_500_RefSB"},
102  { 0, 0, NULL}
103 };
104 static sdslist SDSLIST_QKM[] = {// QKM file
105  { RSB_250, 1, "EV_250_RefSB"},
106  { 0, 0, NULL}
107 };
108 
115 };
116 static const char* l1b_coeff[] = {/* index with enum L1B_COEFFS */
117  "reflectance_scales",
118  "reflectance_offsets",
119  "radiance_scales",
120  "radiance_offsets"
121 };
122 
123 /***********************************************************************/
124 
125 /* Structure holding info for each MODIS input file */
126 typedef struct {
127  int32_t id;
128  char file[FILENAME_MAX];
129  char shortname[FILENAME_MAX];
130  double starttime;
131  int16_t resolution;
132  int32_t nlines;
133  int32_t npixls;
134  int32_t n_sds;
135  sds_struct *sds;
136 } modis_file;
137 
138 /***********************************************************************/
139 
140 /* Structure holding info needed for each swath-based SDS */
141 typedef struct {
142  sds_struct sds;
143  int16_t resolution;
144  int32_t nbands;
145  int32_t nscans;
146  int32_t ndets;
147  int32_t nframes;
148  int32_t nsf;
149  int32_t scandim;
150 } modis_sds;
151 
152 /*
153  Notes:
154  All GEO arrays are dimensioned [nscans*ndets, nframes]
155  Most L1B arrays are dimensioned [nbands, nscans*ndets, nframes*nsf]
156  where:
157  nscans and nframes are at 1KM resolution
158  nsf = #subframes = 1000/resolution
159  ndets = 10*nsubframes
160 
161  Band 26 drops the 1st dimension.
162  Geolocation is inherently 1KM resolution: nbands==1, ndets==10, nsf==1.
163 
164  ndets and scandim are needed for reading data one scan at a time.
165 
166  There is some redundancy between resolution, ndets, and nsf,
167  but it's helpful to have them all defined.
168  */
169 
170 /***********************************************************************/
171 /* Global variables */
172 
173 static modis_file file_geo, file_1km, file_hkm, file_qkm;
174 static modis_sds ref[REF_NUM_SDS], geo[GEO_NUM_SDS], l1b[L1B_NUM_SDS];
175 
176 /***********************************************************************/
177 
178 /* Substitute provided string for LAC/1KM part of filename */
179 int modpath_1km(const char *oldpath, const char* newchars, char* newpath) {
180  int32_t status;
181  char *ptr = NULL;
182  char *tmpbuf = NULL;
183  char *filename = NULL;
184 
185  /* directory name */
186  tmpbuf = strdup(oldpath);
187  strcpy(newpath, dirname(tmpbuf));
188  strcat(newpath, "/");
189 
190  /* filename up to substring */
191  tmpbuf = strdup(oldpath);
192  filename = basename(tmpbuf);
193  status = (((ptr = strstr(filename, "LAC")) == NULL) &&
194  ((ptr = strstr(filename, "1KM")) == NULL));
195  if (status) {
196  fprintf(stderr,
197  "Input file %s does not have standard name; must include "
198  "string \"1KM\" or \"LAC\" in order to find \"%s\" file.\n",
199  oldpath, newchars);
200  newpath[0] = '\0';
201  return status;
202  }
203  strncat(newpath, filename, ptr - filename);
204 
205  /* add new substring and the rest of the filename */
206  strcat(newpath, newchars);
207  strcat(newpath, ptr + strlen(newchars));
208  free(tmpbuf);
209 
210  return status;
211 }
212 
213 /***********************************************************************/
214 
216  int32_t status = FAIL;
217  int32_t i;
218  char result[FILENAME_MAX];
219  sdslist *file_content;
220 
221  /* Open input file */
222  status = ((l1bfile->id = SDstart(l1bfile->file, DFACC_RDONLY)) == FAIL);
223  if (status) {
224  fopen_warn(l1bfile->file, __FILE__, __LINE__);
225  return status;
226  }
227 
228  /*----- Read Global Attributes to determine file type -----*/
229  get_hdfeos_meta(l1bfile->id, "CoreMetadata.0", "SHORTNAME", result);
230  if ((strcmp(result, "MYD021KM") == 0) ||
231  (strcmp(result, "MOD021KM") == 0)) { // 1KM
232  l1bfile->resolution = 1000;
233  file_content = SDSLIST_1KM;
234  } else if ((strcmp(result, "MYD02HKM") == 0) ||
235  (strcmp(result, "MOD02HKM") == 0)) { // HKM
236  l1bfile->resolution = 500;
237  file_content = SDSLIST_HKM;
238  } else if ((strcmp(result, "MYD02QKM") == 0) ||
239  (strcmp(result, "MOD02QKM") == 0)) { // QKM
240  l1bfile->resolution = 250;
241  file_content = SDSLIST_QKM;
242  } else {
243  fprintf(stderr,
244  "Input file %s has type %s; "
245  "does not contain MODIS calibrated radiances.\n",
246  l1bfile->file, result);
247  return FAIL;
248  }
249  strcpy(l1bfile->shortname, result);
250 
251  /*----- Data set info -----*/
252  l1bfile->n_sds = 0;
253  while (file_content[l1bfile->n_sds].name)
254  l1bfile->n_sds++;
255  TRYMEM(__FILE__, __LINE__,
256  (l1bfile->sds = calloc(l1bfile->n_sds, sizeof (sds_struct))));
257  for (i = 0; i < l1bfile->n_sds; i++)
258  init_sds_byname(l1bfile->id, file_content[i].name, &(l1bfile->sds[i]));
259 
260  /*----- Native dimensions -----*/
261  l1bfile->nlines = l1bfile->sds[RSB_250].dimlen[1]; //(bands,NSCANS*NDETS,nframes*nsf)
262  l1bfile->npixls = l1bfile->sds[RSB_250].dimlen[2]; //(bands,nscans*ndets,NFRAMES*NSF)
263 
264  /* Note: Every MODIS L1B file contains a 250m-band SDS, aggregated for
265  coarser resolutions, so we can use [RSB_250] to reference native
266  dimensions for every file. Later we'll use [RSB_250] to index the
267  highest-resolution SDS, to specify output dimensions. */
268 
269  return SUCCESS;
270 }
271 
272 int init_l1b(const char filename[FILENAME_MAX], int32_t *max_resolution) {
273  int32_t i, j;
274  int32_t status = FAIL;
275 
276  /* Start fresh */
277  memset(&file_1km, 0, sizeof (modis_file));
278  memset(&file_hkm, 0, sizeof (modis_file));
279  memset(&file_qkm, 0, sizeof (modis_file));
280 
281  /* Initialize L1B 1KM file info */
282  strcpy(file_1km.file, filename);
283  status = open_modis_l1bfile(&file_1km);
284  if (status) {
285  fprintf(stderr,
286  "Error reading %s; please specify a valid Level 1B file.\n",
287  filename);
288  return status;
289  }
290 
291  /* Input L1B file must have 1KM resolution */
292  status = (file_1km.resolution != 1000);
293  if (status) {
294  fprintf(stderr,
295  "Input file %s has %dm resolution; "
296  "please specify a 1km-resolution Level 1B file.\n",
297  filename, file_1km.resolution);
298  return status;
299  }
300 
301  /* Open other L1B files as needed for specified resolution */
302  switch (*max_resolution) {
303 
304  case 250: // QKM
305  status = modpath_1km(filename, "QKM", file_qkm.file);
306  if (status)
307  return status;
308  status = open_modis_l1bfile(&file_qkm);
309  if (status) {
310  fprintf(stderr, "File not found: %s\n", file_qkm.file);
311  fprintf(stderr,
312  "Processing at %im resolution requires a QKM file to be "
313  "present in the same directory as the 1KM L1B file.\n",
314  *max_resolution);
315  return status;
316  }
317  /* no break: also load HKM */
318 
319  case 500: // HKM
320  status = modpath_1km(filename, "HKM", file_hkm.file);
321  if (status)
322  return status;
323  status = open_modis_l1bfile(&file_hkm);
324  if (status) {
325  fprintf(stderr, "File not found: %s\n", file_hkm.file);
326  fprintf(stderr,
327  "Processing at %im resolution requires a HKM file to be "
328  "present in the same directory as the 1KM L1B file.\n",
329  *max_resolution);
330  return status;
331  }
332  break;
333 
334  case -1: // default to 1KM
335  *max_resolution = 1000;
336  /* no break: proceed as for 1KM */
337 
338  case 1000: // 1KM
339  /* file_1km has already been initialized */
340  break;
341 
342  default:
343  fprintf(stderr, "Invalid resolution %i; ", *max_resolution);
344  fprintf(stderr, "must be 1000, 500 or 250 m.\n");
345  status = FAIL;
346  exit(status);
347  break;
348  }
349 
350  /* Populate global info for each L1B SDS */
351  for (i = 0; i < L1B_NUM_SDS; i++) {
352  l1b[i].sds = file_1km.sds[i]; // initialize to 1KM inputs
353  l1b[i].resolution = file_1km.resolution;
354  l1b[i].nscans = file_1km.nlines / 10;
355  l1b[i].nframes = file_1km.npixls;
356  l1b[i].scandim = SDSLIST_1KM[i].scandim;
357  }
358  if (*max_resolution < 1000) {
359  l1b[RSB_500].sds = file_hkm.sds[RSB_500]; // HKM or QKM
360  l1b[RSB_500].resolution = file_hkm.resolution;
361 
362  l1b[RSB_250].sds = file_hkm.sds[RSB_250]; // HKM only
363  l1b[RSB_250].resolution = file_hkm.resolution;
364  }
365  if (*max_resolution < 500) {
366  l1b[RSB_250].sds = file_qkm.sds[RSB_250]; // QKM only
367  l1b[RSB_250].resolution = file_qkm.resolution;
368  }
369 
370  /* Derived values */
371  for (i = 0; i < L1B_NUM_SDS; i++) {
372  l1b[i].nbands =
373  (l1b[i].scandim == 0) ? 1 // Band 26 only
374  : l1b[i].sds.dimlen[0];
375  l1b[i].nsf = 1000 / l1b[i].resolution;
376  l1b[i].ndets = 10 * l1b[i].nsf;
377 
378  /* Load scaling coefficient attribute(s) */
379  TRYMEM(__FILE__, __LINE__,
380  (l1b[i].sds.atts = calloc(L1B_NUM_COEFFS, sizeof (att_struct))));
381  l1b[i].sds.natts = L1B_NUM_COEFFS;
382  for (j = 0; j < L1B_NUM_COEFFS; j++)
383  load_att_byname(l1b[i].sds.id, l1b_coeff[j], &l1b[i].sds.atts[j]);
384  }
385 
386  return status;
387 }
388 
389 /***********************************************************************/
390 
392  int32_t status = FAIL;
393  int32_t i;
394  char result[FILENAME_MAX];
395 
396  /* Open input file */
397  status = ((geofile->id = SDstart(geofile->file, DFACC_RDONLY)) == FAIL);
398  if (status) {
399  fopen_warn(geofile->file, __FILE__, __LINE__);
400  return status;
401  }
402 
403  /*----- Read Global Attributes to determine file type -----*/
404  get_hdfeos_meta(geofile->id, "CoreMetadata.0", "SHORTNAME", result);
405  if ((strcmp(result, "MYD03") == 0) ||
406  (strcmp(result, "MOD03") == 0)) { // GEO
407  } else {
408  fprintf(stderr, "Input file %s has type %s; "
409  "does not contain MODIS geolocation.\n",
410  geofile->file, result);
411  return FAIL;
412  }
413  strcpy(geofile->shortname, result);
414 
415  /*----- Data set info -----*/
416  geofile->n_sds = GEO_NUM_SDS;
417  TRYMEM(__FILE__, __LINE__,
418  (geofile->sds = calloc(GEO_NUM_SDS, sizeof (sds_struct))));
419  for (i = 0; i < GEO_NUM_SDS; i++)
420  init_sds_byname(geofile->id, SDSLIST_GEO[i].name, &(geofile->sds[i]));
421 
422  /*----- Native dimensions -----*/
423  geofile->nlines = geofile->sds[GEO_LAT].dimlen[0];
424  geofile->npixls = geofile->sds[GEO_LAT].dimlen[1];
425 
426  return SUCCESS;
427 }
428 
429 int init_geo(const char filename[FILENAME_MAX]) {
430  int32_t i, j;
431  int32_t status = FAIL;
432 
433  /* Start fresh */
434  memset(&file_geo, 0, sizeof (modis_file));
435 
436  /* Initialize GEO file info */
437  strcpy(file_geo.file, filename);
438  status = open_modis_geofile(&file_geo);
439  if (status) {
440  fprintf(stderr,
441  "Error reading %s; please specify a valid geolocation file.\n",
442  filename);
443  return status;
444  }
445 
446  /* Populate global info for each GEO SDS */
447  for (i = 0; i < GEO_NUM_SDS; i++) {
448  geo[i].sds = file_geo.sds[i];
449  geo[i].resolution = 1000;
450  geo[i].nbands = 1;
451  geo[i].nscans = file_geo.nlines / 10;
452  geo[i].ndets = 10;
453  geo[i].nframes = file_geo.npixls;
454  geo[i].nsf = 1;
455  geo[i].scandim = SDSLIST_GEO[i].scandim;
456 
457  /* Load scaling coefficient attribute(s) */
458  TRYMEM(__FILE__, __LINE__,
459  (geo[i].sds.atts = calloc(GEO_NUM_COEFFS, sizeof (att_struct))));
460  geo[i].sds.natts = GEO_NUM_COEFFS;
461  for (j = 0; j < GEO_NUM_COEFFS; j++)
462  load_att_byname(geo[i].sds.id, geo_coeff[j], &geo[i].sds.atts[j]);
463  }
464 
465  /*----- Read preloaded variables -----*/
466  for (i = 0; i < REF_NUM_SDS; i++) {
467  init_sds_byname(file_geo.id, SDSLIST_REF[i].name, &(ref[i].sds));
468  readall_sds(&(ref[i].sds));
469  }
470 
471  return status;
472 }
473 
474 int read_sds_1scan(sds_struct *sds,
475  int32_t iscan,
476  int32_t scandim,
477  int32_t per_scan) {
478  int32_t i;
479  int32_t status = FAIL;
480  int32_t nvals = 1;
481  int32_t *start, *edges;
482 
483  /* don't bother if SDS handle is invalid */
484  if (sds->id != FAIL) {
485 
486  /* set array indices */
487  TRYMEM(__FILE__, __LINE__,
488  (start = malloc(sds->ndims * sizeof (int32_t))));
489  TRYMEM(__FILE__, __LINE__,
490  (edges = malloc(sds->ndims * sizeof (int32_t))));
491  for (i = 0; i < sds->ndims; i++) {
492  if (i == scandim) {
493  start[i] = iscan * per_scan;
494  edges[i] = per_scan;
495  } else {
496  start[i] = 0;
497  edges[i] = sds->dimlen[i];
498  }
499  nvals *= edges[i];
500  }
501 
502  /* allocate one scan as needed (first call) */
503  if (sds->data == NULL)
504  TRYMEM(__FILE__, __LINE__,
505  (sds->data = malloc(nvals * DFKNTsize(sds->ntype))));
506 
507  /* read data */
508  status = SDreaddata(sds->id, start, NULL, edges, sds->data);
509  free(start);
510  free(edges);
511  }
512  if (status != FAIL)
513  return nvals;
514  else
515  return status;
516 }
517 
518 /***********************************************************************/
519 /* Everything below this line is specific to OCSSW implementation. */
520 /***********************************************************************/
521 #define MIN_BAD_SI 65500
522 #define MAX_ATTERR 0.017453293 /* radians, ~= 1 degree */
523 #define LOOPI for (i = 0; i < nvals; i++)
524 
525 /* Global variables from l1_hmodis_hdf.c */
526 static int32_t spixl = 0; /* first valid pixel of extracted granule */
527 static float *Fobar;
528 
529 /* RSMAS corrections to Terra IR radiances (20,22,23,27,28,29,31,32) */
530 /* latband v6 sst coefficients were estimated from Bt's from which these */
531 /* corrections have been applied, so we have to keep using them. */
532 static float radoff[][10] = {
533  {-0.000972, 0.014200, 0.000022, 0.000238, -0.000003,
534  0.003760, 0.002337, 0.000084, 0.008848, 0.005050},
535  {-0.002300, 0.000600, 0.000236, -0.000280, -0.000003,
536  0.002798, 0.003496, 0.018035, 0.002942, 0.001787},
537  {-0.003833, 0.003657, 0.002118, 0.000798, -0.000003,
538  0.000208, 0.000399, 0.000553, 0.000258, 0.021128},
539  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
540  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
541  {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
542  {-0.000423, -0.000242, -0.000330, -0.000065, -0.000001,
543  -0.000006, 0.000064, 0.000398, 0.000362, 0.000322},
544  {-0.000433, -0.000246, -0.000222, -0.000148, -0.000001,
545  0.000068, 0.000221, 0.000221, 0.000437, 0.000303}
546 };
547 static float mfact[] = {0.0, 0.0, 0.0005, 0.0, 0.0, 0.0, 0.0028, 0.00278};
548 
552 typedef struct {
553  int16_t sdsindex;
554  int16_t bandindex;
555  int16_t flagindex;
556  int16_t wavelength;
557  char *bandname;
558 } bandinfo;
559 
560 /* Bands used by l2gen, in expected order */
561 static const bandinfo BANDLIST[] = {
562  /*** Reflective Solar Bands ***/
563  { RSB_1KM, 0, 180, 412, "8"}, // 0
564  { RSB_1KM, 1, 190, 443, "9"}, // 1
565  { RSB_500, 0, 80, 469, "3"}, // 2
566  { RSB_1KM, 2, 200, 488, "10"}, // 3
567  { RSB_1KM, 3, 210, 531, "11"}, // 4
568  { RSB_1KM, 4, 220, 547, "12"}, // 5
569  { RSB_500, 1, 100, 555, "4"}, // 6
570  { RSB_250, 0, 0, 645, "1"}, // 7
571  { RSB_1KM, 5, 230, 667, "13lo"}, // 8
572  { RSB_1KM, 7, 250, 678, "14lo"}, // 9
573  { RSB_1KM, 9, 270, 748, "15"}, // 10
574  { RSB_250, 1, 40, 859, "2"}, // 11
575  { RSB_1KM, 10, 280, 869, "16"}, // 12
576  { RSB_500, 2, 120, 1240, "5"}, // 13
577  { RSB_500, 3, 140, 1640, "6"}, // 14
578  { RSB_500, 4, 160, 2130, "7"}, // 15
579  /*** Thermal Emissive Bands ***/
580  { TEB_1KM, 0, 330, 3750, "20"}, // 16
581  { TEB_1KM, 2, 350, 3959, "22"}, // 17
582  { TEB_1KM, 3, 360, 4050, "23"}, // 18
583  { TEB_1KM, 6, 390, 6715, "27"}, // 19
584  { TEB_1KM, 7, 400, 7325, "28"}, // 20
585  { TEB_1KM, 8, 410, 8550, "29"}, // 21
586  { TEB_1KM, 10, 430, 11000, "31"}, // 22
587  { TEB_1KM, 11, 440, 12000, "32"}, // 23
588  /*** Cirrus Band ***/
589  { CIR_1KM, 0, 320, 1375, "26"} // 24
590 };
591 
592 /***********************************************************************/
593 
599 int openl1_hmodis_hdf(filehandle *l1file) {
600  char geometa[EOSMETALEN] = "";
601  char result[EOSMETALEN] = "";
602  char datetime[32] = "";
603  int32_t resolution;
604  int32_t intval;
605  uint32_t gflags[8];
606 
607  /*----- Initialize L1B and GEO info -----*/
608  if (init_l1b(l1file->name, &input->resolution))
609  exit(EXIT_FAILURE);
610  if (init_geo(l1file->geofile))
611  exit(EXIT_FAILURE);
612 
613  /* Make sure this is not an ocean-band subsetted L1A file */
614  int i;
615  for (i = 0; i < file_1km.n_sds; i++) {
616  if (strcmp(file_1km.sds[i].name, "EV_250_Aggr1km_RefSB") == 0) {
617  if (SDfindattr(file_1km.sds[i].id, "Rescaled Ocean R") != FAIL) {
618  fprintf(stderr, "\n"
619  "This L1B file contains only the ocean band subset.\n"
620  "Processing is no longer supported; exiting.\n");
621  exit(EXIT_FAILURE);
622  }
623  break;
624  }
625  }
626 
627  /*----- Populate filehandle structure -----*/
628  l1file->ndets = l1b[RSB_250].ndets;
629  l1file->nscan = l1b[RSB_250].nscans * l1b[RSB_250].ndets;
630  l1file->npix = l1b[RSB_250].nframes * l1b[RSB_250].nsf;
631 
632  /* Capture processing resolution */
634  if (resolution == 1000)
635  strcpy(l1file->spatialResolution, "1 km");
636  else
637  sprintf(l1file->spatialResolution, "%d m", resolution);
638 
639  if (want_verbose) {
640  printf("Processing at %d meter resolution.\n", resolution);
641  printf(" 1000-meter file: %s\n", file_1km.file);
642  if (resolution < 1000) {
643  printf(" 500-meter file: %s\n", file_hkm.file);
644  if (resolution < 500)
645  printf(" 250-meter file: %s\n", file_qkm.file);
646  }
647  printf("\n");
648  }
649 
650  /*----- Read Level 1B metadata -----*/
651  get_hdfeos_meta(file_1km.id, "CoreMetadata.0", "INPUTPOINTER", result);
652  strcpy(l1file->calfile, result);
653 
654  /*----- Read Geolocation metadata -----*/
655  if (read_att(file_geo.id, "CoreMetadata.0", geometa)) {
656  fprintf(stderr,
657  "-E- %s line %d: Error reading CoreMetadata.0 from %s.\n",
658  __FILE__, __LINE__, file_geo.file);
659  exit(EXIT_FAILURE);
660  }
661 
662  /* Orbit info */
663  parse_odl(geometa, "ORBITNUMBER", result);
664  if (strlen(result) > 0)
665  l1file->orbit_number = (int32_t) atoi(result);
666  parse_odl(geometa, "EQUATORCROSSINGLONGITUDE", result);
667  if (strlen(result) > 0)
668  l1file->orbit_node_lon = atof(result);
669 
670  parse_odl(geometa, "EQUATORCROSSINGDATE", result);
671  strcpy(datetime, result);
672  strcat(datetime, "T");
673  parse_odl(geometa, "EQUATORCROSSINGTIME", result);
674  strcat(datetime, result);
675  l1file->node_crossing_time = isodate2unix(datetime);
676 
677  /* Granule start time */
678  parse_odl(geometa, "RANGEBEGINNINGDATE", result);
679  strcpy(datetime, result);
680  strcat(datetime, "T");
681  parse_odl(geometa, "RANGEBEGINNINGTIME", result);
682  strcat(datetime, result);
683  file_geo.starttime = isodate2unix(datetime);
684 
685  /* Terrain height */
686  if (read_att(file_geo.id, "Cumulated gflags", gflags))
687  fprintf(stderr, "-W- %s line %d: Error reading gflags from %s.\n",
688  __FILE__, __LINE__, file_geo.file);
689  else
690  l1file->terrain_corrected = (((int32_t*) gflags)[5] == 0);
691 
692  /*** Adjust dimensions for extracted granules ***/
693  /*
694  When Level 1 MODIS granules are extracted, the full scan
695  width is retained but pixels outside the ROI are set to fill
696  value. We must keep track of pixel offset and count in
697  order to load the ROI correctly. Line offset and count are
698  handled in the main() loop.
699  */
700  if (read_att(file_1km.id, "Extract Pixel Count", &intval) == 0
701  && intval > 0) {
702  int32_t sline, nline, npixl;
703  npixl = intval * l1b[RSB_250].nsf;
704  read_att(file_1km.id, "Extract Pixel Offset", &intval);
705  spixl = intval * l1b[RSB_250].nsf;
706  read_att(file_1km.id, "Extract Line Offset", &intval);
707  sline = intval * l1b[RSB_250].ndets;
708  if (want_verbose) {
709  nline = l1file->nscan;
710  fprintf(stdout,
711  "File was generated from L1A extract of size %d x %d\n",
712  npixl, nline);
713  fprintf(stdout, " Pixels %d - %d\n", spixl + 1, spixl + npixl);
714  fprintf(stdout, " Lines %d - %d\n", sline + 1, sline + nline);
715  }
716  l1file->npix = npixl;
717  }
718 
719  /* Read reflectance to radiance conversion factors */
720  rdsensorinfo(l1file->sensorID, input->evalmask, "Fobar", (void **) &Fobar);
721 
722  return SUCCESS;
723 }
724 
725 /***********************************************************************/
726 
729 static struct {
730  int32_t iscan;
731  int32_t nvals;
732  int32_t npix;
733  int32_t nbands;
735  /* read from input files */
736  double taisec;
737  int32_t mside;
738  double angles[3];
739  double mnorm[3];
740  float *lon;
741  float *lat;
742  float *hgt;
743  float *solz;
744  float *sola;
745  float *senz;
746  float *sena;
747  float *allbands;
748 } scan;
749 
755 int alloc_scan(const int32_t nvals) {
756 
757  /* set stack variables */
758  memset(&scan, 0, sizeof (scan));
759  scan.mside = -1;
760  scan.nvals = nvals;
761  scan.npix = nvals / l1b[RSB_250].ndets;
762  scan.nbands = sizeof (BANDLIST) / sizeof (BANDLIST[0]); // works on stack vars only
763 
764  /* allocate memory for dynamic variables */
765  TRYMEM(__FILE__, __LINE__, (scan.lon = malloc(nvals * sizeof (*scan.lon))));
766  TRYMEM(__FILE__, __LINE__, (scan.lat = malloc(nvals * sizeof (*scan.lat))));
767  TRYMEM(__FILE__, __LINE__, (scan.hgt = malloc(nvals * sizeof (*scan.hgt))));
768  TRYMEM(__FILE__, __LINE__, (scan.solz = malloc(nvals * sizeof (*scan.solz))));
769  TRYMEM(__FILE__, __LINE__, (scan.sola = malloc(nvals * sizeof (*scan.sola))));
770  TRYMEM(__FILE__, __LINE__, (scan.senz = malloc(nvals * sizeof (*scan.senz))));
771  TRYMEM(__FILE__, __LINE__, (scan.sena = malloc(nvals * sizeof (*scan.sena))));
772  TRYMEM(__FILE__, __LINE__,
773  (scan.allbands = malloc(scan.nbands * nvals * sizeof (*scan.allbands))));
774 
775  return SUCCESS;
776 }
777 
778 /***********************************************************************/
779 
780 typedef struct {
781  int16_t nrad;
782  int16_t nsf;
783  double *SI_val;
784  double *corr;
785 } sfcorr_table;
786 
787 void read_sfcorrtables(sfcorr_table *sfcorr, int32_t sensorID) {
788  char *dataroot;
789  char file[FILENAME_MAX] = "";
790  char sensorName[32];
791  FILE *fp = NULL;
792  char line[81];
793  int32_t iband, isf, nsf, irad, nrad;
794  size_t isds, ical;
795  att_struct *atts;
796  double scale, offset;
797  double radiance, SI, radiance2, SI2, rad_corr[4];
798 
799  /* calibration directory path */
800  if ((dataroot = getenv("OCDATAROOT")) == NULL) {
801  printf("OCDATAROOT environment variable is not defined.\n");
802  exit(1);
803  }
804 
805  /* for each band, */
806  /* get max possible subframes */
807  if (want_verbose)
808  printf("\nLoading subframe correction tables:\n");
809  for (iband = 0; iband < scan.nbands; iband++) {
810  nrad = 0;
811  isds = BANDLIST[iband].sdsindex;
812  if (isds == RSB_250)
813  nsf = 4;
814  else if (isds == RSB_500)
815  nsf = 2;
816  else
817  nsf = 1;
818 
819  /* if a hi-res band, */
820  if (nsf > 1) {
821 
822  /* get radiance factors */
823  atts = l1b[isds].sds.atts;
824  ical = BANDLIST[iband].bandindex;
825  scale = (double) ((float*) atts[RAD_SCALE].data)[ical];
826  offset = (double) ((float*) atts[RAD_OFFSET].data)[ical];
827 
828  strcpy(sensorName, sensorId2SensorName(sensorID));
829  lowcase(sensorName);
830 
831  /* open subframe correction file */
832  sprintf(file, "%s/%s/%s/cal/subframecor_%s_%d.dat",
833  dataroot,
836  sensorName,
837  BANDLIST[iband].wavelength);
838  if (want_verbose)
839  printf(" %s\n", file);
840  if ((fp = fopen(file, "r")) == NULL) {
841  fprintf(stderr,
842  "-E- %s line %d: unable to open %s for reading\n",
843  __FILE__, __LINE__, file);
844  exit(1);
845  }
846 
847  /* find number of radiance levels and allocate space */
848  while (fgets(line, 80, fp))
849  if (strncmp(line, "Number of radiance levels", 25) == 0) {
850  sscanf(line, "Number of radiance levels: %d", &nrad);
851  TRYMEM(__FILE__, __LINE__,
852  (sfcorr[iband].SI_val =
853  malloc(nrad * sizeof (*sfcorr[iband].SI_val))));
854  TRYMEM(__FILE__, __LINE__,
855  (sfcorr[iband].corr =
856  malloc(nrad * 4 * sizeof (*sfcorr[iband].corr))));
857  break;
858  }
859  while (fgets(line, 80, fp)) /* skip over rest of header */
860  if (strncmp(line, "/end_header", 11) == 0)
861  break;
862 
863  /* read radiances and correction factors */
864  for (irad = 0; irad < nrad; irad++) {
865  if (fgets(line, 80, fp) == NULL) {
866  fprintf(stderr,
867  "-E- %s line %d: File %s contains only %d data records\n",
868  __FILE__, __LINE__, file, irad);
869  exit(1);
870  }
871  sscanf(line, "%lf %lf %lf %lf %lf", &radiance,
872  &rad_corr[0], &rad_corr[1], &rad_corr[2], &rad_corr[3]);
873 
874  /* Scale reference radiances and corrections to SI
875  space so we don't have to scale each pixel later.
876  Also invert factor so we later multiply instead of
877  divide.
878  */
879  SI = radiance / scale + offset;
880  sfcorr[iband].SI_val[irad] = SI;
881  for (isf = 0; isf < 4; isf++) {
882  radiance2 = radiance / rad_corr[isf];
883  SI2 = radiance2 / scale + offset;
884  sfcorr[iband].corr[4 * irad + isf] = SI2 / SI;
885  }
886 
887  }
888  fclose(fp);
889 
890  } // hi-res band
891  sfcorr[iband].nrad = nrad;
892  sfcorr[iband].nsf = nsf;
893 
894  } // iband
895 
896  if (want_verbose)
897  printf("Subframe destriping corrections enabled.\n\n");
898 }
899 
900 void subframe_correction(const int32_t sensorID,
901  const int32_t iband,
902  const modis_sds mds,
903  double *data) {
904  /* Applies subframe destriping to an entire full-width scan at once,
905  so we don't need to correct for any extract offset. */
906  sfcorr_table lut;
907  int32_t i;
908  int32_t idet, iframe, isf, ip;
909  double x; // input value (scaled radiance)
910  double y; // interpolated subframe correction
911 
912  static sfcorr_table *sfcorr = NULL;
913 
914  /* Initialize subframe correction tables */
915  if (sfcorr == NULL) {
916  TRYMEM(__FILE__, __LINE__,
917  (sfcorr = malloc(scan.nbands * sizeof (*sfcorr))));
918  read_sfcorrtables(sfcorr, sensorID);
919  }
920  lut = sfcorr[iband];
921 
922  /* Destripe hi-res bands at native resolution only */
923  if ((lut.nsf > 1) && (lut.nrad > 0) && (lut.nsf == mds.nsf)) {
924  for (idet = 0; idet < mds.ndets; idet++)
925  for (iframe = 0; iframe < mds.nframes; iframe++)
926  for (isf = 0; isf < mds.nsf; isf++) {
927  ip = (idet * mds.nframes + iframe) * mds.nsf + isf;
928 
929  if (data[ip] < MIN_BAD_SI) {
930  x = data[ip];
931 
932  /* find closest lookup table indices */
933  /* TO DO: implement more efficient lookup */
934  for (i = 0; i < lut.nrad - 1; i++) {
935  if (x < lut.SI_val[i])
936  break;
937  }
938 
939  /* use end value if out-of-bounds */
940  if (i == 0 || i == (lut.nrad - 1))
941  y = lut.corr[4 * i + isf];
942 
943  /* otherwise, interpolate between two closest values */
944  else {
945  double x0 = lut.SI_val[i - 1];
946  double x1 = lut.SI_val[i];
947  double y0 = lut.corr[4 * (i - 1) + isf];
948  double y1 = lut.corr[4 * i + isf];
949  y = y0 + (x - x0) * (y1 - y0) / (x1 - x0);
950  }
951 
952  /* Apply the destriping correction */
953  data[ip] *= (float) y;
954  }
955  } // isf, iframe, idet loops
956  } // eligible band
957 }
958 
959 /***********************************************************************/
960 
961 typedef struct {
962  int32_t ndets;
963  int32_t *prev;
964  int32_t *next;
965 } fill_table;
966 
968  int32_t ndets; // number of detectors
969  int32_t isds, iband, idet, i0, x0, x1;
970  char dead[EOSMETALEN] = "";
971 
972  /* Read "Dead Detector List" global attribute */
973  if (read_att(file_1km.id, "Dead Detector List", dead)) {
974  fprintf(stderr,
975  "-E- %s line %d: Error reading \"Dead Detector List\" from %s.\n",
976  __FILE__, __LINE__, file_1km.file);
977  exit(EXIT_FAILURE);
978  }
979 
980  /* Load values into dead detector table */
981  for (iband = 0; iband < scan.nbands; iband++) {
982  i0 = BANDLIST[iband].flagindex;
983 
984  /* Find number of detectors at native resolution */
985  isds = BANDLIST[iband].sdsindex;
986  if (isds == RSB_250)
987  ndets = 40;
988  else if (isds == RSB_500)
989  ndets = 20;
990  else
991  ndets = 10;
992  fill[iband].ndets = ndets;
993 
994  /* Allocate space for arrays */
995  TRYMEM(__FILE__, __LINE__,
996  (fill[iband].prev = malloc(ndets * sizeof (*fill->prev))));
997  TRYMEM(__FILE__, __LINE__,
998  (fill[iband].next = malloc(ndets * sizeof (*fill->next))));
999 
1000  /* Find previous and next good detectors */
1001  /*
1002  These arrays hold indices for the nearest previous and next
1003  valid detector for each detector. Since detector arrays are
1004  small, opt for code clarity over efficiency.
1005 
1006  Special values:
1007  prev=idet=next: detector is alive
1008  prev = -1 : no previous good detector
1009  next = ndets: no next good detector
1010  */
1011  for (idet = 0; idet < ndets; idet++) {
1012  x0 = idet; // previous
1013  while ((x0 > -1) && dead[i0 + x0]) x0--;
1014  fill[iband].prev[idet] = x0;
1015  x1 = idet; // next
1016  while ((x1 < ndets) && dead[i0 + x1]) x1++;
1017  fill[iband].next[idet] = x1;
1018  }
1019 
1020  } // iband loop
1021 }
1022 
1023 void fill_dead_detectors(const int32_t iband, const modis_sds mds, double *data) {
1024 
1025  int32_t idet, ndets, ipix, npix;
1026  size_t nbytes; // bytes per line (detector)
1027 
1028  int32_t *prev, *next; // previous and next valid detector for each det
1029  int32_t x, x0, x1; // index into data array for current, prev, next det
1030  double y0, y1; // previous and next valid value for current detector
1031  double dx, dy; // change in det, value over valid interval
1032 
1033  static fill_table *fill = NULL;
1034 
1035  /* Initialize previous and next valid detector array table */
1036  if (fill == NULL) {
1037  TRYMEM(__FILE__, __LINE__, (fill = malloc(scan.nbands * sizeof (*fill))));
1038  init_detfill(fill);
1039  }
1040 
1041  /* Fill dead detectors at native resolution only */
1042  ndets = fill[iband].ndets;
1043  if (ndets == mds.ndets) {
1044 
1045  prev = fill[iband].prev;
1046  next = fill[iband].next;
1047  npix = mds.nframes * mds.nsf; // # values per line (detector)
1048  nbytes = npix * sizeof (*data);
1049 
1050  /* Substitute data values as appropriate for each detector: */
1051  for (idet = 0; idet < ndets; idet++) {
1052  if (prev[idet] != next[idet]) { // continue for bad dets only
1053 
1054  /* Between two good detectors: interpolate */
1055  if ((prev[idet] > -1) && (next[idet] < ndets)) {
1056  dx = next[idet] - prev[idet];
1057 
1058  for (ipix = 0; ipix < npix; ipix++) {
1059  x = idet * npix + ipix;
1060  x0 = prev[idet] * npix + ipix;
1061  x1 = next[idet] * npix + ipix;
1062  y0 = data[x0];
1063  y1 = data[x1];
1064 
1065  /* if both values are valid, interpolate */
1066  if ((y0 < MIN_BAD_SI) && (y1 < MIN_BAD_SI)) {
1067  dy = y1 - y0;
1068  data[x] = y0 + (idet - prev[idet]) * (dy / dx);
1069  }
1070  /* otherwise copy adjacent valid value */
1071  else if (y0 < MIN_BAD_SI) data[x] = data[x0];
1072  else if (y1 < MIN_BAD_SI) data[x] = data[x1];
1073 
1074  /* if none, copy maximum adjacent flag */
1075  else
1076  data[x] = MAX(data[x0], data[x1]);
1077  }
1078  }
1079  /* No previous good detector: copy next value or flag */
1080  else if ((prev[idet] == -1) && (next[idet] < ndets))
1081  memcpy(&data[idet * npix], &data[next[idet] * npix], nbytes);
1082 
1083  /* No next good detector: copy previous value or flag */
1084  else if ((prev[idet] > -1) && (next[idet] == ndets))
1085  memcpy(&data[idet * npix], &data[prev[idet] * npix], nbytes);
1086 
1087  /* Default case: no good detectors = no change. */
1088 
1089  } // dead detector
1090  } // idet loop
1091  } // band is at native resolution
1092 }
1093 
1094 /***********************************************************************/
1095 
1096 double interp_bilin(double val[4], double dx, double dy) {
1097  return (1 - dy) * (val[0] * (1 - dx) + val[1] * dx)
1098  + dy * (val[2] * (1 - dx) + val[3] * dx);
1099 }
1100 
1101 double interp_modis_EVdata(double val[4], double dx, double dy) {
1102  int32_t i, ngood = 0;
1103  double sum_val = 0;
1104  double max_val = 0;
1105  double result;
1106 
1107  /* interpolation: account for flagged values */
1108  for (i = 0; i < 4; i++) {
1109  if (val[i] < MIN_BAD_SI) {
1110  sum_val = sum_val + val[i];
1111  ngood++;
1112  }
1113  if (val[i] > max_val)
1114  max_val = val[i];
1115  }
1116 
1117  /* If we have four non-flagged samples, bilinearly interpolate. */
1118  /* If all samples are flagged, set to the max flag value. */
1119  /* If we have 1-3 samples, use average of unflagged values. */
1120  switch (ngood) {
1121  case 4:
1122  result = interp_bilin(val, dx, dy);
1123  break;
1124  case 0:
1125  result = max_val;
1126  break;
1127  default:
1128  /* TO DO: should be WEIGHTED avg of good values */
1129  result = sum_val / ngood;
1130  break;
1131  }
1132  return result;
1133 }
1134 
1135 double interp_modis_Longitude(double val[4], double dx, double dy) {
1136  int32_t i;
1137  double min_val, max_val;
1138 
1139  /* longitude interpolation: account for dateline */
1140  min_val = max_val = val[0];
1141  for (i = 1; i < 4; i++) {
1142  if (val[i] < min_val)
1143  min_val = val[i];
1144  else if (val[i] > max_val)
1145  max_val = val[i];
1146  }
1147  if (max_val > 180.0 + min_val) {
1148  for (i = 0; i < 4; i++) {
1149  if (val[i] < 0.0)
1150  val[i] += 360.0;
1151  }
1152  }
1153  return interp_bilin(val, dx, dy);
1154 }
1155 
1156 double (*interp_modis_var)(double val[4], double dx, double dy);
1157 
1159  const int32_t sdsband,
1160  const int32_t newres,
1161  double *olddata,
1162  double *newdata) {
1163 
1164  int32_t ny = s.ndets;
1165  int32_t nx = s.nframes * s.nsf;
1166  int32_t nvals = ny * nx;
1167  int32_t oldres = s.resolution;
1168 
1169  int32_t nxnew, nynew;
1170  int32_t ixnew, iynew;
1171 
1172  int32_t x0, y0, x1, y1; // 4 points in old array
1173  double fx, fy;
1174  double dx, dy;
1175  double val[4];
1176  int32_t i;
1177  double factor, xoffset, yoffset;
1178  int32_t maxi = 2; /* 1 duplicates end pixel; 2 extrapolates past end */
1179 
1180  int32_t sdsoffset = nvals * sdsband; // offset within SDS
1181  short alloc_old;
1182 
1183  /* Cast input data to double precision */
1184  alloc_old = (olddata == NULL);
1185  if (alloc_old) {
1186  TRYMEM(__FILE__, __LINE__, (olddata = malloc(nvals * sizeof (*olddata))));
1187 
1188  switch (s.sds.ntype) {
1189  case DFNT_FLOAT32:
1190  LOOPI olddata[i] = (double) ((float*) s.sds.data)[sdsoffset + i];
1191  break;
1192  case DFNT_INT16:
1193  LOOPI olddata[i] = (double) ((int16_t*) s.sds.data)[sdsoffset + i];
1194  break;
1195  case DFNT_UINT16:
1196  LOOPI olddata[i] = (double) ((uint16_t*) s.sds.data)[sdsoffset + i];
1197  break;
1198  default:
1199  fprintf(stderr,
1200  "-E- %s line %d: Cannot interpolate type code: %d\n",
1201  __FILE__, __LINE__, s.sds.ntype);
1202  exit(EXIT_FAILURE);
1203  break;
1204  }
1205  } // alloc_old
1206 
1207  /* Set interp function pointer according to SDS */
1208  if (strncmp(s.sds.name, "EV_", 3) == 0)
1210  else if (strcmp(s.sds.name, "Longitude") == 0)
1212  else
1214 
1215  /* Calculate offset between different pixel resolutions */
1216  factor = (double) newres / (double) oldres;
1217  xoffset = 0; // frame centers are aligned
1218  yoffset = 0.5 * (1.0 - factor); // detector edges are aligned
1219  /* yoffset = 0; // for testing only */
1220 
1221  /* Size of new array */
1222  nxnew = nx * oldres / newres;
1223  nynew = ny * oldres / newres;
1224 
1225  /* Calculate index into original data */
1226  for (iynew = 0; iynew < nynew; iynew++) {
1227  fy = factor * iynew - yoffset;
1228  y0 = MIN(MAX((int32_t) floor(fy), 0), ny - maxi);
1229  y1 = MIN(y0 + 1, ny - 1);
1230  dy = fy - (double) y0;
1231 
1232  for (ixnew = 0; ixnew < nxnew; ixnew++) {
1233  fx = factor * ixnew - xoffset;
1234  x0 = MIN(MAX((int32_t) floor(fx), 0), nx - maxi);
1235  x1 = MIN(x0 + 1, nx - 1);
1236  dx = fx - (double) x0;
1237 
1238  /* Do bilinear interpolation */
1239  val[0] = olddata[nx * y0 + x0];
1240  val[1] = olddata[nx * y0 + x1];
1241  val[2] = olddata[nx * y1 + x0];
1242  val[3] = olddata[nx * y1 + x1];
1243  newdata[nxnew * iynew + ixnew] = interp_modis_var(val, dx, dy);
1244  }
1245  }
1246 
1247  if (alloc_old)
1248  free(olddata);
1249 
1250  return SUCCESS;
1251 }
1252 
1253 /***********************************************************************/
1254 
1260 int load_modis_scan(const int32_t iscan,
1261  const int32_t resolution,
1262  int32_t sensorID) {
1263  modis_sds s;
1264  int32_t isds, iband;
1265  int32_t i, i1, i2, nold;
1266  static int32_t nvals = 0; /* pixels per variable at output resolution */
1267  double *olddata = NULL; /* at native resolution */
1268  double *newdata = NULL; /* at output resolution */
1269 
1270  /* Allocate space to hold all data for one scan */
1271  if (nvals == 0) {
1272  nvals = l1b[RSB_250].ndets * l1b[RSB_250].nframes * l1b[RSB_250].nsf;
1273  alloc_scan(nvals);
1274  }
1275 
1276  /* Allocate space to hold interpolation buffer for one variable */
1277  TRYMEM(__FILE__, __LINE__, (newdata = malloc(nvals * sizeof (*newdata))));
1278 
1279  /* ----- Geolocation values ----- */
1280  for (isds = 0; isds < GEO_NUM_SDS; isds++)
1281  read_sds_1scan(&geo[isds].sds,
1282  iscan,
1283  geo[isds].scandim,
1284  geo[isds].ndets);
1285 
1286  if (geo[GEO_LON].resolution == resolution) { /* native resolution */
1287  /* memcpy works only when src and dst have same data type */
1288  memcpy(scan.lon, geo[GEO_LON].sds.data, nvals * sizeof (float));
1289  memcpy(scan.lat, geo[GEO_LAT].sds.data, nvals * sizeof (float));
1290 
1291  /* must copy each element when datatypes differ */
1292  LOOPI{
1293  scan.hgt[i] = (float) ((int16_t *) geo[GEO_HGT].sds.data)[i];
1294  scan.solz[i] = (float) ((int16_t *) geo[GEO_SOLZ].sds.data)[i] * 0.01;
1295  scan.sola[i] = (float) ((int16_t *) geo[GEO_SOLA].sds.data)[i] * 0.01;
1296  scan.senz[i] = (float) ((int16_t *) geo[GEO_SENZ].sds.data)[i] * 0.01;
1297  scan.sena[i] = (float) ((int16_t *) geo[GEO_SENA].sds.data)[i] * 0.01;}
1298 
1299  } else { /* must interpolate */
1300  modis_interp(geo[GEO_LON], 0, resolution, NULL, newdata);
1301  LOOPI scan.lon[i] = (float) newdata[i];
1302  modis_interp(geo[GEO_LAT], 0, resolution, NULL, newdata);
1303  LOOPI scan.lat[i] = (float) newdata[i];
1304  modis_interp(geo[GEO_HGT], 0, resolution, NULL, newdata);
1305  LOOPI scan.hgt[i] = (float) newdata[i];
1306  modis_interp(geo[GEO_SOLZ], 0, resolution, NULL, newdata);
1307  LOOPI scan.solz[i] = (float) newdata[i] * 0.01;
1308  modis_interp(geo[GEO_SOLA], 0, resolution, NULL, newdata);
1309  LOOPI scan.sola[i] = (float) newdata[i] * 0.01;
1310  modis_interp(geo[GEO_SENZ], 0, resolution, NULL, newdata);
1311  LOOPI scan.senz[i] = (float) newdata[i] * 0.01;
1312  modis_interp(geo[GEO_SENA], 0, resolution, NULL, newdata);
1313  LOOPI scan.sena[i] = (float) newdata[i] * 0.01;
1314  }
1315 
1316  /* ----- Radiometric values ----- */
1317  for (isds = 0; isds < L1B_NUM_SDS; isds++)
1318  read_sds_1scan(&l1b[isds].sds,
1319  iscan,
1320  l1b[isds].scandim,
1321  l1b[isds].ndets);
1322 
1323  for (iband = 0; iband < scan.nbands; iband++) {
1324  s = l1b[BANDLIST[iband].sdsindex];
1325  nold = s.ndets * s.nframes * s.nsf;
1326  i1 = nold * BANDLIST[iband].bandindex; // offset within SDS
1327  i2 = nvals * iband; // offset within scan.allbands
1328 
1329  /* Cast SI values to double precision, in native spatial resolution */
1330  TRYMEM(__FILE__, __LINE__, (olddata = malloc(nold * sizeof (*olddata))));
1331  for (i = 0; i < nold; i++)
1332  olddata[i] = (double) ((uint16_t*) s.sds.data)[i1 + i];
1333 
1334  /* Apply subframe correction to destripe BEFORE interpolation */
1335  if ((sensorID == MODISA) && (resolution < 1000))
1336  subframe_correction(sensorID, iband, s, olddata);
1337 
1338  /* Fill any dead detectors with neighboring values */
1339  fill_dead_detectors(iband, s, olddata);
1340 
1341  /* Interpolate to higher resolution, as needed */
1342  if (s.resolution != resolution) {
1343  modis_interp(l1b[BANDLIST[iband].sdsindex],
1344  BANDLIST[iband].bandindex,
1345  resolution, olddata, newdata);
1346 
1347  /* Copy adjusted data into scan structure */
1348  LOOPI scan.allbands[i2++] = (float) newdata[i]; // interpolated
1349  } else
1350  LOOPI scan.allbands[i2++] = (float) olddata[i]; // native resolution
1351 
1352  free(olddata);
1353  }
1354 
1355  /* ----- Store preloaded values for later ----- */
1356  scan.iscan = iscan;
1357  scan.taisec = ((double*) ref[GEO_TAISEC].sds.data)[iscan];
1358  scan.mside = (int32_t) ((int16_t*) ref[GEO_MSIDE].sds.data)[iscan];
1359  for (i = 0; i < 3; i++) {
1360  scan.angles[i] = ((double*) ref[GEO_ANGLES].sds.data)[3 * iscan + i];
1361  scan.mnorm[i] = ((double*) ref[GEO_MNORM].sds.data)[3 * (3 * iscan + i)];
1362  }
1363  free(newdata);
1364  return SUCCESS;
1365 }
1366 
1373 int set_l1rec_scanvals(const int32_t iscan, l1str *l1rec) {
1374  double esdist;
1375  static double lastTAIsec = -1;
1376  int16_t badatt, badmside;
1377  double TAIsec, dsec;
1378  int16_t year, day;
1379  int32_t i;
1380 
1381  /* Mirror Side */
1382  l1rec->mside = scan.mside;
1383  badmside = (l1rec->mside < 0);
1384  if (badmside) /* Fix mirror side if set to -1 */
1385  l1rec->mside = 0;
1386 
1387  /* Check for non-nominal roll, pitch, or yaw */
1388  badatt =
1389  (fabs(scan.angles[0]) > MAX_ATTERR) ||
1390  (fabs(scan.angles[1]) > MAX_ATTERR) ||
1391  (fabs(scan.angles[2]) > MAX_ATTERR);
1392 
1393  /* Flag each pixel of scan for geolocation errors */
1394  for (i = 0; i < l1rec->npix; i++) {
1395  l1rec->pixnum[i] = spixl + i;
1396  l1rec->navfail[i] = badmside;
1397  l1rec->navwarn[i] = badatt;
1398  }
1399 
1400  /* Scan Start Time */
1401  TAIsec = scan.taisec;
1402  if (TAIsec < 0) { /* Workaround for scan time errors */
1403  fprintf(stderr, "-W- %s: Bad time for scan %d; ", __FILE__, iscan);
1404  if (iscan == 0) {
1405  fprintf(stderr, "using granule start time.\n");
1406  TAIsec = unix_to_tai93(file_geo.starttime);
1407  } else {
1408  fprintf(stderr, "extrapolating from last good scan.\n");
1409  TAIsec = lastTAIsec + SCAN_TIME_INTERVAL;
1410  }
1411  }
1412  l1rec->scantime = tai93_to_unix(TAIsec);
1413  lastTAIsec = TAIsec;
1414 
1415  /* Earth-sun distance correction for this scan */
1416  unix2yds(l1rec->scantime, &year, &day, &dsec);
1417  int32_t msec = (int32_t) (dsec * 1000.0);
1418  int32_t yr = (int32_t) year;
1419  int32_t dy = (int32_t) day;
1420  esdist = esdist_(&yr, &dy, &msec);
1421  l1rec->fsol = pow(1.0 / esdist, 2);
1422 
1423  return SUCCESS;
1424 }
1425 
1433 int readl1_hmodis_hdf(filehandle *l1file,
1434  const int32_t line,
1435  l1str *l1rec) {
1436  static int32_t prev_scan = -1;
1437  int32_t iband, iscan, idet, ndets;
1438  size_t ip, ipb;
1439  size_t detoffset, dataptr;
1440  size_t nbytes;
1441  size_t irsb, iteb, ical;
1442  att_struct *atts;
1443  float scale, offset, SI_val, value;
1444 
1445  /* Scan and detector numbers for this line */
1446  ndets = l1b[RSB_250].ndets;
1447  iscan = line / ndets;
1448  idet = line % ndets;
1449 
1450  /* Load all values for next scan as needed */
1451  if (iscan != prev_scan) {
1452  prev_scan = iscan;
1453  load_modis_scan(iscan, input->resolution, l1file->sensorID);
1454  }
1455 
1456  /* The l1rec structure is reinitialized */
1457  /* each line, so always reload basic info. */
1458  l1rec->npix = l1file->npix;
1459  l1rec->detnum = (int32_t) idet;
1461  unix_to_tai93(file_geo.starttime);
1462 
1463  /* "scan" structure data pointer offset due */
1464  /* to line number and extract start pixel */
1465  detoffset = idet * scan.npix + spixl;
1466 
1467  /* ----- Geolocation values ----- */
1468  nbytes = l1rec->npix * sizeof (float);
1469  memcpy(l1rec->lon, &scan.lon[detoffset], nbytes);
1470  memcpy(l1rec->lat, &scan.lat[detoffset], nbytes);
1471  memcpy(l1rec->height, &scan.hgt[detoffset], nbytes);
1472  memcpy(l1rec->solz, &scan.solz[detoffset], nbytes);
1473  memcpy(l1rec->sola, &scan.sola[detoffset], nbytes);
1474  memcpy(l1rec->senz, &scan.senz[detoffset], nbytes);
1475  memcpy(l1rec->sena, &scan.sena[detoffset], nbytes);
1476 
1477  /* Don't need to check for navigation errors
1478  - it's done by readl1() */
1479 
1480  /* Compute polarization frame rotation angles */
1481  compute_alpha(l1rec->lon, l1rec->lat,
1482  l1rec->senz, l1rec->sena,
1483  scan.mnorm, l1rec->npix, l1rec->alpha);
1484 
1485  /* ----- Radiometric values ----- */
1486  irsb = iteb = 0;
1487  for (iband = 0; iband < scan.nbands; iband++) {
1488 
1489  /* calibration info */
1490  atts = l1b[BANDLIST[iband].sdsindex].sds.atts;
1491  ical = BANDLIST[iband].bandindex;
1492 
1493  /* offset within scan.allbands */
1494  dataptr = iband * scan.nvals + detoffset;
1495 
1496  /*** Thermal Emissive Bands ***/
1497  if (BANDLIST[iband].sdsindex == TEB_1KM) {
1498  scale = ((float*) atts[RAD_SCALE].data)[ical];
1499  offset = ((float*) atts[RAD_OFFSET].data)[ical];
1500 
1501  for (ip = 0; ip < l1rec->npix; ip++) {
1502  value = 0.0; /* init to fill value */
1503 
1504  /* Apply calibration for valid values */
1505  SI_val = scan.allbands[dataptr + ip];
1506  if (SI_val < MIN_BAD_SI) {
1507 
1508  /* Convert radiance from Watts/m^2/micrometer/steradian */
1509  /* to mW/cm2/um/sr, as expected by l2gen */
1510  value = (SI_val - offset) * scale / 10.0;
1511 
1512  /* Terra adjustments */
1513  if (l1rec->l1file->sensorID == MODIST) {
1514  value *= (1.0 - radoff[iteb][idet / (ndets / 10)]);
1515  if (l1rec->mside == 1)
1516  value *= (1.0 + mfact[iteb]);
1517  } /* Aqua adjustments */
1518  else {
1519  /* channel 22 detector balancing */
1520  if (ical == 1)
1521  value /= input->ch22detcor[line % 10];
1522  /* channel 23 detector balancing */
1523  if (ical == 2)
1524  value /= input->ch23detcor[line % 10];
1525  }
1526  }
1527 
1528  /* populate l1rec->Ltir */
1529  ipb = ip * l1rec->l1file->nbandsir + iteb; /* band varies fastest */
1530  l1rec->Ltir[ipb] = value;
1531 
1532  } // ip
1533  iteb++;
1534  }// TEB
1535 
1536  /*** Cirrus Band ***/
1537  else if (BANDLIST[iband].sdsindex == CIR_1KM) {
1538  scale = ((float*) atts[REFL_SCALE].data)[ical];
1539  offset = ((float*) atts[REFL_OFFSET].data)[ical];
1540 
1541  for (ip = 0; ip < l1rec->npix; ip++) {
1542  value = BAD_FLT; /* init to fill value */
1543 
1544  /* Apply calibration for valid values */
1545  SI_val = scan.allbands[dataptr + ip];
1546  if (SI_val < MIN_BAD_SI)
1547 
1548  /* Normalize reflectance by solar zenith angle */
1549  value = (SI_val - offset) * scale
1550  / cos(l1rec->solz[ip] / RADEG);
1551 
1552  /* populate l1rec->rho_cirrus */
1553  l1rec->rho_cirrus[ip] = value;
1554 
1555  } // ip
1556  }// CIR
1557 
1558  /*** Reflective Solar Bands ***/
1559  else {
1560  l1rec->Fo[irsb] = Fobar[irsb] * l1rec->fsol;
1561  scale = ((float*) atts[REFL_SCALE].data)[ical];
1562  offset = ((float*) atts[REFL_OFFSET].data)[ical];
1563 
1564  for (ip = 0; ip < l1rec->npix; ip++) {
1565  value = BAD_FLT; /* init to fill value */
1566 
1567  /* skip night data for visible bands */
1568  if (SOLZNIGHT < l1rec->solz[ip])
1569  continue;
1570 
1571  /* Apply calibration for valid values */
1572  SI_val = scan.allbands[dataptr + ip];
1573  if (SI_val < MIN_BAD_SI)
1574 
1575  /* Convert from reflectance to radiance */
1576  value = (SI_val - offset) * scale * l1rec->Fo[irsb] / PI;
1577 
1578  /* Flag any sentinel values */
1579  else
1580  switch ((int32_t) SI_val) {
1581  case 65527: /* not in earth view: SECTOR_ROTATION_SI */
1582  l1rec->navfail[ip] = 1;
1583  break;
1584  case 65529: /* saturation: TEB_OR_RSB_GT_MAX_SI */
1585  case 65533: /* saturation: SATURATED_DETECTOR_SI */
1586  value = 1000.0;
1587  if (irsb < 13)
1588  l1rec->hilt[ip] = 1;
1589  break;
1590  default: /* all others */
1591  break;
1592  }
1593 
1594  /* populate l1rec->Lt */
1595  ipb = ip * l1rec->l1file->nbands + irsb; /* band varies fastest */
1596  l1rec->Lt[ipb] = value;
1597 
1598  } // ip
1599  irsb++;
1600  } // RSB
1601 
1602  } // iband
1603 
1604  /* Convert IR bands to brightness temperature */
1605  radiance2bt(l1rec, input->resolution);
1606 
1607  return SUCCESS;
1608 }
1609 
1610 /***********************************************************************/
1611 
1617  int32_t i;
1618 
1619  /* Free memory allocated for global variables */
1620  for (i = 0; i < REF_NUM_SDS; i++)
1621  free_sds(&(ref[i].sds));
1622  for (i = 0; i < GEO_NUM_SDS; i++)
1623  free_sds(&(geo[i].sds));
1624  for (i = 0; i < L1B_NUM_SDS; i++)
1625  free_sds(&(l1b[i].sds));
1626 
1627  /* Free any memory allocated for input files */
1628  if (file_geo.sds)
1629  for (i = 0; i < file_geo.n_sds; i++)
1630  free_sds(&(file_geo.sds[i]));
1631  if (file_1km.sds)
1632  for (i = 0; i < file_1km.n_sds; i++)
1633  free_sds(&(file_1km.sds[i]));
1634  if (file_hkm.sds)
1635  for (i = 0; i < file_hkm.n_sds; i++)
1636  free_sds(&(file_hkm.sds[i]));
1637  if (file_qkm.sds)
1638  for (i = 0; i < file_qkm.n_sds; i++)
1639  free_sds(&(file_qkm.sds[i]));
1640 
1641  /* End HDF access */
1642  if (file_geo.id) SDend(file_geo.id);
1643  if (file_1km.id) SDend(file_1km.id);
1644  if (file_hkm.id) SDend(file_hkm.id);
1645  if (file_qkm.id) SDend(file_qkm.id);
1646  file_geo.id = file_1km.id = file_hkm.id = file_qkm.id = 0;
1647 
1648  /* TO DO: free scan structure */
1649 
1650  return SUCCESS;
1651 }
1652 
1653 /***********************************************************************/
1654 
1655 int readl1_lonlat_hmodis_hdf(filehandle *l1file, int32_t line, l1str *l1rec) {
1656  static float *lonscan = NULL; /* storage buffer arrays */
1657  static float *latscan = NULL; /* (at output resolution) */
1658 
1659  int32_t ndets = l1b[RSB_250].ndets; /* output resolution dimensions*/
1660  int32_t npixl = l1b[RSB_250].nframes * l1b[RSB_250].nsf;
1661  int32_t nvals = ndets * npixl;
1662  int32_t resolution = input->resolution;
1663 
1664  double *newdata = NULL; /* interpolation buffer */
1665  static int32_t prev_scan = -1;
1666  int32_t i, iscan, idet;
1667  size_t detoffset, nbytes;
1668 
1669  /* Allocate space to hold lon/lat data for one scan at output resolution */
1670  if (lonscan == NULL) {
1671  TRYMEM(__FILE__, __LINE__, (lonscan = malloc(nvals * sizeof (*lonscan))));
1672  TRYMEM(__FILE__, __LINE__, (latscan = malloc(nvals * sizeof (*latscan))));
1673  }
1674 
1675  /* Scan and detector indices for this line */
1676  iscan = line / ndets;
1677  idet = line % ndets;
1678 
1679  /* Load next scan as needed */
1680  if (iscan != prev_scan) {
1681  prev_scan = iscan;
1682 
1683  /* Read values */
1684  read_sds_1scan(&geo[GEO_LON].sds,
1685  iscan,
1686  geo[GEO_LON].scandim,
1687  geo[GEO_LON].ndets);
1688  read_sds_1scan(&geo[GEO_LAT].sds,
1689  iscan,
1690  geo[GEO_LAT].scandim,
1691  geo[GEO_LAT].ndets);
1692 
1693  /* If native resolution, copy values into storage buffers */
1694  if (geo[GEO_LON].resolution == resolution) {
1695  memcpy(lonscan, geo[GEO_LON].sds.data, nvals * sizeof (*lonscan));
1696  memcpy(latscan, geo[GEO_LAT].sds.data, nvals * sizeof (*latscan));
1697  }
1698  /* Otherwise, interpolate first */
1699  else {
1700  TRYMEM(__FILE__, __LINE__,
1701  (newdata = malloc(nvals * sizeof (*newdata))));
1702  modis_interp(geo[GEO_LON], 0, resolution, NULL, newdata);
1703  LOOPI lonscan[i] = (float) newdata[i];
1704  modis_interp(geo[GEO_LAT], 0, resolution, NULL, newdata);
1705  LOOPI latscan[i] = (float) newdata[i];
1706  free(newdata);
1707  }
1708 
1709  } // loaded new scan
1710 
1711  /* Copy specified line from storage buffers into l1rec structure */
1712  detoffset = idet * npixl + spixl;
1713  nbytes = l1rec->npix * sizeof (float);
1714  memcpy(l1rec->lon, &lonscan[detoffset], nbytes);
1715  memcpy(l1rec->lat, &latscan[detoffset], nbytes);
1716 
1717  return SUCCESS;
1718 }
GEO_SDS
Definition: l1_hmodis_hdf.c:53
int open_modis_geofile(modis_file *geofile)
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
#define MAX(A, B)
Definition: swl0_utils.h:26
int32 value
Definition: Granule.c:1235
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:217
#define MIN(x, y)
Definition: rice.h:169
int32_t nscans
int modpath_1km(const char *oldpath, const char *newchars, char *newpath)
double esdist_(int32_t *year, int32_t *day, int32_t *msec)
int16_t sdsindex
int16_t resolution
float * hgt
#define SUCCESS
Definition: ObpgReadGrid.h:15
int j
Definition: decode_rs.h:73
float * allbands
int32_t day
int status
Definition: l1_czcs_hdf.c:31
int32_t ndets
int set_l1rec_scanvals(const int32_t iscan, l1str *l1rec)
int32_t * prev
char * lowcase(char *instr)
Definition: lowcase.c:10
constexpr double RADEG
double mnorm[3]
void read_sfcorrtables(sfcorr_table *sfcorr, int32_t sensorID)
@ RSB_500
Definition: l1_hmodis_hdf.c:85
@ GEO_SENA
Definition: l1_hmodis_hdf.c:60
@ RAD_SCALE
int16_t resolution
@ L1B_NUM_SDS
Definition: l1_hmodis_hdf.c:89
@ GEO_HGT
Definition: l1_hmodis_hdf.c:56
#define FAIL
Definition: ObpgReadGrid.h:18
#define NULL
Definition: decode_rs.h:63
@ RSB_1KM
Definition: l1_hmodis_hdf.c:86
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
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
int32_t scandim
@ GEO_LAT
Definition: l1_hmodis_hdf.c:55
int16_t wavelength
int readl1_hmodis_hdf(filehandle *l1file, const int32_t line, l1str *l1rec)
int init_l1b(const char filename[FILENAME_MAX], int32_t *max_resolution)
read l1rec
int32_t nbands
char file[FILENAME_MAX]
@ RSB_250
Definition: l1_hmodis_hdf.c:84
L1B_COEFFS
sds_struct sds
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending wavelength
int init_sds_byname(int32_t fileid, const char *sdsname, sds_struct *sds)
Definition: hdf4_utils.c:89
double angles[3]
float * lat
int32 * msec
Definition: l1_czcs_hdf.c:30
#define MODIST
Definition: sensorDefs.h:18
void init_detfill(fill_table *fill)
void subframe_correction(const int32_t sensorID, const int32_t iband, const modis_sds mds, double *data)
float * senz
int open_modis_l1bfile(modis_file *l1bfile)
@ TEB_1KM
Definition: l1_hmodis_hdf.c:88
int32_t npix
@ GEO_SOLA
Definition: l1_hmodis_hdf.c:58
@ GEO_MSIDE
Definition: l1_hmodis_hdf.c:39
int sensorId2SubsensorId(int sensorId)
Definition: sensorInfo.c:299
int32_t n_sds
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
instr * input
@ GEO_NUM_COEFFS
Definition: l1_hmodis_hdf.c:76
L1B_SDS
Definition: l1_hmodis_hdf.c:83
#define BAD_FLT
Definition: jplaeriallib.h:19
#define PI
Definition: l3_get_org.c:6
char * bandname
char * name
Definition: l1_hmodis_hdf.c:33
int32_t nsf
int modis_interp(const modis_sds s, const int32_t sdsband, const int32_t newres, double *olddata, double *newdata)
int32_t index
Definition: l1_hmodis_hdf.c:31
void compute_alpha(float lon[], float lat[], float senz[], float sena[], double mnorm[3], int npix, float alpha[])
Definition: compute_alpha.c:8
int readl1_lonlat_hmodis_hdf(filehandle *l1file, int32_t line, l1str *l1rec)
double unix_to_tai93(double unixtime)
Definition: yds2tai.c:11
int32_t mside
float * solz
double(* interp_modis_var)(double val[4], double dx, double dy)
if(BUILD_JPL) add_executable(interp_hycom interp_hycom.f) add_executable(interp_hycom_ascii interp_hycom_ascii.f) add_executable(tec tec.c) add_executable(swh swh.c) target_link_libraries(interp_hycom netcdff dfutils $
Definition: CMakeLists.txt:127
int openl1_hmodis_hdf(filehandle *l1file)
double tai93_to_unix(double tai93)
Definition: yds2tai.c:16
@ GEO_TAISEC
Definition: l1_hmodis_hdf.c:38
double interp_bilin(double val[4], double dx, double dy)
@ GEO_SCALE
Definition: l1_hmodis_hdf.c:75
char * strdup(const char *)
int32_t rdsensorinfo(int32_t sensorID, int32_t evalmask, const char *pname, void **pval)
Definition: rdsensorinfo.c:69
int load_modis_scan(const int32_t iscan, const int32_t resolution, int32_t sensorID)
#define LOOPI
int32_t nbands
int32_t * next
#define MAX_ATTERR
@ GEO_NUM_SDS
Definition: l1_hmodis_hdf.c:61
double * corr
int16_t flagindex
int32_t id
int want_verbose
void unix2yds(double usec, short *year, short *day, double *secs)
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
double interp_modis_EVdata(double val[4], double dx, double dy)
@ CIR_1KM
Definition: l1_hmodis_hdf.c:87
@ REF_NUM_SDS
Definition: l1_hmodis_hdf.c:42
double interp_modis_Longitude(double val[4], double dx, double dy)
HDF4 data type of the output SDS Default is DFNT_FLOAT32 Common types used DFNT_INT16
sds_struct * sds
int get_hdfeos_meta(int32_t sd_id, char *attribute, char *name, char *data)
Definition: hdf4_utils.c:323
@ GEO_SOLZ
Definition: l1_hmodis_hdf.c:57
instead the metadata field ProcessingEnvinronment is filled in from the output of a call to the POSIX compliant function uname from within the L1B code A small bug in L1B_Tables an incorrect comparison of RVS coefficients for TEBs to RVS coefficients for RSBs was being made This was replaced with a comparison between TEB coefficients This error never resulted in an incorrect RVS correction but did lead to recalculating the coefficients for each detector in a thermal band even if the coefficients were the same for all detectors To reduce to overall size of the reflective LUT HDF fill values were eliminated from all LUTs previously dimensioned where and where NUM_TIMES is the number of time dependent table pieces In Preprocess a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed This counter is internal only and is not yet used for any chiefly to casting of were added to make it LINUX compatible Output of code run on LINUX machines displays differences of at most scaled sector incalculable values of the Emissive calibration factor and incalculable values of SV or BB averages was moved outside the loop over frames in Emissive_Cal c since none of these quantities are frame dependent Initialization of b1 and XMS values in Preprocess c routine Process_OBCENG_Emiss was moved inside the detector loops The code was altered so that if up to five scans are dropped between the leading middle or middle trailing the leading or trailing granule will still be used in emissive calibration to form a cross granule average QA bits and are set for a gap between the leading middle and middle trailing granules respectively This may in rare instances lead to a change in emissive calibration coefficients for scans at the beginning or end of a granule A small bug in the Band correction algorithm was corrected an uncertainty value was being checked against an upper bound whereas the proper quantity to be checked was the corresponding SI
Definition: HISTORY.txt:595
double starttime
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
int32_t nlines
@ GEO_LON
Definition: l1_hmodis_hdf.c:54
int32_t scandim
Definition: l1_hmodis_hdf.c:32
#define basename(s)
Definition: l0chunk_modis.c:29
#define SCAN_TIME_INTERVAL
Definition: l1_hmodis_hdf.c:17
const char * sensorId2SensorName(int sensorId)
Definition: sensorInfo.c:175
void radiance2bt(l1str *l1rec, int resolution)
Definition: brightness.c:168
#define MIN_BAD_SI
#define fabs(a)
Definition: misc.h:93
int32_t iscan
float * sena
void fopen_warn(const char *filename, const char *file, const int32_t line)
Definition: hdf4_utils.c:258
int read_sds_1scan(sds_struct *sds, int32_t iscan, int32_t scandim, int32_t per_scan)
int16_t bandindex
void free_sds(sds_struct *sds)
Definition: hdf4_utils.c:158
int32_t npixls
int32_t nvals
int init_geo(const char filename[FILENAME_MAX])
@ L1B_NUM_COEFFS
float * lon
int readall_sds(sds_struct *sds)
Definition: hdf4_utils.c:135
data_t s[NROOTS]
Definition: decode_rs.h:75
real *8 function esdist(iyr, iday, msec)
Definition: esdist.f:3
double taisec
#define read_att(obj_id, attname, valptr)
Definition: hdf4utils.h:177
@ RAD_OFFSET
const char * subsensorId2SubsensorDir(int subsensorId)
Definition: sensorInfo.c:231
l2prod offset
#define TRYMEM(file, line, memstat)
Definition: l1_hmodis_hdf.c:20
HDF4 data type of the output SDS Default is DFNT_FLOAT32 Common types used DFNT_FLOAT32
int32_t ndets
@ GEO_SENZ
Definition: l1_hmodis_hdf.c:59
@ GEO_MNORM
Definition: l1_hmodis_hdf.c:41
void fill_dead_detectors(const int32_t iband, const modis_sds mds, double *data)
REF_SDS
Definition: l1_hmodis_hdf.c:37
double * SI_val
int i
Definition: decode_rs.h:71
int alloc_scan(const int32_t nvals)
msiBandIdx val
Definition: l1c_msi.cpp:36
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int load_att_byname(int32_t obj_id, const char *attname, att_struct *att)
Definition: hdf4_utils.c:28
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:95
#define MODISA
Definition: sensorDefs.h:19
#define EOSMETALEN
Definition: filetype.c:20
@ REFL_OFFSET
int closel1_hmodis_hdf()
int32_t nframes
@ GEO_ANGLES
Definition: l1_hmodis_hdf.c:40
float * sola
char shortname[FILENAME_MAX]
GEO_COEFFS
Definition: l1_hmodis_hdf.c:74
double isodate2unix(const char *isodate)
Definition: unix2isodate.c:61
@ REFL_SCALE