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
ancgen.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <math.h>
3 #include <libgen.h>
4 
5 #include <iostream>
6 #include <fstream>
7 #include <vector>
8 #include <netcdf>
9 
10 #include <cgal_interp.h>
11 #include <allocate2d.h>
12 #include <allocate3d.h>
13 #include <timeutils.h>
14 #include <clo.h>
15 #include <genutils.h>
16 
17 #include <gsl/gsl_matrix_float.h>
18 #include <gsl/gsl_blas.h>
19 
20 #define NCACHE 20
21 const double PI = 3.14159265358979323846;
22 const double RADEG = 180 / PI;
23 
24 using namespace std;
25 using namespace netCDF;
26 using namespace netCDF::exceptions;
27 
28 // Modification history:
29 // Programmer Organization Date Ver Description of change
30 // ---------- ------------ ---- --- ---------------------
31 // Joel Gales SAIC 01/14/24 1.50 Fix type bug in the
32 // calculation of the water cloud
33 // average values for the L1C bins
34 //
35 // Joel Gales SAIC 02/01/24 1.51 Make CLD processing optional
36 // Martin Montes SSAI 07/22/24 1.52 Fix water/ice cloud mask
37 
38 int accum( uint32_t totlines, uint32_t npixels,
39  uint32_t nalong, uint32_t nacross,
40  short *brow, short *bcol,
41  float **cldprod, int8_t **cld_phase,
42  float **binval_ic, short **nobs_ic,
43  float **binval_wc, short **nobs_wc) {
44 
45  // Clear accumulation arrays
46  for (size_t i=0; i<nalong; i++) {
47  for (size_t j=0; j<nacross; j++) {
48  binval_ic[i][j] = 0.0;
49  binval_wc[i][j] = 0.0;
50  nobs_ic[i][j] = 0;
51  nobs_wc[i][j] = 0;
52  }
53  }
54 
55  // Accumulate bin values and counts
56  short *brptr=brow;
57  short *bcptr=bcol;
58  for (size_t i=0; i<totlines; i++) {
59  for (size_t j=0; j<npixels; j++) {
60  if (*brptr != -1) {
61 
62  if (cldprod[i][j] == -32767) {
63  brptr++;
64  bcptr++;
65  continue;
66  }
67 
68  if (cld_phase[i][j] == 1) {
69  binval_ic[*brptr][*bcptr] += cldprod[i][j];
70  nobs_ic[*brptr][*bcptr]++;
71  } else {
72  binval_wc[*brptr][*bcptr] += cldprod[i][j];
73  nobs_wc[*brptr][*bcptr]++;
74  }
75  }
76  brptr++;
77  bcptr++;
78  }
79  }
80 
81  return 0;
82 }
83 
84 // accumulate the ice and water cloud fraction
85 int accum_frac( uint32_t totlines, uint32_t npixels,
86  uint32_t nalong, uint32_t nacross,
87  short *brow, short *bcol,
88  float **cldprod, int8_t **cld_phase,
89  float **binval_ic, short **nobs_ic,
90  float **binval_wc, short **nobs_wc) {
91 
92  // Clear accumulation arrays
93  for (size_t i=0; i<nalong; i++) {
94  for (size_t j=0; j<nacross; j++) {
95  binval_ic[i][j] = 0.0;
96  binval_wc[i][j] = 0.0;
97  nobs_ic[i][j] = 0;
98  nobs_wc[i][j] = 0;
99  }
100  }
101 
102  // Accumulate bin values and counts
103  short *brptr=brow;
104  short *bcptr=bcol;
105  for (size_t i=0; i<totlines; i++) {
106  for (size_t j=0; j<npixels; j++) {
107  if (*brptr != -1) {
108 
109  if(cldprod[i][j] != -32767) {
110 
111  nobs_ic[*brptr][*bcptr]++;
112  nobs_wc[*brptr][*bcptr]++;
113 
114  if (cldprod[i][j] == 1) {
115  if (cld_phase[i][j] == 1) {
116  binval_ic[*brptr][*bcptr]++;
117  } else {
118  binval_wc[*brptr][*bcptr]++;
119  }
120  }
121  }
122 
123  }
124  brptr++;
125  bcptr++;
126  }
127  }
128 
129  return 0;
130 }
131 
132 
133 int accum_wm( uint32_t nwm, short *brow, short *bcol, float *wlval,
134  float **binval, short **nobs) {
135 
136  // Accumulate bin values and counts
137  short *brptr=brow;
138  short *bcptr=bcol;
139  for (size_t i=0; i<nwm; i++) {
140  if (*brptr != -1) {
141 
142  if (wlval[i] == -32767) {
143  brptr++;
144  bcptr++;
145  continue;
146  }
147 
148  binval[*brptr][*bcptr] += wlval[i];
149  nobs[*brptr][*bcptr]++;
150  }
151  brptr++;
152  bcptr++;
153  }
154 
155  return 0;
156 }
157 
158 
159 int lonlat2rowcol( uint32_t nalong, uint32_t nacross, uint32_t ncm,
160  float gridres, float *lonL1C, float *latL1C,
161  float *lonCM, float *latCM, short *brow, short *bcol);
162 
163 int copyVarAtts( NcVar *varin, NcVar *varout, string override[],
164  string overridetype[]);
165 
166 int copyVarAtts( NcVar *varin, NcVar *varout);
167 
168 template<typename T>
169 int readCLD( NcGroup ncGrp[], const char *fldname, T *array,
170  uint32_t npixels, uint32_t nlines[]) {
171 
172  vector<size_t> start, count;
173  NcVar var;
174 
175  uint32_t nlines0 = 0;
176 
177  // Read from trailing granule if specified
178  if ( !ncGrp[0].isNull()) {
179  if(nlines[0] < 250) {
180  nlines0 = nlines[0];
181  } else {
182  nlines0 = 250;
183  }
184 
185  start.clear();
186  start.push_back(nlines[0]-nlines0);
187  start.push_back(0);
188 
189  count.clear();
190  count.push_back(nlines0);
191  count.push_back(npixels);
192 
193  var = ncGrp[0].getVar( fldname);
194  var.getVar( start, count, &array[0]);
195 
196  }
197 
198  var = ncGrp[1].getVar( fldname);
199  var.getVar( &array[nlines0*npixels]);
200 
201  // Read from following granule if specified
202  if ( !ncGrp[2].isNull()) {
203  start.clear();
204  start.push_back(0);
205  start.push_back(0);
206 
207  count.clear();
208  if(nlines[2] < 250)
209  count.push_back(nlines[2]);
210  else
211  count.push_back(250);
212  count.push_back(npixels);
213 
214  var = ncGrp[2].getVar( fldname);
215  var.getVar( start, count, &array[(nlines0+nlines[1])*npixels]);
216  }
217 
218  return 0;
219 }
220 
221 inline
222 int expandEnvVar( std::string *sValue) {
223  if ( (*sValue).find_first_of( "$" ) == std::string::npos) return 0;
224  std::string::size_type posEndIdx = (*sValue).find_first_of( "/" );
225  if ( posEndIdx == std::string::npos) return 0;
226  const std::string envVar = sValue->substr (1, posEndIdx - 1);
227  char *envVar_str = getenv(envVar.c_str());
228  if (envVar_str == 0x0) {
229  printf("Environment variable: %s not defined.\n", sValue->c_str());
230  exit(1);
231  }
232  *sValue = envVar_str + (*sValue).substr( posEndIdx);
233 
234  return 0;
235 }
236 
237 
238 int main (int argc, char* argv[]) {
239 
240  NcVar var, varin, varout;
241 
242  clo_optionList_t *list;
243  clo_option_t *option;
244  char *strVal;
245  char keyword [50];
246  char targetfilename[FILENAME_MAX];
247  char chlfilename[FILENAME_MAX];
248  char profilename[FILENAME_MAX];
249  char metfilename[FILENAME_MAX];
250  char aerfilename[FILENAME_MAX];
251  char geoscffilename[FILENAME_MAX];
252  char cldmask1[FILENAME_MAX];
253  char cldmask2[FILENAME_MAX];
254  char cldmask3[FILENAME_MAX];
255  char cldprod1[FILENAME_MAX];
256  char cldprod2[FILENAME_MAX];
257  char cldprod3[FILENAME_MAX];
258  char albedofilename[FILENAME_MAX];
259  char camsch4filename[FILENAME_MAX];
260  char camsco2filename[FILENAME_MAX];
261  char camsn2ofilename[FILENAME_MAX];
262  char gebcofilename[FILENAME_MAX];
263  char ancoutfilename[FILENAME_MAX];
264 
265  list = clo_createList();
266 
267  option = clo_addOption(list, "targetfile", CLO_TYPE_IFILE, NULL,
268  "Input target file");
269  option = clo_addOption(list, "ancfile", CLO_TYPE_IFILE, NULL,
270  "Output ancillary file");
271  option = clo_addOption(list, "merraprofile", CLO_TYPE_IFILE, NULL,
272  "Input MERRA2 profile ancillary file");
273  option = clo_addOption(list, "merrametfile", CLO_TYPE_IFILE, NULL,
274  "Input MERRA2 MET ancillary file");
275  option = clo_addOption(list, "merraaerfile", CLO_TYPE_IFILE, NULL,
276  "Input MERRA2 AER ancillary file (optional)");
277  option = clo_addOption(list, "geoscffile", CLO_TYPE_IFILE, NULL,
278  "Input GEOS CF ancillary file");
279 
280  option = clo_addOption(list, "cldmask1", CLO_TYPE_IFILE, NULL,
281  "Input trailing cloud mask L2 file");
282  option = clo_addOption(list, "cldmask2", CLO_TYPE_IFILE, NULL,
283  "Input current cloud mask L2 file");
284  option = clo_addOption(list, "cldmask3", CLO_TYPE_IFILE, NULL,
285  "Input following cloud mask L2 file");
286 
287  option = clo_addOption(list, "cldprod1", CLO_TYPE_IFILE, NULL,
288  "Input trailing cloud product L2 file");
289  option = clo_addOption(list, "cldprod2", CLO_TYPE_IFILE, NULL,
290  "Input current cloud product L2 file");
291  option = clo_addOption(list, "cldprod3", CLO_TYPE_IFILE, NULL,
292  "Input following cloud product L2 file");
293 
294 
295  option = clo_addOption(list, "albedofile", CLO_TYPE_IFILE, NULL,
296  "Input albedo ancillary file");
297  option = clo_addOption(list, "chlfile", CLO_TYPE_IFILE, NULL,
298  "Input chlor_a L3 map file (optional)");
299  option = clo_addOption(list, "ch4file", CLO_TYPE_IFILE, NULL,
300  "Input CAMS CH4 ancillary file");
301  option = clo_addOption(list, "co2file", CLO_TYPE_IFILE, NULL,
302  "Input CAMS CO2 ancillary file");
303  option = clo_addOption(list, "n2ofile", CLO_TYPE_IFILE, NULL,
304  "Input CAMS N2O ancillary file");
305 
306  option = clo_addOption(list, "gebcofile", CLO_TYPE_IFILE, NULL,
307  "Input GEBCO landmask file");
308 
309  clo_setVersion2("ancgen", "1.52");
310 
311  if (argc == 1) {
312  clo_printUsage(list);
313  exit(1);
314  }
315 
316  clo_readArgs(list, argc, argv);
317 
318  strVal = clo_getString(list, "targetfile");
319  strcpy(targetfilename, strVal);
320 
321  strVal = clo_getString(list, "ancfile");
322  strcpy(ancoutfilename, strVal);
323 
324  strVal = clo_getString(list, "merraprofile");
325  strcpy(profilename, strVal);
326 
327  strVal = clo_getString(list, "merrametfile");
328  strcpy(metfilename, strVal);
329 
330  strVal = clo_getString(list, "geoscffile");
331  strcpy(geoscffilename, strVal);
332 
333  //strVal = clo_getString(list, "cldmask2");
334  //strcpy(cldmask2, strVal);
335 
336  //strVal = clo_getString(list, "cldprod2");
337  //strcpy(cldprod2, strVal);
338 
339  strVal = clo_getString(list, "albedofile");
340  strcpy(albedofilename, strVal);
341 
342  strVal = clo_getString(list, "ch4file");
343  strcpy(camsch4filename, strVal);
344 
345  strVal = clo_getString(list, "co2file");
346  strcpy(camsco2filename, strVal);
347 
348  strVal = clo_getString(list, "n2ofile");
349  strcpy(camsn2ofilename, strVal);
350 
351  int numOptions, optionId;
352 
353  numOptions = clo_getNumOptions(list);
354  for (optionId = 0; optionId < numOptions; optionId++) {
355  option = clo_getOption(list, optionId);
356 
357  strcpy(keyword, option->key);
358 
359  if (strcmp(keyword, "chlfile") == 0) {
360  if (clo_isOptionSet(option))
361  parse_file_name(clo_getOptionString(option), chlfilename);
362  else
363  strcpy(chlfilename, "");
364 
365  } else if (strcmp(keyword, "merraaerfile") == 0) {
366  if (clo_isOptionSet(option))
367  parse_file_name(clo_getOptionString(option), aerfilename);
368  else
369  strcpy(aerfilename, "");
370 
371  } else if (strcmp(keyword, "cldmask1") == 0) {
372  if (clo_isOptionSet(option))
373  parse_file_name(clo_getOptionString(option), cldmask1);
374  else
375  strcpy(cldmask1, "");
376 
377  } else if (strcmp(keyword, "cldmask2") == 0) {
378  if (clo_isOptionSet(option))
379  parse_file_name(clo_getOptionString(option), cldmask2);
380  else
381  strcpy(cldmask2, "");
382 
383  } else if (strcmp(keyword, "cldmask3") == 0) {
384  if (clo_isOptionSet(option))
385  parse_file_name(clo_getOptionString(option), cldmask3);
386  else
387  strcpy(cldmask3, "");
388 
389  } else if (strcmp(keyword, "cldprod1") == 0) {
390  if (clo_isOptionSet(option))
391  parse_file_name(clo_getOptionString(option), cldprod1);
392  else
393  strcpy(cldprod1, "");
394 
395  } else if (strcmp(keyword, "cldprod2") == 0) {
396  if (clo_isOptionSet(option))
397  parse_file_name(clo_getOptionString(option), cldprod2);
398  else
399  strcpy(cldprod2, "");
400 
401  } else if (strcmp(keyword, "cldprod3") == 0) {
402  if (clo_isOptionSet(option))
403  parse_file_name(clo_getOptionString(option), cldprod3);
404  else
405  strcpy(cldprod3, "");
406  } else if (strcmp(keyword, "gebcofile") == 0) {
407  if (clo_isOptionSet(option))
408  parse_file_name(clo_getOptionString(option), gebcofilename);
409  else
410  strcpy(gebcofilename, "$OCDATAROOT/common/gebco_ocssw_v2020.nc");
411  }
412 
413  }
414 
418 
419  // Read lon/lat of l1c
420  cout << endl << "Opening L1C file" << endl;
421  NcFile *L1Cfile = new NcFile( targetfilename, NcFile::read);
422  // get the start time
423  NcGroupAtt att = L1Cfile->getAtt("time_coverage_start");
424  if(att.isNull()) {
425  cout << "Error - Could not find time_coverage_start global attribute.\n";
426  exit(1);
427  }
428  string time_coverage_start;
429  att.getValues(time_coverage_start);
430 
431  // get the end time
432  att = L1Cfile->getAtt("time_coverage_end");
433  if(att.isNull()) {
434  cout << "Error - Could not find time_coverage_end global attribute.\n";
435  exit(1);
436  }
437  string time_coverage_end;
438  att.getValues(time_coverage_end);
439 
440  // get instrument
441  att = L1Cfile->getAtt("instrument");
442  if(att.isNull()) {
443  cout << "Error - Could not find instrument global attribute.\n";
444  exit(1);
445  }
446  string instrument;
447  att.getValues(instrument);
448 
449  NcDim across_dim = L1Cfile->getDim("bins_across_track");
450  uint32_t nacross = across_dim.getSize();
451  NcDim along_dim = L1Cfile->getDim("bins_along_track");
452  uint32_t nalong = along_dim.getSize();
453  // cout << "lines=" << nalong << ", pixels=" << nacross << endl;
454 
455  float **out_l1c_2d = allocate2d_float(nalong, nacross);
456  float **latL1C = allocate2d_float(nalong, nacross);
457  float **lonL1C = allocate2d_float(nalong, nacross);
458 
459  int bins_along_track = nalong;
460  int bins_across_track = nacross;
461 
462  size_t nl1c = bins_along_track * bins_across_track;
463 
464  vector<NcDim> dims;
465 
466  float gridres=5.2;
467 
471 
472  // Create output ancillary file
473  cout << "Creating ancillary file" << endl;
474  NcFile* nc_output;
475  nc_output = new NcFile( ancoutfilename, NcFile::replace);
476  // "PACE_SPEXONE_ANC.20170115T225222.L1C.5km.nc",
477 
478  // Global metadata
479  string date_created = string(unix2isodate(now(), 'G'));
480 
481  nc_output->putAtt( "date_created", date_created);
482  nc_output->putAtt( "time_coverage_start", time_coverage_start);
483  nc_output->putAtt( "time_coverage_end", time_coverage_end);
484  nc_output->putAtt( "title", instrument + " L1C ancillary file");
485  nc_output->putAtt( "product_name", basename(ancoutfilename));
486  nc_output->putAtt( "creator_name", "NASA/GSFC/OBPG");
487  nc_output->putAtt( "creator_url", "http://oceandata.sci.gsfc.nasa.gov");
488  nc_output->putAtt( "creator_email", "data@oceancolor.gsfc.nasa.gov");
489  nc_output->putAtt( "project",
490  "Ocean Biology Processing Group (NASA/GSFC/OBPG)");
491  nc_output->putAtt( "publisher_name", "NASA/GSFC/OBPG");
492  nc_output->putAtt( "publisher_url", "http://oceandata.sci.gsfc.nasa.gov");
493  nc_output->putAtt( "publisher_email", "data@oceancolor.gsfc.nasa.gov");
494 
495 
496  // Create netCDF dimensions
497  NcDim rows = nc_output->addDim("bins_along_track", bins_along_track);
498  NcDim cols = nc_output->addDim("bins_across_track", bins_across_track);
499 
500  // Add target lat/lon fields
501  dims.push_back(rows);
502  dims.push_back(cols);
503 
504  // Get L1C lon/lat
505  NcGroup gidGEO = L1Cfile->getGroup( "geolocation_data");
506 
507  float *fptr;
508 
509  var = gidGEO.getVar( "latitude");
510  var.getVar( &latL1C[0][0]);
511  varout = nc_output->addVar("latitude", ncFloat, dims);
512 
513  string override_ll[] = {"-32767", "=", "=", "=", "="};
514  string overridetype_ll[] = {"F", "=", "=", "=", "="};
515  copyVarAtts( &var, &varout, override_ll, overridetype_ll);
516 
517  // Convert fill values to -32767
518  fptr = &latL1C[0][0];
519  for (size_t i=0; i<nl1c; i++) {
520  if (*fptr < -900) *fptr = -32767;
521  fptr++;
522  }
523  varout.putVar( &latL1C[0][0]);
524 
525  var = gidGEO.getVar( "longitude");
526  var.getVar( &lonL1C[0][0]);
527  varout = nc_output->addVar("longitude", ncFloat, dims);
528 
529  copyVarAtts( &var, &varout, override_ll, overridetype_ll);
530 
531  // Convert fill values to -32767
532  fptr = &lonL1C[0][0];
533  for (size_t i=0; i<nl1c; i++) {
534  if (*fptr < -900) *fptr = -32767;
535  fptr++;
536  }
537  varout.putVar( &lonL1C[0][0]);
538 
539 
540  // Get lon/lat min/max for L1C granule
541  float lonmax=-180, lonmin=+180;
542  fptr = &lonL1C[0][0];
543  for (size_t i=0; i< (size_t) nl1c; i++) {
544  if (*fptr < lonmin) lonmin = *fptr;
545  if (*fptr > lonmax) lonmax = *fptr;
546  fptr++;
547  }
548 
549  float latmax=-90, latmin=+90;
550  fptr = &latL1C[0][0];
551  for (size_t i=0; i< (size_t) nl1c; i++) {
552  if (*fptr < latmin) latmin = *fptr;
553  if (*fptr > latmax) latmax = *fptr;
554  fptr++;
555  }
556 
557  vector<size_t> start, count;
558  vector<ptrdiff_t> stride;
559 
563  size_t albedo_npix = 43200;
564  size_t albedo_nlines = 21600;
565  size_t albedo_nfields = 5;
566  float albedo_pixel_size = 0.00833333;
567 
568  cout << endl << "Opening ALBEDO file" << endl;
569  NcFile *albedofile = new NcFile( albedofilename, NcFile::read);
570 
571  uint8_t **albedorow = new uint8_t*[albedo_nlines];
572 
573  std::vector<string> ALBEDOfields
574  = {"Albedo_Map_0.659", "Albedo_Map_0.858",
575  "Albedo_Map_1.24", "Albedo_Map_1.64",
576  "Albedo_Map_2.13"};
577 
578  int32_t syear;
579  int32_t sday;
580  int32_t smsec;
581  isodate2ydmsec((char*)time_coverage_start.c_str(), &syear, &sday, &smsec);
582  int32_t period_8day = sday / 8;
583  start.clear();
584  start.push_back(period_8day);
585  start.push_back(0);
586  start.push_back(0);
587 
588  count.clear();
589  count.push_back(1);
590  count.push_back(NCACHE);
591  count.push_back(albedo_npix);
592 
593  // Loop through albedo fields
594  for (size_t k=0; k<albedo_nfields; k++) {
595  var = albedofile->getVar( ALBEDOfields[k].c_str());
596 
597  for (size_t i=0; i<albedo_nlines; i++)
598  albedorow[i] = NULL;
599 
600  int first=-1;
601  for (size_t i=0; i<nalong; i++) {
602  for (size_t j=0; j<nacross; j++) {
603  // cout << "i,j: " << i << " " << j << endl;
604  int latrow = (latL1C[i][j] + 90) / albedo_pixel_size;
605 
606  if (albedorow[latrow] == NULL) {
607 
608  // align the row to the cache size boundry
609  size_t cache_boundary = ((size_t)(latrow / NCACHE)) * NCACHE;
610  if (first == -1) first = cache_boundary;
611  albedorow[cache_boundary] = new uint8_t[albedo_npix*NCACHE];
612  for (size_t k=1; k<NCACHE; k++)
613  albedorow[cache_boundary+k]
614  = albedorow[cache_boundary] + k*albedo_npix;
615 
616  start[1] = cache_boundary;
617  // cout << "Reading: " << latrow << " " << start[1] << " " << count[1]
618  // << endl;
619  var.getVar( start, count, albedorow[cache_boundary]);
620  }
621  int loncol = (lonL1C[i][j] + 180) / albedo_pixel_size;
622  // cout << "latrow, loncol: " << latrow << " " << loncol << endl;
623  uint8_t albedo = albedorow[latrow][loncol];
624  if (albedo == 255)
625  out_l1c_2d[i][j] = -32767;
626  else
627  out_l1c_2d[i][j] = albedo * 0.004;
628  }
629  }
630 
631  cout << "Writing " << ALBEDOfields[k].c_str() << endl;
632  varout = nc_output->addVar(ALBEDOfields[k].c_str(), ncFloat, dims);
633 
634  string override[] = {"-32767", "", "=", "-32767", "", "=", "="};
635  string overridetype[] = {"F", "", "=", "F", "", "=", "="};
636  copyVarAtts( &var, &varout, override, overridetype);
637 
638  varout.putVar(&out_l1c_2d[0][0]);
639 
640  for (size_t i=first; i<albedo_nlines; i+=NCACHE)
641  if (albedorow[i] != NULL) delete [] albedorow[i];
642 
643  } // albedo field loop
644 
645  delete[] albedorow;
646 
647  // SW to NE (landmask)
648 
649 
650  // Determine if L1C granule crosses the dateline
651  // Adjust L1C lon to 0-360 if dateline crossed
652  bool datelinecross=false;
653  if (lonmax - lonmin > 180) {
654  cout << endl << "Correcting longitude for dataline crossing" << endl;
655  datelinecross = true;
656  lonmax=0; lonmin=360;
657  fptr = &lonL1C[0][0];
658  for (size_t i=0; i< (size_t) nl1c; i++) {
659  if (*fptr < 0) *fptr += 360;
660  if (*fptr < lonmin) lonmin = *fptr;
661  if (*fptr > lonmax) lonmax = *fptr;
662  fptr++;
663  }
664  }
665 
666 
670 
671  // Read MERRA2 ancillary file
672  int nlon=576;
673  int nlat=361;
674 
675  cout << endl << "Opening MERRA2 file" << endl;
676  NcFile *MERRAfile = new NcFile( profilename, NcFile::read);
677  // "N202200903_PROFILE_GMAOFP_3h.nc", NcFile::read);
678 
679  float** lat_merra = allocate2d_float(361, 576);
680  float** lon_merra = allocate2d_float(361, 576);
681 
682  for (size_t i=0; i<361; i++) {
683  for (size_t j=0; j<576; j++) {
684  lat_merra[i][j] = i*0.50 - 90;
685  lon_merra[i][j] = j*0.625 - 180;
686 
687  if (datelinecross && lon_merra[i][j] < 0) lon_merra[i][j] += 360;
688  }
689  }
690  int nmerra = nlon*nlat;
691 
692  // Get nearest neighbor lon/lat bin numbers for MERRA
693  short **latbinL1C = allocate2d_short(nalong, nacross);
694  short **lonbinL1C = allocate2d_short(nalong, nacross);
695 
696  fptr = &latL1C[0][0];
697  short *slatptr = &latbinL1C[0][0];
698  for (size_t i=0; i< (size_t) nl1c; i++) {
699  *slatptr = (short) ((*fptr + 90) / 0.50);
700  fptr++;
701  slatptr++;
702  }
703 
704  fptr = &lonL1C[0][0];
705  short *slonptr = &lonbinL1C[0][0];
706  for (size_t i=0; i< (size_t) nl1c; i++) {
707  *slonptr = (short) ((*fptr + 180) / 0.625);
708  fptr++;
709  slonptr++;
710  }
711 
712  cout << "Computing MERRA to L1C interpolation mapping" << endl;
713  int nnc_tag_merra = cgal_nnc( nmerra, &lat_merra[0][0], &lon_merra[0][0],
714  nl1c, &latL1C[0][0], &lonL1C[0][0]);
715 
716  size_t nlev = 42;
717 
718  NcDim levs = nc_output->addDim("levels", nlev);
719 
720  dims.clear();
721  dims.push_back(levs);
722  dims.push_back(rows);
723  dims.push_back(cols);
724 
725  float ***out_l1c_3d = allocate3d_float(nlev, nalong, nacross);
726  float ***merra_data = allocate3d_float(42, 361, 576);
727 
728  // Geo-potential height profile
729  // Temperature vertical profile
730  // Relative humidity profile
731  // Specific humidity profile
732  // Ozone profile
733 
734  string MERRAfields[5]={"H", "T", "RH", "QV", "O3"};
735 
736  // Product loop
737  for (size_t k=0; k<5; k++) {
738 
739  varin = MERRAfile->getVar( MERRAfields[k].c_str());
740  varin.getVar(&merra_data[0][0][0]);
741 
742  // Convert MERRA2 fill values to -1e14
743  float *fptr = &merra_data[0][0][0];
744  for (size_t i=0; i<42*361*576; i++) {
745  if (*fptr > 1e14) *fptr = -1e14;
746  fptr++;
747  }
748 
749  // Level loop
750  for (size_t i=0; i<nlev; i++) {
751  cgal_interp2( nmerra,
752  &lat_merra[0][0], &lon_merra[0][0], &merra_data[i][0][0],
753  nl1c, &out_l1c_3d[i][0][0], nnc_tag_merra);
754 
755  // Convert interpolation failures to -32767
756  fptr = &out_l1c_3d[0][0][0];
757  for (size_t i=0; i<nlev*nalong*nacross; i++) {
758  if (*fptr < 0) *fptr = -32767;
759  fptr++;
760  }
761 
762  // Attempt to correct fill values with nearest-neighbor values
763  fptr = &out_l1c_3d[i][0][0];
764  slatptr = &latbinL1C[0][0];
765  slonptr = &lonbinL1C[0][0];
766 
767  for (size_t j=0; j< (size_t) nl1c; j++) {
768  if (*fptr == -32767)
769  *fptr = merra_data[i][*slatptr][*slonptr];
770 
771  fptr++;
772  slatptr++;
773  slonptr++;
774  }
775  }
776 
777  cout << "Writing " << MERRAfields[k] << endl;
778 
779  varout = nc_output->addVar(MERRAfields[k].c_str(), ncFloat, dims);
780 
781  string override[] = {"-32767", "=", "-32767", "=", "="};
782  string overridetype[] = {"F", "=", "F", "=", "="};
783  copyVarAtts( &varin, &varout, override, overridetype);
784 
785  varout.putVar(&out_l1c_3d[0][0][0]);
786 
787  } // fields loop
788 
789  delete MERRAfile;
790 
791 
795 
796  cout << endl << "Opening MERRA2 MET file" << endl;
797  MERRAfile = new NcFile( metfilename, NcFile::read);
798  // "GMAO_MERRA2.20220321T000000.MET.nc"
799 
800  dims.clear();
801  dims.push_back(rows);
802  dims.push_back(cols);
803 
804  string MERRA_MET_fields[10]={"PS", "QV10M", "SLP", "T10M", "TO3",
805  "TQV", "U10M", "V10M", "FRSNO", "FRSEAICE"};
806 
807  for (size_t k=0; k<10; k++) {
808 
809  varin = MERRAfile->getVar( MERRA_MET_fields[k].c_str());
810  varin.getVar(&merra_data[0][0][0]);
811 
812  // Convert MERRA2 fill values to -1e14
813  float *fptr = &merra_data[0][0][0];
814  for (size_t i=0; i<361*576; i++) {
815  if (*fptr > 1e14) *fptr = -1e14;
816  fptr++;
817  }
818 
819  cgal_interp2( nmerra,
820  &lat_merra[0][0], &lon_merra[0][0], &merra_data[0][0][0],
821  nl1c, &out_l1c_3d[0][0][0], nnc_tag_merra);
822 
823  // Convert interpolation failures to -32767
824  fptr = &out_l1c_3d[0][0][0];
825  for (size_t i=0; i<nalong*nacross; i++) {
826  if (*fptr < 0) *fptr = -32767;
827  fptr++;
828  }
829 
830 
831  // Attempt to correct fill values with nearest-neighbor values
832  fptr = &out_l1c_3d[0][0][0];
833  slatptr = &latbinL1C[0][0];
834  slonptr = &lonbinL1C[0][0];
835 
836  for (size_t j=0; j< (size_t) nl1c; j++) {
837  if (*fptr == -32767)
838  *fptr = merra_data[0][*slatptr][*slonptr];
839 
840  fptr++;
841  slatptr++;
842  slonptr++;
843  }
844 
845  cout << "Writing " << MERRA_MET_fields[k] << endl;
846 
847  varout = nc_output->addVar(MERRA_MET_fields[k].c_str(), ncFloat, dims);
848 
849  string override[] = {"-32767", "=", "-32767", "=", "="};
850  string overridetype[] = {"F", "=", "F", "=", "="};
851  copyVarAtts( &varin, &varout, override, overridetype);
852  varout.putVar(&out_l1c_3d[0][0][0]);
853 
854  } // fields loop
855 
856 
860 
861  if (strcmp(aerfilename, "") != 0) {
862  cout << endl << "Opening MERRA2 AER_file" << endl;
863  MERRAfile = new NcFile( aerfilename, NcFile::read);
864  // "GMAO_MERRA2.20220321T000000.AER.nc"
865 
866  dims.clear();
867  dims.push_back(rows);
868  dims.push_back(cols);
869 
870  string MERRA_AER_fields[13]={"BCEXTTAU", "BCSCATAU", "DUEXTTAU",
871  "DUSCATAU", "SSEXTTAU", "SSSCATAU",
872  "SUEXTTAU", "SUSCATAU", "OCEXTTAU",
873  "OCSCATAU", "TOTEXTTAU","TOTSCATAU",
874  "TOTANGSTR"};
875 
876  for (size_t k=0; k<13; k++) {
877 
878  varin = MERRAfile->getVar( MERRA_AER_fields[k].c_str());
879  varin.getVar(&merra_data[0][0][0]);
880 
881  // Convert MERRA2 fill values to -1e14
882  float *fptr = &merra_data[0][0][0];
883  for (size_t i=0; i<361*576; i++) {
884  if (*fptr > 1e14) *fptr = -1e14;
885  fptr++;
886  }
887 
888  cgal_interp2( nmerra,
889  &lat_merra[0][0], &lon_merra[0][0], &merra_data[0][0][0],
890  nl1c, &out_l1c_3d[0][0][0], nnc_tag_merra);
891 
892  // Convert interpolation failures to -32767
893  fptr = &out_l1c_3d[0][0][0];
894  for (size_t i=0; i<nalong*nacross; i++) {
895  if (*fptr < 0) *fptr = -32767;
896  fptr++;
897  }
898 
899 
900  // Attempt to correct fill values with nearest-neighbor values
901  fptr = &out_l1c_3d[0][0][0];
902  slatptr = &latbinL1C[0][0];
903  slonptr = &lonbinL1C[0][0];
904 
905  for (size_t j=0; j< (size_t) nl1c; j++) {
906  if (*fptr == -32767)
907  *fptr = merra_data[0][*slatptr][*slonptr];
908 
909  fptr++;
910  slatptr++;
911  slonptr++;
912  }
913 
914  cout << "Writing " << MERRA_AER_fields[k] << endl;
915 
916  varout = nc_output->addVar(MERRA_AER_fields[k].c_str(), ncFloat, dims);
917 
918  string override[] = {"-32767", "=", "-32767", "=", "="};
919  string overridetype[] = {"F", "=", "F", "=", "="};
920  copyVarAtts( &varin, &varout, override, overridetype);
921  varout.putVar(&out_l1c_3d[0][0][0]);
922 
923  } // fields loop
924 
925 
926  // cgal_release_tag( nnc_tag_merra);
927  // free2d_float(lon_merra);
928  //free2d_float(lat_merra);
929  free2d_short(lonbinL1C);
930  free2d_short(latbinL1C);
931  free3d_float(merra_data);
932  }
933 
937 
938  NcFile *CAMSfile = new NcFile( camsch4filename, NcFile::read);
939 
940  NcDim camslat_dim = CAMSfile->getDim("latitude_bins");
941  uint32_t camslat = camslat_dim.getSize();
942 
943  NcDim camslon_dim = CAMSfile->getDim("longitude_bins");
944  uint32_t camslon = camslon_dim.getSize();
945 
946  float **lat_cams = allocate2d_float(camslat, camslon);
947  float **lon_cams = allocate2d_float(camslat, camslon);
948 
949  for (size_t i=0; i<camslat; i++) {
950  for (size_t j=0; j<camslon; j++) {
951  lat_cams[i][j] = i*2 - 89.0;
952  lon_cams[i][j] = j*3 - 178.5;
953 
954  if (datelinecross && lon_cams[i][j] < 0) lon_cams[i][j] += 360;
955  }
956  }
957  int ncams = camslon*camslat;
958 
959  cout << endl <<"Computing CAMS SATELLITE to L1C interpolation mapping"
960  << endl;
961  int nnc_tag_cams = cgal_nnc( ncams, &lat_cams[0][0], &lon_cams[0][0],
962  nl1c, &latL1C[0][0], &lonL1C[0][0]);
963 
964  dims.clear();
965  dims.push_back(levs);
966  dims.push_back(rows);
967  dims.push_back(cols);
968 
969  float ***cams_data = allocate3d_float(nlev, camslat, camslon);
970 
971  int month = atoi(time_coverage_start.substr(5,2).c_str());
972 
973  string s = to_string(month);
974  unsigned int number_of_zeros = 2 - s.length();
975  s.insert(0, number_of_zeros, '0');
976  s.insert(0, "CH4_");
977 
978  varin = CAMSfile->getVar( s.c_str());
979  varin.getVar(&cams_data[0][0][0]);
980 
981  // Convert CAMS fill values to -32767
982  fptr = &cams_data[0][0][0];
983  for (size_t i=0; i<nlev*camslat*camslon; i++) {
984  if (*fptr == -999) *fptr = -32767;
985  fptr++;
986  }
987 
988  for (size_t l=0; l<42; l++) {
989  cgal_interp2( ncams,
990  &lat_cams[0][0], &lon_cams[0][0], &cams_data[l][0][0],
991  nl1c, &out_l1c_3d[l][0][0], nnc_tag_cams);
992  }
993 
994  // Convert interpolation failures to -32767
995  fptr = &out_l1c_3d[0][0][0];
996  for (size_t i=0; i<nlev*nalong*nacross; i++) {
997  if (*fptr < 0) *fptr = -32767;
998  fptr++;
999  }
1000 
1001  s.assign("CH4");
1002  cout << "Writing " << s << endl;
1003 
1004  varout = nc_output->addVar(s.c_str(), ncFloat, dims);
1005 
1006  string override[] = {"-32767", "=", "=", "=", "=", "=", "=", "="};
1007  string overridetype[] = {"F", "=", "=", "=", "=", "=", "=", "="};
1008  copyVarAtts( &varin, &varout, override, overridetype);
1009 
1010  varout.putVar(&out_l1c_3d[0][0][0]);
1011 
1012  delete CAMSfile;
1013 
1014  free2d_float(lat_cams);
1015  free2d_float(lon_cams);
1016  free3d_float(cams_data);
1017 
1018  cgal_release_tag( nnc_tag_cams);
1019 
1020 
1024 
1025  CAMSfile = new NcFile( camsco2filename, NcFile::read);
1026 
1027  camslat_dim = CAMSfile->getDim("latitude_bins");
1028  camslat = camslat_dim.getSize();
1029 
1030  camslon_dim = CAMSfile->getDim("longitude_bins");
1031  camslon = camslon_dim.getSize();
1032 
1033  lat_cams = allocate2d_float(camslat, camslon);
1034  lon_cams = allocate2d_float(camslat, camslon);
1035 
1036  for (size_t i=0; i<camslat; i++) {
1037  for (size_t j=0; j<camslon; j++) {
1038  lat_cams[i][j] = i*1.8947368421 - 90.0;
1039  lon_cams[i][j] = j*3.75 - 180.0;
1040 
1041  if (datelinecross && lon_cams[i][j] < 0) lon_cams[i][j] += 360;
1042  }
1043  }
1044  ncams = camslon*camslat;
1045 
1046  cout << endl << "Computing CAMS INST to L1C interpolation mapping"
1047  << endl;
1048  nnc_tag_cams = cgal_nnc( ncams, &lat_cams[0][0], &lon_cams[0][0],
1049  nl1c, &latL1C[0][0], &lonL1C[0][0]);
1050 
1051  cams_data = allocate3d_float(nlev, camslat, camslon);
1052 
1053  s = to_string(month);
1054  number_of_zeros = 2 - s.length();
1055  s.insert(0, number_of_zeros, '0');
1056  s.insert(0, "CO2_");
1057 
1058  varin = CAMSfile->getVar( s.c_str());
1059  varin.getVar(&cams_data[0][0][0]);
1060 
1061  // Convert CAMS fill values to -32767
1062  fptr = &cams_data[0][0][0];
1063  for (size_t i=0; i<nlev*camslat*camslon; i++) {
1064  if (*fptr == -999) *fptr = -32767;
1065  fptr++;
1066  }
1067 
1068  for (size_t l=0; l<42; l++) {
1069  cgal_interp2( ncams,
1070  &lat_cams[0][0], &lon_cams[0][0], &cams_data[l][0][0],
1071  nl1c, &out_l1c_3d[l][0][0], nnc_tag_cams);
1072  }
1073 
1074  // Convert interpolation failures to -32767
1075  fptr = &out_l1c_3d[0][0][0];
1076  for (size_t i=0; i<nlev*nalong*nacross; i++) {
1077  if (*fptr < 0) *fptr = -32767;
1078  fptr++;
1079  }
1080 
1081  s.assign("CO2");
1082  cout << "Writing " << s << endl;
1083 
1084  varout = nc_output->addVar(s.c_str(), ncFloat, dims);
1085  copyVarAtts( &varin, &varout, override, overridetype);
1086 
1087  varout.putVar(&out_l1c_3d[0][0][0]);
1088 
1089  delete CAMSfile;
1090 
1091  CAMSfile = new NcFile( camsn2ofilename, NcFile::read);
1092 
1093  s = to_string(month);
1094  number_of_zeros = 2 - s.length();
1095  s.insert(0, number_of_zeros, '0');
1096  s.insert(0, "N2O_");
1097 
1098  varin = CAMSfile->getVar( s.c_str());
1099  varin.getVar(&cams_data[0][0][0]);
1100 
1101  // Convert CAMS fill values to -32767
1102  fptr = &cams_data[0][0][0];
1103  for (size_t i=0; i<nlev*camslat*camslon; i++) {
1104  if (*fptr == -999) *fptr = -32767;
1105  fptr++;
1106  }
1107 
1108  for (size_t l=0; l<42; l++) {
1109  cgal_interp2( ncams,
1110  &lat_cams[0][0], &lon_cams[0][0], &cams_data[l][0][0],
1111  nl1c, &out_l1c_3d[l][0][0], nnc_tag_cams);
1112  }
1113 
1114  // Convert interpolation failures to -32767
1115  fptr = &out_l1c_3d[0][0][0];
1116  for (size_t i=0; i<nlev*nalong*nacross; i++) {
1117  if (*fptr < 0) *fptr = -32767;
1118  fptr++;
1119  }
1120 
1121  s.assign("N2O");
1122  cout << "Writing " << s << endl;
1123 
1124  varout = nc_output->addVar(s.c_str(), ncFloat, dims);
1125  copyVarAtts( &varin, &varout, override, overridetype);
1126 
1127  varout.putVar(&out_l1c_3d[0][0][0]);
1128 
1129  delete CAMSfile;
1130 
1131  free2d_float(lat_cams);
1132  free2d_float(lon_cams);
1133  free3d_float(cams_data);
1134 
1135  cgal_release_tag( nnc_tag_cams);
1136 
1137 
1141 
1142  // Read and process GEOS-CF data
1143 
1144  cout << endl << "Opening GEOS-CF file" << endl;
1145  NcFile *GEOSfile = new NcFile( geoscffilename, NcFile::read);
1146  //"GMAO_GEOS-CF.20210428T020000.PROFILE.nc"
1147 
1148  nlon=1440;
1149  nlat=721;
1150  int ngeos = nlon*nlat;
1151 
1152  float **lat_geos = allocate2d_float(721, 1440);
1153  float **lon_geos = allocate2d_float(721, 1440);
1154 
1155  // Compute lon/lat geos
1156  for (size_t i=0; i<721; i++) {
1157  for (size_t j=0; j<1440; j++) {
1158  lat_geos[i][j] = i*0.25 - 90;
1159  lon_geos[i][j] = j*0.25 - 180;
1160 
1161  if (datelinecross && lon_geos[i][j] < 0) lon_geos[i][j] += 360;
1162  }
1163  }
1164 
1165  // Write level field
1166  dims.clear();
1167  dims.push_back(levs);
1168  varout = nc_output->addVar("levels", ncFloat, dims);
1169 
1170  float levels[nlev];
1171  varin = GEOSfile->getVar( "lev");
1172  varin.getVar(levels);
1173  copyVarAtts( &varin, &varout);
1174  varout.putVar(levels);
1175 
1176  float ***geos_data = allocate3d_float(42, 721, 1440);
1177 
1178  cout << "Computing GEOS to L1C interpolation mapping" << endl;
1179  int nnc_tag_geos = cgal_nnc( ngeos, &lat_geos[0][0], &lon_geos[0][0],
1180  nl1c, &latL1C[0][0], &lonL1C[0][0]);
1181 
1182  string GEOSfields[5]={"NO2", "SO2",
1183  "TOTCOL_NO2", "STRATCOL_NO2", "TROPCOL_NO2"};
1184 
1185  for (size_t k=0; k<5; k++) {
1186 
1187  varin = GEOSfile->getVar( GEOSfields[k].c_str());
1188  varin.getVar(&geos_data[0][0][0]);
1189 
1190  if (GEOSfields[k].compare("NO2") == 0 ||
1191  GEOSfields[k].compare("SO2") == 0)
1192  nlev = 42; else nlev = 1;
1193 
1194  // Convert GOES fill values to -1e14
1195  float *fptr = &geos_data[0][0][0];
1196  for (size_t i=0; i<nlev*721*1440; i++) {
1197  if (*fptr > 1e14) *fptr = -1e14;
1198  fptr++;
1199  }
1200 
1201  for (size_t i=0; i<nlev; i++) {
1202  cgal_interp2( ngeos,
1203  &lat_geos[0][0], &lon_geos[0][0], &geos_data[i][0][0],
1204  nl1c, &out_l1c_3d[i][0][0], nnc_tag_geos);
1205  }
1206 
1207  // Convert interpolation failures to -32767
1208  fptr = &out_l1c_3d[0][0][0];
1209  for (size_t i=0; i<nlev*nalong*nacross; i++) {
1210  if (*fptr < 0) *fptr = -32767;
1211  fptr++;
1212  }
1213 
1214  dims.clear();
1215  if (GEOSfields[k].compare("NO2") == 0 ||
1216  GEOSfields[k].compare("SO2") == 0) {
1217  dims.push_back(levs);
1218  dims.push_back(rows);
1219  dims.push_back(cols);
1220  } else {
1221  dims.push_back(rows);
1222  dims.push_back(cols);
1223  }
1224 
1225  cout << "Writing " << GEOSfields[k].c_str() << endl;
1226 
1227  varout = nc_output->addVar(GEOSfields[k].c_str(), ncFloat, dims);
1228 
1229  string overrideGEOS[] = {"-32767", "=", "-32767", "=", "="};
1230  string overridetypeGEOS[] = {"F", "=", "F", "=", "="};
1231  copyVarAtts( &varin, &varout, overrideGEOS, overridetypeGEOS);
1232 
1233  varout.putVar(&out_l1c_3d[0][0][0]);
1234  } // fields loop
1235 
1236  cgal_release_tag( nnc_tag_geos);
1237  free2d_float(lat_geos);
1238  free2d_float(lon_geos);
1239  free3d_float(geos_data);
1240  free3d_float(out_l1c_3d);
1241 
1242 
1246 
1247  if ( strcmp(cldmask2, "") != 0 &&
1248  strcmp(cldprod2, "") != 0) {
1249 
1250  cout << endl << "Opening CLDMSK files" << endl;
1251 
1252  NcFile *CLDMSKfile[3]={NULL, NULL, NULL};
1253 
1254  uint32_t nlines[3]={0,0,0};
1255  uint32_t npixels;
1256  NcDim lines_dim;
1257  NcDim pixel_dim;
1258 
1259  uint32_t nlines0 = 0;
1260 
1261  // Open trailing/following L2 CLDMASK datasets
1262  if ( strcmp(cldmask1, "") != 0) {
1263  CLDMSKfile[0] = new NcFile( cldmask1, NcFile::read);
1264  lines_dim = CLDMSKfile[0]->getDim("number_of_lines");
1265  nlines[0] = lines_dim.getSize();
1266  }
1267 
1268  if ( strcmp(cldmask3, "") != 0) {
1269  CLDMSKfile[2] = new NcFile( cldmask3, NcFile::read);
1270  lines_dim = CLDMSKfile[2]->getDim("number_of_lines");
1271  nlines[2] = lines_dim.getSize();
1272  }
1273 
1274  // Open current L2 LDMASK dataset
1275  CLDMSKfile[1] = new NcFile( cldmask2, NcFile::read);
1276  lines_dim = CLDMSKfile[1]->getDim("number_of_lines");
1277  nlines[1] = lines_dim.getSize();
1278 
1279  pixel_dim = CLDMSKfile[1]->getDim("pixels_per_line");
1280  npixels = pixel_dim.getSize();
1281 
1282  // Use only 250 (or less) lines from trailing/following datasets
1283  uint32_t totlines = nlines[1];
1284  if (nlines[0] != 0) {
1285  if (nlines[0] < 250)
1286  totlines += nlines[0];
1287  else
1288  totlines += 250;
1289  }
1290 
1291  if (nlines[2] != 0) {
1292  if (nlines[2] < 250)
1293  totlines += nlines[2];
1294  else
1295  totlines += 250;
1296  }
1297 
1298  float **latCM = allocate2d_float(totlines, npixels);
1299  float **lonCM = allocate2d_float(totlines, npixels);
1300  int8_t **adj_mask = allocate2d_schar(totlines, npixels);
1301 
1302  for (size_t i=0; i<totlines; i++) {
1303  for (size_t j=0; j<npixels; j++) {
1304  latCM[i][j] = 0.0;
1305  lonCM[i][j] = 0.0;
1306  adj_mask[i][j] = 0;
1307  }
1308  }
1309 
1310  NcGroup ncGroup;
1311  // Read from trailing granule if specified
1312  if ( CLDMSKfile[0] != NULL) {
1313  if(nlines[0] < 250) {
1314  nlines0 = nlines[0];
1315  start[0] = 0;
1316  } else {
1317  nlines0 = 250;
1318  start[0] = nlines[0] - 250;
1319  }
1320  count[0] = nlines0;
1321  start[1] = 0;
1322  count[1] = npixels;
1323 
1324  ncGroup = CLDMSKfile[0]->getGroup("navigation_data");
1325  if(ncGroup.isNull()) {
1326  var = CLDMSKfile[0]->getVar( "latitude");
1327  var.getVar( start, count, &latCM[0][0]);
1328 
1329  var = CLDMSKfile[0]->getVar( "longitude");
1330  var.getVar( start, count, &lonCM[0][0]);
1331 
1332  var = CLDMSKfile[0]->getVar( "ADJ_MASK");
1333  var.getVar( start, count, &adj_mask[0][0]);
1334  } else {
1335  var = ncGroup.getVar( "latitude");
1336  var.getVar( start, count, &latCM[0][0]);
1337 
1338  var = ncGroup.getVar( "longitude");
1339  var.getVar( start, count, &lonCM[0][0]);
1340 
1341  ncGroup = CLDMSKfile[0]->getGroup("geophysical_data");
1342  var = ncGroup.getVar("cloud_flag");
1343  var.getVar( start, count, &adj_mask[0][0]);
1344  }
1345  }
1346 
1347  ncGroup = CLDMSKfile[1]->getGroup("navigation_data");
1348  if(ncGroup.isNull()) {
1349  var = CLDMSKfile[1]->getVar( "latitude");
1350  var.getVar( &latCM[nlines0][0]);
1351 
1352  var = CLDMSKfile[1]->getVar( "longitude");
1353  var.getVar( &lonCM[nlines0][0]);
1354 
1355  var = CLDMSKfile[1]->getVar( "ADJ_MASK");
1356  var.getVar( &adj_mask[nlines0][0]);
1357  } else {
1358  var = ncGroup.getVar( "latitude");
1359  var.getVar( &latCM[nlines0][0]);
1360 
1361  var = ncGroup.getVar( "longitude");
1362  var.getVar( &lonCM[nlines0][0]);
1363 
1364  ncGroup = CLDMSKfile[1]->getGroup("geophysical_data");
1365  var = ncGroup.getVar( "cloud_flag");
1366  var.getVar( &adj_mask[nlines0][0]);
1367  }
1368 
1369  // Read from following granule if specified
1370  if ( CLDMSKfile[2] != NULL) {
1371  start[0] = 0;
1372  if(nlines[2] < 250) {
1373  count[0] = nlines[2];
1374  } else {
1375  count[0] = 250;
1376  }
1377  start[1] = 0;
1378  count[1] = npixels;
1379 
1380  ncGroup = CLDMSKfile[2]->getGroup("navigation_data");
1381  if(ncGroup.isNull()) {
1382  var = CLDMSKfile[2]->getVar( "latitude");
1383  var.getVar( start, count, &latCM[nlines0+nlines[1]][0]);
1384 
1385  var = CLDMSKfile[2]->getVar( "longitude");
1386  var.getVar( start, count, &lonCM[nlines0+nlines[1]][0]);
1387 
1388  var = CLDMSKfile[2]->getVar( "ADJ_MASK");
1389  var.getVar( start, count, &adj_mask[nlines0+nlines[1]][0]);
1390  } else {
1391  var = ncGroup.getVar( "latitude");
1392  var.getVar( start, count, &latCM[nlines0+nlines[1]][0]);
1393 
1394  var = ncGroup.getVar( "longitude");
1395  var.getVar( start, count, &lonCM[nlines0+nlines[1]][0]);
1396 
1397  ncGroup = CLDMSKfile[2]->getGroup("geophysical_data");
1398  var = ncGroup.getVar( "cloud_flag");
1399  var.getVar( start, count, &adj_mask[nlines0+nlines[1]][0]);
1400  }
1401  }
1402 
1403  uint32_t ncm = totlines*npixels;
1404 
1405 
1406  // Determine L1C row/col of L2 CLD lon/lat
1407  // Ported from code developed by F.Patt
1408 
1409  short *brow = new short[ncm];
1410  short *bcol = new short[ncm];
1411  lonlat2rowcol( nalong, nacross, ncm, gridres, &lonL1C[0][0], &latL1C[0][0],
1412  &lonCM[0][0], &latCM[0][0], brow, bcol);
1413 
1414  ofstream tempOut;
1415 
1416  // tempOut.open ("rc.bin", ios::out | ios::trunc | ios::binary);
1417  //tempOut.write((char *) brow, sizeof(short)*ncm);
1418  //tempOut.write((char *) bcol, sizeof(short)*ncm);
1419  //tempOut.close();
1420 
1421 
1422  // Open CLD products file
1423 
1424  // Note: CLD files uses same lon/lat as CLDMSK
1425 
1426  NcFile *CLDfile[3];
1427  NcGroup ncGrp[3];
1428 
1429  if ( strcmp(cldprod1, "") != 0) {
1430  CLDfile[0] = new NcFile( cldprod1, NcFile::read);
1431  ncGrp[0] = CLDfile[0]->getGroup( "geophysical_data");
1432  }
1433 
1434  if ( strcmp(cldprod3, "") != 0) {
1435  CLDfile[2] = new NcFile( cldprod3, NcFile::read);
1436  ncGrp[2] = CLDfile[2]->getGroup( "geophysical_data");
1437  }
1438 
1439  CLDfile[1] = new NcFile( cldprod2, NcFile::read);
1440  ncGrp[1] = CLDfile[1]->getGroup( "geophysical_data");
1441 
1442  int8_t **cld_phase = allocate2d_schar(totlines, npixels);
1443 
1444  // Read cloud phase
1445  readCLD<int8_t>( ncGrp, "cld_phase_21", &cld_phase[0][0], npixels, nlines);
1446 
1447  // Set ice cloud pixels to 1 otherwise 0
1448  for (size_t i=0; i<totlines; i++)
1449  for (size_t j=0; j<npixels; j++)
1450  if (cld_phase[i][j] == 3)
1451  cld_phase[i][j] = 1;
1452  else
1453  cld_phase[i][j] = 0;
1454 
1455  // Allocate bin variables
1456  float **binval_ic = allocate2d_float(nalong, nacross);
1457  float **binval_wc = allocate2d_float(nalong, nacross);
1458  short **nobs_ic = allocate2d_short(nalong, nacross);
1459  short **nobs_wc = allocate2d_short(nalong, nacross);
1460 
1461  float **cldprod = allocate2d_float(totlines, npixels);
1462 
1463  // Average ADJ_MASK for water/ice cloud pixels
1464 
1465  // Convert BYTE to FLOAT
1466  for (size_t i=0; i<totlines; i++)
1467  for (size_t j=0; j<npixels; j++)
1468  if(adj_mask[i][j] == -128)
1469  cldprod[i][j] = -32767;
1470  else
1471  cldprod[i][j] = (float) adj_mask[i][j];
1472 
1473  free2d_schar(adj_mask);
1474 
1475  // tempOut.open ("adj_mask.bin", ios::out | ios::trunc | ios::binary);
1476  //tempOut.write((char *) &cldprod[0][0], sizeof(float)*totlines*npixels);
1477  //tempOut.close();
1478 
1479  accum_frac( totlines, npixels, nalong, nacross, brow, bcol, cldprod, cld_phase,
1480  binval_ic, nobs_ic, binval_wc, nobs_wc);
1481 
1482  // Compute average bin values for ice clouds
1483  for (size_t i=0; i<nalong; i++) {
1484  for (size_t j=0; j<nacross; j++) {
1485  if (nobs_ic[i][j] != 0)
1486  out_l1c_2d[i][j] = binval_ic[i][j] / nobs_ic[i][j];
1487  else
1488  out_l1c_2d[i][j] = -32767;
1489  }
1490  }
1491 
1492  dims.clear();
1493  dims.push_back(rows);
1494  dims.push_back(cols);
1495 
1496  cout << "Writing ice cloud fraction" << endl;
1497  varout = nc_output->addVar("ice_cloud_fraction", ncFloat, dims);
1498 
1499  string override_icf[] =
1500  {"-32767", "", "", "", "", "Ice cloud fraction", "", "1.0", "0.0"};
1501  string overridetype_icf[] = {"F", "", "", "", "", "=", "", "F", "F"};
1502  copyVarAtts( &var, &varout, override_icf, overridetype_icf);
1503 
1504  varout.putVar(&out_l1c_2d[0][0]);
1505 
1506  // Compute average bin values for water clouds
1507  for (size_t i=0; i<nalong; i++) {
1508  for (size_t j=0; j<nacross; j++) {
1509  if (nobs_wc[i][j] != 0)
1510  out_l1c_2d[i][j] = binval_wc[i][j] / nobs_wc[i][j];
1511  else
1512  out_l1c_2d[i][j] = -32767;
1513  }
1514  }
1515 
1516  cout << "Writing water cloud fraction" << endl;
1517  varout = nc_output->addVar("water_cloud_fraction", ncFloat, dims);
1518 
1519  string override_wcf[] =
1520  {"-32767", "", "", "", "", "Water cloud fraction", "", "1.0", "0.0"};
1521  string overridetype_wcf[] = {"F", "", "", "", "", "=", "", "F", "F"};
1522  copyVarAtts( &var, &varout, override_wcf, overridetype_wcf);
1523 
1524  varout.putVar(&out_l1c_2d[0][0]);
1525 
1526 
1527  // Cloud Top Pressure (CTP)
1528  // Cloud Top Temperature (CTT)
1529  // Cloud Top Height (CTH)
1530  // Particle effective radius cer_21 (micron)
1531  // COT cot_21
1532 
1533  string CLDfields[5]={"ctp", "ctt", "cth", "cer_21", "cot_21"};
1534 
1535  for (size_t k=0; k<5; k++) {
1536  cout << "Writing " << CLDfields[k].c_str() << endl;
1537 
1538  var = ncGrp[1].getVar( CLDfields[k].c_str());
1539 
1540  for (size_t i=0; i<totlines; i++)
1541  for (size_t j=0; j<npixels; j++) cldprod[i][j] = 0.0;
1542 
1543  readCLD<float>( ncGrp, CLDfields[k].c_str(), &cldprod[0][0],
1544  npixels, nlines);
1545 
1546  accum( totlines, npixels, nalong, nacross, brow, bcol, cldprod, cld_phase,
1547  binval_ic, nobs_ic, binval_wc, nobs_wc);
1548 
1549  // Compute average bin values for ice clouds
1550  for (size_t i=0; i<nalong; i++) {
1551  for (size_t j=0; j<nacross; j++) {
1552  if (nobs_ic[i][j] != 0)
1553  out_l1c_2d[i][j] = binval_ic[i][j] / nobs_ic[i][j];
1554  else
1555  out_l1c_2d[i][j] = -32767;
1556  }
1557  }
1558 
1559  string outname;
1560 
1561  outname.assign(CLDfields[k]);
1562  outname.append("_ice_cloud");
1563 
1564  varout = nc_output->addVar(outname.c_str(), ncFloat, dims);
1565  copyVarAtts( &var, &varout);
1566  varout.putVar(&out_l1c_2d[0][0]);
1567 
1568  // Compute average bin values for water clouds
1569  for (size_t i=0; i<nalong; i++) {
1570  for (size_t j=0; j<nacross; j++) {
1571  if (nobs_wc[i][j] != 0)
1572  out_l1c_2d[i][j] = binval_wc[i][j] / nobs_wc[i][j];
1573  else
1574  out_l1c_2d[i][j] = -32767;
1575  }
1576  }
1577 
1578  outname.assign(CLDfields[k]);
1579  outname.append("_water_cloud");
1580 
1581  varout = nc_output->addVar(outname.c_str(), ncFloat, dims);
1582  copyVarAtts( &var, &varout);
1583  varout.putVar(&out_l1c_2d[0][0]);
1584  } // CLD field loop
1585 
1586  delete [] brow;
1587  delete [] bcol;
1588 
1589  free2d_float(latCM);
1590  free2d_float(lonCM);
1591  free2d_float(cldprod);
1592 
1593  free2d_float(binval_ic);
1594  free2d_float(binval_wc);
1595  free2d_short(nobs_ic);
1596  free2d_short(nobs_wc);
1597  } // if CLD processing
1598 
1599 
1603 
1604  float wm_pixel_size = 360.0 / 86400;
1605  int wm_latrow = (latmin + 90) / wm_pixel_size;
1606  int wm_loncol = (lonmin + 180) / wm_pixel_size;
1607 
1608  //cout << latmin << " " << latmax << " " << lonmin << " " << lonmax << endl;
1609  //cout << wm_latrow << " " << wm_loncol << endl;
1610 
1611  start.clear();
1612  start.push_back(wm_latrow);
1613  start.push_back(wm_loncol);
1614 
1615  size_t n_wm_latrow = ((latmax + 90) / wm_pixel_size - wm_latrow + 1)/4;
1616  size_t n_wm_loncol = ((lonmax + 180) / wm_pixel_size - wm_loncol + 1)/4;
1617 
1618  count.clear();
1619  count.push_back(n_wm_latrow);
1620  count.push_back(n_wm_loncol);
1621 
1622  stride.clear();
1623  stride.push_back(4);
1624  stride.push_back(4);
1625 
1626  float *lon_wm = new float[n_wm_latrow*n_wm_loncol];
1627  float *lat_wm = new float[n_wm_latrow*n_wm_loncol];
1628 
1629  //cout << n_wm_latrow << " " << n_wm_loncol << endl;
1630 
1631  int k=0;
1632  for (size_t i=0; i<(size_t) n_wm_latrow; i++) {
1633  for (size_t j=0; j<(size_t) n_wm_loncol; j++) {
1634  lat_wm[k] = (4*i+wm_latrow)*wm_pixel_size - 90;
1635  lon_wm[k] = (4*j+wm_loncol)*wm_pixel_size - 180;
1636 
1637  if (datelinecross && lon_wm[k] < 0) lon_wm[k] += 360;
1638  k++;
1639  }
1640  }
1641 
1642  uint32_t nwm = n_wm_latrow * n_wm_loncol;
1643 
1644  // Reading landmask
1645  string gebcofilestring = gebcofilename;
1646  expandEnvVar(&gebcofilestring);
1647 
1648  NcFile *WMfile = new NcFile( gebcofilestring.c_str(), NcFile::read);
1649 
1650  float *wm = new float[n_wm_latrow*n_wm_loncol];
1651  uint8_t *wmbyte = new uint8_t[n_wm_latrow*n_wm_loncol];
1652  var = WMfile->getVar( "watermask");
1653  if (datelinecross) {
1654  uint8_t *wm0 = new uint8_t[n_wm_latrow*n_wm_loncol];
1655  count[1] = (86400 - wm_loncol) / 4;
1656  var.getVar( start, count, stride, &wm0[0]);
1657  for (size_t i=0; i<(size_t) n_wm_latrow; i++)
1658  memcpy(&wmbyte[i*n_wm_loncol], &wm0[i*count[1]],
1659  count[1]*sizeof(uint8_t));
1660 
1661  start[1] = 0;
1662  count[1] = n_wm_loncol - count[1];
1663  //cout << count[1] << endl;
1664  var.getVar( start, count, stride, &wm0[0]);
1665  for (size_t i=0; i<(size_t) n_wm_latrow; i++)
1666  memcpy(&wmbyte[i*n_wm_loncol+(86400-wm_loncol)/4],
1667  &wm0[i*count[1]], count[1]*sizeof(uint8_t));
1668  delete [] wm0;
1669  } else {
1670  var.getVar( start, count, stride, &wmbyte[0]);
1671  }
1672 
1673  // Convert byte to float
1674  k=0;
1675  for (size_t i=0; i<(size_t) n_wm_latrow; i++) {
1676  for (size_t j=0; j<(size_t) n_wm_loncol; j++) {
1677  wm[k] = (float) wmbyte[k];
1678  k++;
1679  }
1680  }
1681 
1682 
1683  // Determine L1C row/col of L2 CLD lon/lat
1684 
1685  short *brow = new short[nwm];
1686  short *bcol = new short[nwm];
1687 
1688  float **binval_wm = allocate2d_float(nalong, nacross);
1689  short **nobs_wm = allocate2d_short(nalong, nacross);
1690 
1691  // Clear accumulation arrays
1692  for (size_t i=0; i<nalong; i++) {
1693  for (size_t j=0; j<nacross; j++) {
1694  binval_wm[i][j] = 0.0;
1695  nobs_wm[i][j] = 0;
1696  }
1697  }
1698 
1699  size_t N = 2147483648 / (n_wm_loncol*nalong);
1700  size_t M = n_wm_latrow/N;
1701  if (n_wm_latrow % N != 0) M++;
1702  for (size_t i=0; i<M; i++) {
1703  // cout << i*N << " " << n_wm_latrow << endl;
1704  size_t L=N;
1705  if (n_wm_latrow - i*N < N) L = n_wm_latrow - i*N;
1706  // cout << "L: " << L << endl;
1707  lonlat2rowcol( nalong, nacross, L*n_wm_loncol, gridres,
1708  &lonL1C[0][0], &latL1C[0][0],
1709  &lon_wm[i*N*n_wm_loncol], &lat_wm[i*N*n_wm_loncol],
1710  brow, bcol);
1711 
1712  accum_wm( L*n_wm_loncol, brow, bcol, &wm[i*N*n_wm_loncol],
1713  binval_wm, nobs_wm);
1714  }
1715 
1716  for (size_t i=0; i<nalong; i++) {
1717  for (size_t j=0; j<nacross; j++) {
1718 
1719  if (nobs_wm[i][j] != 0)
1720  out_l1c_2d[i][j] = binval_wm[i][j] / nobs_wm[i][j];
1721  else
1722  out_l1c_2d[i][j] = -32767;
1723  }
1724  }
1725 
1726  cout << "Writing waterfraction" << endl;
1727  varout = nc_output->addVar("waterfraction", ncFloat, dims);
1728 
1729  string override_wf[] =
1730  {"-32767", "=", "Water fraction", "waterfraction", "1.0", "0.0"};
1731  string overridetype_wf[] = {"F", "", "", "=", "F", "F"};
1732  copyVarAtts( &var, &varout, override_wf, overridetype_wf);
1733  varout.putVar(&out_l1c_2d[0][0]);
1734 
1735  delete [] lat_wm;
1736  delete [] lon_wm;
1737  delete [] wm;
1738 
1739  delete [] brow;
1740  delete [] bcol;
1741 
1742  free2d_float(binval_wm);
1743  free2d_short(nobs_wm);
1744 
1745 
1749 
1750  k=0;
1751 
1752  if (strcmp(chlfilename, "") != 0) {
1753  float chl_pixel_size = 360.0 / 8640;
1754  int chl_latrow = (90 - latmax) / chl_pixel_size;
1755  int chl_loncol = (lonmin + 180) / chl_pixel_size;
1756 
1757  cout << latmin << " " << latmax << " " << lonmin << " " << lonmax << endl;
1758  cout << chl_latrow << " " << chl_loncol << endl;
1759 
1760  start.clear();
1761  start.push_back(chl_latrow);
1762  start.push_back(chl_loncol);
1763 
1764  int n_chl_latrow = (90 - latmin) / chl_pixel_size - chl_latrow + 1;
1765  int n_chl_loncol = (lonmax + 180) / chl_pixel_size - chl_loncol + 1;
1766 
1767  count.clear();
1768  count.push_back(n_chl_latrow);
1769  count.push_back(n_chl_loncol);
1770 
1771  float *lon_chl = new float[n_chl_latrow*n_chl_loncol];
1772  float *lat_chl = new float[n_chl_latrow*n_chl_loncol];
1773 
1774  cout << n_chl_latrow << " " << n_chl_loncol << endl;
1775 
1776  k=0;
1777  for (size_t i=0; i<(size_t) n_chl_latrow; i++) {
1778  for (size_t j=0; j<(size_t) n_chl_loncol; j++) {
1779  lat_chl[k] = 90 - (i+chl_latrow)*chl_pixel_size;
1780  lon_chl[k] = (j+chl_loncol)*chl_pixel_size - 180;
1781 
1782  if (datelinecross && lon_chl[k] < 0) lon_chl[k] += 360;
1783  k++;
1784  }
1785  }
1786 
1787  int nchl = n_chl_latrow * n_chl_loncol;
1788 
1789  cout << "Computing CHLOR_A map to L1C interpolation mapping" << endl;
1790  int nnc_tag_chl = cgal_nnc( nchl, &lat_chl[0], &lon_chl[0],
1791  nl1c, &latL1C[0][0], &lonL1C[0][0]);
1792 
1793  // Reading chlor_a from L3M file
1794  NcFile *CHLfile = new NcFile( chlfilename, NcFile::read);
1795 
1796  float *chl = new float[n_chl_latrow*n_chl_loncol];
1797  var = CHLfile->getVar( "chlor_a");
1798 
1799  if (datelinecross) {
1800  float *chl0 = new float[n_chl_latrow*n_chl_loncol];
1801  count[1] = 8640 - chl_loncol;
1802  var.getVar( start, count, &chl0[0]);
1803  for (size_t i=0; i<(size_t) n_chl_latrow; i++)
1804  memcpy(&chl[i*n_chl_loncol], &chl0[i*count[1]], count[1]*sizeof(float));
1805 
1806  start[1] = 0;
1807  count[1] = n_chl_loncol + chl_loncol - 8640;
1808  cout << count[1] << endl;
1809  var.getVar( start, count, &chl0[0]);
1810  for (size_t i=0; i<(size_t) n_chl_latrow; i++)
1811  memcpy(&chl[i*n_chl_loncol+(8640-chl_loncol)],
1812  &chl0[i*count[1]], count[1]*sizeof(float));
1813  delete [] chl0;
1814  } else {
1815  var.getVar( start, count, &chl[0]);
1816  }
1817 
1818  cgal_interp2( nchl, &lat_chl[0], &lon_chl[0], &chl[0],
1819  nl1c, &out_l1c_2d[0][0], nnc_tag_chl);
1820 
1821  // Convert interpolation failures to -32767
1822  fptr = &out_l1c_2d[0][0];
1823  for (size_t i=0; i<nalong*nacross; i++) {
1824  if (*fptr < 0) *fptr = -32767;
1825  fptr++;
1826  }
1827 
1828  cout << "Writing chlor_a" << endl;
1829  varout = nc_output->addVar("chlor_a", ncFloat, dims);
1830  copyVarAtts( &var, &varout);
1831  varout.putVar(&out_l1c_2d[0][0]);
1832 
1833  cgal_release_tag(nnc_tag_chl);
1834  delete [] lat_chl;
1835  delete [] lon_chl;
1836  delete [] chl;
1837  }
1838 
1839  free2d_float(out_l1c_2d);
1840  free2d_float(latL1C);
1841  free2d_float(lonL1C);
1842 
1843  return 0;
1844 }
1845 
1846 int lonlat2rowcol( uint32_t nalong, uint32_t nacross, uint32_t ncm,
1847  float gridres, float *lonL1C, float *latL1C,
1848  float *lonCM, float *latCM, short *brow, short *bcol) {
1849 
1850  ofstream tempOut;
1851 
1852  // Determine L1C row/col of L2 CLD lon/lat
1853  // Ported from code developed by F.Patt
1854 
1855  float *fptrlon, *fptrlat;
1856 
1857  // Convert lon/lat to unit vectors (L1C)
1858  float ***gvec = allocate3d_float(nalong, nacross, 3);
1859  fptrlon = lonL1C;
1860  fptrlat = latL1C;
1861  for (size_t i=0; i<nalong; i++) {
1862  for (size_t j=0; j<nacross; j++) {
1863  gvec[i][j][0] = cos(*fptrlon/RADEG) * cos(*fptrlat/RADEG);
1864  gvec[i][j][1] = sin(*fptrlon/RADEG) * cos(*fptrlat/RADEG);
1865  gvec[i][j][2] = sin(*fptrlat/RADEG);
1866 
1867  fptrlon++;
1868  fptrlat++;
1869  }
1870  }
1871 
1872  // Convert lon/lat to unit vectors (L2 CLD)
1873  float **bvec = allocate2d_float(ncm, 3);
1874  fptrlon = lonCM;
1875  fptrlat = latCM;
1876  for (size_t i=0; i<ncm; i++) {
1877  bvec[i][0] = cos(*fptrlon/RADEG)*cos(*fptrlat/RADEG);
1878  bvec[i][1] = sin(*fptrlon/RADEG)*cos(*fptrlat/RADEG);
1879  bvec[i][2] = sin(*fptrlat/RADEG);
1880 
1881  fptrlon++;
1882  fptrlat++;
1883  }
1884 
1885  //tempOut.open ("gb.bin", ios::out | ios::trunc | ios::binary);
1886  //tempOut.write((char *) &gvec[0][0][0], sizeof(float)*nalong*nacross*3);
1887  //tempOut.write((char *) &bvec[0][0], sizeof(float)*ncm*3);
1888  //tempOut.close();
1889 
1890  // Compute normal vectors for L1C rows
1891  // (cross product of first and last vector in each row)
1892  float **gnvec = allocate2d_float(nalong, 3);
1893  float *gnvm = new float[nalong];
1894  float vecm[3];
1895  for (size_t i=0; i<nalong; i++) {
1896  gnvec[i][0] =
1897  gvec[i][nacross-1][1]*gvec[i][0][2] -
1898  gvec[i][nacross-1][2]*gvec[i][0][1];
1899 
1900  gnvec[i][1] =
1901  gvec[i][nacross-1][2]*gvec[i][0][0] -
1902  gvec[i][nacross-1][0]*gvec[i][0][2];
1903 
1904  gnvec[i][2] =
1905  gvec[i][nacross-1][0]*gvec[i][0][1] -
1906  gvec[i][nacross-1][1]*gvec[i][0][0];
1907 
1908  vecm[0] = gnvec[i][0];
1909  vecm[1] = gnvec[i][1];
1910  vecm[2] = gnvec[i][2];
1911 
1912  gnvm[i] = sqrt(vecm[0]*vecm[0] + vecm[1]*vecm[1] + vecm[2]*vecm[2]);
1913  }
1914  for (size_t i=0; i<nalong; i++) {
1915  gnvec[i][0] /= gnvm[i];
1916  gnvec[i][1] /= gnvm[i];
1917  gnvec[i][2] /= gnvm[i];
1918  }
1919  delete [] gnvm;
1920 
1921  //tempOut.open ("gnvec.bin", ios::out | ios::trunc | ios::binary);
1922  //tempOut.write((char *) &gnvec[0][0], sizeof(float)*nalong*3);
1923  //tempOut.close();
1924 
1925  // Compute dot products of L1B vectors with normal vectors
1926 
1927  gsl_matrix_float_view G =
1928  gsl_matrix_float_view_array(&gnvec[0][0], nalong, 3);
1929  gsl_matrix_float_view B =
1930  gsl_matrix_float_view_array(&bvec[0][0], ncm, 3);
1931 
1932  float **bdotgn = allocate2d_float(ncm, nalong);
1933  gsl_matrix_float_view bdotgn_mat =
1934  gsl_matrix_float_view_array(&bdotgn[0][0], ncm, nalong);
1935 
1936  //cout << "matrix multiplication" << endl;
1937  gsl_blas_sgemm(CblasNoTrans, CblasTrans, 1.0, &B.matrix, &G.matrix,
1938  0.0, &bdotgn_mat.matrix);
1939 
1940 
1941  // tempOut.open ("bd.bin", ios::out | ios::trunc | ios::binary);
1942  //tempOut.write((char *) &bdotgn[0][0], sizeof(float)*ncm*nalong);
1943  //tempOut.close();
1944 
1945  // Determine row and column for each pixel within grid
1946  // Compute normals to center columns
1947  uint32_t ic = nacross / 2;
1948  float **cnvec = allocate2d_float(nalong, 3);
1949  for (size_t i=0; i<nalong; i++) {
1950  cnvec[i][0] = gvec[i][ic][1]*gnvec[i][2] - gvec[i][ic][2]*gnvec[i][1];
1951  cnvec[i][1] = gvec[i][ic][2]*gnvec[i][0] - gvec[i][ic][0]*gnvec[i][2];
1952  cnvec[i][2] = gvec[i][ic][0]*gnvec[i][1] - gvec[i][ic][1]*gnvec[i][0];
1953  }
1954 
1955  // Compute grid row nadir resolution
1956  float *dcm = new float[nalong];
1957  for (size_t i=0; i<nalong; i++) {
1958  vecm[0] = gvec[i][ic+1][0] - gvec[i][ic][0];
1959  vecm[1] = gvec[i][ic+1][1] - gvec[i][ic][1];
1960  vecm[2] = gvec[i][ic+1][2] - gvec[i][ic][2];
1961 
1962  dcm[i] = sqrt(vecm[0]*vecm[0] + vecm[1]*vecm[1] + vecm[2]*vecm[2]);
1963  }
1964 
1965 
1966  float db = gridres/6311/2;
1967 
1968  // tempOut.open ("ibad.bin", ios::out | ios::trunc | ios::binary);
1969 
1970  for (size_t i=0; i<ncm; i++) {
1971  if (bdotgn[i][0] > db || bdotgn[i][nalong-1] < -db) {
1972  brow[i] = -1;
1973  bcol[i] = -1;
1974  // tempOut.write((char *) &i, sizeof(size_t));
1975 
1976  } else {
1977  brow[i] = -1;
1978  bcol[i] = -1;
1979  //if (i % 100000 == 0) cout << i << endl;
1980  for (size_t j=0; j<nalong; j++) {
1981  if (bdotgn[i][j] > -db && bdotgn[i][j] <= db) {
1982  brow[i] = j;
1983 
1984  float bdotcn = 0.0;
1985  for (size_t k=0; k<3; k++)
1986  bdotcn += (bvec[i][k] * cnvec[j][k]);
1987 
1988  bcol[i] = (short) round(ic + bdotcn/dcm[j]);
1989 
1990  if (bcol[i] < 0 || bcol[i] >= (short) nacross) {
1991  bcol[i] = -1;
1992  brow[i] = -1;
1993  }
1994  break;
1995  }
1996  }
1997  }
1998  }
1999 
2000  // tempOut.close();
2001 
2002  free3d_float(gvec);
2003  free2d_float(bvec);
2004  free2d_float(gnvec);
2005  free2d_float(cnvec);
2006  free2d_float(bdotgn);
2007 
2008  delete [] dcm;
2009 
2010  return 0;
2011 }
2012 
cgal_interp.h
clo_addOption
clo_option_t * clo_addOption(clo_optionList_t *list, const char *key, enum clo_dataType_t dataType, const char *defaultVal, const char *desc)
Definition: clo.c:684
cgal_nnc
int cgal_nnc(int nxy, float *x, float *y, int nxyq, float *xq, float *yq)
Definition: cgal_interp.cpp:27
allocate3d.h
Utility functions for allocating and freeing three-dimensional arrays of various types.
j
int j
Definition: decode_rs.h:73
clo_getString
char * clo_getString(clo_optionList_t *list, const char *key)
Definition: clo.c:1357
L
#define L(lambda, T)
Definition: PreprocessP.h:185
clo_readArgs
void clo_readArgs(clo_optionList_t *list, int argc, char *argv[])
Definition: clo.c:2103
free2d_schar
void free2d_schar(signed char **p)
Free a two-dimensional array created by allocate2d_schar.
Definition: allocate2d.c:74
mapgen.levels
list levels
Definition: mapgen.py:182
expandEnvVar
int expandEnvVar(std::string *sValue)
Definition: ancgen.cpp:222
lonlat2rowcol
int lonlat2rowcol(uint32_t nalong, uint32_t nacross, uint32_t ncm, float gridres, float *lonL1C, float *latL1C, float *lonCM, float *latCM, short *brow, short *bcol)
Definition: ancgen.cpp:1846
N
int N
Definition: Usds.c:60
NULL
#define NULL
Definition: decode_rs.h:63
clo_option_t::key
char * key
Definition: clo.h:107
syear
int syear
Definition: l1_czcs_hdf.c:15
smsec
int32 smsec
Definition: l1_czcs_hdf.c:16
string
@ string
Definition: GEO_read_param_file.c:64
clo_isOptionSet
int clo_isOptionSet(clo_option_t *option)
Definition: clo.c:2257
readCLD
int readCLD(NcGroup ncGrp[], const char *fldname, T *array, uint32_t npixels, uint32_t nlines[])
Definition: ancgen.cpp:169
free2d_short
void free2d_short(short **p)
Free a two-dimensional array created by allocate2d_short.
Definition: allocate2d.c:96
nobs
int32_t nobs
Definition: atrem_cor.h:93
clo_setVersion2
void clo_setVersion2(const char *programName, const char *versionStr)
Definition: clo.c:464
sday
int sday
Definition: l1_czcs_hdf.c:15
RADEG
const double RADEG
Definition: ancgen.cpp:22
clo_optionList_t
Definition: clo.h:126
NCACHE
#define NCACHE
Definition: ancgen.cpp:20
tempOut
ofstream tempOut
Definition: l1agen_oci.cpp:181
clo_createList
clo_optionList_t * clo_createList()
Definition: clo.c:532
list
list(APPEND LIBS ${NETCDF_LIBRARIES}) find_package(GSL REQUIRED) include_directories($
Definition: CMakeLists.txt:8
allocate2d_schar
signed char ** allocate2d_schar(size_t h, size_t w)
Allocate a two-dimensional array of type signed char of a given size.
Definition: allocate2d.c:57
clo_option_t
Definition: clo.h:106
clo_getOptionString
char * clo_getOptionString(clo_option_t *option)
Definition: clo.c:1050
free2d_float
void free2d_float(float **p)
Free a two-dimensional array created by allocate2d_float.
Definition: allocate2d.c:140
hktgen_pace.att
int att
Definition: hktgen_pace.py:28
main
int main(int argc, char *argv[])
Definition: ancgen.cpp:238
clo_printUsage
void clo_printUsage(clo_optionList_t *list)
Definition: clo.c:1988
oci_hdr_strip.argc
argc
Definition: oci_hdr_strip.py:7
allocate3d_float
float *** allocate3d_float(size_t nz, size_t ny, size_t nx)
Allocate a three-dimensional array of type float of a given size.
Definition: allocate3d.c:77
cgal_interp2
int cgal_interp2(int nxy, float *x, float *y, float *v, int nxyq, float *vq, int nnc_tag)
Definition: cgal_interp.cpp:74
CLO_TYPE_IFILE
@ CLO_TYPE_IFILE
Definition: clo.h:87
clo_getOption
clo_option_t * clo_getOption(clo_optionList_t *list, int i)
Definition: clo.c:908
M
int M[]
Definition: Usds.c:107
isodate2ydmsec
void isodate2ydmsec(char *date, int32_t *year, int32_t *day, int32_t *msec)
Definition: date2ydmsec.c:20
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
allocate2d.h
Utility functions for allocating and freeing two-dimensional arrays of various types.
clo_getNumOptions
int clo_getNumOptions(clo_optionList_t *list)
Definition: clo.c:1017
parse_file_name
void parse_file_name(const char *inpath, char *outpath)
Definition: parse_file_name.c:23
free3d_float
void free3d_float(float ***p)
Free a three-dimensional array created by allocate3d_float.
Definition: allocate3d.c:103
basename
#define basename(s)
Definition: l0chunk_modis.c:29
B
const float B
Definition: calc_par.c:101
accum_wm
int accum_wm(uint32_t nwm, short *brow, short *bcol, float *wlval, float **binval, short **nobs)
Definition: ancgen.cpp:133
now
u5 which has been done in the LOCALGRANULEID metadata should have an extension NRT It is requested to identify the NRT production Changes from v6 which may affect scientific the sector rotation may actually occur during one of the scans earlier than the one where it is first reported As a the b1 values are about the LOCALGRANULEID metadata should have an extension NRT It is requested to identify the NRT to fill pixels affected by dead subframes with a special value Output the metadata of noisy and dead subframe Dead Subframe EV and Detector Quality Flag2 Removed the function call of Fill_Dead_Detector_SI to stop interpolating SI values for dead but also for all downstream products for science test only Changes from v5 which will affect scientific to conform to MODIS requirements Removed the Mixed option from the ScanType in the code because the L1A Scan Type is never Mixed Changed for ANSI C compliance and comments to better document the fact that when the HDF_EOS metadata is stricly the and products are off by and in the track respectively Corrected some misspelling of RCS swir_oob_sending_detector to the Reflective LUTs to enable the SWIR OOB correction detector so that if any of the sending detectors becomes noisy or non near by good detectors from the same sending band can be specified as the substitute in the new look up table Code change for adding an additional dimension of mirror side to the Band_21_b1 LUT to separate the coefficient of the two mirror sides for just like other thermal emissive so that the L1B code can calibrate Band scan to scan with mirror side dependency which leads better calibration result Changes which do not affect scientific when the EV data are not provided in this Crosstalk Correction will not be performed to the Band calibration data Changes which do not affect scientific and BB_500m in L1A Logic was added to turn off the or to spatial aggregation processes and the EV_250m_Aggr1km_RefSB and EV_500m_Aggr1km_RefSB fields were set to fill values when SDSs EV_250m and EV_500m are absent in L1A file Logic was added to skip the processing and turn off the output of the L1B QKM and HKM EV data when EV_250m and EV_500m are absent from L1A In this the new process avoids accessing and reading the and L1A EV skips and writing to the L1B and EV omits reading and subsampling SDSs from geolocation file and writing them to the L1B and omits writing metadata to L1B and EV and skips closing the L1A and L1B EV and SDSs Logic was added to turn off the L1B OBC output when the high resolution OBC SDSs are absent from L1A This is accomplished by skipping the openning the writing of metadata and the closing of the L1B OBC hdf which is Bit in the scan by scan bit QA has been changed Until now
Definition: HISTORY.txt:361
unix2isodate
char * unix2isodate(double dtime, char zone)
Definition: unix2isodate.c:10
cgal_release_tag
int cgal_release_tag(int nnc_tag)
Definition: cgal_interp.cpp:111
s
data_t s[NROOTS]
Definition: decode_rs.h:75
allocate2d_short
short ** allocate2d_short(size_t h, size_t w)
Allocate a two-dimensional array of type short of a given size.
Definition: allocate2d.c:79
allocate2d_float
float ** allocate2d_float(size_t h, size_t w)
Allocate a two-dimensional array of type float of a given size.
Definition: allocate2d.c:123
PI
const double PI
Definition: ancgen.cpp:21
copyVarAtts
int copyVarAtts(NcVar *varin, NcVar *varout, string override[], string overridetype[])
Definition: copyvaratts.cpp:10
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")
accum_frac
int accum_frac(uint32_t totlines, uint32_t npixels, uint32_t nalong, uint32_t nacross, short *brow, short *bcol, float **cldprod, int8_t **cld_phase, float **binval_ic, short **nobs_ic, float **binval_wc, short **nobs_wc)
Definition: ancgen.cpp:85
k
int k
Definition: decode_rs.h:73
accum
int accum(uint32_t totlines, uint32_t npixels, uint32_t nalong, uint32_t nacross, short *brow, short *bcol, float **cldprod, int8_t **cld_phase, float **binval_ic, short **nobs_ic, float **binval_wc, short **nobs_wc)
Definition: ancgen.cpp:38
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