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
l1c_latlongrid.cpp
Go to the documentation of this file.
1 //**********************************************************
2 // library for searching row/col from lat/lon of L1C file
3 // Originally created by Martin Montes on 11/4/2022
4 // last version 3/13/2023
5 //**********************************************************
6 #include <allocate2d.h>
7 #include <allocate3d.h>
8 #include "allocate4d.h"
9 #include <iostream>
10 #include <chrono>
11 #include <sys/stat.h>
12 #include <genutils.h>
13 #include <iostream>
14 #include "l1c_latlongrid.h"
15 #include <netcdf>
16 #include <GeographicLib/Geodesic.hpp>
17 #include <GeographicLib/Constants.hpp>
18 #include <GeographicLib/Rhumb.hpp>
19 
20 #include <fstream>
21 #include <string>
22 
23 using namespace std;
24 using namespace netCDF;
25 using namespace netCDF::exceptions;
26 using namespace GeographicLib;
27 
28 bin_str::bin_str() {
29  // global attrs
30  version = "";
31  history = "";
32  pversion = "";
33  doi = "";
34  sensor = -1;
35 
36  // dimensions
37  nviews = -1;
38  nbands = -1;
39  num_gridlines = -1;
40  nbinx = -1;
41  num_pixels = -1;
42  nscans = -1;
43  verbose = 0;
44  // geo
45  gdshift = -1;
46  lat_gd = nullptr;
47  lon_gd = nullptr;
48  time_l1b = nullptr;
49  time_gd = nullptr;
50  alt = nullptr;
51 
52  // binned vars
53 
54  nrec_2D = nullptr; // row/col
55  nrec_3D = nullptr;
56  nrec_3D_view = nullptr;
57  nrec_4D_band = nullptr; // row/col/view
58  alt_mean = nullptr;
59  alt_rmse = nullptr;
60  alt_2D = nullptr;
61  alt_diff2 = nullptr; // row/col
62  i_diff2 = nullptr;
63  sca_3D = nullptr; // row/col/view
64  rot_angle = nullptr;
65  QC_bitwise_4D = nullptr; // row/col/view/bands
66  QC_4D = nullptr; // row/col/view
67  I_4D = nullptr; // row/col/view/bands
68  I_noise_4D = nullptr; // #pixels
69 
70  // OCIS line by line
71  obs_per_view = nullptr;
72  QC = nullptr;
73  I = nullptr;
74  I_noise = nullptr;
75 
76  l1file = nullptr;
77  full_l1cgrid = "";
78  l1cgrid = "";
79 
80  fillval1 = BAD_FLT;
81  fillval2 = float(BAD_FLT);
82  fillval3 = double(BAD_FLT);
83 
84  inpix = 0;
85  outpix = 0;
86  badgeo = 0;
87 
88  tini_l1c = -1;
89  tend_l1c = -1;
90  tini_l1b = -1;
91  tend_l1b = -1;
92  date_mid_grid = "";
93 
94  bintype = 0;
95  outlist[0] = '\0';
96  l1c_anc[0] = '\0';
97  cloudem_flag = 0;
98  cloud_type = 0; // 0 water, 1 ice
99  dem_flag=0;//dem 0 geoid height or L1C grid height between ellipsoid and geoid wgs84 in m, 1 orthometric height or dem
100 }
101 
102 bin_str::~bin_str() {
103 }
104 
105 int bin_str::close_bin(bin_str *binl1c) {
106  if (binl1c->time_l1b != nullptr)
107  free3d_double(binl1c->time_l1b);
108  binl1c->time_l1b = nullptr;
109  if (binl1c->time_gd != nullptr)
110  free(binl1c->time_gd);
111  binl1c->time_gd = nullptr;
112  if (binl1c->alt != nullptr)
113  free2d_float(binl1c->alt);
114  binl1c->alt = nullptr;
115  if (binl1c->lat_gd != nullptr)
116  free2d_float(binl1c->lat_gd);
117  binl1c->lat_gd = nullptr;
118  if (binl1c->lon_gd != nullptr)
119  free2d_float(binl1c->lon_gd);
120  binl1c->lon_gd = nullptr;
121  if (binl1c->alt_rmse != nullptr)
122  free2d_float(binl1c->alt_rmse);
123  binl1c->alt_rmse = nullptr;
124  if (binl1c->alt_mean != nullptr)
125  free2d_float(binl1c->alt_mean);
126  binl1c->alt_mean = nullptr;
127  if (binl1c->alt_2D != nullptr)
128  free2d_short(binl1c->alt_2D);
129  binl1c->alt_2D = nullptr;
130  if (binl1c->alt_diff2 != nullptr)
131  free2d_short(binl1c->alt_diff2);
132  binl1c->alt_diff2 = nullptr;
133  if (binl1c->i_diff2 != nullptr)
134  free4d_float(binl1c->i_diff2);
135  binl1c->i_diff2 = nullptr;
136  if (binl1c->suna_3D != nullptr)
137  free3d_float(binl1c->suna_3D);
138  binl1c->suna_3D = nullptr;
139  if (binl1c->sunz_3D != nullptr)
140  free3d_float(binl1c->sunz_3D);
141  binl1c->sunz_3D = nullptr;
142  if (binl1c->sena_3D != nullptr)
143  free3d_float(binl1c->sena_3D);
144  binl1c->sena_3D = nullptr;
145  if (binl1c->senz_3D != nullptr)
146  free3d_float(binl1c->senz_3D);
147  binl1c->senz_3D = nullptr;
148  if (binl1c->sca_3D != nullptr)
149  free3d_float(binl1c->sca_3D);
150  binl1c->sca_3D = nullptr;
151  if (binl1c->rot_angle != nullptr)
152  free3d_float(binl1c->rot_angle);
153  binl1c->rot_angle = nullptr;
154  if (binl1c->nrec_2D != nullptr)
155  free2d_short(binl1c->nrec_2D);
156  binl1c->nrec_2D = nullptr;
157  if (binl1c->nrec_3D != nullptr)
158  free3d_short(binl1c->nrec_3D);
159  binl1c->nrec_3D = nullptr;
160  if (binl1c->nrec_3D_view != nullptr)
161  free3d_short(binl1c->nrec_3D_view);
162  binl1c->nrec_3D_view = nullptr;
163  if (binl1c->nrec_4D_band != nullptr)
164  free4d_float(binl1c->nrec_4D_band);
165  binl1c->nrec_4D_band = nullptr;
166  if (binl1c->I_4D != nullptr)
167  free4d_float(binl1c->I_4D);
168  binl1c->I_4D = nullptr;
169  if (binl1c->I_noise_4D != nullptr)
170  free4d_float(binl1c->I_noise_4D);
171  binl1c->I_noise_4D = nullptr;
172 
173  return 0;
174 }
175 
176 int bin_str::alloc_bin(bin_str *binl1c) {
177  if (binl1c->verbose) {
178  cout << "allocating/init high order dimensional arrays for L1C binning----------------------------"
179  << endl;
180  cout << "nviews.." << binl1c->nviews << "nbands.." << binl1c->nbands << "#gridlines.."
181  << binl1c->num_gridlines << "nbinx.." << binl1c->nbinx << endl;
182  }
183 
184  binl1c->nrec_2D = allocate2d_short(binl1c->num_gridlines, binl1c->nbinx);
185  binl1c->alt = allocate2d_float(binl1c->num_gridlines, binl1c->nbinx);
186  binl1c->alt_rmse = allocate2d_float(binl1c->num_gridlines, binl1c->nbinx);
187  binl1c->alt_mean = allocate2d_float(binl1c->num_gridlines, binl1c->nbinx);
188  binl1c->alt_2D = allocate2d_short(binl1c->num_gridlines, binl1c->nbinx);
189  binl1c->alt_diff2 = allocate2d_short(binl1c->num_gridlines, binl1c->nbinx);
190  if (binl1c->bintype == 0) {
191  binl1c->lat_gd = allocate2d_float(binl1c->num_gridlines, binl1c->nbinx);
192  binl1c->lon_gd = allocate2d_float(binl1c->num_gridlines, binl1c->nbinx);
193  }
194 
195  binl1c->time_gd = (double *)calloc(binl1c->num_gridlines, sizeof(double));
196  binl1c->time_l1b = allocate3d_double(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
197  binl1c->suna_3D = allocate3d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
198  binl1c->sunz_3D = allocate3d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
199  binl1c->sena_3D = allocate3d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
200  binl1c->senz_3D = allocate3d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
201  binl1c->sca_3D = allocate3d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
202  binl1c->rot_angle = allocate3d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
203  binl1c->nrec_3D = allocate3d_short(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
204  binl1c->nrec_3D_view = allocate3d_short(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
205 
206  binl1c->nrec_4D_band =
207  allocate4d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews, binl1c->nbands);
208  binl1c->I_4D = allocate4d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews, binl1c->nbands);
209  binl1c->I_noise_4D =
210  allocate4d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews, binl1c->nbands);
211  binl1c->i_diff2 = allocate4d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews, binl1c->nbands);
212 
213  // init
214  for (int i = 0; i < binl1c->num_gridlines; i++) {
215  binl1c->time_gd[i] = 0;
216  for (int j = 0; j < binl1c->nbinx; j++) {
217  binl1c->nrec_2D[i][j] = 0;
218  binl1c->alt[i][j] = 0.;
219  binl1c->alt_rmse[i][j] = 0.;
220  binl1c->alt_mean[i][j] = 0.;
221  binl1c->alt_2D[i][j] = 0;
222  binl1c->alt_diff2[i][j] = 0;
223  if (binl1c->bintype == 0) {
224  binl1c->lat_gd[i][j] = 0.;
225  binl1c->lon_gd[i][j] = 0.;
226  }
227 
228  for (int v = 0; v < binl1c->nviews; v++) {
229  binl1c->time_l1b[i][j][v] = 0.;
230  binl1c->suna_3D[i][j][v] = 0.;
231  binl1c->sunz_3D[i][j][v] = 0.;
232  binl1c->sena_3D[i][j][v] = 0.;
233  binl1c->senz_3D[i][j][v] = 0.;
234  binl1c->sca_3D[i][j][v] = 0.;
235  binl1c->rot_angle[i][j][v] = 0.;
236  binl1c->nrec_3D[i][j][v] = 0;
237  binl1c->nrec_3D_view[i][j][v] = 0;
238  for (int sb = 0; sb < binl1c->nbands; sb++) {
239  binl1c->nrec_4D_band[i][j][v][sb] = 0.;
240  binl1c->I_4D[i][j][v][sb] = 0.;
241  binl1c->I_noise_4D[i][j][v][sb] = 0.;
242  binl1c->i_diff2[i][j][v][sb] = 0;
243  }
244  }
245  }
246  }
247 
248  return 0;
249 }
250 
251 int rmse_l1c_alt(filehandle *l1file, bin_str *binl1c, l1str *l1rec, short **gdindex) {
252  short gd_row = -1, gd_col = -1;
253  float fill_tilt = BAD_FLT;
254  short fill_height = BAD_FLT;
255  int npix = l1file->npix;
256 
257  for (int pix = 0; pix < npix; pix++)
258  {
259  gd_row = gdindex[pix][0] - 1;
260  gd_col = gdindex[pix][1] - 1;
261 
262  if (l1rec->height[pix] != fill_height && l1rec->tilt != fill_tilt && gd_row >= 0 && gd_col >= 0 &&
263  l1rec->lat[pix] != binl1c->fillval2 && l1rec->lon[pix] != binl1c->fillval2)
264  {
265 
266  binl1c->alt_diff2[gd_row][gd_col] =
267  binl1c->alt_diff2[gd_row][gd_col] + (l1rec->height[pix] - binl1c->alt[gd_row][gd_col]) *
268  (l1rec->height[pix] - binl1c->alt[gd_row][gd_col]);
269  }
270  }
271 
272  return 0;
273 }
274 
275 int meta_l1c_altvar(bin_str *binl1c, NcFile *nc_output) {
276  float term;
277 
278  if (binl1c->verbose)
279  cout << "adding alt_rmse to L1C file..." << endl;
280 
281  NcGroup geo_grp = nc_output->getGroup("geolocation_data");
282  NcVar v1 = geo_grp.getVar("height_stdev");
283 
284  for (int i = 0; i < binl1c->num_gridlines; i++) {
285  for (int j = 0; j < binl1c->nbinx; j++) {
286  term = binl1c->alt_diff2[i][j];
287  if (term >= 0 && binl1c->nrec_2D[i][j] > 0) {
288  binl1c->alt_rmse[i][j] = sqrt(binl1c->alt_diff2[i][j] / binl1c->nrec_2D[i][j]);
289  } else
290  binl1c->alt_rmse[i][j] = binl1c->fillval2;
291  }
292  }
293  v1.putVar(&binl1c->alt_rmse[0][0]);
294 
295 
296  //i_stdev
297  float ****temp2 = allocate4d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews, binl1c->nbands);
298 
299  for (int i = 0; i < binl1c->num_gridlines; i++) {
300  for (int j = 0; j < binl1c->nbinx; j++) {
301  for (int v = 0; v < binl1c->nviews; v++) {
302  for (int sb = 0; sb < binl1c->nbands; sb++) {
303  if (binl1c->nrec_4D_band[i][j][v][sb] > 0) {
304  binl1c->i_diff2[i][j][v][sb]-=0.01;
305  temp2[i][j][v][sb] = sqrt(binl1c->i_diff2[i][j][v][sb]/binl1c->nrec_4D_band[i][j][v][sb]);
306  } else {
307  temp2[i][j][v][sb] = binl1c->fillval2;
308  }
309  }
310  }
311  }
312  }
313 
314  NcGroup od_grp = nc_output->getGroup("observation_data");
315  v1 = od_grp.getVar("i_stdev");
316  if (binl1c->bintype == 0) {
317  v1 = od_grp.getVar("i_stdev");
318  v1.putVar(&temp2[0][0][0][0]);
319  }
320 
321  free4d_float(temp2);
322 
323  return 0;
324 }
325 
326 int meta_l1c_bin(filehandle *l1file, bin_str *binl1c, NcFile *nc_output) {
327  if (binl1c->verbose)
328  cout << "adding final binned vars to ...." << binl1c->full_l1cgrid << endl;
329 
330  NcDim ydim = nc_output->getDim("bins_along_track");
331  NcDim xdim = nc_output->getDim("bins_across_track");
332  NcDim vdim = nc_output->getDim("number_of_views");
333 
334  std::vector<NcDim> dimvec3;
335  dimvec3.push_back(ydim);
336  dimvec3.push_back(xdim);
337  dimvec3.push_back(vdim);
338 
339  NcGroup geo_grp = nc_output->getGroup("geolocation_data");
340 
341  // time_offsets
342 
343  double ***time_offsets = nullptr, toff, toff_fill, toff_min, toff_max;
344 
345 
346  time_offsets = allocate3d_double(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
347 
348  for (int i = 0; i < binl1c->num_gridlines; i++) {
349  for (int j = 0; j < binl1c->nbinx; j++) {
350  for (int v = 0; v < binl1c->nviews; v++) {
351  time_offsets[i][j][v] = 0;
352  }
353  }
354  }
355 
356  NcGroup ba_grp = nc_output->getGroup("bin_attributes");
357  NcVar v1 = ba_grp.getVar("nadir_view_time");
358  v1.getVar(&binl1c->time_gd[0]);
359 
360  v1 = ba_grp.getVar("view_time_offsets");
361  NcVarAtt a1 = v1.getAtt("_FillValue"); // root group
362  a1.getValues(&toff_fill);
363  a1 = v1.getAtt("valid_min"); // root group
364  a1.getValues(&toff_min);
365  a1 = v1.getAtt("valid_max"); // root group
366  a1.getValues(&toff_max);
367 
368  // sensor azimuth---
369  for (int i = 0; i < binl1c->num_gridlines; i++) {
370  for (int j = 0; j < binl1c->nbinx; j++) {
371  for (int v = 0; v < binl1c->nviews; v++) {
372  if (binl1c->nrec_3D[i][j][v] > 0) {
373  toff = binl1c->time_gd[i] - (binl1c->time_l1b[i][j][v] / binl1c->nrec_3D[i][j][v]);
374 
375  if (toff >= toff_min && toff <= toff_max) {
376  time_offsets[i][j][v] = toff;
377  } else
378  time_offsets[i][j][v] = toff_fill;
379  } else {
380  time_offsets[i][j][v] = toff_fill;
381  }
382  }
383  }
384  }
385 
386  v1 = ba_grp.getVar("view_time_offsets");
387  v1.putVar(&time_offsets[0][0][0]);
388 
389  free3d_double(time_offsets);
390 
391  v1 = geo_grp.getVar("sensor_azimuth_angle");
392 
393  short ***temp;
394  float angleScale = 100.0;
395 
396  temp = allocate3d_short(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews);
397 
398  // sensor azimuth---
399  for (int i = 0; i < binl1c->num_gridlines; i++) {
400  for (int j = 0; j < binl1c->nbinx; j++) {
401  for (int v = 0; v < binl1c->nviews; v++) {
402  if (binl1c->nrec_3D[i][j][v] > 0) {
403  temp[i][j][v] = (binl1c->sena_3D[i][j][v] / binl1c->nrec_3D[i][j][v]) * angleScale;
404  } else {
405  temp[i][j][v] = binl1c->fillval1;
406  }
407  }
408  }
409  }
410 
411  v1.putVar(&temp[0][0][0]);
412 
413  // sensor zenith---
414  for (int i = 0; i < binl1c->num_gridlines; i++) {
415  for (int j = 0; j < binl1c->nbinx; j++) {
416  for (int v = 0; v < binl1c->nviews; v++) {
417  if (binl1c->nrec_3D[i][j][v] > 0) {
418  temp[i][j][v] = (binl1c->senz_3D[i][j][v] / binl1c->nrec_3D[i][j][v] )* angleScale;
419  } else {
420  temp[i][j][v] = binl1c->fillval1;
421  }
422  }
423  }
424  }
425 
426  v1 = geo_grp.getVar("sensor_zenith_angle");
427  v1.putVar(&temp[0][0][0]);
428 
429  // solar azimuth---
430  for (int i = 0; i < binl1c->num_gridlines; i++) {
431  for (int j = 0; j < binl1c->nbinx; j++) {
432  for (int v = 0; v < binl1c->nviews; v++) {
433  if (binl1c->nrec_3D[i][j][v] > 0) {
434  temp[i][j][v] = (binl1c->suna_3D[i][j][v] / binl1c->nrec_3D[i][j][v] )* angleScale;
435  } else {
436  temp[i][j][v] = binl1c->fillval1;
437  }
438  }
439  }
440  }
441 
442  v1 = geo_grp.getVar("solar_azimuth_angle");
443  v1.putVar(&temp[0][0][0]);
444 
445  // sensor zenith---
446  for (int i = 0; i < binl1c->num_gridlines; i++) {
447  for (int j = 0; j < binl1c->nbinx; j++) {
448  for (int v = 0; v < binl1c->nviews; v++) {
449  if (binl1c->nrec_3D[i][j][v] > 0) {
450  temp[i][j][v] = (binl1c->sunz_3D[i][j][v] / binl1c->nrec_3D[i][j][v] )* angleScale;
451  } else {
452  temp[i][j][v] = binl1c->fillval1;
453  }
454  }
455  }
456  }
457 
458  v1 = geo_grp.getVar("solar_zenith_angle");
459  v1.putVar(&temp[0][0][0]);
460 
461  // scattering angle---
462  for (int i = 0; i < binl1c->num_gridlines; i++) {
463  for (int j = 0; j < binl1c->nbinx; j++) {
464  for (int v = 0; v < binl1c->nviews; v++) {
465  if (binl1c->nrec_3D[i][j][v] > 0) {
466  temp[i][j][v] = (binl1c->sca_3D[i][j][v] / binl1c->nrec_3D[i][j][v]) * angleScale;
467  } else {
468  temp[i][j][v] = binl1c->fillval1;
469  }
470  }
471  }
472  }
473 
474  v1 = geo_grp.getVar("scattering_angle");
475  v1.putVar(&temp[0][0][0]);
476 
477  // rotation angle---
478  // for (int i = 0; i < binl1c->num_gridlines; i++) {
479  // for (int j = 0; j < binl1c->nbinx; j++) {
480  // for (int v = 0; v < binl1c->nviews; v++) {
481  // if (binl1c->nrec_3D[i][j][v] > 0) {
482  // temp[i][j][v] = (binl1c->rot_angle[i][j][v] / binl1c->nrec_3D[i][j][v]) * angleScale;
483  // } else {
484  // temp[i][j][v] = binl1c->fillval1;
485  // }
486  // }
487  // }
488  // }
489 
490  // v1 = geo_grp.getVar("rotation_angle");
491  // v1.putVar(&temp[0][0][0]);
492  free3d_short(temp);
493 
494  NcGroup od_grp = nc_output->getGroup("observation_data");
495  v1 = od_grp.getVar("number_of_observations"); // ONLY 1 BAND!!!!!! NO QUALITY CONTROL,constrained no fillvalues,
496  // and withing min/max Lt
497 
498  if (binl1c->bintype == 0) {
499  v1.putVar(&binl1c->nrec_3D[0][0][0]);
500  }
501 
502  // I/reflectance
503  float ****temp2 = allocate4d_float(binl1c->num_gridlines, binl1c->nbinx, binl1c->nviews, binl1c->nbands);
504  for (int i = 0; i < binl1c->num_gridlines; i++) {
505  for (int j = 0; j < binl1c->nbinx; j++) {
506  for (int v = 0; v < binl1c->nviews; v++) {
507  for (int sb = 0; sb < binl1c->nbands; sb++) {
508  if (binl1c->nrec_4D_band[i][j][v][sb] > 0) {
509  temp2[i][j][v][sb] = binl1c->I_4D[i][j][v][sb] / binl1c->nrec_4D_band[i][j][v][sb];
510 // cout<<"v "<<v+1<<"sb "<<sb+1<<"row "<<i+1<<"col.."<<j+1<<"MEAN radiance.. OCI = "<<temp2[i][j][v][sb]<<endl;
511  } else {
512  temp2[i][j][v][sb] = binl1c->fillval2;
513  }
514  }
515  }
516  }
517  }
518 
519  if (binl1c->bintype == 0) {
520  v1 = od_grp.getVar("i");
521  v1.putVar(&temp2[0][0][0][0]);
522  }
523 
524  free4d_float(temp2);
525 
526  return 0;
527 }
528 
529 int meta_l1c_full(filehandle *l1file, bin_str *binl1c, const char *l1c_grid, NcFile *nc_output) {
530  string senstr, titlestr, prodstr, ifile_str;
531  int NVIEWS, NBANDS; // NBANDS_POL;
532  int32_t xbins; //,ybins;
533  char *ifile_char = l1file->name;
534  file_format format = getFormat(ifile_char);
535  int16_t syear, smon, sday, syear2, smon2, sday2;
536  double secs, secs2;
537  double tgran_ini, tgran_ini_sec, tgran_end, tgran_end_sec;
538  string gfull;
539  int16_t gtime;
540  int16_t y_ini, mo_ini, d_ini, h_ini, mi_ini, y_end, mo_end, d_end, h_end, mi_end;
541  double sec_ini, sec_end;
542  int32_t gransize = 5;
543  int logoff = -1;
544 
545  if (binl1c->outlist[0] == '\0')
546  strcpy(binl1c->outlist, "l1c.tmp");
547  string outxt(binl1c->outlist);
548 
549  std::ofstream outf;
550 
551  outf.open(outxt, std::ofstream::out | std::ofstream::trunc);
552 
553  if (outf) {
554  if (binl1c->verbose)
555  cout << "writing L1C granules to outfile..." << outxt << endl;
556  } else {
557  std::cerr << "output file.." << outxt << " could not be opened for writing!\n";
558  return 1;
559  }
560 
561  if (format.type == FT_HKT || format.type == FT_L1C)
562  format.type = FT_OCIS; // if HKT then FT_OCIS----
563 
564  if (format.type == FT_SPEXONE) {
565  senstr = "SPEXONE";
566  titlestr = "PACE SPEXone Level-1C Data";
567  xbins = 29;
568  NVIEWS = 5;
569  NBANDS = 400;
570  // NBANDS_POL=50;
571  } else if (format.type == FT_HARP2) {
572  senstr = "HARP2";
573  titlestr = "PACE HARP2 Level-1C Data";
574  xbins = 457;
575  NVIEWS = 90;
576  NBANDS = 1;
577  // NBANDS_POL=1;
578  } else if (format.type == FT_OCIS) {
579  senstr = "OCI";
580  titlestr = "PACE OCI Level-1C Data";
581  xbins = 519;
582  NVIEWS = 2;
583  NBANDS = 239; // 249 originally no polarization bands
584  } else if (format.type == FT_OCIL1B) {
585  senstr = "OCI";
586  titlestr = "PACE OCI Level-1C Data";
587  xbins = 519;
588  NVIEWS = 2;
589  NBANDS = 286;
590  } else // OCI
591  {
592  senstr = "OCI";
593  titlestr = "PACE OCI Level-1C Data";
594  xbins = 519;
595  NVIEWS = 2;
596  NBANDS = 286;
597  }
598 
599  // l1c grid
600  string l1c_str = l1c_grid;
601 
602  NcFile *nc_l1cgrid;
603  try {
604  nc_l1cgrid = new NcFile(l1c_grid, NcFile::read);
605  } catch (NcException &e) {
606  e.what();
607  cerr << "l1cgen l1c_pflag= 8:: Failure reading L1C grid: " + l1c_str << endl;
608  exit(1);
609  }
610 
611  // l1c_grid--------
612  NcDim yd = nc_l1cgrid->getDim("bins_along_track");
613  int32_t num_gridlines = yd.getSize();
614  binl1c->nviews = NVIEWS;
615  binl1c->nbands = NBANDS;
616  binl1c->num_gridlines = num_gridlines;
617  binl1c->nbinx = xbins;
618 
619  // alloc binned vars in binl1c class
620  binl1c->alloc_bin(binl1c);
621 
622  NcDim vdim = nc_l1cgrid->getDim("number_of_views");
623 
624  // NcDim idim=nc_l1cgrid->getDim("intensity_bands_per_view");
625 
626  meta_l1c_global(ifile_char, binl1c,num_gridlines, nc_output);
627 
628  NcDim idim = nc_output->getDim("intensity_bands_per_view");
629 
630  // adding fill values of binned variables----
631  // sena,senz,sola,solz,sca, I
632  NcGroup geo_grp = nc_output->getGroup("geolocation_data");
633  NcVar v1 = geo_grp.getVar("sensor_azimuth_angle");
634  short value;
635  NcVarAtt a1 = v1.getAtt("_FillValue"); // root group
636  a1.getValues(&value);
637 
638  binl1c->fillval1 = value; //-32767s
639 
640  NcGroup od_grp = nc_output->getGroup("observation_data");
641  v1 = od_grp.getVar("i");
642  float value2;
643  a1 = v1.getAtt("_FillValue"); // root group
644  a1.getValues(&value2);
645 
646  binl1c->fillval2 = value2;
647 
648  // global attributes
649  if(!binl1c->pversion.empty())
650  nc_output->putAtt("processing_version", binl1c->pversion);
651  if(!binl1c->doi.empty()) {
652  nc_output->putAtt("identifier_product_doi", binl1c->doi);
653  nc_output->putAtt("identifier_product_doi_authority", "http://dx.doi.org");
654  }
655 
656  nc_output->putAtt("history", binl1c->history);
657  string name;
658  nc_output->putAtt("product_name", binl1c->full_l1cgrid);
659  NcGroupAtt i1 = nc_l1cgrid->getAtt("startDirection"); // NcGroupAtt is a global attr!!
660  i1.getValues(name);
661  nc_output->putAtt("startDirection", name);
662  i1 = nc_l1cgrid->getAtt("endDirection");
663  i1.getValues(name);
664  nc_output->putAtt("endDirection", name);
665  i1 = nc_l1cgrid->getAtt("time_coverage_start");
666  i1.getValues(name);
667  nc_output->putAtt("time_coverage_start", name);
668  tgran_ini = isodate2unix(name.c_str());
669  i1 = nc_l1cgrid->getAtt("time_coverage_end");
670  i1.getValues(name);
671  tgran_end = isodate2unix(name.c_str());
672  nc_output->putAtt("time_coverage_end", name);
673  i1 = nc_l1cgrid->getAtt("sun_earth_distance");
674  i1.getValues(&value2);
675  nc_output->putAtt("sun_earth_distance", ncFloat, value2);
676 
677  unix2ymds(tgran_ini, &syear, &smon, &sday, &secs);
678  unix2ymds(tgran_end, &syear2, &smon2, &sday2, &secs2);
679  tgran_ini_sec = ymds2unix(syear, smon, sday, secs);
680  tgran_end_sec = ymds2unix(syear2, smon2, sday2, secs2);
681 
682  unix2ymdhms(tgran_ini_sec, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
683  unix2ymdhms(tgran_end_sec, &y_end, &mo_end, &d_end, &h_end, &mi_end, &sec_end);
684 
685  if (mi_end * 60 == 0) {
686  gtime = ((60 * 60 + round(sec_end)) - mi_ini * 60 - round(sec_ini)) / 60;
687  } else {
688  gtime = ((mi_end * 60 + round(sec_end)) - mi_ini * 60 - round(sec_ini)) / 60;
689  }
690 
691  if (gtime > gransize) {
692  cout << "gtime = " << gtime << " is greater than granule size = " << gransize << endl;
693  exit(1);
694  }
695  if ((gtime % gransize) == 0)
696  gfull = "1";
697  else
698  gfull = "0";
699 
700  // time coverage start----
701  string secstr = std::to_string(sec_ini);
702  string mistr = std::to_string(mi_ini);
703  string hstr = std::to_string(h_ini);
704  string daystr = std::to_string(d_ini);
705  string monstr = std::to_string(mo_ini);
706  string yearstr = std::to_string(y_ini);
707 
708  int length = (int)floor(log10(mo_ini)) + 1;
709  if (length == 1)
710  monstr = "0" + monstr;
711  length = (int)floor(log10(d_ini)) + 1;
712  if (length == 1)
713  daystr = "0" + daystr;
714 
715  if (h_ini == 0)
716  logoff = 1;
717  else
718  logoff = 0;
719  length = (int)floor(log10(h_ini + logoff)) + 1;
720  if (length == 1)
721  hstr = "0" + hstr;
722  if (mi_ini == 0)
723  logoff = 1;
724  else
725  logoff = 0;
726  length = (int)floor(log10(mi_ini + logoff)) + 1;
727  if (length == 1)
728  mistr = "0" + mistr;
729  if (sec_ini == 0)
730  logoff = 1;
731  else
732  logoff = 0;
733  length = (int)floor(log10(round(sec_ini + logoff))) + 1;
734  if (length == 1)
735  secstr = "0" + secstr;
736 
737  string fdatetimestr1 = yearstr + monstr + daystr + "T" + hstr + mistr + secstr.substr(0, 2);
738  string datetimestr1 =
739  yearstr + "-" + monstr + "-" + daystr + "T" + hstr + ":" + mistr + ":" + secstr.substr(0, 2);
740  // time coevrage end----
741 
742  secstr = std::to_string(sec_end);
743  mistr = std::to_string(mi_end);
744  hstr = std::to_string(h_end);
745  daystr = std::to_string(d_end);
746  monstr = std::to_string(mo_end);
747  yearstr = std::to_string(y_end);
748 
749  length = (int)floor(log10(mo_end)) + 1;
750  if (length == 1)
751  monstr = "0" + monstr;
752  length = (int)floor(log10(d_end)) + 1;
753  if (length == 1)
754  daystr = "0" + daystr;
755 
756  if (h_end == 0)
757  logoff = 1;
758  else
759  logoff = 0;
760  length = (int)floor(log10(h_end + logoff)) + 1;
761  if (length == 1)
762  hstr = "0" + hstr;
763  if (mi_end == 0)
764  logoff = 1;
765  else
766  logoff = 0;
767  length = (int)floor(log10(mi_end + logoff)) + 1;
768  if (length == 1)
769  mistr = "0" + mistr;
770  if (sec_end == 0)
771  logoff = 1;
772  else
773  logoff = 0;
774  length = (int)floor(log10(round(sec_end + logoff))) + 1;
775  if (length == 1)
776  secstr = "0" + secstr;
777 
778  string datetimestr2 =
779  yearstr + "-" + monstr + "-" + daystr + "T" + hstr + ":" + mistr + ":" + secstr.substr(0, 2);
780  string extstr = ".nc";
781  string fname_out_nopath = binl1c->full_l1cgrid;
782  outf << fname_out_nopath << "," << datetimestr1 << "," << datetimestr2 << "," << gfull << "\n";
783  outf.close();
784 
785  double *time_nad = (double *)calloc(num_gridlines, sizeof(double));
786  NcGroup ba_grp = nc_l1cgrid->getGroup("bin_attributes");
787  v1 = ba_grp.getVar("nadir_view_time");
788  v1.getVar(time_nad);
789  string tmpUnits;
790  v1.getAtt("units").getValues(tmpUnits);
791  ba_grp = nc_output->getGroup("bin_attributes");
792  v1 = ba_grp.getVar("nadir_view_time");
793  v1.putVar(&time_nad[0]);
794  if(!tmpUnits.empty())
795  v1.putAtt("units", tmpUnits);
796  if (time_nad != nullptr)
797  free(time_nad);
798  time_nad = nullptr;
799 
800  geo_grp = nc_l1cgrid->getGroup("geolocation_data");
801 
802  if (binl1c->bintype == 0) {
803  v1 = geo_grp.getVar("latitude");
804  v1.getVar(&binl1c->lat_gd[0][0]);
805  geo_grp = nc_output->getGroup("geolocation_data");
806  v1 = geo_grp.getVar("latitude");
807  v1.putVar(&binl1c->lat_gd[0][0]);
808 
809  geo_grp = nc_l1cgrid->getGroup("geolocation_data");
810  v1 = geo_grp.getVar("longitude");
811  v1.getVar(&binl1c->lon_gd[0][0]);
812  geo_grp = nc_output->getGroup("geolocation_data");
813  v1 = geo_grp.getVar("longitude");
814  v1.putVar(&binl1c->lon_gd[0][0]);
815  }
816 
817  geo_grp = nc_l1cgrid->getGroup("geolocation_data");
818  v1 = geo_grp.getVar("height");
819  // v1=geo_grp.getVar("altitude");
820  v1.getVar(&binl1c->alt[0][0]);
821  geo_grp = nc_output->getGroup("geolocation_data");
822  v1 = geo_grp.getVar("height");
823  v1.putVar(&binl1c->alt[0][0]);
824 
825  // GRING----
826  string gs_bounds, crs;
827  float gs_lat_max=-999, gs_lat_min=-999, gs_lon_max=-999, gs_lon_min=-999;
828  i1 = nc_l1cgrid->getAtt("geospatial_bounds");
829  i1.getValues(gs_bounds);
830  nc_output->putAtt("geospatial_bounds", gs_bounds);
831 
832  i1 = nc_l1cgrid->getAtt("geospatial_bounds_crs");
833  i1.getValues(crs);
834  nc_output->putAtt("geospatial_bounds_crs", crs);
835 
836  i1 = nc_l1cgrid->getAtt("geospatial_lat_max");
837  i1.getValues(&gs_lat_max);
838  nc_output->putAtt("geospatial_lat_max", ncFloat, gs_lat_max);
839 
840  i1 = nc_l1cgrid->getAtt("geospatial_lat_min");
841  i1.getValues(&gs_lat_min);
842  nc_output->putAtt("geospatial_lat_min", ncFloat, gs_lat_min);
843 
844  i1 = nc_l1cgrid->getAtt("geospatial_lon_max");
845  i1.getValues(&gs_lon_max);
846  nc_output->putAtt("geospatial_lon_max", ncFloat, gs_lon_max);
847 
848  i1 = nc_l1cgrid->getAtt("geospatial_lon_min");
849  i1.getValues(&gs_lon_min);
850  nc_output->putAtt("geospatial_lon_min", ncFloat, gs_lon_min);
851 
852  std::vector<NcDim> dimvec2_rad;
853  dimvec2_rad.push_back(vdim);
854  dimvec2_rad.push_back(idim);
855 
856  float *Fobar = l1file->Fobar;
857  float *fwave = l1file->fwave;
858 
859  float fwave_view[NVIEWS][NBANDS], fwhm_view[NVIEWS][NBANDS], Fobar_view[NVIEWS][NBANDS];
860 
861  for (int i = 0; i < NVIEWS; i++) {
862  for (int j = 0; j < NBANDS; j++) {
863  fwave_view[i][j] = fwave[j];
864  fwhm_view[i][j] = 5.0;
865  Fobar_view[i][j] = 10 * Fobar[j];
866  }
867  }
868 
869  // intensity_wavelengths, bandpasses and F0---
870  NcGroup svb_grp = nc_output->getGroup("sensor_views_bands");
871  v1 = svb_grp.addVar("intensity_wavelength", ncFloat, dimvec2_rad); // this is infull meta
872  string longName = "Intensity field center wavelengths at each view";
873  v1.putAtt("long_name", longName);
874  string units = "nm";
875  v1.putAtt("units", units);
876  v1.putAtt("_FillValue", ncFloat, binl1c->fillval2);
877  float valid_min = 320;
878  v1.putAtt("valid_min", ncFloat, valid_min);
879  float valid_max = 2260;
880  v1.putAtt("valid_max", ncFloat, valid_max);
881  v1.putVar(fwave_view);
882 
883  // intensity_wavelengths, bandpasses and F0---
884  v1 = svb_grp.addVar("intensity_bandpass", ncFloat, dimvec2_rad);
885  longName = "Intensity field bandpass at each view";
886  v1.putAtt("long_name", longName);
887  units = "nm";
888  v1.putAtt("units", units);
889  v1.putAtt("_FillValue", ncFloat, binl1c->fillval2);
890  valid_min = 2.5;
891  v1.putAtt("valid_min", ncFloat, valid_min);
892  valid_max = 100;
893  v1.putAtt("valid_max", ncFloat, valid_max);
894  v1.putVar(fwhm_view);
895 
896  v1 = svb_grp.addVar("intensity_f0", ncFloat, dimvec2_rad);
897  longName = "Intensity band solar irradiance";
898  v1.putAtt("long_name", longName);
899  units = "W m^-2 um^-1";
900  v1.putAtt("units", units);
901  v1.putAtt("_FillValue", ncFloat,binl1c->fillval2);
902  valid_min = 0;
903  v1.putAtt("valid_min", ncFloat, valid_min);
904  valid_max = 4000;
905  v1.putAtt("valid_max", ncFloat, valid_max);
906  v1.putVar(Fobar_view);
907 
908  nc_l1cgrid->close();
909 
910  return 0;
911 }
912 
913 int meta_l1c_grid(char *gridname, bin_str *binl1c, int16_t num_gridlines, NcFile *nc_output) {
914  string senstr, titlestr, prodstr, ifile_str;
915  string date_created;
916  int NVIEWS, NBANDS, NBANDS_POL;
917  int32_t xbins, ybins;
918  int nadir_bin_index;
919  ybins = num_gridlines;
920 
921  if (binl1c->sensor == 34) {
922  senstr = "SPEXONE";
923  titlestr = "PACE SPEXone Level-1C Data";
924  nadir_bin_index = 14;
925  xbins = 29;
926  NVIEWS = 5;
927  NBANDS = 400;
928  NBANDS_POL = 50;
929  } else if (binl1c->sensor == 35) {
930  senstr = "HARP2";
931  titlestr = "PACE HARP2 Level-1C Data";
932  nadir_bin_index = 228;
933  xbins = 457;
934  NVIEWS = 90;
935  NBANDS = 1;
936  NBANDS_POL = 1;
937  } else if (binl1c->sensor == 30 || binl1c->sensor == 31) {
938  senstr = "OCI";
939  titlestr = "PACE OCI Level-1C Data";
940  nadir_bin_index = 259;
941  xbins = 519;
942  NVIEWS = 2;
943  NBANDS = 239; // 249 originally no polarization bands
944  } else // OCIS
945  {
946  senstr = "OCI";
947  titlestr = "PACE OCI Level-1C Data";
948  nadir_bin_index =259;
949  xbins = 519;
950  NVIEWS = 2;
951  NBANDS = 239; // 249 ORIGINALLY no polarization bands
952  }
953 
954  prodstr = string(gridname);
955  // get rid of sensor--
956  if (binl1c->verbose)
957  cout << "sensor.." << senstr << "xbins.." << xbins << endl;
958 
959  // creation date---
960  date_created = unix2isodate(now(), 'G');
961 
962  // global attributes---
963 
964  nc_output->putAtt("title", titlestr);
965  nc_output->putAtt("instrument", senstr);
966  nc_output->putAtt("Conventions", "CF-1.8 ACDD-1.3");
967  nc_output->putAtt("institution", "NASA Goddard Space Flight Center, Ocean Biology Processing Group");
968  nc_output->putAtt("license",
969  "https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/");
970  nc_output->putAtt("naming_authority", "gov.nasa.gsfc.sci.oceancolor");
971  nc_output->putAtt("keywords_vocabulary", "NASA Global Change Master Directory (GCMD) Science Keywords");
972  nc_output->putAtt("stdname_vocabulary", "NetCDF Climate and Forecast (CF) Metadata Convention");
973  nc_output->putAtt("creator_name", "NASA/GSFC");
974  nc_output->putAtt("creator_email", "data@oceancolor.gsfc.nasa.gov");
975  nc_output->putAtt("creator_url", "https://oceancolor.gsfc.nasa.gov");
976  nc_output->putAtt("project", "PACE Project");
977  nc_output->putAtt("publisher_name", "NASA/GSFC");
978  nc_output->putAtt("publisher_email", "data@oceancolor.gsfc.nasa.gov");
979  nc_output->putAtt("publisher_url", "https://oceancolor.gsfc.nasa.gov");
980  nc_output->putAtt("processing_level", "L1C");
981  nc_output->putAtt("cdm_data_type", "swath");
982  nc_output->putAtt("product_name", prodstr);
983  nc_output->putAtt("date_created", date_created);
984  nc_output->putAtt("sun_earth_distance", ncFloat, BAD_FLT);
985  nc_output->putAtt("nadir_bin", NC_INT,nadir_bin_index);
986  nc_output->putAtt("bin_size_at_nadir", "5.2km2");
987 
988  NcDim ydim = nc_output->addDim("bins_along_track", ybins);
989  NcDim xdim = nc_output->addDim("bins_across_track", xbins);
990  NcDim vdim = nc_output->addDim("number_of_views", NVIEWS);
991  NcDim idim = nc_output->addDim("intensity_bands_per_view", NBANDS);
992 
993  // dims
994  std::vector<NcDim> dimvec2_geo;
995  dimvec2_geo.push_back(ydim);
996  dimvec2_geo.push_back(xdim);
997  std::vector<NcDim> dimvec2_rad;
998  dimvec2_rad.push_back(vdim);
999  dimvec2_rad.push_back(idim);
1000  std::vector<NcDim> dimvec3;
1001  dimvec3.push_back(ydim);
1002  dimvec3.push_back(xdim);
1003  dimvec3.push_back(vdim);
1004  std::vector<NcDim> dimvec4;
1005  dimvec4.push_back(ydim);
1006  dimvec4.push_back(xdim);
1007  dimvec4.push_back(vdim);
1008  dimvec4.push_back(idim);
1009 
1010  // groups
1011  NcGroup ba_grp = nc_output->addGroup("bin_attributes");
1012  NcGroup geo_grp = nc_output->addGroup("geolocation_data");
1013 
1014  if (binl1c->sensor == 34 || binl1c->sensor == 35) { // SPEX 34, HARP2 35
1015 
1016  NcDim pdim = nc_output->addDim("polarization_bands_per_view", NBANDS_POL);
1017  std::vector<NcDim> dimvec2b_rad;
1018  dimvec2b_rad.push_back(vdim);
1019  dimvec2b_rad.push_back(pdim);
1020  std::vector<NcDim> dimvec4b;
1021  dimvec4b.push_back(ydim);
1022  dimvec4b.push_back(xdim);
1023  dimvec4b.push_back(vdim);
1024  dimvec4b.push_back(pdim);
1025  }
1026 
1027  NcVar v1 = ba_grp.addVar("nadir_view_time", ncDouble, ydim);
1028  string longName = "Time bin was viewed at nadir view";
1029  v1.putAtt("long_name", longName);
1030  string units = "seconds since date";
1031  v1.putAtt("units", units);
1032  v1.putAtt("_FillValue", ncDouble, binl1c->fillval3);
1033  double valid_min_d = 0;
1034  v1.putAtt("valid_min", ncDouble, valid_min_d);
1035  double valid_max_d = 172800;
1036  v1.putAtt("valid_max", ncDouble, valid_max_d);
1037 
1038  v1 = geo_grp.addVar("latitude", ncFloat, dimvec2_geo);
1039  longName = "Latitudes of bin locations";
1040  v1.putAtt("long_name", longName);
1041  units = "degrees_north";
1042  v1.putAtt("units", units);
1043  v1.putAtt("_FillValue", ncFloat, binl1c->fillval2);
1044  float valid_min = -90;
1045  v1.putAtt("valid_min", ncFloat, valid_min);
1046  float valid_max = 90;
1047  v1.putAtt("valid_max", ncFloat, valid_max);
1048 
1049  v1 = geo_grp.addVar("longitude", ncFloat, dimvec2_geo);
1050  longName = "Longitudes of bin locations";
1051  v1.putAtt("long_name", longName);
1052  units = "degrees_east";
1053  v1.putAtt("units", units);
1054  v1.putAtt("_FillValue", ncFloat, binl1c->fillval2);
1055  valid_min = -180;
1056  v1.putAtt("valid_min", ncFloat, valid_min);
1057  valid_max = 180;
1058  v1.putAtt("valid_max", ncFloat, valid_max);
1059 
1060  v1 = geo_grp.addVar("height", ncShort, dimvec2_geo);
1061  longName = "Altitude at bin locations";
1062  v1.putAtt("long_name", longName);
1063  units = "meters";
1064  v1.putAtt("units", units);
1065  v1.putAtt("_FillValue", ncShort, binl1c->fillval1);
1066  // use min/max from gebco file
1067  // short valid_min_s = -10952;
1068  short valid_min_s = -10000;
1069  v1.putAtt("valid_min", ncShort, valid_min_s);
1070 // short valid_max_s = 8627;
1071  short valid_max_s = 10000;
1072  v1.putAtt("valid_max", ncShort, valid_max_s);
1073 
1074  return 0;
1075 }
1076 
1077 int meta_l1c_global(char *gridname, bin_str *binl1c,int16_t num_gridlines, NcFile *nc_output) {
1078  string senstr, GATT_VAL1, prodstr,ifile_str;
1079  string date_created;
1080  int NVIEWS, NBANDS, NBANDS_POL;
1081  int32_t xbins, ybins;
1082  int nadir_bin_index;
1083 
1084  ybins = num_gridlines;
1085  prodstr = string(gridname);
1086 
1087  file_format format = getFormat(gridname);
1088  if (format.type == FT_HKT || format.type == FT_L1C)
1089  format.type = FT_OCIS;
1090 
1091  if (format.type == FT_SPEXONE) {
1092  senstr = "SPEXONE";
1093  GATT_VAL1 = "PACE SPEXone Level-1C Data";
1094  nadir_bin_index =14;
1095  xbins = 29;
1096  NVIEWS = 5;
1097  NBANDS = 400;
1098  NBANDS_POL = 50;
1099  } else if (format.type == FT_HARP2) {
1100  senstr = "HARP2";
1101  GATT_VAL1 = "PACE HARP2 Level-1C Data";
1102  nadir_bin_index =228;
1103  xbins = 457;
1104  NVIEWS = 90;
1105  NBANDS = 1;
1106  NBANDS_POL = 1;
1107  } else if (format.type == FT_OCIS) {
1108  senstr = "OCIS";
1109  GATT_VAL1 = "PACE OCIS Level-1C Data";
1110  nadir_bin_index =259;
1111  xbins = 519;
1112  NVIEWS = 2;
1113  NBANDS = 239; // 249 originally no polarization bands
1114  } else if (format.type == FT_OCIL1B) {
1115  senstr = "OCI";
1116  GATT_VAL1 = "PACE OCI Level-1C Data";
1117  nadir_bin_index =259;
1118  xbins = 519;
1119  NVIEWS = 2;
1120  NBANDS = 286; // 249 originally no polarization bands
1121  } else // OCI
1122  {
1123  senstr = "OCI";
1124  GATT_VAL1 = "PACE OCI Level-1C Data";
1125  nadir_bin_index =259;
1126  xbins = 519;
1127  NVIEWS = 2;
1128  NBANDS = 286; // 249 ORIGINALLY no polarization bands
1129  }
1130 
1131  // views
1132  float views[NVIEWS];
1133  if (format.type == FT_SPEXONE) {
1134  views[0] = -58;
1135  views[1] = -20;
1136  views[2] = 0;
1137  views[3] = 20;
1138  views[4] = 58;
1139  } else if (format.type == FT_OCIS) {
1140  views[0] = -20; // OCI
1141  views[1] = 20;
1142  } else if (format.type == FT_HARP2) {
1143  cout << "# views TBD.....ERROR.." << endl;
1144  exit(1);
1145  } else { // OCI default in HKT---
1146  views[0] = -20; // OCI
1147  views[1] = 20;
1148  }
1149 
1150  // creation date---
1151  date_created = unix2isodate(now(), 'G');
1152 
1153  // global attributes---
1154 
1155  nc_output->putAtt("title", GATT_VAL1);
1156  nc_output->putAtt("instrument", senstr);
1157  nc_output->putAtt("Conventions", "CF-1.8 ACDD-1.3");
1158  nc_output->putAtt("institution", "NASA Goddard Space Flight Center, Ocean Biology Processing Group");
1159  nc_output->putAtt("license",
1160  "https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/");
1161  nc_output->putAtt("naming_authority", "gov.nasa.gsfc.sci.oceancolor");
1162  nc_output->putAtt("keywords_vocabulary", "NASA Global Change Master Directory (GCMD) Science Keywords");
1163  nc_output->putAtt("stdname_vocabulary", "NetCDF Climate and Forecast (CF) Metadata Convention");
1164  nc_output->putAtt("creator_name", "NASA/GSFC");
1165  nc_output->putAtt("creator_email", "data@oceancolor.gsfc.nasa.gov");
1166  nc_output->putAtt("creator_url", "https://oceancolor.gsfc.nasa.gov");
1167  nc_output->putAtt("project", "PACE Project");
1168  nc_output->putAtt("publisher_name", "NASA/GSFC");
1169  nc_output->putAtt("publisher_email", "data@oceancolor.gsfc.nasa.gov");
1170  nc_output->putAtt("publisher_url", "https://oceancolor.gsfc.nasa.gov");
1171  nc_output->putAtt("processing_level", "L1C");
1172  nc_output->putAtt("cdm_data_type", "swath");
1173  nc_output->putAtt("product_name", prodstr);
1174  nc_output->putAtt("date_created", date_created);
1175  nc_output->putAtt("sun_earth_distance", ncFloat, BAD_FLT);
1176  nc_output->putAtt("nadir_bin",NC_INT,nadir_bin_index);
1177  nc_output->putAtt("bin_size_at_nadir", "5.2km2");
1178 
1179  NcDim ydim = nc_output->addDim("bins_along_track", ybins);
1180  NcDim xdim = nc_output->addDim("bins_across_track", xbins);
1181  NcDim vdim = nc_output->addDim("number_of_views", NVIEWS);
1182  NcDim idim = nc_output->addDim("intensity_bands_per_view", NBANDS);
1183 
1184  // dims
1185  std::vector<NcDim> dimvec2_geo;
1186  dimvec2_geo.push_back(ydim);
1187  dimvec2_geo.push_back(xdim);
1188  std::vector<NcDim> dimvec2_rad;
1189  dimvec2_rad.push_back(vdim);
1190  dimvec2_rad.push_back(idim);
1191  std::vector<NcDim> dimvec3;
1192  dimvec3.push_back(ydim);
1193  dimvec3.push_back(xdim);
1194  dimvec3.push_back(vdim);
1195  std::vector<NcDim> dimvec4;
1196  dimvec4.push_back(ydim);
1197  dimvec4.push_back(xdim);
1198  dimvec4.push_back(vdim);
1199  dimvec4.push_back(idim);
1200 
1201  // groups
1202  NcGroup svb_grp = nc_output->addGroup("sensor_views_bands");
1203  NcGroup ba_grp = nc_output->addGroup("bin_attributes");
1204  NcGroup geo_grp = nc_output->addGroup("geolocation_data");
1205  NcGroup od_grp = nc_output->addGroup("observation_data");
1206 
1207 
1208  // vars
1209  NcVar v1 = svb_grp.addVar("sensor_view_angle", ncFloat, vdim);
1210  string longName = "Along-track view angles for sensor";
1211  v1.putAtt("long_name", longName);
1212  string units = "degrees";
1213  v1.putAtt("units", units);
1214  v1.putAtt("_FillValue", ncFloat,binl1c->fillval2);
1215  float valid_min = -89;
1216  v1.putAtt("valid_min", ncFloat, valid_min);
1217  float valid_max = 89;
1218  v1.putAtt("valid_max", ncFloat, valid_max);
1219  v1.putVar(views);
1220 
1221 
1222  if (format.type == FT_HARP2 || format.type == FT_SPEXONE) {
1223  NcDim pdim = nc_output->addDim("polarization_bands_per_view", NBANDS_POL);
1224  std::vector<NcDim> dimvec2b_rad;
1225  dimvec2b_rad.push_back(vdim);
1226  dimvec2b_rad.push_back(pdim);
1227  std::vector<NcDim> dimvec4b;
1228  dimvec4b.push_back(ydim);
1229  dimvec4b.push_back(xdim);
1230  dimvec4b.push_back(vdim);
1231  dimvec4b.push_back(pdim);
1232  }
1233 
1234  v1 = ba_grp.addVar("nadir_view_time", ncDouble, ydim);
1235  longName = "Time bin was viewed at nadir view";
1236  v1.putAtt("long_name", longName);
1237  units = "seconds since date";
1238  v1.putAtt("units", units);
1239  v1.putAtt("_FillValue", ncDouble, binl1c->fillval3);
1240  double valid_min_d = 0;
1241  v1.putAtt("valid_min", ncDouble, valid_min_d);
1242  double valid_max_d = 172800;
1243  v1.putAtt("valid_max", ncDouble, valid_max_d);
1244 
1245  v1 = ba_grp.addVar("view_time_offsets", ncDouble, dimvec3);
1246  longName = "Time offsets of views from nadir view";
1247  v1.putAtt("long_name", longName);
1248  units = "seconds";
1249  v1.putAtt("units", units);
1250  v1.putAtt("_FillValue", ncDouble, binl1c->fillval3);
1251  valid_min_d = -200;
1252  v1.putAtt("valid_min", ncDouble, valid_min_d);
1253  valid_max_d = 200;
1254  v1.putAtt("valid_max", ncDouble, valid_max_d);
1255 
1256  v1 = geo_grp.addVar("latitude", ncFloat, dimvec2_geo);
1257  longName = "Latitudes of bin locations";
1258  v1.putAtt("long_name", longName);
1259  units = "degrees_north";
1260  v1.putAtt("units", units);
1261  v1.putAtt("_FillValue", ncFloat,binl1c->fillval2);
1262  valid_min = -90;
1263  v1.putAtt("valid_min", ncFloat, valid_min);
1264  valid_max = 90;
1265  v1.putAtt("valid_max", ncFloat, valid_max);
1266 
1267  v1 = geo_grp.addVar("longitude", ncFloat, dimvec2_geo);
1268  longName = "Longitudes of bin locations";
1269  v1.putAtt("long_name", longName);
1270  units = "degrees_east";
1271  v1.putAtt("units", units);
1272  v1.putAtt("_FillValue", ncFloat, binl1c->fillval2);
1273  valid_min = -180;
1274  v1.putAtt("valid_min", ncFloat, valid_min);
1275  valid_max = 180;
1276  v1.putAtt("valid_max", ncFloat, valid_max);
1277 
1278  v1 = geo_grp.addVar("height", ncShort, dimvec2_geo);
1279  longName = "Altitude at bin locations";
1280  v1.putAtt("long_name", longName);
1281  units = "meters";
1282  v1.putAtt("units", units);
1283  v1.putAtt("_FillValue", ncShort,binl1c->fillval1);
1284  // use min/max from gebco file
1285  // short valid_min_s = -10952;
1286  short valid_min_s = -10000;
1287  v1.putAtt("valid_min", ncShort, valid_min_s);
1288 // short valid_max_s = 8627;
1289  short valid_max_s=10000;
1290  v1.putAtt("valid_max", ncShort, valid_max_s);
1291 
1292  v1 = geo_grp.addVar("height_stdev", ncShort, dimvec2_geo);
1293  longName = "Standard deviation of terrain altitude within bin";
1294  v1.putAtt("long_name", longName);
1295  units = "meters";
1296  v1.putAtt("units", units);
1297  v1.putAtt("_FillValue", ncShort, binl1c->fillval1);
1298  valid_min_s = 0;
1299  v1.putAtt("valid_min", ncShort, valid_min_s);
1300  valid_max_s = 1000;
1301  v1.putAtt("valid_max", ncShort, valid_max_s);
1302 
1303  v1 = geo_grp.addVar("sensor_azimuth_angle", ncShort, dimvec3);
1304  longName = "Sensor azimuth angles at bin locations";
1305  v1.putAtt("long_name", longName);
1306  units = "degrees";
1307  v1.putAtt("units", units);
1308  v1.putAtt("_FillValue", ncShort, binl1c->fillval1);
1309  valid_min_s = -18000;
1310  v1.putAtt("valid_min", ncShort, valid_min_s);
1311  valid_max_s = 18000;
1312  v1.putAtt("valid_max", ncShort, valid_max_s);
1313  float scale_factor = 0.01;
1314  v1.putAtt("scale_factor", ncFloat, scale_factor);
1315  float add_offset = 0;
1316  v1.putAtt("add_offset", ncFloat, add_offset);
1317 
1318  v1 = geo_grp.addVar("sensor_zenith_angle", ncShort, dimvec3);
1319  longName = "Sensor zenith angles at bin locations";
1320  v1.putAtt("long_name", longName);
1321  units = "degrees";
1322  v1.putAtt("units", units);
1323  v1.putAtt("_FillValue", ncShort,binl1c->fillval1);
1324  valid_min_s = 0;
1325  v1.putAtt("valid_min", ncShort, valid_min_s);
1326  valid_max_s = 18000;
1327  v1.putAtt("valid_max", ncShort, valid_max_s);
1328  scale_factor = 0.01;
1329  v1.putAtt("scale_factor", ncFloat, scale_factor);
1330  add_offset = 0;
1331  v1.putAtt("add_offset", ncFloat, add_offset);
1332 
1333  v1 = geo_grp.addVar("solar_azimuth_angle", ncShort, dimvec3);
1334  longName = "Solar azimuth angle at bin locations";
1335  v1.putAtt("long_name", longName);
1336  units = "degrees";
1337  v1.putAtt("units", units);
1338  v1.putAtt("_FillValue", ncShort,binl1c->fillval1);
1339  valid_min_s = -18000;
1340  v1.putAtt("valid_min", ncShort, valid_min_s);
1341  valid_max_s = 18000;
1342  v1.putAtt("valid_max", ncShort, valid_max_s);
1343  scale_factor = 0.01;
1344  v1.putAtt("scale_factor", ncFloat, scale_factor);
1345  add_offset = 0;
1346  v1.putAtt("add_offset", ncFloat, add_offset);
1347 
1348  v1 = geo_grp.addVar("solar_zenith_angle", ncShort, dimvec3);
1349  longName = "Solar zenith angle at bin locations";
1350  v1.putAtt("long_name", longName);
1351  units = "degrees";
1352  v1.putAtt("units", units);
1353  v1.putAtt("_FillValue", ncShort,binl1c->fillval1);
1354  valid_min_s = 0;
1355  v1.putAtt("valid_min", ncShort, valid_min_s);
1356  valid_max_s = 18000;
1357  v1.putAtt("valid_max", ncShort, valid_max_s);
1358  scale_factor = 0.01;
1359  v1.putAtt("scale_factor", ncFloat, scale_factor);
1360  add_offset = 0;
1361  v1.putAtt("add_offset", ncFloat, add_offset);
1362 
1363  v1 = geo_grp.addVar("scattering_angle", ncShort, dimvec3);
1364  longName = "Scattering angle at bin locations";
1365  v1.putAtt("long_name", longName);
1366  units = "degrees";
1367  v1.putAtt("units", units);
1368  v1.putAtt("_FillValue", ncShort,binl1c->fillval1);
1369  valid_min_s = 0;
1370  v1.putAtt("valid_min", ncShort, valid_min_s);
1371  valid_max_s = 18000;
1372  v1.putAtt("valid_max", ncShort, valid_max_s);
1373  scale_factor = 0.01;
1374  v1.putAtt("scale_factor", ncFloat, scale_factor);
1375  add_offset = 0;
1376  v1.putAtt("add_offset", ncFloat, add_offset);
1377 
1378  // v1 = geo_grp.addVar("rotation_angle", ncShort, dimvec3);
1379  // longName = "Rotation angle at bin locations";
1380  // v1.putAtt("long_name", longName);
1381  // units = "degrees";
1382  // v1.putAtt("units", units);
1383  // v1.putAtt("_FillValue", ncShort, binl1c->fillval1);
1384  // valid_min_s = 0;
1385  // v1.putAtt("valid_min", ncShort, valid_min_s);
1386  // valid_max_s = 18000;
1387  // v1.putAtt("valid_max", ncShort, valid_max_s);
1388  // scale_factor = 0.01;
1389  // v1.putAtt("scale_factor", ncFloat, scale_factor);
1390  // add_offset = 0;
1391  // v1.putAtt("add_offset", ncFloat, add_offset);
1392 
1393  v1 = od_grp.addVar("number_of_observations", ncShort, dimvec3);
1394  longName = "Observations contributing to bin from each view";
1395  v1.putAtt("long_name", longName);
1396  valid_min_s = 0;
1397  v1.putAtt("valid_min", ncShort, valid_min_s);
1398  valid_max_s = 999;
1399  v1.putAtt("valid_max", ncShort, valid_max_s);
1400  v1.putAtt("coordinates", "geolocation_data/longitude geolocation_data/latitude");
1401 
1402  if (format.type == FT_OCIL1B || format.type == FT_OCIS || format.type == FT_HARP2) {
1403  uint8_t FillValue3, valid_min2, valid_max2;
1404 
1405  v1 = od_grp.addVar("qc", ncUbyte, dimvec4);
1406  longName = "quality indicator";
1407  v1.putAtt("long_name", longName);
1408  FillValue3 = 255;
1409  v1.putAtt("_FillValue", ncUbyte, FillValue3);
1410  valid_min2 = 0;
1411  v1.putAtt("valid_min", ncUbyte, valid_min2);
1412  valid_max2 = 10;
1413  v1.putAtt("valid_max", ncUbyte, valid_max2);
1414  v1.putAtt("coordinates", "geolocation_data/longitude geolocation_data/latitude");
1415 
1416  v1 = od_grp.addVar("i", ncFloat, dimvec4);
1417  longName = "I Stokes vector component";
1418  v1.putAtt("long_name", longName);
1419  units = "W m^-2 sr^-1 um^-1";
1420  v1.putAtt("units", units);
1421  v1.putAtt("_FillValue", ncFloat, binl1c->fillval2);
1422  valid_min = 0;
1423  v1.putAtt("valid_min", ncFloat, valid_min);
1424  valid_max = 999;
1425  v1.putAtt("valid_max", ncFloat, valid_max);
1426  v1.putAtt("coordinates", "geolocation_data/longitude geolocation_data/latitude");
1427 
1428  v1 = od_grp.addVar("i_stdev", ncFloat, dimvec4);
1429  longName = "Random stdev of i in bin";
1430  v1.putAtt("long_name", longName);
1431  units = "W m^-2 sr^-1 um^-1";
1432  v1.putAtt("units", units);
1433  v1.putAtt("_FillValue", ncFloat,binl1c->fillval2);
1434  valid_min = 0;
1435  v1.putAtt("valid_min", ncFloat, valid_min);
1436  valid_max = 800;
1437  v1.putAtt("valid_max", ncFloat, valid_max);
1438  v1.putAtt("coordinates", "geolocation_data/longitude geolocation_data/latitude");
1439 
1440  } else if (format.type == FT_SPEXONE) // spexone
1441  {
1442  }
1443 
1444  return 0;
1445 }
1446 
1447 int bintime_l1c(filehandle *l1file, l1str *l1rec, bin_str *binstr, short **gdindex, double scantime,
1448  NcFile *nc_output) {
1449  short gd_row = -1, gd_col = -1;
1450  int npix = l1file->npix;
1451  int view;
1452  double timefill, mintime, maxtime;
1453  float fill_tilt = binstr->fillval2;
1454  float tilt = binstr->fillval2;
1455 
1456  NcGroup ba_grp = nc_output->getGroup("bin_attributes");
1457  NcVar v1 = ba_grp.getVar("nadir_view_time");
1458  NcVarAtt a1 = v1.getAtt("_FillValue"); // root group
1459  a1.getValues(&timefill);
1460  a1 = v1.getAtt("valid_min"); // root group
1461  a1.getValues(&mintime);
1462  a1 = v1.getAtt("valid_max"); // root group
1463  a1.getValues(&maxtime);
1464 
1465  // ASSUMING same time for each line!!
1466  for (int pix = 0; pix < npix; pix++) {
1467  gd_row = gdindex[pix][0] - 1;
1468  gd_col = gdindex[pix][1] - 1;
1469  tilt = l1rec->tilt;
1470  view = BAD_FLT;
1471 
1472  // 1 view at the time for OCI!!!!
1473  if (l1rec->tilt != fill_tilt && l1rec->tilt < -18) { //-19.9 previous
1474  view = 0;
1475  tilt = l1rec->tilt;
1476  } else if (l1rec->tilt != fill_tilt && l1rec->tilt > 18) { // 19.9 previous
1477  view = 1;
1478  tilt = l1rec->tilt;
1479  }
1480 
1481  if (gd_row >= 2000) {
1482  tilt = 22;
1483  view = 1;
1484  // if(binstr->verbose) cout<<"pix # "<<pix+1<<"pos tilt beyond L1C
1485  // grid"<<tilt<<"gd_row"<<gd_row<<endl;
1486  tilt = fill_tilt;
1487  }
1488  if (gd_row <= -1) {
1489  tilt = -22;
1490  view = 0;
1491  // if(binstr->verbose) cout<<"pix # "<<pix+1<<"neg tilt beyond L1C
1492  // grid"<<tilt<<"gd_row"<<gd_row<<endl;
1493  tilt = fill_tilt;
1494  }
1495 
1496  // if(binstr->verbose) cout<<"pix # "<<pix+1<<"scantime "<<scantime<<"min scantime
1497  // "<<mintime<<"max scantime "<<maxtime<<"gd_row "<<gd_row<<"gd_col "<<gd_col<<"tilt
1498  // "<<tilt<<endl;
1499  // Cumulative values---
1500  if (view != BAD_FLT && tilt != fill_tilt && gd_row >= 0 && gd_col >= 0 &&
1501  gd_row <= binstr->num_gridlines - 1 && gd_col <= binstr->nbinx - 1 && scantime != timefill &&
1502  scantime >= mintime && scantime <= maxtime) { // and row<=ybin and col<xbin, I need time counter
1503  binstr->time_l1b[gd_row][gd_col][view] = binstr->time_l1b[gd_row][gd_col][view] + scantime;
1504  // if(binstr->verbose) cout<<"pix # "<<pix+1<<"binning time at row
1505  // "<<gd_row+1<<"col "<<gd_col+1<<endl;
1506  }
1507 
1508  } // end pixel
1509 
1510  return 0;
1511 }
1512 
1513 int parallax(filehandle *l1file, const char *l1c_anc, const char *l1c_grid, l1str *l1rec, bin_str *binl1c,
1514  short **gdindex, NcFile *nc_output,int32_t sline,int firstcall) {
1515  short gd_row = -1, gd_col = -1;
1516  string l1c_str = l1c_grid;
1517  string l1c_str2 = l1c_anc;
1518  float *lat_new = nullptr, *lon_new = nullptr;
1519  float dv, Re = 6378.137, Hsat = 676.5;
1520  float **cth = nullptr, **height = nullptr;
1521  float scale_factor;
1522  short sena_min, sena_max, senz_min, senz_max;
1523  int result = 0,status;
1524  float look_angle,*scan_angle=nullptr;//in addition to ccd bands also there is scan angle swir
1525  int navGrp,scangId;
1526  size_t start[] = { 0, 0,};
1527  size_t count[] = { 1, 1 };
1528  start[0]=sline;
1529  start[1]=0;
1530  count[0]=1;
1531  count[1]=l1file->npix;
1532 
1533  if(firstcall==0){ scan_angle = (float*)calloc(l1file->npix,sizeof(float));firstcall=1;}
1534 
1535 
1536 
1537  status = nc_inq_grp_ncid(l1file->sd_id, "navigation_data", &navGrp);
1538  check_err(status, __LINE__, __FILE__);
1539  status = nc_inq_varid(navGrp, "CCD_scan_angles", &scangId);
1540  check_err(status, __LINE__, __FILE__);
1541  status = nc_get_vara_float(navGrp, scangId, start, count, scan_angle);
1542  check_err(status, __LINE__, __FILE__);
1543 
1544  NcFile *nc_l1cgrid;
1545  try {
1546  nc_l1cgrid = new NcFile(l1c_grid, NcFile::read);
1547  } catch (NcException &e) {
1548  e.what();
1549  cerr << "l1cgen l1c_pflag= 8:: Failure reading L1C grid: " + l1c_str << endl;
1550  exit(1);
1551  }
1552 
1553  NcDim yd1 = nc_l1cgrid->getDim("bins_along_track");
1554  NcDim xd1 = nc_l1cgrid->getDim("bins_across_track");
1555  // get height
1556  NcGroup geo_grp = nc_l1cgrid->getGroup("geolocation_data");
1557  NcVar v1 = geo_grp.getVar("height");
1558  height = allocate2d_float(yd1.getSize(), xd1.getSize());
1559  v1.getVar(&height[0][0]);
1560 
1561  // anc file
1562  NcFile *nc_l1canc;
1563  if (binl1c->cloudem_flag == 1){
1564  try {
1565  nc_l1canc = new NcFile(l1c_anc, NcFile::read);
1566  } catch (NcException &e) {
1567  e.what();
1568  cerr << "l1cgen l1c_pflag= 8:: Failure reading L1C-ANCfile: " + l1c_str2 << endl;
1569  exit(1);
1570  }
1571 
1572  // 2-d as L1C grid but different nlines, get cth_ice_cloud, cth_water_cloud vars
1573  NcDim yd2 = nc_l1canc->getDim("bins_along_track");
1574  NcDim xd2 = nc_l1canc->getDim("bins_across_track");
1575  cth = allocate2d_float(yd2.getSize(), xd2.getSize());
1576  if (binl1c->cloud_type == 0) {
1577  v1 = nc_l1canc->getVar("cth_water_cloud");
1578  v1.getVar(&cth[0][0]);
1579  }
1580  if (binl1c->cloud_type == 1) {
1581  v1 = nc_l1canc->getVar("cth_ice_cloud");
1582  v1.getVar(&cth[0][0]);
1583  }
1584  }
1585  // compute displacement and new lat/lon for each pixel and 1 LINE!
1586  // compute new lat and lon
1587 
1588  NcGroup geo_grp2 = nc_output->getGroup("geolocation_data");
1589  v1 = geo_grp2.getVar("sensor_azimuth_angle");
1590  NcVarAtt a1 = v1.getAtt("scale_factor"); // root group
1591  a1.getValues(&scale_factor);
1592  a1 = v1.getAtt("valid_min"); // root group
1593  a1.getValues(&sena_min);
1594  sena_min *= scale_factor;
1595  a1 = v1.getAtt("valid_max"); // root group
1596  a1.getValues(&sena_max);
1597  sena_max *= scale_factor;
1598  v1 = geo_grp2.getVar("sensor_zenith_angle");
1599  a1 = v1.getAtt("scale_factor"); // root group
1600  a1.getValues(&scale_factor);
1601  a1 = v1.getAtt("valid_min"); // root group
1602  a1.getValues(&senz_min);
1603  senz_min *= scale_factor;
1604  a1 = v1.getAtt("valid_max"); // root group
1605  a1.getValues(&senz_max);
1606  senz_max *= scale_factor;
1607 
1608  lat_new = (float *)calloc(l1file->npix, sizeof(float));
1609  lon_new = (float *)calloc(l1file->npix, sizeof(float));
1610 
1611  for (int pix = 0; pix < l1file->npix; pix++) {
1612  gd_row = gdindex[pix][0] - 1;
1613  gd_col = gdindex[pix][1] - 1;
1614  // Displacement vector -- wang et al. 2011 ----
1615  if (l1rec->senz[pix]!= binl1c->fillval1 && l1rec->sena[pix]!= binl1c->fillval1 && l1rec->senz[pix] >= senz_min && l1rec->senz[pix] <= senz_max && l1rec->sena[pix] >= sena_min &&
1616  l1rec->sena[pix] <= sena_max) {
1617  look_angle=acos(cos(l1rec->tilt * M_PI / 180)*cos(scan_angle[pix]*M_PI/180));
1618  if (binl1c->cloudem_flag == 1) // CTH constant in km
1619  {
1620  dv = Hsat * cth[gd_row][gd_col] * tan(look_angle) /
1621  (Hsat - cth[gd_row][gd_col]);
1622  } else if (binl1c->cloudem_flag == 2) // DEM constant in km
1623  {
1624  dv = Hsat * 0.001 * l1rec->height[pix] * tan(look_angle) /
1625  (Hsat - 0.001 * l1rec->height[pix]);
1626  }
1627 
1628  lat_new[pix] = l1rec->lat[pix] * M_PI / 180 - dv * cos(l1rec->sena[pix] * M_PI / 180 + M_PI) / Re;
1629  lon_new[pix] = l1rec->lon[pix] * M_PI / 180 - dv * sin((l1rec->sena[pix] * M_PI / 180 + M_PI)) /
1630  (Re * cos(lat_new[pix] * M_PI / 180));
1631 
1632  if (lon_new[pix] * 180 / M_PI > 180) {
1633  lon_new[pix] -= 2 * M_PI;
1634  }
1635  if (lon_new[pix] * 180 / M_PI < -180) {
1636  lon_new[pix] += 2 * M_PI;
1637  }
1638  lat_new[pix] *= 180 / M_PI;
1639  lon_new[pix] *= 180 / M_PI;
1640  } else {
1641  lat_new[pix] = binl1c->fillval2;
1642  lon_new[pix] = binl1c->fillval2;
1643  }
1644  if(pix==635 && binl1c->verbose) cout<<"sline # "<<sline+1<<"look_angle in m "<<look_angle*180/M_PI<<"dv "<<dv<<"scan_angle "<<180/M_PI*scan_angle[pix]<<"lat_new "<<lat_new[pix]<<"lon_new "<<lon_new[pix]<<"height "<<l1rec->height[pix]<<endl;
1645  }
1646 
1647  // recompute row/col
1648  result = search_l1c_parallax(l1file, lat_new, lon_new, binl1c, gdindex);
1649 
1650  nc_l1cgrid->close();
1651  if (binl1c->cloudem_flag == 1) nc_l1canc->close();
1652  if (cth != nullptr)
1653  free2d_float(cth);
1654  if (height != nullptr)
1655  free2d_float(height);
1656  if (lat_new != nullptr)
1657  free(lat_new);
1658  if (lon_new != nullptr)
1659  free(lon_new);
1660 
1661  if(sline==l1file->nscan-1){ if (scan_angle != nullptr) free(scan_angle);}
1662 
1663  return result;
1664 }
1665 
1666 int bin_l1c(filehandle *l1file, l1str *l1rec, bin_str *binl1c, short **gdindex, NcFile *nc_output) {
1667  int32_t ibp = 0, sb = 0;
1668  short gd_row = -1, gd_col = -1;
1669  int npix = l1file->npix;
1670  int nbands = l1file->nbands;
1671  int view;
1672 
1673  // int nadpix=(npix-1)/2;
1674  float term1, term2, term3, sca_pix, cosu, cose;
1675  float scale_all[7], offset_all[7];
1676  short minval1_all[7], maxval1_all[7];
1677  float minval2_all[7], maxval2_all[7];
1678  float fill_tilt = binl1c->fillval2;
1679  short fill_height = binl1c->fillval2;
1680  float tilt = binl1c->fillval2;
1681  double rotangle;
1682  float seed_mean;
1683  NcGroup geo_grp = nc_output->getGroup("geolocation_data");
1684  NcGroup od_grp = nc_output->getGroup("observation_data");
1685 
1686  // altitude
1687  NcVar v1 = geo_grp.getVar("height");
1688  scale_all[0] = 1; // no scale and offset for height
1689  offset_all[0] = 0;
1690  NcVarAtt a1 = v1.getAtt("valid_min"); // root group
1691  a1.getValues(&minval1_all[0]);
1692  a1 = v1.getAtt("valid_max"); // root group
1693  a1.getValues(&maxval1_all[0]);
1694 
1695  // sensor azimuth
1696  v1 = geo_grp.getVar("sensor_azimuth_angle");
1697  a1 = v1.getAtt("scale_factor"); // root group
1698  a1.getValues(&scale_all[1]);
1699  a1 = v1.getAtt("add_offset"); // root group
1700  a1.getValues(&offset_all[1]);
1701  a1 = v1.getAtt("valid_min"); // root group
1702  a1.getValues(&minval1_all[1]);
1703  a1 = v1.getAtt("valid_max"); // root group
1704  a1.getValues(&maxval1_all[1]);
1705 
1706  // sensor zenith
1707  v1 = geo_grp.getVar("sensor_zenith_angle");
1708  a1 = v1.getAtt("scale_factor"); // root group
1709  a1.getValues(&scale_all[2]);
1710  a1 = v1.getAtt("add_offset"); // root group
1711  a1.getValues(&offset_all[2]);
1712  a1 = v1.getAtt("valid_min"); // root group
1713  a1.getValues(&minval1_all[2]);
1714  a1 = v1.getAtt("valid_max"); // root group
1715  a1.getValues(&maxval1_all[2]);
1716 
1717  minval1_all[2] = scale_all[2] * minval1_all[2] + offset_all[2];
1718  maxval1_all[2] = scale_all[2] * maxval1_all[2] + offset_all[2];
1719 
1720  // solar zenith
1721  v1 = geo_grp.getVar("solar_azimuth_angle");
1722  a1 = v1.getAtt("scale_factor"); // root group
1723  a1.getValues(&scale_all[3]);
1724  a1 = v1.getAtt("add_offset"); // root group
1725  a1.getValues(&offset_all[3]);
1726  a1 = v1.getAtt("valid_min"); // root group
1727  a1.getValues(&minval1_all[3]);
1728  a1 = v1.getAtt("valid_max"); // root group
1729  a1.getValues(&maxval1_all[3]);
1730 
1731  // solar zenith
1732  v1 = geo_grp.getVar("solar_zenith_angle");
1733  a1 = v1.getAtt("scale_factor"); // root group
1734  a1.getValues(&scale_all[4]);
1735  a1 = v1.getAtt("add_offset"); // root group
1736  a1.getValues(&offset_all[4]);
1737  a1 = v1.getAtt("valid_min"); // root group
1738  a1.getValues(&minval1_all[4]);
1739  a1 = v1.getAtt("valid_max"); // root group
1740  a1.getValues(&maxval1_all[4]);
1741 
1742  // scattering angle
1743  v1 = geo_grp.getVar("scattering_angle");
1744  a1 = v1.getAtt("scale_factor"); // root group
1745  a1.getValues(&scale_all[5]);
1746  a1 = v1.getAtt("add_offset"); // root group
1747  a1.getValues(&offset_all[5]);
1748  a1 = v1.getAtt("valid_min"); // root group
1749  a1.getValues(&minval1_all[5]);
1750  a1 = v1.getAtt("valid_max"); // root group
1751  a1.getValues(&maxval1_all[5]);
1752 
1753  // rotation angle
1754  // v1 = geo_grp.getVar("rotation_angle");
1755  // a1 = v1.getAtt("scale_factor"); // root group
1756  // a1.getValues(&scale_all[6]);
1757  // a1 = v1.getAtt("add_offset"); // root group
1758  // a1.getValues(&offset_all[6]);
1759  // a1 = v1.getAtt("valid_min"); // root group
1760  // a1.getValues(&minval1_all[6]);
1761  // a1 = v1.getAtt("valid_max"); // root group
1762  // a1.getValues(&maxval1_all[6]);
1763 
1764  // I
1765  v1 = od_grp.getVar("i");
1766  a1 = v1.getAtt("valid_min"); // root group
1767  a1.getValues(&minval2_all[0]);
1768  a1 = v1.getAtt("valid_max"); // root group
1769  a1.getValues(&maxval2_all[0]);
1770 
1771  for (int pix = 0; pix < npix; pix++) {
1772  gd_row = gdindex[pix][0] - 1;
1773  gd_col = gdindex[pix][1] - 1;
1774  tilt = l1rec->tilt;
1775  view = BAD_FLT;
1776 
1777  if (l1rec->tilt != fill_tilt && l1rec->tilt < -18) { //-19.9 previous
1778  view = 0;
1779  tilt = l1rec->tilt;
1780  } else if (l1rec->tilt != fill_tilt && l1rec->tilt > 18) { // 19.9 previous
1781  view = 1;
1782  tilt = l1rec->tilt;
1783  }
1784  if (gd_row >= 2000) {
1785  tilt = 22;
1786  view = 1;
1787  // if(binl1c->verbose) cout<<"pix # "<<pix+1<<"pos tilt beyond L1C
1788  // grid"<<tilt<<"gd_row"<<gd_row<<endl;
1789  tilt = fill_tilt;
1790  }
1791  if (gd_row <= -1) {
1792  tilt = -22;
1793  view = 0;
1794  // if(binl1c->verbose) cout<<"pix # "<<pix+1<<"neg tilt beyond L1C
1795  // grid"<<tilt<<"gd_row"<<gd_row<<endl;
1796  tilt = fill_tilt;
1797  }
1798 
1799  // 2-D arrays---needing to add fillvalue != for l1rec->alt[pix]
1800  if (view != BAD_FLT && tilt != fill_tilt && gd_row >= 0 && gd_col >= 0 &&
1801  gd_row <= binl1c->num_gridlines - 1 && gd_col <= binl1c->nbinx - 1 &&
1802  l1rec->lat[pix] != binl1c->fillval2 && l1rec->lon[pix] != binl1c->fillval2) {
1803  if (l1rec->height[pix] != fill_height) // ADD MIN/MAX CONSTRAINT
1804  {
1805  binl1c->nrec_2D[gd_row][gd_col] += 1;
1806  binl1c->alt_2D[gd_row][gd_col] =
1807  binl1c->alt_2D[gd_row][gd_col] + l1rec->height[pix];
1808  }
1809  binl1c->inpix++;
1810 
1811  } else {
1812  binl1c->outpix++;
1813  }
1814 
1815  // Cumulative values---
1816  // geometry
1817  if (view != BAD_FLT && tilt != fill_tilt && gd_row >= 0 && gd_col >= 0 &&
1818  gd_row <= binl1c->num_gridlines - 1 && gd_col <= binl1c->nbinx - 1 &&
1819  l1rec->senz[pix] != binl1c->fillval1 && l1rec->senz[pix] >= minval1_all[2] &&
1820  l1rec->senz[pix] <= maxval1_all[2]) {
1821  binl1c->nrec_3D[gd_row][gd_col][view] += 1;
1822  binl1c->sena_3D[gd_row][gd_col][view] = binl1c->sena_3D[gd_row][gd_col][view] + l1rec->sena[pix];
1823  binl1c->senz_3D[gd_row][gd_col][view] =
1824  binl1c->senz_3D[gd_row][gd_col][view] + l1rec->senz[pix]; // sensor zenith goes from 0 to 180
1825  binl1c->suna_3D[gd_row][gd_col][view] = binl1c->suna_3D[gd_row][gd_col][view] + l1rec->sola[pix];
1826  binl1c->sunz_3D[gd_row][gd_col][view] = binl1c->sunz_3D[gd_row][gd_col][view] + l1rec->solz[pix];
1827 
1828  cose = cos((l1rec->senz[pix] + 180) * M_PI / 180);
1829  cosu = cos((l1rec->solz[pix]) * M_PI / 180);
1830  term1 = cose * cosu;
1831  term2 = sqrt((1 - cose * cose)) * sqrt((1 - cosu * cosu));
1832  term3 = cos((l1rec->sena[pix] + 180 - l1rec->sola[pix]) * M_PI / 180);
1833  sca_pix = acos(term1 + term2 * term3) * 180 / M_PI;
1834  binl1c->sca_3D[gd_row][gd_col][view] = binl1c->sca_3D[gd_row][gd_col][view] + sca_pix;
1835  rotangle = rot_angle(l1rec->senz[pix], l1rec->solz[pix], l1rec->sena[pix],
1836  l1rec->sola[pix]); // in degrees
1837  binl1c->rot_angle[gd_row][gd_col][view] = binl1c->rot_angle[gd_row][gd_col][view] + rotangle;
1838  }
1839 
1840  ibp = pix * nbands;
1841 
1842  for (sb = 0; sb < nbands; sb++) {
1843  if (view != BAD_FLT && tilt != fill_tilt && gd_row >= 0 && gd_col >= 0 &&
1844  gd_row <= binl1c->num_gridlines - 1 && gd_col <= binl1c->nbinx - 1 && l1rec->Lt[ibp] != binl1c->fillval2 && l1rec->Lt[ibp] >= minval2_all[0] && 10 * l1rec->Lt[ibp] <= maxval2_all[0]){
1845 /*
1846  if (view != BAD_FLT && tilt != fill_tilt && gd_row >= 0 && gd_col >= 0 &&
1847  gd_row <= binl1c->num_gridlines - 1 && gd_col <= binl1c->nbinx - 1 && l1rec->Lt[ibp] != binl1c->fillval2 && 10 * l1rec->Lt[ibp] >= minval2_all[0] &&
1848  10 * l1rec->Lt[ibp] <= maxval2_all[0] && l1rec->senz[pix] != binl1c->fillval1 && l1rec->senz[pix] >= minval1_all[2]&& l1rec->senz[pix] <= maxval1_all[2]) {
1849 */
1850  binl1c->I_4D[gd_row][gd_col][view][sb] =
1851  binl1c->I_4D[gd_row][gd_col][view][sb] + 10 * l1rec->Lt[ibp];
1852  binl1c->nrec_4D_band[gd_row][gd_col][view][sb] += 1;
1853  binl1c->nrec_3D_view[gd_row][gd_col][view] += 1;
1854 
1855  if(binl1c->i_diff2[gd_row][gd_col][view][sb]==0)
1856  {
1857  seed_mean=10 * l1rec->Lt[ibp];
1858  binl1c->i_diff2[gd_row][gd_col][view][sb]=0.01;//offset substracted before stdev
1859  }
1860  else{
1861  binl1c->i_diff2[gd_row][gd_col][view][sb]+=(10 * l1rec->Lt[ibp]-seed_mean)*(10 * l1rec->Lt[ibp]-seed_mean);
1862  }
1863 
1864  if (gd_row > binl1c->num_gridlines - 1 || gd_col > binl1c->nbinx - 1) {
1865  if (binl1c->verbose) {
1866  cout << "ERROR IN BINNING Lt/rhot --gd_row/gd-col outside the grid limits" << endl;
1867  cout << "row.." << gd_row << "gd_col.." << gd_col << endl;
1868  }
1869  }
1870  // missing QC an QC_bitwise variables
1871  }
1872  ibp++;
1873  } // bands
1874 
1875  } // pixels
1876 
1877  return 0;
1878 }
1879 
1880 int check_l1c_time(const char *l1b_file, const char *l1c_grid, bin_str *binl1c) {
1881  string l1c_str = l1c_grid;
1882  string l1b_str = l1b_file;
1883 
1884  NcFile *nc_l1cgrid;
1885  try {
1886  nc_l1cgrid = new NcFile(l1c_grid, NcFile::read);
1887  } catch (NcException &e) {
1888  e.what();
1889  cerr << "l1cgen l1c_pflag= 8:: Failure reading L1C grid: " + l1c_str << endl;
1890  exit(1);
1891  }
1892 
1893  NcDim yd = nc_l1cgrid->getDim("bins_along_track");
1894  int32_t num_gridlines = yd.getSize();
1895  double *time_nad_l1c = (double *)calloc(num_gridlines, sizeof(double));
1896  NcGroup ba_grp = nc_l1cgrid->getGroup("bin_attributes");
1897  NcVar v1 = ba_grp.getVar("nadir_view_time");
1898  v1.getVar(time_nad_l1c);
1899 
1900  NcFile *nc_l1b;
1901 
1902  try {
1903  nc_l1b = new NcFile(l1b_file, NcFile::read);
1904  } catch (NcException &e) {
1905  e.what();
1906  cerr << "l1cgen l1c_pflag= 8:: Failure reading L1B or L1C-ANCfile: " + l1b_str << endl;
1907  exit(1);
1908  }
1909 
1910  // only if L1C ANC file
1911  int32_t nscans;
1912 
1913  if (binl1c->l1c_anc[0] != '\0' && strcmp(binl1c->l1c_anc, l1b_file) == 0) {
1914  if (binl1c->verbose)
1915  cout << "L1C ANC file provided" << endl;
1916  yd = nc_l1b->getDim("bins_along_track");
1917  int32_t num_gridlines2 = yd.getSize();
1918  nscans = num_gridlines2;
1919 
1920  if (num_gridlines != num_gridlines2) {
1921  cerr << "WARNING number of gridlines ARE NOT THE SAME IN L1C grid and ANC files"
1922  << "L1C grid # " << num_gridlines << "ANC file # " << num_gridlines2 << endl;
1923  // exit(1);
1924  }
1925  } else {
1926  yd = nc_l1b->getDim("number_of_scans");
1927  nscans = yd.getSize();
1928  }
1929 
1930  double tend_l1b, *time_l1b = nullptr;
1931 
1932  if (binl1c->l1c_anc[0] == '\0' || strcmp(binl1c->l1c_anc, l1b_file) != 0) {
1933  time_l1b = (double *)calloc(nscans, sizeof(double));
1934  NcGroup sla_grp = nc_l1b->getGroup("scan_line_attributes");
1935  v1 = sla_grp.getVar("time");
1936  v1.getVar(time_l1b);
1937  tend_l1b = time_l1b[nscans - 1];
1938  } else {
1939  time_l1b = time_nad_l1c;
1940  tend_l1b = time_l1b[nscans - 2];
1941  }
1942 
1943  // check time limits
1944  double tini_l1c = time_nad_l1c[0];
1945  double tend_l1c = time_nad_l1c[num_gridlines - 1];
1946  double tini_l1b = time_l1b[0];
1947 
1948  binl1c->tini_l1c = tini_l1c;
1949  binl1c->tend_l1c = tend_l1c;
1950  binl1c->tini_l1b = tini_l1b;
1951  binl1c->tend_l1b = tend_l1b;
1952 
1953  int time_flag = -1;
1954  if (tini_l1b > (tini_l1c - 5 * 60 - 1) && tini_l1b < (tend_l1c + 5 * 60 + 1)) {
1955  time_flag = 0;
1956  if (binl1c->verbose) {
1957  cout << "OK--L1B granule is INSIDE of the L1C grid--" << endl;
1958  cout << "tini_l1c.." << tini_l1c << "tend_l1c.." << tend_l1c << "tini_l1b.." << tini_l1b
1959  << "tend_l1b.." << tend_l1b << endl;
1960  cout << "tini_l1c-5*60" << tini_l1c - 5 * 60 << "tend_l1c+5*60" << tend_l1c + 5 * 60 << endl;
1961  }
1962  } else {
1963  time_flag = 1;
1964  if (binl1c->verbose) {
1965  cout << "ERROR--L1B granule is outside of the L1C grid--" << endl;
1966  cout << "tini_l1c.." << tini_l1c << "tend_l1c.." << tend_l1c << "tini_l1b.." << tini_l1b
1967  << "tend_l1b.." << tend_l1b << endl;
1968  cout << "tini_l1c-5*60" << tini_l1c - 5 * 60 << "tend_l1c+5*60" << tend_l1c + 5 * 60 << endl;
1969  }
1970  }
1971 
1972  if (time_nad_l1c != nullptr)
1973  free(time_nad_l1c);
1974  time_nad_l1c = nullptr;
1975 
1976  if (binl1c->l1c_anc[0] == '\0' || strcmp(binl1c->l1c_anc, l1b_file) != 0) {
1977  if (time_l1b != nullptr)
1978  free(time_l1b);
1979  time_l1b = nullptr;
1980  }
1981 
1982  nc_l1cgrid->close();
1983  nc_l1b->close();
1984 
1985  return time_flag;
1986 }
1987 
1988 int open_l1c(const char *ifile_l1c, size_t *ybins, size_t *xbins, float **lat_gd, float **lon_gd) {
1989  std::string str;
1990  std::string ifile_str;
1991  string gridname, azeast_name;
1992  std::string fname_out, pathstr, senstr, monstr, daystr, yearstr, prodstr, gdstr, swtstr, swtnum, extstr,
1993  granstr, timestr, azstr, missionstr, ofilestr;
1994 
1995  // if(binl1c->verbose)cout<<"Opening L1C
1996  // grid........................................................................."<<endl;
1997  string l1c_str(ifile_l1c);
1998 
1999  NcFile *nc_l1cgrid;
2000 
2001  try {
2002  nc_l1cgrid = new NcFile(ifile_l1c, NcFile::read);
2003  } catch (NcException &e) {
2004  e.what();
2005  cerr << "l1cgen :: Failure reading L1C grid: " + l1c_str << endl;
2006  exit(1);
2007  }
2008  NcGroup geo_grp = nc_l1cgrid->getGroup("geolocation_data");
2009  NcVar v1 = geo_grp.getVar("latitude");
2010  v1.getVar(&lat_gd[0][0]);
2011  v1 = geo_grp.getVar("longitude");
2012  v1.getVar(&lon_gd[0][0]);
2013 
2014  nc_l1cgrid->close();
2015 
2016  return 0;
2017 }
2018 
2019 int open_l1c_grid(const char *ifile_l1c, bin_str *binl1c, float **lat_gd, float **lon_gd, float **alt_gd) {
2020  std::string str;
2021  std::string ifile_str;
2022  string gridname, azeast_name;
2023  std::string fname_out, pathstr, senstr, monstr, daystr, yearstr, prodstr, gdstr, swtstr, swtnum, extstr,
2024  granstr, timestr, azstr, missionstr, ofilestr;
2025 
2026  // if(binl1c->verbose)cout<<"Opening L1C
2027  // grid........................................................................."<<endl;
2028  string l1c_str(ifile_l1c);
2029 
2030  NcFile *nc_l1cgrid;
2031 
2032  try {
2033  nc_l1cgrid = new NcFile(ifile_l1c, NcFile::read);
2034  } catch (NcException &e) {
2035  e.what();
2036  cerr << "l1cgen :: Failure reading L1C grid: " + l1c_str << endl;
2037  exit(1);
2038  }
2039  NcGroup geo_grp = nc_l1cgrid->getGroup("geolocation_data");
2040  NcVar v1 = geo_grp.getVar("latitude");
2041  v1.getVar(&lat_gd[0][0]);
2042  v1 = geo_grp.getVar("longitude");
2043  v1.getVar(&lon_gd[0][0]);
2044  v1 = geo_grp.getVar("height");
2045  v1.getVar(&alt_gd[0][0]);
2046 
2047  NcDim ydim = nc_l1cgrid->getDim("bins_along_track");
2048  NcDim xdim = nc_l1cgrid->getDim("bins_across_track");
2049  binl1c->num_gridlines = ydim.getSize();
2050  binl1c->nbinx = xdim.getSize();
2051 
2052  nc_l1cgrid->close();
2053 
2054  return 0;
2055 }
2056 
2057 
2058 double search_calc_dotprod(bin_str *binl1c, float bvec[3], int32_t gridline) {
2059  int32_t nbinx = binl1c->nbinx;
2060  float gnvec[3];
2061  float gnvm;
2062 
2063  // normal vectors for L1C rows
2064  if (binl1c->lat_gd[gridline][nbinx - 1] > 90 || binl1c->lat_gd[gridline][nbinx - 1] < -90 ||
2065  binl1c->lon_gd[gridline][nbinx - 1] < -180 || binl1c->lon_gd[gridline][nbinx - 1] > 180) {
2066  if (binl1c->verbose)
2067  cout << "lat lon for L1C GRID out of the boundaries.." << endl;
2068  exit(1);
2069  }
2070  if (binl1c->lat_gd[gridline][0] > 90 || binl1c->lat_gd[gridline][0] < -90 || binl1c->lon_gd[gridline][0] < -180 ||
2071  binl1c->lon_gd[gridline][0] > 180) {
2072  if (binl1c->verbose)
2073  cout << "lat lon for L1C GRID out of the boundaries.." << endl;
2074  exit(1);
2075  }
2076 
2077  gnvec[0] = sin(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2078  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2079  sin(binl1c->lat_gd[gridline][0] * M_PI / 180) -
2080  sin(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2081  sin(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180);
2082  gnvec[1] = sin(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2083  cos(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180) -
2084  cos(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2085  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2086  sin(binl1c->lat_gd[gridline][0] * M_PI / 180);
2087  gnvec[2] = cos(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2088  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2089  sin(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180) -
2090  sin(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2091  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2092  cos(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180);
2093 
2094  // vector norm
2095  gnvm = sqrt(gnvec[0] * gnvec[0] + gnvec[1] * gnvec[1] + gnvec[2] * gnvec[2]);
2096  if (isnan(gnvm) == 1) {
2097  if (binl1c->verbose)
2098  cout << "NAN value for gnvm.." << endl;
2099  exit(1);
2100  }
2101  if (gnvm == 0) {
2102  if (binl1c->verbose)
2103  cout << "ERROR gnvm == 0--- WE CANT NORMALIZE..." << endl;
2104  exit(1);
2105  }
2106  // normalization
2107  gnvec[0] = gnvec[0] / gnvm;
2108  gnvec[1] = gnvec[1] / gnvm;
2109  gnvec[2] = gnvec[2] / gnvm;
2110 
2111  // for each pixels
2112  // dot prod, orbital normaliz and transposed by pixel vector
2113  return gnvec[0] * bvec[0] + gnvec[1] * bvec[1] + gnvec[2] * bvec[2];
2114 }
2115 
2116 
2117 
2118 
2119 int search_l1c_parallax(filehandle *l1file, float *lat_new, float *lon_new, bin_str *binl1c,
2120  short **gdindex) {
2121  int32_t num_gridlines, nbinx;
2122  int flag_out = -1;
2123  size_t pix,badpix;
2124  int32_t i;
2125  size_t num_pixels;
2126  int irow = -1, col = -1;
2127  float bmcm, bm = 100;
2128  double db;
2129  float c1, c2, c3;
2130  double dotprod, dot_firstline, dot_lastline;
2131  int32_t last_dotprod_index = 0;
2132  double last_dotprod;
2133  float gvec[3], bvec[3];
2134 
2135  num_gridlines = binl1c->num_gridlines;
2136  nbinx = binl1c->nbinx;
2137  num_pixels = l1file->npix;
2138  db = (5.2) / 6371 / 2; // Half of bin size in radians, resolution in km
2139  double db_fudge = db * 1.5;
2140 
2141  // big loop
2142  // dot product
2143  for (pix = 0; pix < num_pixels; pix++) {
2144  gdindex[pix][0] = -1; // actually these are row/col not irow/icol like in my searching algos!!!
2145  gdindex[pix][1] = -1;
2146 
2147  if (lat_new[pix] > 90 || lat_new[pix] < -90 || lon_new[pix] < -180 || lon_new[pix] > 180 ||
2148  lat_new[pix] == BAD_FLT || lon_new[pix] == BAD_FLT) {
2149  if (binl1c->verbose)
2150  cout
2151  << "ERROR in search_l1cgen --latitude longitude pixel out of the boundaries.OR FILLVALUE."
2152  << "latpix>90 or <-90.." << lat_new[pix] << "lonpix>180 or <-180.." << lon_new[pix]
2153  << endl;
2154  flag_out = 110;
2155  return flag_out;
2156  }
2157  bvec[0] = cos(lon_new[pix] * M_PI / 180) * cos(lat_new[pix] * M_PI / 180);
2158  bvec[1] = sin(lon_new[pix] * M_PI / 180) * cos(lat_new[pix] * M_PI / 180);
2159  bvec[2] = sin(lat_new[pix] * M_PI / 180);
2160 
2161  last_dotprod = search_calc_dotprod(binl1c, bvec, last_dotprod_index);
2162  if(last_dotprod < 0)
2163  last_dotprod = 0 - last_dotprod;
2164  double min_dotprod = last_dotprod;
2165  int min_dotprod_index = last_dotprod_index;
2166  bool found_new = false;
2167 
2168  // look right
2169  for (i = last_dotprod_index+1; i < num_gridlines; i++) {
2170  dotprod = search_calc_dotprod(binl1c, bvec, i);
2171  if(dotprod < 0)
2172  dotprod = 0 - dotprod;
2173  if(dotprod > last_dotprod)
2174  break;
2175  if (dotprod < min_dotprod) {
2176  found_new = true;
2177  min_dotprod = dotprod;
2178  min_dotprod_index = i;
2179  last_dotprod = dotprod;
2180  last_dotprod_index = i;
2181  }
2182  } // end lines
2183 
2184  // look left
2185  if(!found_new && last_dotprod_index > 0) {
2186  for (i = last_dotprod_index-1; i >= 0; i--) {
2187  dotprod = search_calc_dotprod(binl1c, bvec, i);
2188  if(dotprod < 0)
2189  dotprod = 0 - dotprod;
2190  if(dotprod > last_dotprod)
2191  break;
2192  if (dotprod < min_dotprod) {
2193  min_dotprod = dotprod;
2194  min_dotprod_index = i;
2195  last_dotprod = dotprod;
2196  last_dotprod_index = i;
2197  }
2198  } // end lines
2199  } // not found_new
2200 
2201  if (min_dotprod <= db_fudge) {
2202  gdindex[pix][0] = min_dotprod_index + 1;
2203  }
2204 
2205  // calc first gridline
2206  dot_firstline = search_calc_dotprod(binl1c, bvec, 0);
2207 
2208  // calc last gridline
2209  dot_lastline = search_calc_dotprod(binl1c, bvec, num_gridlines - 1);
2210 
2211  // for each pixels
2212  if (dot_firstline <= db && dot_lastline > -db) {
2213  // find column
2214  irow = gdindex[pix][0] - 1;
2215 
2216  if (irow < 0) {
2217  if (binl1c->verbose)cout << "ERROR irow<0 in search_l1c..." << irow << "at pix#.." << pix + 1 << endl;
2218  flag_out = 110;
2219  return flag_out;
2220  }
2221 
2222  for (int j = 0; j < nbinx; j++) {
2223  gvec[0] =
2224  cos(binl1c->lon_gd[irow][j] * M_PI / 180) * cos(binl1c->lat_gd[irow][j] * M_PI / 180);
2225  gvec[1] =
2226  sin(binl1c->lon_gd[irow][j] * M_PI / 180) * cos(binl1c->lat_gd[irow][j] * M_PI / 180);
2227  gvec[2] = sin(binl1c->lat_gd[irow][j] * M_PI / 180);
2228 
2229  c1 = bvec[0] - gvec[0];
2230  c2 = bvec[1] - gvec[1];
2231  c3 = bvec[2] - gvec[2];
2232 
2233  bmcm = sqrt(c1 * c1 + c2 * c2 + c3 * c3);
2234  if (bmcm < bm) {
2235  bm = bmcm;
2236  col = j + 1;
2237  }
2238  }
2239 
2240  if (col < 1) {
2241  if (binl1c->verbose)
2242  cout << "ERROR col<1 in search_l1c..." << col << "at pix#.." << pix + 1 << endl;
2243  flag_out = 110;
2244  return flag_out;
2245  }
2246  gdindex[pix][1] = col;
2247 
2248  if (irow == num_gridlines - 1 && lat_new[pix] < 85.0) {
2249  gdindex[pix][0] = -1;
2250  gdindex[pix][1] = -1;
2251  }
2252 
2253  bm = 100;
2254  col = -1;
2255  } else {
2256  gdindex[pix][0] = -1; // actually these are row/col not irow/icol like in my searching algos!!!
2257  gdindex[pix][1] = -1;
2258  }
2259  } // end pixels
2260 
2261  for (pix = 0; pix < num_pixels; pix++) {
2262  if (gdindex[pix][0] < 1 && gdindex[pix][1] < 1) {
2263  badpix++;
2264  }
2265  }
2266 
2267  if (badpix == num_pixels) {
2268  flag_out = 110;
2269  } else
2270  flag_out = 0;
2271 
2272  if (binl1c->verbose && flag_out == 110)
2273  cout << "THIS LINE WILL BE SKIPPED -- NO PIXELS BINNED.............." << endl;
2274 
2275  return flag_out;
2276 
2277 }
2278 
2279 
2280 /*
2281 double search_calc_dotprod(bin_str *binl1c, float bvec[3], int32_t gridline) {
2282  int32_t nbinx = binl1c->nbinx;
2283  float gnvec[3];
2284  float gnvm;
2285 
2286  // normal vectors for L1C rows
2287  if (binl1c->lat_gd[gridline][nbinx - 1] > 90 || binl1c->lat_gd[gridline][nbinx - 1] < -90 ||
2288  binl1c->lon_gd[gridline][nbinx - 1] < -180 || binl1c->lon_gd[gridline][nbinx - 1] > 180) {
2289  if (binl1c->verbose)
2290  cout << "lat lon for L1C GRID out of the boundaries.." << endl;
2291  exit(1);
2292  }
2293  if (binl1c->lat_gd[gridline][0] > 90 || binl1c->lat_gd[gridline][0] < -90 || binl1c->lon_gd[gridline][0] < -180 ||
2294  binl1c->lon_gd[gridline][0] > 180) {
2295  if (binl1c->verbose)
2296  cout << "lat lon for L1C GRID out of the boundaries.." << endl;
2297  exit(1);
2298  }
2299 
2300  gnvec[0] = sin(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2301  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2302  sin(binl1c->lat_gd[gridline][0] * M_PI / 180) -
2303  sin(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2304  sin(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180);
2305  gnvec[1] = sin(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2306  cos(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180) -
2307  cos(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2308  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2309  sin(binl1c->lat_gd[gridline][0] * M_PI / 180);
2310  gnvec[2] = cos(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2311  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2312  sin(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180) -
2313  sin(binl1c->lon_gd[gridline][nbinx - 1] * M_PI / 180) *
2314  cos(binl1c->lat_gd[gridline][nbinx - 1] * M_PI / 180) *
2315  cos(binl1c->lon_gd[gridline][0] * M_PI / 180) * cos(binl1c->lat_gd[gridline][0] * M_PI / 180);
2316 
2317  // vector norm
2318  gnvm = sqrt(gnvec[0] * gnvec[0] + gnvec[1] * gnvec[1] + gnvec[2] * gnvec[2]);
2319  if (isnan(gnvm) == 1) {
2320  if (binl1c->verbose)
2321  cout << "NAN value for gnvm.." << endl;
2322  exit(1);
2323  }
2324  if (gnvm == 0) {
2325  if (binl1c->verbose)
2326  cout << "ERROR gnvm == 0--- WE CANT NORMALIZE..." << endl;
2327  exit(1);
2328  }
2329  // normalization
2330  gnvec[0] = gnvec[0] / gnvm;
2331  gnvec[1] = gnvec[1] / gnvm;
2332  gnvec[2] = gnvec[2] / gnvm;
2333 
2334  // for each pixels
2335  // dot prod, orbital normaliz and transposed by pixel vector
2336  return gnvec[0] * bvec[0] + gnvec[1] * bvec[1] + gnvec[2] * bvec[2];
2337 }
2338 */
2339 
2340 
2341 int search_l1c(filehandle *l1file, l1str *l1rec, bin_str *binl1c, short **gdindex) {
2342  int32_t num_gridlines, nbinx;
2343  float gvec[3], bvec[3];
2344  int flag_out = 0;
2345  size_t pix;
2346  size_t badpix = 0;
2347  int32_t i;
2348  size_t num_pixels;
2349  int irow = -1, col = -1;
2350  float bmcm, bm = 100;
2351  double db;
2352  float c1, c2, c3;
2353  double dotprod, dot_firstline, dot_lastline;
2354 
2355  int32_t last_dotprod_index = 0;
2356  double last_dotprod;
2357 
2358  num_gridlines = binl1c->num_gridlines;
2359  nbinx = binl1c->nbinx;
2360  num_pixels = l1file->npix;
2361 
2362  db = (5.2) / 6371 / 2; // Half of bin size in radians, resolution in km
2363  double db_fudge = db * 1.5;
2364 
2365  // big loop
2366  // dot product
2367  for (pix = 0; pix < num_pixels; pix++) {
2368  gdindex[pix][0] = -1; // actually these are row/col not irow/icol like in my searching algos!!!
2369  gdindex[pix][1] = -1;
2370  flag_out = 0;
2371 
2372  if (l1rec->lat[pix] > 90 || l1rec->lat[pix] < -90 || l1rec->lon[pix] < -180 ||
2373  l1rec->lon[pix] > 180 || l1rec->lat[pix] == BAD_FLT || l1rec->lon[pix] == BAD_FLT) {
2374  if (binl1c->verbose)
2375  cout
2376  << "ERROR in search_l1cgen --latitude longitude pixel out of the boundaries.OR FILLVALUE."
2377  << "latpix>90 or <-90.." << l1rec->lat[pix] << "lonpix>180 or <-180.." << l1rec->lon[pix]
2378  << endl;
2379  binl1c->badgeo++;
2380  flag_out = 110;
2381  return flag_out;
2382  }
2383 
2384  bvec[0] = cos(l1rec->lon[pix] * M_PI / 180) * cos(l1rec->lat[pix] * M_PI / 180);
2385  bvec[1] = sin(l1rec->lon[pix] * M_PI / 180) * cos(l1rec->lat[pix] * M_PI / 180);
2386  bvec[2] = sin(l1rec->lat[pix] * M_PI / 180);
2387 
2388  last_dotprod = search_calc_dotprod(binl1c, bvec, last_dotprod_index);
2389  if(last_dotprod < 0)
2390  last_dotprod = 0 - last_dotprod;
2391  double min_dotprod = last_dotprod;
2392  int min_dotprod_index = last_dotprod_index;
2393  bool found_new = false;
2394 
2395  // look right
2396  for (i = last_dotprod_index+1; i < num_gridlines; i++) {
2397  dotprod = search_calc_dotprod(binl1c, bvec, i);
2398  if(dotprod < 0)
2399  dotprod = 0 - dotprod;
2400  if(dotprod > last_dotprod)
2401  break;
2402  if (dotprod < min_dotprod) {
2403  found_new = true;
2404  min_dotprod = dotprod;
2405  min_dotprod_index = i;
2406  last_dotprod = dotprod;
2407  last_dotprod_index = i;
2408  }
2409  } // end lines
2410 
2411  // look left
2412  if(!found_new && last_dotprod_index > 0) {
2413  for (i = last_dotprod_index-1; i >= 0; i--) {
2414  dotprod = search_calc_dotprod(binl1c, bvec, i);
2415  if(dotprod < 0)
2416  dotprod = 0 - dotprod;
2417  if(dotprod > last_dotprod)
2418  break;
2419  if (dotprod < min_dotprod) {
2420  min_dotprod = dotprod;
2421  min_dotprod_index = i;
2422  last_dotprod = dotprod;
2423  last_dotprod_index = i;
2424  }
2425  } // end lines
2426  } // not found_new
2427 
2428  if (min_dotprod <= db_fudge) {
2429  gdindex[pix][0] = min_dotprod_index + 1;
2430  }
2431 
2432  // calc first gridline
2433  dot_firstline = search_calc_dotprod(binl1c, bvec, 0);
2434 
2435  // calc last gridline
2436  dot_lastline = search_calc_dotprod(binl1c, bvec, num_gridlines - 1);
2437 
2438  // for each pixels
2439  if (dot_firstline <= db && dot_lastline > -db) {
2440  // find column
2441  irow = gdindex[pix][0] - 1;
2442 
2443  if (irow < 0) {
2444  if (binl1c->verbose)
2445  cout << "ERROR irow<0 in search_l1c..." << irow << "at pix#.." << pix + 1 << endl;
2446  flag_out = 110;
2447  return flag_out;
2448  }
2449 
2450  for (int j = 0; j < nbinx; j++) {
2451  gvec[0] =
2452  cos(binl1c->lon_gd[irow][j] * M_PI / 180) * cos(binl1c->lat_gd[irow][j] * M_PI / 180);
2453  gvec[1] =
2454  sin(binl1c->lon_gd[irow][j] * M_PI / 180) * cos(binl1c->lat_gd[irow][j] * M_PI / 180);
2455  gvec[2] = sin(binl1c->lat_gd[irow][j] * M_PI / 180);
2456 
2457  c1 = bvec[0] - gvec[0];
2458  c2 = bvec[1] - gvec[1];
2459  c3 = bvec[2] - gvec[2];
2460 
2461  bmcm = sqrt(c1 * c1 + c2 * c2 + c3 * c3);
2462  if (bmcm < bm) {
2463  bm = bmcm;
2464  col = j + 1;
2465  }
2466  }
2467 
2468  if (col < 1) {
2469  if (binl1c->verbose)
2470  cout << "ERROR col<1 in search_l1c..." << col << "at pix#.." << pix + 1 << endl;
2471  flag_out = 110;
2472  return flag_out;
2473  }
2474  gdindex[pix][1] = col;
2475 
2476  if (irow == num_gridlines - 1 && l1rec->lat[pix] < 85.0) {
2477  gdindex[pix][0] = -1;
2478  gdindex[pix][1] = -1;
2479  }
2480 
2481  bm = 100;
2482  col = -1;
2483  } else {
2484  gdindex[pix][0] = -1; // actually these are row/col not irow/icol like in my searching algos!!!
2485  gdindex[pix][1] = -1;
2486  }
2487  } // end pixels
2488 
2489  for (pix = 0; pix < num_pixels; pix++) {
2490  if (gdindex[pix][0] < 1 && gdindex[pix][1] < 1) {
2491  badpix++;
2492  }
2493  }
2494 
2495  if (badpix == num_pixels) {
2496  flag_out = 110;
2497  } else
2498  flag_out = 0;
2499 
2500  if (binl1c->verbose && flag_out == 110)
2501  cout << "THIS LINE WILL BE SKIPPED -- NO PIXELS BINNED.............." << endl;
2502 
2503  return flag_out;
2504 }
2505 
2506 // fred/richard implemtation
2507 // one pixel at the time----
2508 int search2_l1c(size_t ybins, size_t xbins, float lat, float lon, float **lat_gd, float **lon_gd, short *row,
2509  short *col) {
2510  int32_t num_gridlines, nbinx, ic;
2511  short gdrow = -1, gdcol = -1;
2512  float pvec[3], gnvm, db, pdotgn_first, pdotgn_last;
2513  float **gnvec = nullptr, **cnvec = nullptr, **dc = nullptr;
2514  int flag_out = -1;
2515  int32_t i;
2516  float dotprod, dotprod2, dcm;
2517  flag_out = -1;
2518 
2519  num_gridlines = ybins;
2520  nbinx = xbins;
2521 
2522  db = (5.2) / 6371 / 2; // Half of bin size in radians, resolution in km
2523 
2524  // fred/richard implemtation
2525  gnvec = allocate2d_float(3, num_gridlines);
2526  cnvec = allocate2d_float(3, num_gridlines);
2527  dc = allocate2d_float(3, num_gridlines);
2528 
2529  ic = nbinx / 2;
2530 
2531  for (i = 0; i < num_gridlines; i++) {
2532  // normal vectors for L1C rows
2533  if (lat_gd[i][nbinx - 1] > 90 || lat_gd[i][nbinx - 1] < -90 || lon_gd[i][nbinx - 1] < -180 ||
2534  lon_gd[i][nbinx - 1] > 180) {
2535  // if(binl1c->verbose)cout<<"lat lon out of the boundaries.."<<endl;
2536  exit(1);
2537  }
2538  if (lat_gd[i][0] > 90 || lat_gd[i][0] < -90 || lon_gd[i][0] < -180 || lon_gd[i][0] > 180) {
2539  // if(binl1c->verbose)cout<<"lat lon out of the boundaries.."<<endl;
2540  exit(1);
2541  }
2542  // compute normal vector for each rowgrid----
2543  gnvec[0][i] = sin(lon_gd[i][nbinx - 1] * M_PI / 180) * cos(lat_gd[i][nbinx - 1] * M_PI / 180) *
2544  sin(lat_gd[i][0] * M_PI / 180) -
2545  sin(lat_gd[i][nbinx - 1] * M_PI / 180) * sin(lon_gd[i][0] * M_PI / 180) *
2546  cos(lat_gd[i][0] * M_PI / 180);
2547  gnvec[1][i] = sin(lat_gd[i][nbinx - 1] * M_PI / 180) * cos(lon_gd[i][0] * M_PI / 180) *
2548  cos(lat_gd[i][0] * M_PI / 180) -
2549  cos(lon_gd[i][nbinx - 1] * M_PI / 180) * cos(lat_gd[i][nbinx - 1] * M_PI / 180) *
2550  sin(lat_gd[i][0] * M_PI / 180);
2551  gnvec[2][i] = cos(lon_gd[i][nbinx - 1] * M_PI / 180) * cos(lat_gd[i][nbinx - 1] * M_PI / 180) *
2552  sin(lon_gd[i][0] * M_PI / 180) * cos(lat_gd[i][0] * M_PI / 180) -
2553  sin(lon_gd[i][nbinx - 1] * M_PI / 180) * cos(lat_gd[i][nbinx - 1] * M_PI / 180) *
2554  cos(lon_gd[i][0] * M_PI / 180) * cos(lat_gd[i][0] * M_PI / 180);
2555 
2556  gnvm = sqrt(gnvec[0][i] * gnvec[0][i] + gnvec[1][i] * gnvec[1][i] + gnvec[2][i] * gnvec[2][i]);
2557  if (isnan(gnvm) == 1) {
2558  // if(binl1c->verbose)cout<<"NAN value for gnvm.."<<endl;
2559  exit(1);
2560  }
2561  if (gnvm == 0) {
2562  // if(binl1c->verbose)cout<<"ERROR gnvm == 0--- WE CANT NORMALIZE..."<<endl;
2563  exit(1);
2564  }
2565  // normalization
2566  gnvec[0][i] = gnvec[0][i] / gnvm;
2567  gnvec[1][i] = gnvec[1][i] / gnvm;
2568  gnvec[2][i] = gnvec[2][i] / gnvm;
2569 
2570  // Compute normals to center columns
2571  cnvec[0][i] = sin(lon_gd[i][ic] * M_PI / 180) * cos(lat_gd[i][ic] * M_PI / 180) * gnvec[2][i] -
2572  sin(lat_gd[i][ic] * M_PI / 180) * gnvec[1][i];
2573  cnvec[1][i] = sin(lat_gd[i][ic] * M_PI / 180) * gnvec[0][i] -
2574  cos(lon_gd[i][ic] * M_PI / 180) * cos(lat_gd[i][ic] * M_PI / 180) * gnvec[2][i];
2575  cnvec[2][i] = cos(lon_gd[i][ic] * M_PI / 180) * cos(lat_gd[i][ic] * M_PI / 180) * gnvec[1][i] -
2576  sin(lon_gd[i][ic] * M_PI / 180) * cos(lat_gd[i][ic] * M_PI / 180) * gnvec[0][i];
2577 
2578  //; Compute grid row nadir resolution
2579  dc[0][i] = cos(lon_gd[i][ic + 1] * M_PI / 180) * cos(lat_gd[i][ic + 1] * M_PI / 180) -
2580  cos(lon_gd[i][ic] * M_PI / 180) * cos(lat_gd[i][ic] * M_PI / 180);
2581  dc[1][i] = sin(lon_gd[i][ic + 1] * M_PI / 180) * cos(lat_gd[i][ic + 1] * M_PI / 180) -
2582  sin(lon_gd[i][ic] * M_PI / 180) * cos(lat_gd[i][ic] * M_PI / 180);
2583  dc[2][i] = sin(lat_gd[i][ic + 1] * M_PI / 180) - sin(lat_gd[i][ic] * M_PI / 180);
2584  }
2585 
2586  // dot product
2587  // for(pix=0;pix<num_pixels;pix++){
2588  for (i = 0; i < num_gridlines; i++) {
2589  if (lat > 90 || lat < -90 || lon < -180 || lon > 180) {
2590  // if(binl1c->verbose)cout<<"latitude longitude pixel out of the boundaries.."<<endl;
2591  exit(1);
2592  }
2593 
2594  pvec[0] = cos(lon * M_PI / 180) * cos(lat * M_PI / 180);
2595  pvec[1] = sin(lon * M_PI / 180) * cos(lat * M_PI / 180);
2596  pvec[2] = sin(lat * M_PI / 180);
2597 
2598  dotprod = pvec[0] * gnvec[0][i] + pvec[1] * gnvec[1][i] + pvec[2] * gnvec[2][i];
2599 
2600  if (i == 0) {
2601  pdotgn_first = dotprod; // first gridline
2602  }
2603  if (i == num_gridlines - 1) {
2604  pdotgn_last = dotprod; // last line
2605  }
2606 
2607  if (pdotgn_first <= db && pdotgn_last > -db && gdrow < 0 && gdcol < 0) {
2608  // if(binl1c->verbose) cout<<"this pix is inside the L1C grid.."<<endl;
2609  }
2610 
2611  // check row for pixel j
2612  if (dotprod <= db && dotprod > -db && gdrow < 0) {
2613  gdrow = i;
2614  }
2615 
2616  if (gdrow >= 0) {
2617  dotprod2 = pvec[0] * cnvec[0][gdrow] + pvec[1] * cnvec[1][gdrow] + pvec[2] * cnvec[2][gdrow];
2618  dcm =
2619  sqrt(dc[0][gdrow] * dc[0][gdrow] + dc[1][gdrow] * dc[1][gdrow] + dc[2][gdrow] * dc[2][gdrow]);
2620  gdcol = ic + dotprod2 / dcm + .5;
2621  if (gdrow >= 0 && gdcol >= 0) {
2622  flag_out = 0;
2623  } else {
2624  flag_out = 1;
2625  }
2626  }
2627  }
2628  *row = gdrow;
2629  *col = gdcol;
2630 
2631  gdrow = -1;
2632  gdcol = -1;
2633 
2634  // bmcm only checked one pixel!!I
2635  if (dc != nullptr)
2636  free2d_float(dc);
2637  if (gnvec != nullptr)
2638  free2d_float(gnvec);
2639  if (cnvec != nullptr)
2640  free2d_float(cnvec);
2641 
2642  if (flag_out == 0)
2643  return 0;
2644  else
2645  return 1;
2646 }
2647 
2648 double rot_angle(double senz, double solz, double sena, double suna) {
2649  double rotangle = double(BAD_FLT), term2;
2650  if (senz != BAD_FLT && solz != BAD_FLT && sena != BAD_FLT && suna != BAD_FLT) {
2651  double cosunz = cos(solz * M_PI / 180.);
2652  double cosenz = cos(senz * M_PI / 180.);
2653  double term1 = -1 * cosunz + cosenz * cos(senz * M_PI / 180. + M_PI / 180.);
2654  double cos2 = 0.5 * (cos(2 * senz * M_PI / 180.) + 1);
2655  double term3 =
2656  sqrt(1 - cos(senz * M_PI / 180. + M_PI / 180.) * cos(senz * M_PI / 180. + M_PI / 180.));
2657 
2658  if ((sena - suna) < 2 * 180. && (sena - suna) > 180.)
2659  term2 = sqrt(1 - cos2) * term3;
2660  if ((sena - suna) < 180. && (sena - suna) > 0.)
2661  term2 = -1 * sqrt(1 - cos2) * term3;
2662 
2663  double cos_alpha = term1 / term2;
2664 
2665  rotangle = acos(cos_alpha) * 180. / M_PI; // in degrees
2666  }
2667 
2668  return rotangle;
2669 }
l1file
int32 l1file(int32 sdfid, int32 *nsamp, int32 *nscans, int16 *dtynum)
Definition: l1stat_chk.c:586
bin_str::cloud_type
int cloud_type
Definition: l1c_latlongrid.h:100
value
int32 value
Definition: Granule.c:1235
bin_str::sca_3D
float *** sca_3D
Definition: l1c_latlongrid.h:60
bin_str::lat_gd
float ** lat_gd
Definition: l1c_latlongrid.h:41
search_calc_dotprod
double search_calc_dotprod(bin_str *binl1c, float bvec[3], int32_t gridline)
Definition: l1c_latlongrid.cpp:2058
bin_str::doi
std::string doi
Definition: l1c_latlongrid.h:27
bin_str::alloc_bin
virtual int alloc_bin(bin_str *binstr)
Definition: l1c_latlongrid.cpp:176
allocate3d.h
Utility functions for allocating and freeing three-dimensional arrays of various types.
meta_l1c_altvar
int meta_l1c_altvar(bin_str *binl1c, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:275
j
int j
Definition: decode_rs.h:73
status
int status
Definition: l1_czcs_hdf.c:32
bin_str::bin_str
bin_str()
Definition: l1c_latlongrid.cpp:28
bin_str::~bin_str
virtual ~bin_str()
Definition: l1c_latlongrid.cpp:102
check_err
void check_err(const int stat, const int line, const char *file)
Definition: nc4utils.c:35
bin_str::nrec_3D
short *** nrec_3D
Definition: l1c_latlongrid.h:49
allocate4d_float
float **** allocate4d_float(size_t nr, size_t nz, size_t ny, size_t nx)
Allocate a four-dimensional array of type float of a given size.
Definition: allocate4d.c:56
bin_str::senz_3D
float *** senz_3D
Definition: l1c_latlongrid.h:59
bin_str::nbands
int16_t nbands
Definition: l1c_latlongrid.h:34
rmse_l1c_alt
int rmse_l1c_alt(filehandle *l1file, bin_str *binl1c, l1str *l1rec, short **gdindex)
Definition: l1c_latlongrid.cpp:251
viirs_ler_luts::sb
real, dimension(260) sb
Definition: viirs_ler_luts.f95:12
out
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 out
Definition: HISTORY.txt:422
create_images.sensor
string sensor
Definition: create_images.py:73
search_l1c_parallax
int search_l1c_parallax(filehandle *l1file, float *lat_new, float *lon_new, bin_str *binl1c, short **gdindex)
Definition: l1c_latlongrid.cpp:2119
bin_str::sensor
int sensor
Definition: l1c_latlongrid.h:28
bin_str::I_noise_4D
float **** I_noise_4D
Definition: l1c_latlongrid.h:66
open_l1c_grid
int open_l1c_grid(const char *ifile_l1c, bin_str *binl1c, float **lat_gd, float **lon_gd, float **alt_gd)
Definition: l1c_latlongrid.cpp:2019
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
color_dtdb.v
v
Definition: color_dtdb.py:430
meta_l1c_global
int meta_l1c_global(char *gridname, bin_str *binl1c, int16_t num_gridlines, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:1077
unix2ymds
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
hico.proc_hico.valid_min
valid_min
Definition: proc_hico.py:246
parallax
int parallax(filehandle *l1file, const char *l1c_anc, const char *l1c_grid, l1str *l1rec, bin_str *binl1c, short **gdindex, NcFile *nc_output, int32_t sline, int firstcall)
Definition: l1c_latlongrid.cpp:1513
syear
int syear
Definition: l1_czcs_hdf.c:15
bin_str::nviews
int16_t nviews
Definition: l1c_latlongrid.h:33
file_format
Definition: filetype.h:69
hico.proc_hico.valid_max
valid_max
Definition: proc_hico.py:247
bin_str::suna_3D
float *** suna_3D
Definition: l1c_latlongrid.h:56
string
@ string
Definition: GEO_read_param_file.c:64
hico.proc_hico.solz
solz
Definition: proc_hico.py:240
free2d_short
void free2d_short(short **p)
Free a two-dimensional array created by allocate2d_short.
Definition: allocate2d.c:96
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
meta_l1c_bin
int meta_l1c_bin(filehandle *l1file, bin_str *binl1c, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:326
bin_str::i_diff2
float **** i_diff2
Definition: l1c_latlongrid.h:65
sday
int sday
Definition: l1_czcs_hdf.c:15
allocate3d_short
short *** allocate3d_short(size_t nz, size_t ny, size_t nx)
Allocate a three-dimensional array of type short of a given size.
Definition: allocate3d.c:13
run_local.lat
lat
Definition: run_local.py:69
rot_angle
double rot_angle(double senz, double solz, double sena, double suna)
Definition: l1c_latlongrid.cpp:2648
check_l1c_time
int check_l1c_time(const char *l1b_file, const char *l1c_grid, bin_str *binl1c)
Definition: l1c_latlongrid.cpp:1880
Re
float Re
Definition: l1c.cpp:57
bin_str::alt_2D
short ** alt_2D
Definition: l1c_latlongrid.h:54
run_local.lon
lon
Definition: run_local.py:69
bin_str::bintype
int bintype
Definition: l1c_latlongrid.h:95
allocate4d.h
Utility functions for allocating and freeing four-dimensional arrays of various types.
bin_str::outlist
char outlist[FILENAME_MAX]
Definition: l1c_latlongrid.h:97
bin_str::badgeo
int32_t badgeo
Definition: l1c_latlongrid.h:87
I
#define I
Definition: outernetDecode.c:21
getFormat
file_format getFormat(char *filename)
Definition: filetype.c:192
hico.proc_hico.sena
sena
Definition: proc_hico.py:241
bin_str::time_l1b
double *** time_l1b
Definition: l1c_latlongrid.h:43
bin_str::fillval1
short fillval1
Definition: l1c_latlongrid.h:81
free3d_short
void free3d_short(short ***p)
Free a three-dimensional array created by allocate3d_short.
Definition: allocate3d.c:39
free2d_float
void free2d_float(float **p)
Free a two-dimensional array created by allocate2d_float.
Definition: allocate2d.c:140
bin_str::tini_l1b
double tini_l1b
Definition: l1c_latlongrid.h:91
FT_OCIS
@ FT_OCIS
Definition: filetype.h:44
term
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to already in place for MODIS TERRA was implemented for MODIS AQUA A time dependent LUT was added which gives coefficients for a detector specific crosstalk correction based on the aggregated Band radiances The Band scaled integers are adjusted by the Band correction term
Definition: HISTORY.txt:439
bin_str::rot_angle
float *** rot_angle
Definition: l1c_latlongrid.h:61
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
bin_str::nrec_3D_view
short *** nrec_3D_view
Definition: l1c_latlongrid.h:50
bin_str::fillval2
float fillval2
Definition: l1c_latlongrid.h:82
free4d_float
void free4d_float(float ****p)
Free a four-dimensional array created by allocate4d_float.
Definition: allocate4d.c:93
generalauxtype::double
integer, parameter double
Definition: GeneralAuxType.f90:27
bin_str
Definition: l1c_latlongrid.h:15
color_dtdb.result
result
Definition: color_dtdb.py:197
ncattredit.history
string history
Definition: ncattredit.py:30
bin_str::close_bin
virtual int close_bin(bin_str *binstr)
Definition: l1c_latlongrid.cpp:105
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
allocate2d.h
Utility functions for allocating and freeing two-dimensional arrays of various types.
bin_str::tend_l1b
double tend_l1b
Definition: l1c_latlongrid.h:92
bin_str::alt_diff2
short ** alt_diff2
Definition: l1c_latlongrid.h:55
FT_SPEXONE
@ FT_SPEXONE
Definition: filetype.h:61
bin_str::num_gridlines
int16_t num_gridlines
Definition: l1c_latlongrid.h:31
bin_str::nrec_4D_band
float **** nrec_4D_band
Definition: l1c_latlongrid.h:51
free3d_float
void free3d_float(float ***p)
Free a three-dimensional array created by allocate3d_float.
Definition: allocate3d.c:103
bin_str::alt_mean
float ** alt_mean
Definition: l1c_latlongrid.h:52
trunc
double trunc(double)
hico.proc_hico.units
units
Definition: proc_hico.py:245
bin_str::nbinx
int16_t nbinx
Definition: l1c_latlongrid.h:32
dotprod
long dotprod(void *dp, signed short a[])
bin_str::lon_gd
float ** lon_gd
Definition: l1c_latlongrid.h:42
FT_L1C
@ FT_L1C
Definition: filetype.h:65
meta_l1c_grid
int meta_l1c_grid(char *gridname, bin_str *binl1c, int16_t num_gridlines, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:913
FT_HARP2
@ FT_HARP2
Definition: filetype.h:62
l1c_latlongrid.h
FT_HKT
@ FT_HKT
Definition: filetype.h:64
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
bin_str::alt
float ** alt
Definition: l1c_latlongrid.h:45
bin_str::time_gd
double * time_gd
Definition: l1c_latlongrid.h:44
NBANDS
@ NBANDS
Definition: make_L3_v1.1.c:53
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
meta_l1c_full
int meta_l1c_full(filehandle *l1file, bin_str *binl1c, const char *l1c_grid, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:529
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
search_l1c
int search_l1c(filehandle *l1file, l1str *l1rec, bin_str *binl1c, short **gdindex)
Definition: l1c_latlongrid.cpp:2341
bin_str::cloudem_flag
int cloudem_flag
Definition: l1c_latlongrid.h:99
bin_str::nrec_2D
short ** nrec_2D
Definition: l1c_latlongrid.h:48
free3d_double
void free3d_double(double ***p)
Free a three-dimensional array created by allocate3d_double.
Definition: allocate3d.c:135
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
ymds2unix
double ymds2unix(short year, short month, short day, double secs)
allocate3d_double
double *** allocate3d_double(size_t nz, size_t ny, size_t nx)
Allocate a three-dimensional array of type double of a given size.
Definition: allocate3d.c:109
open_l1c
int open_l1c(const char *ifile_l1c, size_t *ybins, size_t *xbins, float **lat_gd, float **lon_gd)
Definition: l1c_latlongrid.cpp:1988
bin_str::l1c_anc
char l1c_anc[FILENAME_MAX]
Definition: l1c_latlongrid.h:98
bin_str::fillval3
double fillval3
Definition: l1c_latlongrid.h:83
FT_OCIL1B
@ FT_OCIL1B
Definition: filetype.h:43
unix2ymdhms
void unix2ymdhms(double usec, int16_t *year, int16_t *mon, int16_t *day, int16_t *hour, int16_t *min, double *sec)
Definition: unix2ymdhms.c:8
bin_str::outpix
int32_t outpix
Definition: l1c_latlongrid.h:86
tilt
int16_t * tilt
Definition: l2bin.cpp:80
i
int i
Definition: decode_rs.h:71
bin_str::I_4D
float **** I_4D
Definition: l1c_latlongrid.h:64
bintime_l1c
int bintime_l1c(filehandle *l1file, l1str *l1rec, bin_str *binstr, short **gdindex, double scantime, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:1447
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
bin_str::tend_l1c
double tend_l1c
Definition: l1c_latlongrid.h:90
bin_str::inpix
int32_t inpix
Definition: l1c_latlongrid.h:85
bin_str::verbose
bool verbose
Definition: l1c_latlongrid.h:37
verbose
int verbose
Definition: fmt_check.c:6
npix
int npix
Definition: get_cmp.c:28
bin_str::sena_3D
float *** sena_3D
Definition: l1c_latlongrid.h:58
bin_str::history
std::string history
Definition: l1c_latlongrid.h:25
str
#define str(s)
Definition: leveltwo2csv-main.c:19
MDN.setup.version
version
Definition: setup.py:15
search2_l1c
int search2_l1c(size_t ybins, size_t xbins, float lat, float lon, float **lat_gd, float **lon_gd, short *row, short *col)
Definition: l1c_latlongrid.cpp:2508
bin_str::tini_l1c
double tini_l1c
Definition: l1c_latlongrid.h:89
mapgen_overlay.crs
def crs
Definition: mapgen_overlay.py:176
bin_str::alt_rmse
float ** alt_rmse
Definition: l1c_latlongrid.h:53
bin_l1c
int bin_l1c(filehandle *l1file, l1str *l1rec, bin_str *binl1c, short **gdindex, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:1666
bin_str::full_l1cgrid
std::string full_l1cgrid
Definition: l1c_latlongrid.h:78
isodate2unix
double isodate2unix(const char *isodate)
Definition: unix2isodate.c:61
bin_str::sunz_3D
float *** sunz_3D
Definition: l1c_latlongrid.h:57
bin_str::pversion
std::string pversion
Definition: l1c_latlongrid.h:26
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