A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
NASA Logo
Ocean Color Science Software

ocssw V2022
l1_misr.c
Go to the documentation of this file.
1 #include "l1.h"
2 #include "l1_misr.h"
3 #include <hdf4utils.h>
4 
5 #include <hdf.h>
6 #include <mfhdf.h>
7 
8 #include <gsl/gsl_math.h>
9 #include <gsl/gsl_interp2d.h>
10 #include <gsl/gsl_spline2d.h>
11 #include <gsl/gsl_multifit.h>
12 
13 #include <libgen.h>
14 
15 static gsl_spline2d *spline_misr;
16 static gsl_interp_accel *xacc, *yacc;
17 static double xa[32], ya[2];
18 
19 static int32_t geoFileId;
20 static int32_t lat_id;
21 static int32_t lon_id;
22 
23 static int32_t sd_id[N_CAMERAS];
24 static int32_t red_id[N_CAMERAS];
25 static int32_t grn_id[N_CAMERAS];
26 static int32_t blu_id[N_CAMERAS];
27 
28 static int32_t nir_id[N_CAMERAS];
29 
30 static int icamera;
31 static int single_ifile=0;
32 
33 int getRadScaleFactors( char *file, double rad_scl_factors[4]);
34 int reduce_res(uint16_t rad_data[4][2048]);
35 int interp_values( double *grid_values, double interpolated_values[16][512]);
36 int interp_values_dbl( int32_t diff_offset, double *grid_values,
37  double interpolated_values[16][512]);
38 
39 int openl1_misr(filehandle *file) {
40  int32_t i;
41  int32_t sds_id;
42  int32_t gmp_id;
43  int32_t refn;
44  int32_t start[3]={0,0,0};
45  int32_t count[3]={180,8,32};
46  intn status;
47  char camera_type[N_CAMERAS][3]=
48  {"DF","CF","BF","AF","AN","AA","BA","CA","DA"};
49  char camera_type_mc[N_CAMERAS][3]=
50  {"Df","Cf","Bf","Af","An","Aa","Ba","Ca","Da"};
51  char namebuf[512];
52 
53 
54  misr_t *private_data = file->private_data;
55  //(misr_t *) calloc(1, sizeof(misr_t));
56 
57  //file->private_data = private_data;
58 
59  for (i=0; i<N_CAMERAS; i++) sd_id[i] = -1;
60 
61  if ( private_data->multipleInput == 1) {
62  // Multiple input files
63 
64  for (i=0; i<N_CAMERAS; i++) {
65  strcpy(namebuf, dirname(file->name));
66  strcat(namebuf, "/");
67  strncat(namebuf, basename(file->name), 39);
68  strcat(namebuf, camera_type[i]);
69  strcat(namebuf, basename(file->name)+41);
70  sd_id[i] = SDstart(namebuf, DFACC_RDONLY);
71  if (sd_id[i] == FAIL) {
72  fprintf(stderr, "-E- %s line %d: SDstart(%s, %d) failed.\n",
73  __FILE__, __LINE__, namebuf, DFACC_RDONLY);
74  return (HDF_FUNCTION_ERROR);
75  }
76 
77  // Read radiance scale factors
78  if (i == 0)
79  getRadScaleFactors( namebuf, private_data->radScaleFactors);
80 
81  // Get BlockTime IDs
82  private_data->fileID[i] = Hopen(namebuf, DFACC_RDONLY, 0);
83  Vstart ( private_data->fileID[i]);
84  refn = VSfind( private_data->fileID[i], "PerBlockMetadataTime");
85  private_data->blockTimeID[i] =
86  VSattach( private_data->fileID[i], refn, "r");
87  } // camera loop
88  icamera = 0;
89  // private_data->isSingleFile = 0;
90  } else {
91  // Single input file
92 
93  for (i=0; i<N_CAMERAS; i++) {
94  if ( strncmp(basename(file->name)+39, camera_type[i], 2) == 0) {
95  icamera = i;
96  single_ifile = 1;
97  // private_data->isSingleFile = 1;
98  sd_id[icamera] = SDstart(file->name, DFACC_RDONLY);
99  if (sd_id[icamera] == FAIL) {
100  fprintf(stderr, "-E- %s line %d: SDstart(%s, %d) failed.\n",
101  __FILE__, __LINE__, file->name, DFACC_RDONLY);
102  return (HDF_FUNCTION_ERROR);
103  }
104 
105  // Read radiance scale factors
106  getRadScaleFactors( file->name, private_data->radScaleFactors);
107  // Get BlockTime IDs
108  private_data->fileID[icamera] = Hopen(file->name, DFACC_RDONLY, 0);
109  Vstart ( private_data->fileID[icamera]);
110  refn = VSfind( private_data->fileID[icamera], "PerBlockMetadataTime");
111  private_data->blockTimeID[icamera] =
112  VSattach( private_data->fileID[icamera], refn, "r");
113 
114  break;
115  }
116  } // camera loop
117  }
118 
119  for (i=0; i<N_CAMERAS; i++) {
120  if (sd_id[i] == -1) continue;
121  red_id[i] =
122  SDselect(sd_id[i],SDnametoindex(sd_id[i], "Red Radiance/RDQI"));
123  grn_id[i] =
124  SDselect(sd_id[i], SDnametoindex(sd_id[i], "Green Radiance/RDQI"));
125  blu_id[i] =
126  SDselect(sd_id[i], SDnametoindex(sd_id[i], "Blue Radiance/RDQI"));
127  nir_id[i] =
128  SDselect(sd_id[i], SDnametoindex(sd_id[i], "NIR Radiance/RDQI"));
129  } // camera loop
130 
131  // Open the AGP geofile
132  //
133  // Dimensions: Block 180
134  // Row 128
135  // Col 512
136 
137  // -111 = Fill above data
138  // -222 = Fill below data
139  // -333 = Fill IPI invalid
140  // -444 = Fill to side of data
141  // -555 = Fill not processed
142  // -999 = Fill IPI error
143 
144  geoFileId = SDstart(file->geofile, DFACC_RDONLY);
145  if (geoFileId == FAIL) {
146  fprintf(stderr, "-E- %s line %d: SDstart(%s, %d) failed.\n",
147  __FILE__, __LINE__, file->geofile, DFACC_RDONLY);
148  return (HDF_FUNCTION_ERROR);
149  }
150 
151  lat_id = SDselect(geoFileId, SDnametoindex(geoFileId, "GeoLatitude"));
152  lon_id = SDselect(geoFileId, SDnametoindex(geoFileId, "GeoLongitude"));
153 
154 
155  // Open the GMP file
156  //
157  // Dimensions: Block 180
158  // Row 8
159  // Col 32
160  //
161  // Camera order: DF CF BF AF AN AA BA CA DA
162  //
163 
164  // Azimuth: 0.0 to 360.0
165  // Zenith: 0.0 to 90.0
166 
167  // Fill Value: -555
168  gmp_id = SDstart(file->gmpfile, DFACC_RDONLY);
169  if (gmp_id == FAIL) {
170  fprintf(stderr, "-E- %s line %d: SDstart(%s, %d) failed.\n",
171  __FILE__, __LINE__, file->gmpfile, DFACC_RDONLY);
172  return (HDF_FUNCTION_ERROR);
173  }
174 
175  // Read Solar Azimuth
176  sds_id = SDselect(gmp_id, SDnametoindex(gmp_id, "SolarAzimuth"));
177  status = SDreaddata(sds_id, start, NULL, count,
178  (VOIDP) &private_data->SolAzimuth);
179  if (status != 0) {
180  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
181  __LINE__, "SolAzimuth", file->gmpfile);
182  exit(1);
183  }
184  SDendaccess(sds_id);
185 
186  // Read Solar Zenith
187  sds_id = SDselect(gmp_id, SDnametoindex(gmp_id, "SolarZenith"));
188  status = SDreaddata(sds_id, start, NULL, count,
189  (VOIDP) &private_data->SolZenith);
190  if (status != 0) {
191  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
192  __LINE__, "SolZenith", file->gmpfile);
193  exit(1);
194  }
195  SDendaccess(sds_id);
196 
197  // Read Sensor Azimuth/Zenith
198  for (i=0; i<N_CAMERAS; i++) {
199  if (sd_id[i] == -1) continue;
200 
201  strcpy(namebuf, camera_type_mc[i]);
202  strcat(namebuf, "Azimuth");
203  sds_id = SDselect(gmp_id, SDnametoindex(gmp_id, namebuf));
204  status = SDreaddata(sds_id, start, NULL, count,
205  (VOIDP) &private_data->SenAzimuth[i]);
206  if (status != 0) {
207  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
208  __LINE__, namebuf, file->gmpfile);
209  exit(1);
210  }
211  SDendaccess(sds_id);
212 
213  strcpy(namebuf, camera_type_mc[i]);
214  strcat(namebuf, "Zenith");
215  sds_id = SDselect(gmp_id, SDnametoindex(gmp_id, namebuf));
216  status = SDreaddata(sds_id, start, NULL, count,
217  (VOIDP) &private_data->SenZenith[i]);
218  if (status != 0) {
219  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
220  __LINE__, namebuf, file->gmpfile);
221  exit(1);
222  }
223  SDendaccess(sds_id);
224  } // camera loop
225 
226  SDend(gmp_id);
227 
228 
229  // Read start and end block numbers
230  SDreadattr(sd_id[icamera],
231  SDfindattr(sd_id[icamera], "Start_block"),
232  (VOIDP) &private_data->startBlock);
233  SDreadattr(sd_id[icamera],
234  SDfindattr(sd_id[icamera], "End block"),
235  (VOIDP) &private_data->endBlock);
236 
237  file->nbands = 4;
238  file->npix = 512;
239  file->nscan = 128*(private_data->endBlock);
240  if (single_ifile == 0) file->nscan *= N_CAMERAS;
241 
242 
243  // Read ocean block numbers (1-based)
244  SDreadattr(sd_id[icamera],
245  SDfindattr(sd_id[icamera], "Ocean_blocks.numbers"),
246  (VOIDP) private_data->ocean_block_numbers);
247  memset(&private_data->isOceanBlock, 0, 180);
248  for (i=0; i<180; i++) {
249  int j = private_data->ocean_block_numbers[i];
250  if (j != 0) private_data->isOceanBlock[j-1] = 1;
251  }
252 
253  // Compute block offsets
254  memset(&private_data->offset, 0, 180);
255  for (i=0; i<179; i++) {
256  double p1 = private_data->SolZenith[8*i+6][15];
257  double p2 = private_data->SolZenith[8*i+7][15];
258  double f1 = private_data->SolZenith[8*(i+1)+0][15];
259  //double f2 = private_data->SolZenith[8*(i+1)+1][15];
260  double f1_l = private_data->SolZenith[8*(i+1)+0][14];
261  double f1_r = private_data->SolZenith[8*(i+1)+0][16];
262  double diff1a = fabs((p2-p1)-(f1-p2));
263  double diff1b = fabs((p2-p1)-(f1_l-p2));
264  double diff1c = fabs((p2-p1)-(f1_r-p2));
265  if (diff1a < diff1b && diff1a < diff1c)
266  private_data->offset[i] = 0;
267  if (diff1b < diff1a && diff1b < diff1c)
268  private_data->offset[i] = -1;
269  if (diff1c < diff1a && diff1c < diff1b)
270  private_data->offset[i] = +1;
271  }
272 
273 
274  // Set up interpolation grid
275 
276  // 512 / 32 = 16
277  // (16-1)/2 = 7.5
278  const gsl_interp2d_type *T = gsl_interp2d_bilinear;
279  spline_misr = gsl_spline2d_alloc(T, 32, 2);
280  xacc = gsl_interp_accel_alloc();
281  yacc = gsl_interp_accel_alloc();
282 
283  for (size_t i=0; i<32; i++) xa[i] = 7.5 + 16*i;
284  ya[0] = -0.5;
285  ya[1] = 15.5;
286 
287  return 0;
288 }
289 
290 
291 int getRadScaleFactors( char *file, double rad_scl_factors[4]) {
292 
293  int i,j;
294  int32_t file_id, vg_band_ref, refn, vg_id, vd_id;
295  int32_t nObj, tag, ref;
296  char name[100];
297  char *grpTags[4]={"BlueBand", "GreenBand", "RedBand", "NIRBand"};
298 
299  file_id = Hopen(file, DFACC_RDONLY, 0);
300  Vstart (file_id);
301 
302  // Loop over bands
303  for (i=0; i<4; i++) {
304 
305  // Find "Grid Attributes" group for each band
306  vg_band_ref = Vfind(file_id, grpTags[i]);
307  refn = vg_band_ref;
308  while (1) {
309  refn = Vgetid(file_id, refn);
310  vg_id = Vattach(file_id, refn, "r");
311  Vgetname(vg_id, name);
312  if (strcmp(name, "Grid Attributes") == 0) break;
313  Vdetach(vg_id);
314  }
315  nObj = Vntagrefs(vg_id);
316 
317  // Search within"Grid Attributes" group for "Scale factor" vdata
318  for (j=0; j<nObj; j++) {
319  Vgettagref(vg_id, j, &tag, &ref);
320  vd_id = VSattach(file_id, ref, "r");
321  VSgetname(vd_id, name);
322  if (strcmp(name, "Scale factor") == 0) {
323  VSread(vd_id, (uint8 *) &rad_scl_factors[i], 1, NO_INTERLACE);
324  VSdetach(vd_id);
325  break;
326  }
327  VSdetach(vd_id);
328  }
329  rad_scl_factors[i] /= 10;
330  } // band loop
331  Vdetach(vg_id);
332  Vend(file_id);
333  Hclose(file_id);
334 
335  return 0;
336 }
337 
338 
339 int readl1_misr(filehandle *l1file, l1str *l1rec) {
340 
341  int i,j,k, jcamera;
342 
343  int32_t start[3]={0,0,0};
344  int32_t count[3]={1,1,512};
345 
346  //static int32_t last_recnum=-1;
347  static int32_t last_recnum_gp=-1;
348 
349  //int32_t diff_offset, interp_index;
350  int32_t recnum, recnum_red, recnum_gp, block;
351 
352  ushort rad_data[4][4][2048];
353 
354  double dbl_data[2*32];
355  double dbl_xy[2*32];
356  ushort *usptr;
357  char timestring[28];
358  int32_t year, day, msec;
359  int16_t yr16, day16;
360  double sec;
361 
362  static double sla_interp_data[16][512];
363  static double slz_interp_data[16][512];
364  static double sna_interp_data[N_CAMERAS][16][512];
365  static double snz_interp_data[N_CAMERAS][16][512];
366 
367  static double xy_interp_data[2][16][512];
368  static int32_t gp_index;
369 
370  static int first=1;
371  intn status;
372 
373  misr_t *private_data = l1file->private_data;
374 
375  // Initialize to fill value
376  if (first == 1) {
377  for (i=0; i<16; i++) {
378  for (j=0; j<512; j++) {
379  sla_interp_data[i][j] = -32767;
380  slz_interp_data[i][j] = -32767;
381  }
382  }
383 
384  for (k=0; k<N_CAMERAS; k++) {
385  for (i=0; i<16; i++) {
386  for (j=0; j<512; j++) {
387  // sna_interp_data[k][i][j] = -32767;
388  //snz_interp_data[k][i][j] = -32767;
389 
390  sna_interp_data[k][i][j] = -555;
391  snz_interp_data[k][i][j] = -555;
392  }
393  }
394  }
395  }
396  first = 0;
397 
398  recnum = l1rec->iscan;
399  if (single_ifile == 0) recnum /= N_CAMERAS;
400 
401  block = recnum / 128;
402 
403  // If no ocean data set time fill and bad nav and bail
404  // if(private_data->isOceanBlock[block] == 0) {
405  // l1rec->scantime = BAD_FLT;
406  // for (i=0; i<l1rec->npix; i++) l1rec->navfail[i] = 1;
407  //return 0;
408  //}
409 
410  // Bail if record already processed
411  //if (recnum == last_recnum) last_recnum_gp = -1;
412 
413  // Get corresponding camera number for scan
414  if (single_ifile) jcamera = icamera; else jcamera = l1rec->iscan % N_CAMERAS;
415  // Get scantime (from PerBlockMetadataTime vdata)
416  //i = 0;
417  // while (1) {
418  //if (block+1 == private_data->ocean_block_numbers[i]) break;
419  //i++;
420  //}
421 
422  VSseek(private_data->blockTimeID[jcamera], block);
423  VSread(private_data->blockTimeID[jcamera], (uint8 *) timestring, 1,
424  NO_INTERLACE);
425 
426  // Convert time string (isodate) to seconds
427  isodate2ydmsec(timestring, &year, &day, &msec);
428  yr16 = year;
429  day16 = day;
430  sec = (double) msec / 1000;
431 
432  if (yr16 > 1999)
433  l1rec->scantime = yds2unix(yr16, day16, sec);
434  else
435  return 0;
437  // Read lon/lat data //
439  start[0] = block;
440  start[1] = recnum - start[0]*128;
441 
442  status = SDreaddata(lat_id, start, NULL, count, (VOIDP) l1rec->lat);
443  if (status != 0) {
444  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
445  __LINE__, "GeoLatitude", l1file->geofile);
446  exit(1);
447  }
448 
449  status = SDreaddata(lon_id, start, NULL, count, (VOIDP) l1rec->lon);
450  if (status != 0) {
451  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
452  __LINE__, "GeoLongitude", l1file->geofile);
453  exit(1);
454  }
455 
456 
458  // Generate interpolated Azimuth/Zenith values if necessary //
460  recnum_gp = (recnum - 8) / 16;
461  if (recnum_gp != last_recnum_gp && recnum >= 8) {
462 
463  gp_index = 0;
464 
465 // printf("recnum_gp: %d block: %d %d jcamera: %d\n",
466 // recnum_gp, block, recnum_gp % 8, jcamera);
467 
468 // interp_index = recnum_gp % 8;
469 // if (interp_index == 7) {
470 // diff_offset = private_data->offset[block];
471 // } else {
472 // diff_offset = 0;
473 // }
474 
475  if (single_ifile || (l1rec->iscan % N_CAMERAS) == 0) {
477  // Interpolate Solar Azimuth values //
479  for (i=0; i<32; i++) {
480  dbl_data[i] = private_data->SolAzimuth[recnum_gp+0][i];
481  if (dbl_data[i] >= 0 && dbl_data[i] < 180)
482  dbl_data[i] = dbl_data[i] + 360;
483  dbl_data[i+32] = private_data->SolAzimuth[recnum_gp+1][i];
484  if (dbl_data[i+32] >= 0 && dbl_data[i+32] < 180)
485  dbl_data[i+32] = dbl_data[i+32] + 360;
486  }
487  interp_values( dbl_data, sla_interp_data);
488 
490  // Interpolate Solar Zenith values //
492  for (i=0; i<32; i++) {
493  dbl_data[i] = private_data->SolZenith[recnum_gp+0][i];
494  dbl_data[i+32] = private_data->SolZenith[recnum_gp+1][i];
495  }
496  interp_values( dbl_data, slz_interp_data);
497  }
499  // Interpolate Sensor Azimuth values //
501 
502  // In order to avoid problems interpolating over 0/360 jump the
503  // unit vector xy-values will be used rather than the azimuth angle
504 
505  for (j=0; j<N_CAMERAS; j++) {
506  if (sd_id[j] == -1) continue;
507 
508  for (i=0; i<32; i++) {
509  dbl_xy[i] = -555;
510  if ( private_data->SenAzimuth[j][recnum_gp+0][i] > 0)
511  dbl_xy[i] = cos(private_data->SenAzimuth[j][recnum_gp+0][i]*D2R);
512  dbl_xy[i+32] = -555;
513  if ( private_data->SenAzimuth[j][recnum_gp+1][i] > 0)
514  dbl_xy[i+32] = cos(private_data->SenAzimuth[j][recnum_gp+1][i]*D2R);
515  }
516  interp_values( dbl_xy, xy_interp_data[0]);
517 
518  for (i=0; i<32; i++) {
519  dbl_xy[i] = -555;
520  if ( private_data->SenAzimuth[j][recnum_gp+0][i] > 0)
521  dbl_xy[i] = sin(private_data->SenAzimuth[j][recnum_gp+0][i]*D2R);
522  dbl_xy[i+32] = -555;
523  if ( private_data->SenAzimuth[j][recnum_gp+1][i] > 0)
524  dbl_xy[i+32] = sin(private_data->SenAzimuth[j][recnum_gp+1][i]*D2R);
525  }
526  interp_values( dbl_xy, xy_interp_data[1]);
527 
528  for (i=0; i<16; i++) {
529  for (k=0; k<512; k++) {
530  if ( xy_interp_data[1][i][k] != -555 &&
531  xy_interp_data[0][i][k] != -555)
532  sna_interp_data[j][i][k] = atan2(xy_interp_data[1][i][k],
533  xy_interp_data[0][i][k]) / D2R;
534  else
535  sna_interp_data[j][i][k] = -555;
536  }
537  }
538  } // camera loop
539 
540 
542  // Interpolate Sensor Zenith values //
544  for (j=0; j<N_CAMERAS; j++) {
545  if (sd_id[j] == -1) continue;
546  for (i=0; i<32; i++) {
547  dbl_data[i] = private_data->SenZenith[j][recnum_gp+0][i];
548  dbl_data[i+32] = private_data->SenZenith[j][recnum_gp+1][i];
549  }
550  interp_values( dbl_data, snz_interp_data[j]);
551  }
552 
553  last_recnum_gp = recnum_gp;
554 
555  } // Generate interpolated GMP values ... if (recnum_gp != last_recnum_gp)
556 
557 
559  // Process radiance data //
561  recnum_red = 4*recnum;
562  start[0] = recnum_red / 512;
563  start[1] = recnum_red - start[0]*512;
564  count[1] = 4;
565  count[2] = 2048;
566 
567  // printf("Reading radiance record: %d\n", recnum);
568 
569  // Red
570  status = SDreaddata(red_id[jcamera], start, NULL, count, (VOIDP) rad_data[2]);
571  if (status != 0) {
572  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
573  __LINE__, "Red Radiance/RDQI", l1file->name);
574  exit(1);
575  }
576  usptr = &rad_data[2][0][0];
577  for (size_t i=0; i<count[1]*count[2]; i++) *usptr++ /= 4;
578 
579  reduce_res(rad_data[2]);
580 
581  if ((jcamera+1) != 5) {
582  start[0] = recnum / 128;
583  start[1] = recnum - start[0]*128;
584  count[1] = 1;
585  count[2] = 512;
586  }
587 
588  // Blue
589  status = SDreaddata(blu_id[jcamera], start, NULL, count, (VOIDP) rad_data[0]);
590  if (status != 0) {
591  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
592  __LINE__, "Blue Radiance/RDQI", l1file->name);
593  exit(1);
594  }
595  usptr = &rad_data[0][0][0];
596  for (size_t i=0; i<count[1]*count[2]; i++) *usptr++ /= 4;
597  if ((jcamera+1) == 5) reduce_res(rad_data[0]);
598 
599 
600  // Green
601  status = SDreaddata(grn_id[jcamera], start, NULL, count, (VOIDP) rad_data[1]);
602  if (status != 0) {
603  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
604  __LINE__, "Green Radiance/RDQI", l1file->name);
605  exit(1);
606  }
607  usptr = &rad_data[1][0][0];
608  for (size_t i=0; i<count[1]*count[2]; i++) *usptr++ /= 4;
609  if ((jcamera+1) == 5) reduce_res(rad_data[1]);
610 
611  // NIR
612  status = SDreaddata(nir_id[jcamera], start, NULL, count, (VOIDP) rad_data[3]);
613  if (status != 0) {
614  printf("-E- %s line %d: Error reading SDS %s from %s.\n", __FILE__,
615  __LINE__, "NIR Radiance/RDQI", l1file->name);
616  exit(1);
617  }
618  usptr = &rad_data[3][0][0];
619  for (size_t i=0; i<count[1]*count[2]; i++) *usptr++ /= 4;
620  if ((jcamera+1) == 5) reduce_res(rad_data[3]);
621 
622 
623  // Convert scaled radiance values to floating point and copy to Lt
624  for (size_t ip=8; ip<504; ++ip) {
625  if (rad_data[0][0][ip] < 16378)
626  l1rec->Lt[4*ip+0] = rad_data[0][0][ip]*private_data->radScaleFactors[0];
627  if (rad_data[1][0][ip] < 16378)
628  l1rec->Lt[4*ip+1] = rad_data[1][0][ip]*private_data->radScaleFactors[1];
629  if (rad_data[2][0][ip] < 16378)
630  l1rec->Lt[4*ip+2] = rad_data[2][0][ip]*private_data->radScaleFactors[2];
631  if (rad_data[3][0][ip] < 16378)
632  l1rec->Lt[4*ip+3] = rad_data[3][0][ip]*private_data->radScaleFactors[3];
633  }
634 
635  // Copy Azimuth/Zenith angles to l1rec arrays
636  for (size_t ip=8; ip<504; ++ip) {
637  if (sla_interp_data[gp_index][ip] != -555) {
638  l1rec->sola[ip] = sla_interp_data[gp_index][ip];
639  if (l1rec->sola[ip] > 360) l1rec->sola[ip] = l1rec->sola[ip] - 360;
640  if (l1rec->sola[ip] > 180) l1rec->sola[ip] = l1rec->sola[ip] - 360;
641  }
642 
643  if (slz_interp_data[gp_index][ip] != -555) {
644  l1rec->solz[ip] = slz_interp_data[gp_index][ip];
645  l1rec->csolz[ip] = cos(l1rec->solz[ip] / RADEG);
646  }
647 
648  if (sna_interp_data[jcamera][gp_index][ip] != -555) {
649  l1rec->sena[ip] = sna_interp_data[jcamera][gp_index][ip];
650  if (l1rec->sena[ip] > 360) l1rec->sena[ip] = l1rec->sena[ip] - 360;
651  if (l1rec->sena[ip] > 180) l1rec->sena[ip] = l1rec->sena[ip] - 360;
652  }
653 
654  if (snz_interp_data[jcamera][gp_index][ip] != -555) {
655  l1rec->senz[ip] = snz_interp_data[jcamera][gp_index][ip];
656  l1rec->csenz[ip] = cos(l1rec->senz[ip] / RADEG);
657  }
658 
659  if (sla_interp_data[gp_index][ip] != -555 &&
660  sna_interp_data[jcamera][gp_index][ip] != -555) {
661  l1rec->delphi[ip] = l1rec->sena[ip] - 180.0 - l1rec->sola[ip];
662 
663  if (l1rec->delphi[ip] < -180.)
664  l1rec->delphi[ip] += 360.0;
665  else if (l1rec->delphi[ip] > 180.0)
666  l1rec->delphi[ip] -= 360.0;
667  }
668  }
669 
670  if (single_ifile || (l1rec->iscan % N_CAMERAS) == 8) gp_index++;
671 
672  //last_recnum = recnum;
673 
674  return 0;
675 }
676 
677 
678 int interp_values( double *grid_values, double interpolated_values[16][512]) {
679 
680  const gsl_interp2d_type *T = gsl_interp2d_bilinear;
681  static gsl_spline2d *spline;
682  static gsl_interp_accel *xacc, *yacc;
683  static double xa[32], ya[2];
684 
685  double dbl_data[2*32];
686 
687  static int first=1;
688 
689  // 512 / 32 = 16
690  // (16-1)/2 = 7.5
691 
692  if (first) {
693  // Set up interpolation grid
694  spline = gsl_spline2d_alloc(T, 32, 2);
695  xacc = gsl_interp_accel_alloc();
696  yacc = gsl_interp_accel_alloc();
697 
698  for (size_t i=0; i<32; i++) xa[i] = 7.5 + 16*i;
699  ya[0] = -0.5;
700  ya[1] = 15.5;
701 
702  first = 0;
703  }
704 
705  for (size_t i=0; i<32; i++) {
706  gsl_spline2d_set(spline, dbl_data, i, 0, grid_values[i]);
707  gsl_spline2d_set(spline, dbl_data, i, 1, grid_values[i+32]);
708  }
709 
710  gsl_spline2d_init(spline, xa, ya, dbl_data, 32, 2);
711 
712  // -111 = Fill above data
713  // -222 = Fill below data
714  // -333 = Fill IPI invalid
715  // -444 = Fill to side of data
716  // -555 = Fill not processed
717  // -999 = Fill IPI error
718 
719  for (size_t j=0; j<16; ++j) {
720  double yi = (double) j;
721  for (size_t i=8; i<504; ++i) {
722  double xi = (double) i;
723 
724  if (grid_values[(i-8)/16] == -111 || grid_values[(i-8)/16] == -222 ||
725  grid_values[(i-8)/16] == -333 || grid_values[(i-8)/16] == -444 ||
726  grid_values[(i-8)/16] == -555 || grid_values[(i-8)/16] == -999 ||
727 
728  grid_values[(i-8)/16+1] == -111 || grid_values[(i-8)/16+1] == -222 ||
729  grid_values[(i-8)/16+1] == -333 || grid_values[(i-8)/16+1] == -444 ||
730  grid_values[(i-8)/16+1] == -555 || grid_values[(i-8)/16+1] == -999 ||
731 
732  grid_values[(i-8)/16+32] == -111 || grid_values[(i-8)/16+32] == -222 ||
733  grid_values[(i-8)/16+32] == -333 || grid_values[(i-8)/16+32] == -444 ||
734  grid_values[(i-8)/16+32] == -555 || grid_values[(i-8)/16+32] == -999 ||
735 
736  grid_values[(i-8)/16+32+1] == -111 || grid_values[(i-8)/16+32+1] == -222 ||
737  grid_values[(i-8)/16+32+1] == -333 || grid_values[(i-8)/16+32+1] == -444 ||
738  grid_values[(i-8)/16+32+1] == -555 || grid_values[(i-8)/16+32+1] == -999)
739 
740  interpolated_values[j][i] = -555;
741  else
742  interpolated_values[j][i] =
743  gsl_spline2d_eval(spline, xi, yi, xacc, yacc);
744 
745 
746  }
747  }
748 
749  return 0;
750 }
751 
752 
753 int interp_values_dbl( int32_t diff_offset, double *grid_values,
754  double interpolated_values[16][512]) {
755 
756  double dbl_data[2*32];
757 
758  if (diff_offset == 0) {
759  for (size_t i=0; i<32; i++) {
760  gsl_spline2d_set(spline_misr, dbl_data, i, 0, (double) grid_values[i]);
761  gsl_spline2d_set(spline_misr, dbl_data, i, 1, (double) grid_values[i+32]);
762  }
763 
764  gsl_spline2d_init(spline_misr, xa, ya, dbl_data, 32, 2);
765 
766  for (size_t j=0; j<16; ++j) {
767  double yi = (double) j;
768  for (size_t i=8; i<504; ++i) {
769  double xi = (double) i;
770  //printf("%d %f %d %f\n", j, yi, i, xi);
771  interpolated_values[j][i] =
772  gsl_spline2d_eval(spline_misr, xi, yi, xacc, yacc);
773  }
774  }
775  } else if (diff_offset == -1) {
776 
777  for (size_t i=0; i<32; i++)
778  gsl_spline2d_set(spline_misr, dbl_data, i, 0, (double) grid_values[i]);
779  for (size_t i=0; i<31; i++)
780  gsl_spline2d_set(spline_misr, dbl_data, i+1, 1,
781  (double) grid_values[i+32]);
782 
783  for (size_t j=0; j<8; ++j) {
784  double yi = (double) j;
785  for (size_t i=8; i<504; ++i) {
786  double xi = (double) i;
787  //printf("%d %f %d %f\n", j, yi, i, xi);
788  interpolated_values[j][i] =
789  gsl_spline2d_eval(spline_misr, xi, yi, xacc, yacc);
790  }
791  }
792 
793  for (size_t i=0; i<32; i++)
794  gsl_spline2d_set(spline_misr, dbl_data, i, 0, (double) grid_values[i]);
795  for (size_t i=0; i<32; i++)
796  gsl_spline2d_set(spline_misr, dbl_data, i, 1,
797  (double) grid_values[i+32]);
798 
799  } else {
800  for (size_t i=1; i<32; i++)
801  gsl_spline2d_set(spline_misr, dbl_data, i-1, 1,
802  (double) grid_values[i]);
803 
804  for (size_t j=0; j<16; ++j) {
805  double yi = (double) j;
806  for (size_t i=8; i<504; ++i) {
807  double xi = (double) i;
808  //printf("%d %f %d %f\n", j, yi, i, xi);
809  interpolated_values[j][i] =
810  gsl_spline2d_eval(spline_misr, xi, yi, xacc, yacc);
811  }
812  }
813 
814  }
815 
816  return 0;
817 }
818 
819 
820 int reduce_res(uint16_t rad_data[4][2048]) {
821 
822  static int first=1;
823  static gsl_matrix *X, *cov;
824  static gsl_vector *y, *w, *c;
825 
826  static gsl_multifit_linear_workspace *work;
827 
828  double chisq;
829 
830  if (first) {
831 
832  X = gsl_matrix_alloc (16, 3);
833  y = gsl_vector_alloc (16);
834  w = gsl_vector_alloc (16);
835 
836  c = gsl_vector_alloc (3);
837  cov = gsl_matrix_alloc (3, 3);
838 
839  work = gsl_multifit_linear_alloc (16, 3);
840 
841  for (size_t i=0; i<16; i++) {
842  gsl_matrix_set (X, i, 0, 1.0);
843 
844  int row, col;
845  row = i / 4;
846  col = i - row*4;
847 
848  double x_val = col - 1.5;
849  double y_val = 1.5 - row;
850 
851  gsl_matrix_set (X, i, 1, x_val);
852  gsl_matrix_set (X, i, 2, y_val);
853 
854  first = 0;
855  }
856  }
857 
858  for (size_t ip=0; ip<512; ip++) {
859  int ngood = 0;
860  for (size_t irow=0; irow<4; irow++) {
861  for (size_t icol=0; icol<4; icol++) {
862  int index = irow*4+icol;
863  double d_val = (double) rad_data[irow][4*ip+icol];
864  gsl_vector_set (y, index, d_val);
865  if (rad_data[irow][4*ip+icol] >= 16378) {
866  gsl_vector_set (w, index, 0.0);
867  } else {
868  gsl_vector_set (w, index, 1.0);
869  ngood++;
870  }
871  }
872  }
873  if (ngood >= 5) {
874 
875  /*
876  double *ptr = gsl_matrix_ptr(X,0,0);
877  for (size_t i=0; i<16; i++)
878  printf("%lf %lf %lf\n", ptr[i*3+0],ptr[i*3+1],ptr[i*3+2]);
879  printf("\n");
880 
881  ptr = gsl_vector_ptr(w,0);
882  for (size_t i=0; i<16; i++) printf("%lf\n", ptr[i]);
883  printf("\n");
884 
885  ptr = gsl_vector_ptr(y,0);
886  for (size_t i=0; i<16; i++) printf("%lf\n", ptr[i]);
887  printf("\n");
888  */
889  gsl_multifit_wlinear (X, w, y, c, cov, &chisq, work);
890  double fit_val = gsl_vector_get(c,0);
891  rad_data[0][ip] = (ushort) ( fit_val + 0.5);
892  } else {
893  rad_data[0][ip] = 16378;
894  }
895  } // i loop
896 
897  return 0;
898 }
899 
900 
901 int closel1_misr(filehandle *file) {
902 
903  int i;
904 
905  SDendaccess(lat_id);
906  SDendaccess(lon_id);
907 
908  for (i=0; i<N_CAMERAS; i++) {
909  if (sd_id[i] == -1) continue;
910  SDendaccess(blu_id[i]);
911  SDendaccess(grn_id[i]);
912  SDendaccess(red_id[i]);
913  SDendaccess(nir_id[i]);
914 
915  SDend(sd_id[i]);
916  }
917  // VSdetach(vd_id);
918  //Vend(file_id);
919  //Hclose(file_id);
920 
921  return 0;
922 }
923 
924 
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
index
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
j
int j
Definition: decode_rs.h:73
day
int32_t day
Definition: atrem_corl1v2.h:308
status
int status
Definition: l1_czcs_hdf.c:32
yds2unix
double yds2unix(int16_t year, int16_t day, double secs)
Definition: yds2unix.c:7
f1
float f1(float x)
viirs_ler_luts::ref
real ref
Definition: viirs_ler_luts.f95:60
FAIL
#define FAIL
Definition: ObpgReadGrid.h:18
closel1_misr
int closel1_misr(filehandle *file)
Definition: l1_misr.c:901
NULL
#define NULL
Definition: decode_rs.h:63
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
spline
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
msec
int32 * msec
Definition: l1_czcs_hdf.c:31
N_CAMERAS
#define N_CAMERAS
Definition: l1_misr.h:4
l1_misr.h
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
openl1_misr
int openl1_misr(filehandle *file)
Definition: l1_misr.c:39
interp_values
int interp_values(double *grid_values, double interpolated_values[16][512])
Definition: l1_misr.c:678
reduce_res
int reduce_res(uint16_t rad_data[4][2048])
Definition: l1_misr.c:820
recnum
read recnum
Definition: HOWTO_Add_a_sensor.txt:72
ushort
unsigned short ushort
Definition: l1_viirs_l1b.c:16
color_dtdb.y
y
Definition: color_dtdb.py:430
D2R
#define D2R
Definition: proj_define.h:91
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
isodate2ydmsec
void isodate2ydmsec(char *date, int32_t *year, int32_t *day, int32_t *msec)
Definition: date2ydmsec.c:20
getRadScaleFactors
int getRadScaleFactors(char *file, double rad_scl_factors[4])
Definition: l1_misr.c:291
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
basename
#define basename(s)
Definition: l0chunk_modis.c:29
fabs
#define fabs(a)
Definition: misc.h:93
interp_values_dbl
int interp_values_dbl(int32_t diff_offset, double *grid_values, double interpolated_values[16][512])
Definition: l1_misr.c:753
readl1_misr
int readl1_misr(filehandle *l1file, l1str *l1rec)
Definition: l1_misr.c:339
c
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 as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
i
int i
Definition: decode_rs.h:71
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
viirs_ler_luts::p1
real p1
Definition: viirs_ler_luts.f95:37
k
int k
Definition: decode_rs.h:73
HDF_FUNCTION_ERROR
#define HDF_FUNCTION_ERROR
Definition: passthebuck.h:7
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

How to Use Earthdata Search

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

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster