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.cpp
Go to the documentation of this file.
1 
2 // Created by Martin Montes on 10/26/20.
3 // latest update 6/28/2023
4 //****************************************************8
5 
6 #include <filehandle.h>
7 #include <l1.h>
8 #include "l1c.h"
9 #include "l1c_filehandle.h"
10 #include "l1c_input.h"
11 #include "l1c_str.h"
12 #include "l2_str.h"
13 #include <string>
14 #include <stdlib.h>
15 #include <stdio.h>
16 #include <math.h>
17 #include <list>
18 #include <stack>
19 #include <vector>
20 #include <iostream>
21 #include <fstream>
22 #include <sstream>
23 #include <netcdf.h>
24 #include <nc4utils.h>
25 #include "hawkeye_methods.h"
26 #include "allocate4d.h"
27 #include <allocate3d.h>
28 #include <allocate2d.h>
29 #include <algorithm>
30 // #include <bits/stdc++.h>
31 #include <cstdlib>
32 #include <ctime>
33 #include <chrono>
34 #include <GeographicLib/Geodesic.hpp>
35 #include <GeographicLib/GeodesicLine.hpp>
36 #include <GeographicLib/Constants.hpp>
37 #include <GeographicLib/Rhumb.hpp>
38 #include <timeutils.h>
39 #include <genutils.h>
40 #include <stdint.h>
41 #include "l2prod.h"
42 #include "Strassen.h"
43 #include <libnav.h>
44 
45 static size_t num_scans, num_pixels;
46 static size_t num_blue_bands, num_red_bands, num_SWIR_bands;
47 
48 // Navigation data
49 static short **senz = nullptr, **sena = nullptr;
50 
51 // static float *orbp;
52 // static float *attang;
53 // static double *timeArr;
54 
55 // L1C vars
56 static int16_t nswath;
57 float Re = 6378.137;
58 using namespace std;
59 using namespace boost::assign;
60 using namespace std::chrono;
61 using namespace GeographicLib;
62 
63 using namespace netCDF;
64 using namespace netCDF::exceptions;
65 
66 namespace bg = boost::geometry;
67 typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree>> Point_t;
68 typedef bg::model::polygon<Point_t> Polygon_t;
69 typedef bg::model::box<Point_t> Box_t;
70 
71 namespace l1c {
72 
73 L1C::L1C() {
74 }
75 
76 L1C::~L1C() {
77 }
78 
79 int L1C::add_proc_group_l1c(L1C_input *l1cinput,l1c_filehandle *l1cfile,const char* l1c_filename)
80 {
81 //additional attributes
82 
83  int ncid_out,grp_pc,grp_ip;
84  string str,str2;
85  int status = nc_open(l1c_filename, NC_WRITE, &ncid_out);
86  if ((status = nc_def_grp(ncid_out, "processing_control", &grp_pc))) // netcdf-4
87  check_err(status, __LINE__, __FILE__);
88  str="software_name";
89  str2=l1cfile->progname;
90  status = nc_put_att_text(grp_pc,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
91  str="software_version";
92  str2=l1cfile->version;
93  str2=str2.substr(0,5);
94  status = nc_put_att_text(grp_pc,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
95 
96  if ((status = nc_def_grp(grp_pc, "input_parameters", &grp_ip))) // netcdf-4
97  check_err(status, __LINE__, __FILE__);
98  str="ifile";
99  str2="";
100  int nfiles=l1cfile->ifiles.size();
101  for(int f=0;f<nfiles;f++)
102  {
103  if(f==nfiles-1)str2=str2+l1cfile->ifiles[f];else str2=str2+l1cfile->ifiles[f]+",";
104  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
105  }
106  str="ofile";
107  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),strlen(l1c_filename),l1c_filename);
108  str="outlist";
109  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),strlen(l1cinput->outlist),l1cinput->outlist);
110  str="l1c_grid";
111  str2=string(l1cinput->l1c_grid);
112  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
113  str="l1c_anc";
114  str2=string(l1cinput->l1c_anc);
115  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
116  str="demfile";
117  str2="$OCDATAROOT/common/gebco_ocssw_v2020.nc";
118  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
119 
120  str="verbose";
121  int value=(int)l1cinput->verbose;
122  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
123  str="mode";
124  value=l1cinput->l1c_pflag;
125  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
126  str="north";
127  float value2=l1cinput->north;
128  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
129  str="south";
130  value2=l1cinput->south;
131  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
132  str="east";
133  value2=l1cinput->east;
134  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
135  str="west";
136  value2=l1cinput->west;
137  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
138  str="selgran";
139  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,10,l1cinput->selgran);
140  str="selday";
141  value=l1cinput->selday;
142  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
143  str="selmon";
144  value=l1cinput->selmon;
145  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
146  str="selyear";
147  value=l1cinput->selyear;
148  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
149  str="gransize";
150  value=l1cinput->gransize;
151  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
152  str="grantype";
153  value=l1cinput->grantype;
154  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
155  str="bintype";
156  value=l1cinput->bintype;
157  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
158  str="start_timeflag";
159  value=l1cinput->start_timeflag;
160  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
161  str="swath_num";
162  value=l1cinput->swath_num;
163  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
164  str="grid_resolution";
165  value2=l1cinput->grid_resolution;
166  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
167  str="sensor";
168  value=l1cinput->sensor;
169  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
170  str="l2prod";
171  str2=l1cinput->l2prod;
172  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
173  str="ix_l1cprod";
174  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,3,l1cinput->ix_l1cprod);
175  str="start_time";
176  str2=l1cinput->start_time;
177  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
178  str="end_time";
179  str2=l1cinput->end_time;
180  status = nc_put_att_text(grp_ip,NC_GLOBAL,str.c_str(),str2.size(),str2.c_str());
181  str="projection";
182  value=l1cinput->projection;
183  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
184  str="sort_method";
185  value=l1cinput->sort_method;
186  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
187  str="demcloud_flag";
188  value=l1cinput->demcloud_flag;
189  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
190  str="terrain_correct";
191  value=(int)l1cinput->terrain_correct;
192  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
193  str="cloud_height";
194  value=l1cinput->cloud_height;
195  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
196  str="cloud_correct";
197  value=l1cinput->cloud_correct;
198  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
199  str="cloud_type";
200  value=l1cinput->cloud_type;
201  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
202  str="overlap_vflag";
203  value=(int)l1cinput->overlap_vflag;
204  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
205  str="overlap_pflag";
206  value=(int)l1cinput->overlap_pflag;
207  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
208  str="overlap_bflag";
209  value=(int)l1cinput->overlap_bflag;
210  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
211  str="unc_meth";
212  value=l1cinput->unc_meth;
213  status = nc_put_att_int(grp_ip,NC_GLOBAL,str.c_str(),NC_INT,1,&value);
214  str="unc_thres_v";
215  value2=l1cinput->unc_thres_v;
216  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
217  str="unc_thres_p";
218  value2=l1cinput->unc_thres_p;
219  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
220  str="unc_thres_b";
221  value2=l1cinput->unc_thres_b;
222  status = nc_put_att_float(grp_ip,NC_GLOBAL,str.c_str(),NC_FLOAT,1,&value2);
223  status = nc_close(ncid_out);
224 
225 return status;
226 }
227 
228 // function to calculate cross product of two vectors, 3 components each vector
229 void cross_product_double(double vector_a[], double vector_b[], double temp[]) {
230  temp[0] = vector_a[1] * vector_b[2] - vector_a[2] * vector_b[1];
231  temp[1] = -(vector_a[0] * vector_b[2] - vector_a[2] * vector_b[0]);
232  temp[2] = vector_a[0] * vector_b[1] - vector_a[1] * vector_b[0];
233 }
234 
235 double cross_product_norm_double(double vector_a[], double vector_b[]) {
236  double temp[3], nvec;
237  temp[0] = vector_a[1] * vector_b[2] - vector_a[2] * vector_b[1];
238  temp[1] = -(vector_a[0] * vector_b[2] - vector_a[2] * vector_b[0]);
239  temp[2] = vector_a[0] * vector_b[1] - vector_a[1] * vector_b[0];
240 
241  nvec = sqrt(temp[0] * temp[0] + temp[1] * temp[1] + temp[2] * temp[2]);
242  return nvec;
243 }
244 
245 int32_t L1C::l1_cloud_correct(L1C_input* l1cinput, l1c_filehandle* l1cfile) {
246  double *ptime = nullptr, *ptime_l1c = nullptr;
247  const char* l1c_gran;
248  size_t num_ybin, num_xbin;
249  float **latpix = nullptr, **lonpix = nullptr, **l1clat = nullptr, **l1clon = nullptr, **l1calt = nullptr;
250  float **lat_new = nullptr, **lon_new = nullptr;
251  float Rpole = 6356;
252  float Xcorr,Ycorr, Zcorr;
253  float Rsurf, radius_ratio = Re / Rpole;
254  float latcorr;
255  size_t gran = 0;
256  float dv, Hsat = 676.5; // sensor height in meters above ellipsoid
257  short **senapix = nullptr, **senzpix = nullptr;
258  float **cth = nullptr, temp;
259  NcFile *nc_l1c, *nc_l1c_new, *nc_l12_new, *nc_l12;
260  string name;
261  double FillValue2=BAD_FLT;
262  float FillValue=BAD_FLT;
263 
264  if (l1cinput->cloud_correct == 0) {
265  cout << "cloud_correct is 0,thus no paralallax correction is needed!!..use 1 (constant cloud height) "
266  "or 2 (variable cloud height)....."
267  << endl;
268  exit(1);
269  } else if (l1cinput->cloud_correct == 1) {
270  cout << "l1cinput->cloud_correct is: " << l1cinput->cloud_correct
271  << ".LIC->(L1C paralallax corrected)...cloud height in km is constant throughout the L1C grid "
272  "but different from zero so parallax correction id needed..."
273  << endl;
274  if (l1cinput->cloud_height == 0) {
275  cout << "ERROR cloud_height is NOT greater than 0..or was not provided as input:"
276  << l1cinput->cloud_height << endl;
277  exit(1);
278  }
279 
280  //--- OPEN L1C FILES ---
281  for (unsigned int i = 0; i < l1cinput->files_l1c.size(); i++) {
282  gran++;
283  string l1c_str =
284  l1cinput->files_l1c[i]; // ifile is used here for old grid with no parallax correction, this
285  // is a list of l1c grid granules or full grid L1C
286  l1c_gran = l1c_str.c_str();
287 
288  cout << "Opening L1C granule .#.." << i + 1 << "...." << l1c_str
289  << "for parallax-cloud corrections......." << endl;
290 
291  try {
292  nc_l1c = new NcFile(l1c_gran, NcFile::read);
293  } catch (NcException& e) {
294  e.what();
295  cerr << "l1cgen l1c_pflag= 8:: Failure write FULL L1C grid: " + l1c_str << endl;
296  exit(1);
297  }
298 
299  NcDim yd = nc_l1c->getDim("bins_along_track");
300  num_ybin = yd.getSize();
301  NcDim xd = nc_l1c->getDim("bins_across_track");
302  num_xbin = xd.getSize();
303 
304  ptime_l1c = (double*)calloc(num_ybin, sizeof(double));
305  l1clat = allocate2d_float(num_ybin, num_xbin);
306  l1clon = allocate2d_float(num_ybin, num_xbin);
307  l1calt = allocate2d_float(num_ybin, num_xbin);
308 
309  NcGroup ba_grp = nc_l1c->getGroup("bin_attributes");
310  NcVar v1 = ba_grp.getVar("nadir_view_time");
311  v1.getVar(ptime_l1c);
312 
313  NcGroup geo_grp = nc_l1c->getGroup("geolocation_data");
314  v1 = geo_grp.getVar("latitude");
315  v1.getVar(&l1clat[0][0]);
316  v1 = geo_grp.getVar("longitude");
317  v1.getVar(&l1clon[0][0]);
318  v1 = geo_grp.getVar("altitude");
319  v1.getVar(&l1calt[0][0]);
320 
321  cout << "num along bins L1C..." << num_ybin << "num_across bins L1C...." << num_xbin
322  << "l1c time_ini.." << ptime_l1c[0] << "l1c time end.." << ptime_l1c[num_ybin - 1] << endl;
323 
324  lat_new = allocate2d_float(num_ybin, num_xbin);
325  lon_new = allocate2d_float(num_ybin, num_xbin);
326 
327  for (size_t i = 0; i < num_ybin; i++) {
328  for (size_t j = 0; j < num_xbin; j++) {
329  // Rsurf is Earth radius at pixel positions
330  Rsurf = Re / (sqrt(cos(l1clat[i][j] * M_PI / 180) * cos(l1clat[i][j] * M_PI / 180) +
331  radius_ratio * radius_ratio * sin(l1clat[i][j] * M_PI / 180) *
332  sin(l1clat[i][j] * M_PI / 180)));
333  Xcorr = (l1cinput->cloud_height + Rsurf) * cos(l1clat[i][j] * M_PI / 180) *
334  sin(l1clon[i][j] * M_PI / 180); // lon speherical and geodetic should be equal
335  Ycorr = (l1cinput->cloud_height + Rsurf) * sin(l1clat[i][j] * M_PI / 180);
336  Zcorr = (l1cinput->cloud_height + Rsurf) * cos(l1clat[i][j] * M_PI / 180) *
337  cos(l1clon[i][j] * M_PI / 180);
338  // convert back to latitude and longitude
339  latcorr = atan(Ycorr / sqrt(Xcorr * Xcorr + Zcorr * Zcorr));
340  lat_new[i][j] = atan(tan(latcorr) / radius_ratio * radius_ratio) * 180.0 / M_PI;
341  lon_new[i][j] = atan2(Xcorr, Zcorr) * 180.0 / M_PI;
342  lon_new[i][j] = atan2(Ycorr, Xcorr) * 180 / M_PI;
343 
344  if(l1cinput->verbose) cout<<"l1cinput->cloud_height.."<<l1cinput->cloud_height<<"latnew.."<<lat_new[i][j]<<"lon_new.."<<lon_new[i][j]<<"l1clat.."<<l1clat[i][j]<<"l1clon.."<<l1clon[i][j]<<endl;
345  }
346  }
347 
348  // write the new L1C gridor grids
349  string gran_str = to_string(gran);
350  string cth_str = to_string(l1cinput->cloud_height);
351  string l1c_gran_new = l1c_str.substr(0, 24) + "_CTH" + cth_str + ".gran" + gran_str + ".nc";
352  if(l1cinput->verbose)cout << "corrected L1C grid" << l1c_gran_new << " granule #" << gran_str << endl;
353 
354  const char* l1c_gran = l1c_gran_new.c_str();
355 
356  try {
357  nc_l1c_new = new NcFile(l1c_gran, NcFile::replace);
358  } catch (NcException& e) {
359  e.what();
360  cerr << "l1cgen l1c_pflag= 8:: Failure write FULL L1C grid: " + l1c_gran_new << endl;
361  exit(1);
362  }
363  // global attributes--
364  char* gridchar = strdup(l1c_gran);
365 
366  bin_str binl1c;
367  meta_l1c_grid(gridchar, &binl1c, num_ybin, nc_l1c_new);
368 
369  nc_l1c_new->putAtt("processing_version", l1cinput->pversion);
370  nc_l1c_new->putAtt("history", l1cinput->history);
371  string name;
372  nc_l1c_new->putAtt("product_name", l1c_gran_new);
373  NcGroupAtt i1 = nc_l1c->getAtt("startDirection"); // NcGroupAtt is a global attr!!
374  i1.getValues(name);
375  nc_l1c_new->putAtt("startDirection", name);
376  i1 = nc_l1c->getAtt("endDirection");
377  i1.getValues(name);
378  nc_l1c_new->putAtt("endDirection", name);
379  i1 = nc_l1c->getAtt("time_coverage_start");
380  i1.getValues(name);
381  nc_l1c_new->putAtt("time_coverage_start", name);
382  i1 = nc_l1c->getAtt("time_coverage_end");
383  i1.getValues(name);
384  nc_l1c_new->putAtt("time_coverage_end", name);
385  i1 = nc_l1c->getAtt("date_created");
386  i1.getValues(name);
387  nc_l1c_new->putAtt("date_created", name);
388 
389  ba_grp = nc_l1c_new->getGroup("bin_attributes");
390  v1 = ba_grp.getVar("nadir_view_time");
391  v1.putVar(&ptime_l1c[0]);
392  geo_grp = nc_l1c_new->getGroup("geolocation_data");
393  v1 = geo_grp.getVar("latitude");
394  v1.putVar(&lat_new[0][0]);
395  v1 = geo_grp.getVar("longitude");
396  v1.putVar(&lon_new[0][0]);
397  v1 = geo_grp.getVar("altitude");
398  v1.putVar(&l1calt[0][0]);
399 
400  nc_l1c->close();
401  nc_l1c_new->close();
402 
403  free (lat_new);
404  free (lon_new);
405  free (ptime_l1c);
406  free (l1clat);
407  free (l1clon);
408  free (l1calt);
409  }
410  } else if (l1cinput->cloud_correct == 2) {
411  if(l1cinput->verbose)cout << "l1cinput->cloud_correct :" << l1cinput->cloud_correct
412  << "LIB->(L1B paralallax corrected)----cloud_height is constant between pixels......" << endl;
413 
414  if (l1cinput->cloud_height == 0) {
415  cout << "ERROR cloud_height is NOT greater than 0..or was not provided as input:"
416  << l1cinput->cloud_height << endl;
417  exit(1);
418  }
419 
420  for (unsigned int i = 0; i < l1cinput->files.size(); i++) {
421  gran++;
422  string l12_str = l1cinput->files[0];
423  const char* l_12 = l12_str.c_str();
424 
425  if(l1cinput->verbose)cout << "Opening L1B/L2 granule ..." << l1cinput->files[0]
426  << "....for parallax-cloud corrections......." << endl;
427 
428  try {
429  nc_l12 = new NcFile(l_12, NcFile::read);
430  } catch (NcException& e) {
431  e.what();
432  cerr << "l1cgen l1c_pflag= 8:: Failure write FULL L1C grid: " + l12_str << endl;
433  exit(1);
434  }
435 
436  // DIMENSIONS
437  NcDim yd = nc_l12->getDim("number_of_scans");
438  NcDim xd = nc_l12->getDim("ccd_pixels");
439  num_scans = yd.getSize();
440  num_pixels = xd.getSize();
441 
442  lat_new = allocate2d_float(num_scans, num_pixels);
443  lon_new = allocate2d_float(num_scans, num_pixels);
444 
445  if(l1cinput->verbose)cout << "number of scans.." << num_scans << "num of pixels.." << num_pixels << endl;
446 
447  latpix = allocate2d_float(num_scans, num_pixels);
448  lonpix = allocate2d_float(num_scans, num_pixels);
449  ptime = (double*)calloc(num_scans, sizeof(double));
450  senzpix = allocate2d_short(num_scans, num_pixels);
451  senapix = allocate2d_short(num_scans, num_pixels);
452  cth = allocate2d_float(num_scans, num_pixels);
453 
454  for (unsigned int i = 0; i < num_scans; i++) {
455  for (unsigned int j = 0; j < num_pixels; j++) {
456  cth[i][j] = l1cinput->cloud_height;
457  }
458  }
459 
460  NcGroup ba_grp = nc_l12->getGroup("scan_line_attributes");
461  NcVar v1 = ba_grp.getVar("time");
462  v1.getVar(ptime);
463  NcGroup geo_grp = nc_l12->getGroup("geolocation_data");
464  v1 = geo_grp.getVar("latitude");
465  v1.getVar(&latpix[0][0]);
466  v1 = geo_grp.getVar("longitude");
467  v1.getVar(&lonpix[0][0]);
468  float fillval_geo, scale_factor;
469  short senz_min, senz_max, sena_min, sena_max;
470  NcVarAtt a1 = v1.getAtt("_FillValue"); // root group
471  a1.getValues(&fillval_geo);
472  v1 = geo_grp.getVar("sensor_azimuth_angle");
473  v1.getVar(&senapix[0][0]);
474  a1 = v1.getAtt("valid_min"); // root group
475  a1.getValues(&sena_min);
476  a1 = v1.getAtt("valid_max"); // root group
477  a1.getValues(&sena_max);
478  v1 = geo_grp.getVar("sensor_zenith_angle");
479  v1.getVar(&senzpix[0][0]);
480  a1 = v1.getAtt("valid_min"); // root group
481  a1.getValues(&senz_min);
482  a1 = v1.getAtt("valid_max"); // root group
483  a1.getValues(&senz_max);
484  a1 = v1.getAtt("scale_factor"); // root group
485  a1.getValues(&scale_factor);
486 
487  senz_min *= scale_factor;
488  senz_max *= scale_factor;
489  sena_min *= scale_factor;
490  sena_max *= scale_factor;
491 
492  if(l1cinput->verbose)cout << "computing parallax for each pixel.............." << endl;
493 
494  for (unsigned int i = 0; i < num_scans; i++) {
495  for (unsigned int j = 0; j < num_pixels; j++) {
496  // Displacement vector -- wang et al. 2011 ----
497  temp = senapix[i][j] * scale_factor;
498  if (senzpix[i][j] * scale_factor >= senz_min &&
499  senzpix[i][j] * scale_factor <= senz_max && temp >= sena_min && temp <= sena_max) {
500  dv = Hsat * cth[i][j] * tan(senzpix[i][j] * scale_factor * M_PI / 180) /
501  (Hsat - cth[i][j]);
502  lat_new[i][j] = latpix[i][j] * M_PI / 180 - dv * cos(temp * M_PI / 180 + M_PI) / Re;
503  lon_new[i][j] =
504  lonpix[i][j] * M_PI / 180 -
505  dv * sin((temp * M_PI / 180 + M_PI)) / (Re * cos(lat_new[i][j] * M_PI / 180));
506  if (lon_new[i][j] * 180 / M_PI > 180) {
507  lon_new[i][j] -= 2 * M_PI;
508  }
509  if (lon_new[i][j] * 180 / M_PI < -180) {
510  lon_new[i][j] += 2 * M_PI;
511  }
512  lat_new[i][j] *= 180 / M_PI;
513  lon_new[i][j] *= 180 / M_PI;
514  } else {
515  lat_new[i][j] = fillval_geo;
516  lon_new[i][j] = fillval_geo;
517  }
518  }
519  }
520 
521  // write the new L1B with parallax correction
522  string gran_str = to_string(gran);
523  string cth_str = to_string(l1cinput->cloud_height);
524  string l12_gran_new =
525  l12_str.substr(0, 24) + "_CTH" + cth_str + ".L1B" + ".gran" + gran_str + ".nc";
526 
527  if(l1cinput->verbose)cout << "corrected L1B grid" << l12_gran_new << " granule #" << gran_str << endl;
528 
529  const char* l12_gran = l12_gran_new.c_str();
530 
531  try {
532  nc_l12_new = new NcFile(l12_gran, NcFile::replace);
533  } catch (NcException& e) {
534  e.what();
535  cerr << "l1cgen mode= 9:: Failure write FULL L1B parallax corrected: " + l12_gran_new << endl;
536  exit(1);
537  }
538  // global attributes--
539  nc_l12_new->putAtt("title", "PACE OCIS Level-1B Data");
540  nc_l12_new->putAtt("instrument", "OCIS");
541  nc_l12_new->putAtt("processing_version", l1cinput->pversion);
542  nc_l12_new->putAtt("Conventions", "CF-1.8 ACDD-1.3");
543  nc_l12_new->putAtt("institution",
544  "NASA Goddard Space Flight Center, Ocean Biology Processing Group");
545  nc_l12_new->putAtt(
546  "license",
547  "https://science.nasa.gov/earth-science/earth-science-data/data-information-policy/");
548  nc_l12_new->putAtt("naming_authority", "gov.nasa.gsfc.sci.oceancolor");
549  nc_l12_new->putAtt("keywords_vocabulary",
550  "NASA Global Change Master Directory (GCMD) Science Keywords");
551  nc_l12_new->putAtt("stdname_vocabulary", "NetCDF Climate and Forecast (CF) Metadata Convention");
552  nc_l12_new->putAtt("creator_name", "NASA/GSFC");
553  nc_l12_new->putAtt("creator_email", "data@oceancolor.gsfc.nasa.gov");
554  nc_l12_new->putAtt("creator_url", "https://oceancolor.gsfc.nasa.gov");
555  nc_l12_new->putAtt("project", "PACE Project");
556  nc_l12_new->putAtt("publisher_name", "NASA/GSFC");
557  nc_l12_new->putAtt("publisher_email", "data@oceancolor.gsfc.nasa.gov");
558  nc_l12_new->putAtt("publisher_url", "https://oceancolor.gsfc.nasa.gov");
559  nc_l12_new->putAtt("processing_level", "L1B parallax");
560  nc_l12_new->putAtt("cdm_data_type", "swath");
561  nc_l12_new->putAtt("normalizedLt", "0LL");
562  nc_l12_new->putAtt("history", l1cinput->history);
563  nc_l12_new->putAtt("product_name", l12_gran_new);
564  NcGroupAtt i1 = nc_l12->getAtt("startDirection"); // NcGroupAtt is a global attr!!
565  i1.getValues(name);
566  nc_l12_new->putAtt("startDirection", name);
567  i1 = nc_l12->getAtt("endDirection");
568  i1.getValues(name);
569  nc_l12_new->putAtt("endDirection", name);
570  nc_l12_new->putAtt("sample_offset", "0LL");
571  i1 = nc_l12->getAtt("time_coverage_start");
572  string timeStartStr;
573  i1.getValues(timeStartStr);
574  nc_l12_new->putAtt("time_coverage_start", timeStartStr);
575  i1 = nc_l12->getAtt("time_coverage_end");
576  i1.getValues(name);
577  nc_l12_new->putAtt("time_coverage_end", name);
578  i1 = nc_l12->getAtt("date_created");
579  i1.getValues(name);
580  nc_l12_new->putAtt("date_created", name);
581  nc_l12_new->putAtt("earth_sun_distance_correction",ncFloat,FillValue);
582 
583  NcDim ydim = nc_l12_new->addDim("number_of_scans", num_scans);
584  NcDim xdim = nc_l12_new->addDim("ccd_pixels", num_pixels);
585 
586  std::vector<NcDim> dimvec2_geo;
587  dimvec2_geo.push_back(ydim);
588  dimvec2_geo.push_back(xdim);
589 
590  nc_l12_new->addDim("SWIR_pixels", num_pixels);
591  num_blue_bands = 120;
592  num_red_bands = 120;
593  num_SWIR_bands = 9;
594  size_t vector_elements = 3;
595  nc_l12_new->addDim("blue_bands", num_blue_bands);
596  nc_l12_new->addDim("red_bands", num_red_bands);
597  nc_l12_new->addDim("SWIR_bands", num_SWIR_bands);
598  nc_l12_new->addDim("vector_elements", vector_elements);
599 
600  ba_grp = nc_l12_new->addGroup("bin_attributes");
601  geo_grp = nc_l12_new->addGroup("geolocation_data");
602  v1 = ba_grp.addVar("nadir_view_time", ncDouble, ydim);
603 
604  string longName = "Time bin was viewed at nadir view";
605  v1.putAtt("long_name", longName);
606  string units = "seconds since " + timeStartStr.substr(0, 10); // just add the date, no time
607  v1.putAtt("units", units);
608  v1.putAtt("_FillValue", ncDouble, FillValue2);
609  double valid_min_d = 0;
610  v1.putAtt("valid_min", ncDouble, valid_min_d);
611  double valid_max_d = 172800;
612  v1.putAtt("valid_max", ncDouble, valid_max_d);
613  v1.putVar(&ptime[0]);
614 
615  v1 = geo_grp.addVar("latitude", ncFloat, dimvec2_geo);
616  longName = "Latitudes of bin locations";
617  v1.putAtt("long_name", longName);
618  units = "degrees_north";
619  v1.putAtt("units", units);
620  v1.putAtt("_FillValue", ncFloat, FillValue);
621  float valid_min = -90;
622  v1.putAtt("valid_min", ncFloat, valid_min);
623  float valid_max = 90;
624  v1.putAtt("valid_max", ncFloat, valid_max);
625  v1.putVar(&lat_new[0][0]);
626 
627  v1 = geo_grp.addVar("longitude", ncFloat, dimvec2_geo);
628  longName = "Longitudes of bin locations";
629  v1.putAtt("long_name", longName);
630  units = "degrees_east";
631  v1.putAtt("units", units);
632  v1.putAtt("_FillValue", ncFloat, FillValue);
633  valid_min = -180;
634  v1.putAtt("valid_min", ncFloat, valid_min);
635  valid_max = 180;
636  v1.putAtt("valid_max", ncFloat, valid_max);
637  v1.putVar(&lon_new[0][0]);
638 
639  nc_l12->close();
640  nc_l12_new->close();
641 
642  free (lat_new);
643  free (lon_new);
644  free (ptime);
645  free (latpix);
646  free (lonpix);
647  free (senapix);
648  free (senzpix);
649  free (cth);
650  }
651  } else {
652  cout << "l1cinput->cloud_correct WAS NOT setup!" << endl;
653  exit(1);
654  }
655 
656  return 0;
657 }
658 
659 
660 
661 int32_t L1C::swtime_swt2_segment(int swt, L1C_input* l1cinput, l1c_filehandle* l1cfile, int32_t norbs,
662  double* tswt, double tcross, double mgv, double* tmgv) {
663  int16_t bina = 0, binb = 0, ngridlines, gd = 0;
664  double tg = 0., mot = 0, t_start = -1, t_end = -1;
665 
666  t_start = tswt[0];
667  t_end = tswt[norbs - 1];
668 
669  mot = ((l1cinput->grid_resolution) * 1000) / mgv; // in seconds
670 
671  ngridlines = l1cfile->num_gridlines;
672 
673  int flag_time = -1;
674 
675  if (tcross > 0.0) {
676  flag_time = 0;
677  if(l1cinput->verbose)cout << "computing time series assuming mean swath velocity..for swath#." << swt << endl;
678 
679  // backward section of swath--before equat
680  tg = tcross - mot / 2;
681  gd = ngridlines / 2 - 1;
682  tmgv[gd] = tg;
683  binb = 1;
684  if(tg<0 && l1cinput->verbose){
685  // cout<<"NEGATIVE TIME swtime_swt2_segment for mean vel tmgv calc "<<endl;
686  }
687  while (tg >= t_start) {
688  binb++;
689  tg -= mot;
690  gd--;
691  tmgv[gd] = tg;
692  if(tg<0 && l1cinput->verbose){
693  // cout<<"NEGATIVE TIME swtime_swt2_segment for mean vel tmgv calc "<<endl;
694  }
695  }
696  // forward section of swath--after equator
697  tg = tcross + mot / 2;
698  gd = ngridlines / 2;
699  tmgv[gd] = tg;
700  bina = 1;
701 
702  while (tg <= t_end) {
703  bina++;
704  tg += mot;
705  gd++;
706  tmgv[gd] = tg;
707  }
708 
709  if(l1cinput->verbose)cout << "bina.." << bina << "binb.." << binb << endl;
710 
711  l1cfile->num_gridlines = binb + bina;
712 
713  if(l1cinput->verbose)cout << "number of L1C gridlines along-track..." << l1cfile->num_gridlines << "for swath #.." << swt
714  << "t_start.." << t_start << "t_end..." << t_end << endl;
715  } // end equat crossing
716  else {
717  if(l1cinput->verbose)cout << "time series not possible for swath #.." << swt << "tcross<0...." << endl;
718  flag_time = 1;
719  }
720 
721  return flag_time;
722 }
723 
724 int32_t L1C::swtime_swt2(int swt, L1C_input* l1cinput, l1c_filehandle* l1cfile, int32_t norbs, double* tswt,
725  double tcross, double mgv, double* tmgv) {
726  int16_t bina = 0, binb = 0, ngridlines, gd = 0;
727  double tg = 0., mot = 0;
728 
729  mot = ((l1cinput->grid_resolution) * 1000) / mgv; // in seconds
730 
731  ngridlines = l1cfile->num_gridlines;
732 
733  int flag_time = -1;
734 
735  if (tcross > 0.0) {
736  flag_time = 0;
737  if(l1cinput->verbose)cout << "computing time series assuming mean swath velocity..for swath#." << swt << endl;
738 
739  tg = tcross - mot / 2;
740  gd = ngridlines / 2 - 1;
741  tmgv[gd] = tg;
742  binb = 1;
743  if(tg<0 && l1cinput->verbose){
744  cout<<"NEGATIVE FIRST TIME BEFORE CROSSING"<<tg<<"in swtime_swt2 for mean vel tmgv calc "<<mgv<<"mot/2 "<<mot/2<<"tcross "<<tcross<<endl;
745  }
746  while (binb < ngridlines / 2) {//&& tg>=(0+mot
747  binb++;
748  tg -= mot;
749  gd--;
750  tmgv[gd] = tg;
751  if(tg<0 && l1cinput->verbose){
752  // cout<<"NEGATIVE TIME "<<tg<<"in swtime_swt2 for mean vel tmgv calc "<<mgv<<"mot/2 "<<mot/2<<"tcross "<<tcross<<endl;
753  }
754 
755  }
756  tg = tcross + mot / 2;
757  gd = ngridlines / 2;
758  tmgv[gd] = tg;
759  bina = 1;
760 
761  while (bina < ngridlines / 2 + 1) {// && tg<=(86400-mot)
762  bina++;
763  tg += mot;
764  gd++;
765  tmgv[gd] = tg;
766  }
767 
768  if(l1cinput->verbose)cout << "bina.." << bina << "binb.." << binb << endl;
769 
770  l1cfile->num_gridlines = binb + bina;
771 
772  if(l1cinput->verbose)cout << "number of L1C gridlines along-track..." << l1cfile->num_gridlines << "for swath #.." << swt
773  << endl;
774  } // end equat crossing
775  else {
776  if(l1cinput->verbose)cout << "time series not possible for swath #.." << swt << "tcross<0...." << endl;
777  flag_time = 1;
778  }
779 
780  return flag_time;
781 }
782 
783 
784 int32_t L1C::ect_swt(l1c_filehandle* l1cfile, int ix1, int ix2, double* tswt_tot, double* latswt_tot,
785  double* lonswt_tot, double* ovel_tot, double* gvel_tot, double* tswt, double* latswt,
786  double* lonswt, float* tcross, float* loncross, double* ovel, double* gvel) {
787  // determine equatorial crossing time--------------
788  double t1 = -1, t2 = -1;
789  int asc_mode = -1;
790  size_t tindex = -1;
791  int scan_brack[2] = {-1, -1};
792  float lat1, lat2, lon1, lon2; // lat,lon defined globally
793  int flag_ect = -1;
794  double sum1 = 0, sum3 = 0;
795  size_t c = 0, norbs;
796 
797  norbs = ix2 - ix1 + 1;
798  int kk = ix1;
799 
800  while (kk <= ix2) {
801  latswt[c] = latswt_tot[kk];
802  lonswt[c] = lonswt_tot[kk];
803  tswt[c] = tswt_tot[kk];
804 
805  sum1 += ovel_tot[kk];
806  sum3 += gvel_tot[kk];
807  kk++;
808  c++;
809  }
810 
811  // determine orbit direction--asc/desc
812  if (latswt[0] > latswt[1])
813  asc_mode = 0; // descending
814  else
815  asc_mode = 1; // ascending
816 
817 
818  for (size_t i = 0; i < norbs - 1; i++) { // improve with binary search
819  if (latswt[i] < 0. && latswt[i + 1] > 0.) {//ascending orbit
820  asc_mode=1;
821  scan_brack[0] = i + 1;
822  scan_brack[1] = i + 2;
823  lat1 = latswt[i];
824  lat2 = latswt[i + 1];
825  lon1 = lonswt[i];
826  lon2 = lonswt[i + 1];
827  if(i<norbs - 2 && i>0)
828  {
829  if(lat1==latswt[i-1] || lat2==latswt[i + 2])
830  scan_brack[0] = -1;
831  scan_brack[1] = -1;
832  }
833 
834  tindex = i;
835  i = norbs;
836  break;
837  }
838 
839  else if (latswt[i] > 0. && latswt[i + 1] < 0) { // descending orbit
840  asc_mode=0;
841  scan_brack[0] = i + 2; // negative lat first convention
842  scan_brack[1] = i + 1;
843  lat1 = latswt[i + 1];
844  lat2 = latswt[i];
845  lon1 = lonswt[i + 1];
846  lon2 = lonswt[i];
847  if(i<norbs - 2 && i>0)
848  {
849  if(lat2==latswt[i-1] || lat1==latswt[i + 2])
850  scan_brack[0] = -1;
851  scan_brack[1] = -1;
852  }
853 
854  tindex = i;
855  i = norbs;
856  break;
857  }
858  } // end for
859 
860  if(l1cfile->verbose) cout<<"computing equatorial crossing time and longitude at crossing for a specific swath..orbit direction 0: descending, 1: ascending..."<<asc_mode<<endl;
861  // interpolate time--linear
862  if (scan_brack[0] > -1) { // if file with equat crossing
863  t1 = tswt[tindex];
864  t2 = tswt[tindex + 1];
865  if(l1cfile->verbose)cout << "lat1.." << lat1 << "lat2..." << lat2 << "t1.." << t1 << "t2..." << t2 << endl;
866  *tcross = t1 - lat1 * (t2 - t1) / (lat2 - lat1); // equat crossing time
867 
868  float dtcross = (*tcross - t1) / (t2 - t1);
869  *loncross = lon1 + (lon2 - lon1) * dtcross; // latcross is zero
870 
871  l1cfile->eqt = *tcross;
872  l1cfile->orbit_node_lon = *loncross;
873  l1cfile->orb_dir = asc_mode;
874  flag_ect = 0;
875  } else {
876  l1cfile->eqt = -999.0;
877  l1cfile->orbit_node_lon = -999.0;
878  l1cfile->orb_dir = asc_mode;
879  flag_ect = 1;
880  }
881 
882  // computing mean swath velocity (m/s) --------------------------
883 
884  *ovel = (sum1 / norbs); // mean velocity in m/s
885  *gvel = (sum3 / norbs) * 1000; // in meters/s
886 
887  return flag_ect;
888 }
889 
890 int sunz_swt(int ix1_swt, int ix2_swt, int16_t* hour_swt, int16_t* day_swt, int16_t* year_swt, double* olat,
891  double* olon) {
892  int day_mode = -1; // daytime is 1, nighttime is 0
893  float sunz1, suna1;
894  int year1 = year_swt[ix1_swt];
895  int day1 = day_swt[ix1_swt];
896  float hour1 = hour_swt[ix1_swt];
897  float lat1 = olat[ix1_swt];
898  float lon1 = olon[ix1_swt];
899 
900  sunangs_(&year1, &day1, &hour1, &lon1, &lat1, &sunz1, &suna1);
901  if (sunz1 <= 93)
902  day_mode = 1;
903  else
904  day_mode = 0; // nighttime
905 
906  // if(l1cfile->verbose)cout << "lat1.." << lat1 << "lon1.." << lon1 << "year1" << year1 << "day1" << day1 << "hour1" << hour1<< "sunz1..first index" << sunz1 << "ix1" << ix1_swt << "ix2" << ix2_swt << endl;
907 
908  year1 = year_swt[ix2_swt];
909  day1 = day_swt[ix2_swt];
910  hour1 = hour_swt[ix2_swt];
911  lat1 = olat[ix2_swt];
912  lon1 = olon[ix2_swt];
913 
914  sunangs_(&year1, &day1, &hour1, &lon1, &lat1, &sunz1, &suna1);
915  if (sunz1 <= 93 || day_mode == 1)
916  day_mode = 1;
917  else
918  day_mode = 0;
919 
920  // if(l1cfile->verbose)cout << "lat1.." << lat1 << "lon1.." << lon1 << "year1" << year1 << "day1" << day1 << "hour1" << hour1<< "sunz1..second index" << sunz1 << "ix1" << ix1_swt << "ix2" << ix2_swt << endl;
921 
922  return day_mode;
923 }
924 
925 int32_t L1C::create_time_swt(int num_gridlines, double tfile_ini_sec, double* tmgvf, double tswt_ini_sec,
926  double tswt_end_sec, string* tswt_ini, string* tswt_ini_file, string* tswt_mid,string* tswt_end) {
927  int16_t oneyear, onemon, oneday, onehour, onemin;
928  double onesec,tswt_mid_sec;
929  string y_swt, mo_swt, d_swt, h_swt, mi_swt, s_swt, s_swt2;
930  int logoff = -1;
931 
932  unix2ymdhms(tswt_ini_sec, &oneyear, &onemon, &oneday, &onehour, &onemin, &onesec);
933  y_swt = std::to_string(oneyear);
934  mo_swt = std::to_string(onemon);
935  d_swt = std::to_string(oneday);
936  h_swt = std::to_string(onehour);
937  mi_swt = std::to_string(onemin);
938  s_swt2 = std::to_string(round(onesec));
939 
940  length = (int)floor(log10(onemon)) + 1;
941  if (length == 1)
942  mo_swt = "0" + mo_swt;
943  length = (int)floor(log10(oneday)) + 1;
944  if (length == 1)
945  d_swt = "0" + d_swt;
946  if (onehour == 0)
947  logoff = 1;
948  else
949  logoff = 0;
950  length = (int)floor(log10(onehour + logoff)) + 1;
951  if (length == 1)
952  h_swt = "0" + h_swt;
953  if (onemin == 0)
954  logoff = 1;
955  else
956  logoff = 0;
957  length = (int)floor(log10(onemin + logoff)) + 1;
958  if (length == 1)
959  mi_swt = "0" + mi_swt;
960  if (onesec == 0)
961  logoff = 1;
962  else
963  logoff = 0;
964  length = (int)floor(log10(round(onesec) + logoff)) + 1;
965  if (length == 1)
966  s_swt2 = "0" + s_swt2;
967  if (s_swt2.substr(1, 1) == ".")
968  s_swt2 = "00";
969 
970  *tswt_ini = y_swt + "-" + mo_swt + "-" + d_swt + "T" + h_swt + ":" + mi_swt + ":" + s_swt2.substr(0, 2);
971  *tswt_ini_file = y_swt + mo_swt + d_swt + "T" + h_swt + mi_swt + s_swt2.substr(0, 2);
972 
973  tswt_end_sec = tfile_ini_sec + tmgvf[num_gridlines - 2];
974  unix2ymdhms(tswt_end_sec, &oneyear, &onemon, &oneday, &onehour, &onemin, &onesec);
975 
976  y_swt = std::to_string(oneyear);
977  mo_swt = std::to_string(onemon);
978  d_swt = std::to_string(oneday);
979  h_swt = std::to_string(onehour);
980  mi_swt = std::to_string(onemin);
981  s_swt2 = std::to_string(round(onesec));
982 
983  length = (int)floor(log10(onemon)) + 1;
984  if (length == 1)
985  mo_swt = "0" + mo_swt;
986  length = (int)floor(log10(oneday)) + 1;
987  if (length == 1)
988  d_swt = "0" + d_swt;
989  if (onehour == 0)
990  logoff = 1;
991  else
992  logoff = 0;
993  length = (int)floor(log10(onehour + logoff)) + 1;
994  if (length == 1)
995  h_swt = "0" + h_swt;
996  if (onemin == 0)
997  logoff = 1;
998  else
999  logoff = 0;
1000  length = (int)floor(log10(onemin + logoff)) + 1;
1001  if (length == 1)
1002  mi_swt = "0" + mi_swt;
1003  if (onesec == 0)
1004  logoff = 1;
1005  else
1006  logoff = 0;
1007  length = (int)floor(log10(round(onesec) + logoff)) + 1;
1008  if (length == 1)
1009  s_swt2 = "0" + s_swt2;
1010  if (s_swt2.substr(1, 1) == ".")
1011  s_swt2 = "00";
1012 
1013  *tswt_end = y_swt + "-" + mo_swt + "-" + d_swt + "T" + h_swt + ":" + mi_swt + ":" + s_swt2.substr(0, 2);
1014 
1015  tswt_mid_sec=(tswt_ini_sec+tswt_end_sec)/2;
1016  unix2ymdhms(tswt_mid_sec, &oneyear, &onemon, &oneday, &onehour, &onemin, &onesec);
1017 
1018  y_swt = std::to_string(oneyear);
1019  mo_swt = std::to_string(onemon);
1020  d_swt = std::to_string(oneday);
1021  h_swt = std::to_string(onehour);
1022  mi_swt = std::to_string(onemin);
1023  s_swt2 = std::to_string(round(onesec));
1024 
1025  length = (int)floor(log10(onemon)) + 1;
1026  if (length == 1)
1027  mo_swt = "0" + mo_swt;
1028  length = (int)floor(log10(oneday)) + 1;
1029  if (length == 1)
1030  d_swt = "0" + d_swt;
1031  if (onehour == 0)
1032  logoff = 1;
1033  else
1034  logoff = 0;
1035  length = (int)floor(log10(onehour + logoff)) + 1;
1036  if (length == 1)
1037  h_swt = "0" + h_swt;
1038  if (onemin == 0)
1039  logoff = 1;
1040  else
1041  logoff = 0;
1042  length = (int)floor(log10(onemin + logoff)) + 1;
1043  if (length == 1)
1044  mi_swt = "0" + mi_swt;
1045  if (onesec == 0)
1046  logoff = 1;
1047  else
1048  logoff = 0;
1049  length = (int)floor(log10(round(onesec) + logoff)) + 1;
1050  if (length == 1)
1051  s_swt2 = "0" + s_swt2;
1052  if (s_swt2.substr(1, 1) == ".")
1053  s_swt2 = "00";
1054 
1055  *tswt_mid = y_swt + "-" + mo_swt + "-" + d_swt + "T" + h_swt + ":" + mi_swt + ":" + s_swt2.substr(0, 2);
1056 
1057  return 0;
1058 }
1059 
1060 // open telemetry file parameters for creating L1C grid----
1061 int32_t L1C::open_l1atol1c3(L1C_input* l1cinput, l1c_filehandle* l1cfile) {
1062  int32_t n_files;
1063  const char* ptstr;
1064  char* ifile_char;
1065  string str, ifile_str;
1066  int status = -1, status1 = -1, status2 = -1;
1067  unsigned char** hkpackets = nullptr;
1068  uint8_t* apids = nullptr;
1069  size_t number_hkpack, number_scpack, number_orecords, number_hkpack_tot=0, number_scpack_tot=0,
1070  number_orecords_tot, nr=0;
1071  uint8_t ubnum1, ubnum2;
1072  double *sec = nullptr, *tai58_sec = nullptr;
1073  int16_t *year = nullptr, *mon = nullptr, *day = nullptr, *hour = nullptr, *min = nullptr;
1074  int16_t oneyear, onemon, oneday, onehour, onemin;
1075  double onesec;
1076  double *orb_time = nullptr, *orb_lat = nullptr, *orb_lon = nullptr;
1077  float **lat_gd = nullptr, **lon_gd = nullptr, **alt = nullptr;
1078  double omeg = 7.29211585494e-5;
1079  short n_ephem;
1080  string temp_str, tai_str;
1081  double vxyz = 0, mov1 = 0., mgv1 = 0., mov2 = 0, mgv2, *orb_vel = nullptr, *grn_vel = nullptr;
1082  int* orb_dir = nullptr;
1083  int ix1 = -1, ix2 = -1, ix3 = -1, ix4 = -1, ix5 = -1, ix6 = -1;
1084  double *tmgv1 = nullptr, *tmgv2 = nullptr, *tmgvf = nullptr, *tmgvf2 = nullptr;
1085  int32_t num_gridlines = -1, norbs = -1;
1086  double *tswt = nullptr, *latswt = nullptr, *lonswt = nullptr;
1087  double *tswt2 = nullptr, *latswt2 = nullptr, *lonswt2 = nullptr;
1088  const char* outlist;
1089  string ofile_str, senstr;
1090  string y_swt, mo_swt, d_swt, h_swt, mi_swt, s_swt, tswt_ini, tswt_end,tswt_mid;
1091  int32_t gd_per_gran = -1, numgran = -1;
1092  double deltasec = -1;
1093  string s_swt2;
1094  string tswt_ini_file;
1095  double rl2, pos_norm, clatg2, fe = 1 / 298.257;
1096  int day_mode = -1; // 0 is nighttime 1 is dayttime
1097  double tai58unix;
1098  double tswt_ini_sec, tswt_end_sec;
1099  vector<size_t> ix;
1100  double tfile_ini, tfile_end;
1101  int16_t syear, smon, sday, shour, smin, syear2, smon2, sday2, shour2, smin2;
1102  double secs, second, secs2, second2;
1103  double tfile_ini_sec, tfile_end_sec;
1104  size_t ix_ini=-1,ix_end=-1;
1105 
1106  // TOTAL ARRAYS---ASSUMING 6000 rcords
1107  number_orecords_tot = 6000;
1108  size_t cc = 0, c = 0;
1109 
1110  int16_t* year_tot = (int16_t*)calloc(number_orecords_tot, sizeof(int16_t));
1111  int16_t* day_tot = (int16_t*)calloc(number_orecords_tot, sizeof(int16_t));
1112  int16_t* hour_tot = (int16_t*)calloc(number_orecords_tot, sizeof(int16_t));
1113  double* orb_time_tot = (double*)calloc(number_orecords_tot, sizeof(double));
1114  double* orb_lat_tot = (double*)calloc(number_orecords_tot, sizeof(double));
1115  double* orb_lon_tot = (double*)calloc(number_orecords_tot, sizeof(double));
1116  int* orb_dir_tot = (int*)calloc(number_orecords_tot, sizeof(int));
1117  orb_array2* posr_tot = new orb_array2[number_orecords_tot]();
1118  orb_array2* velr_tot = new orb_array2[number_orecords_tot]();
1119  double* orb_vel_tot = (double*)calloc(number_orecords_tot, sizeof(double));
1120  double* grn_vel_tot = (double*)calloc(number_orecords_tot, sizeof(double));
1121 
1122  //*************************************************************
1123  ifile_str = l1cinput->files[0];
1124  ifile_char = (char*)ifile_str.c_str();
1125  ofile_str = ifile_str.substr(0, 24);
1126  outlist = "list_l1c_granules.txt"; // DEFAULT OFILE
1127 
1128  if (l1cinput->outlist[0] == '\0') { // only when no txt output list file is provided
1129  strcpy(l1cinput->outlist, outlist);
1130  if(l1cinput->verbose)cout << "L1C granules written to DEFAULT file........" << outlist << endl;
1131  } else {
1132  if(l1cinput->verbose)cout << "L1C granules written to "
1133  "file......................................................................................."
1134  << l1cinput->outlist << endl;
1135  }
1136 
1137  l1cfile->gransize = l1cinput->gransize;
1138 
1139  file_format format = getFormat(ifile_char);
1140  if(l1cinput->verbose)cout << "format.type.." << format.type << endl;
1141 
1142  l1cfile->sensor = l1cinput->sensor; // SPEX 1, OCI 2 and HARP 3
1143 
1144  if (l1cinput->sensor == 34) {
1145  senstr = "SPEXONE";
1146  l1cfile->nbinx = 29;
1147  l1cfile->n_views = 10;
1148  } else if (l1cinput->sensor == 30) {
1149  senstr = "OCI";
1150  l1cfile->nbinx = 519;
1151  l1cfile->n_views = 2;
1152  } else if (l1cinput->sensor == 35) {
1153  senstr = "HARP2";
1154  l1cfile->nbinx = 457;
1155  l1cfile->n_views = 90;
1156  } else {
1157  if(l1cinput->verbose)cout << "sensor by default is OCI option 2....." << endl;
1158  senstr = "OCI";
1159  l1cfile->nbinx = 514;
1160  l1cfile->n_views = 2;
1161  }
1162 
1163  if(l1cinput->verbose)cout << "PACE sensor to be griddded....." << senstr << endl;
1164 
1165  n_files = l1cfile->ifiles.size();
1166  if(l1cinput->verbose)cout << "number of files in the list...." << n_files << endl;
1167 
1168  int32_t nfiles = l1cfile->ifiles.size();
1169 
1170  int FirsTerrain=0;
1171 
1172  for (int fi = 0; fi < nfiles; fi++)
1173  {
1174  if(fi==nfiles-1)FirsTerrain=1;
1175  str = l1cfile->ifiles[fi];
1176  ptstr = str.c_str();
1177  if(l1cinput->verbose)
1178  {
1179  cout << "********************************************************************************************"
1180  "*************"
1181  << endl;
1182  cout << "Opening L1A file ..." << ptstr << "for telemetry......." << endl;
1183  cout << "********************************************************************************************"
1184  "*************"
1185  << endl;
1186  }
1187 
1188  NcFile* nc_l1a;
1189  try {
1190  nc_l1a = new NcFile(ptstr, NcFile::read);
1191  } catch (NcException& e) {
1192  e.what();
1193  cerr << "l1cgen l1c_pflag= 8:: Failure write FULL L1C grid: " + str << endl;
1194  exit(1);
1195  }
1196 
1197  string name;
1198  NcGroupAtt i1 = nc_l1a->getAtt("time_coverage_start"); // NcGroupAtt is a global attr!!
1199  i1.getValues(name);
1200 
1201  tfile_ini = isodate2unix(name.c_str());
1202 
1203  i1 = nc_l1a->getAtt("time_coverage_end"); // NcGroupAtt is a global attr!!
1204  i1.getValues(name);
1205  tfile_end = isodate2unix(name.c_str());
1206 
1207  // this is default start time, start_timeflag=1 or using start time of swath file HKT------
1208  unix2ymds(tfile_ini, &syear, &smon, &sday, &secs);
1209  unix2ymds(tfile_end, &syear2, &smon2, &sday2, &secs2);
1210 
1211  if(l1cinput->verbose)cout << "secs elapsed.." << secs << "initial granule #..." << round(secs / (l1cfile->gransize * 60))
1212  << endl;
1213  unix2ymdhms(tfile_ini, &syear, &smon, &sday, &shour, &smin, &second);
1214  if(l1cinput->verbose)cout << "HKT file start time................."
1215  << "year.." << syear << "month..." << smon << "day..." << sday << "hour.." << shour << "min.."
1216  << smin << "sec..." << second << endl;
1217  unix2ymdhms(tfile_end, &syear2, &smon2, &sday2, &shour2, &smin2, &second2);
1218  if(l1cinput->verbose)cout << "HKT file end time................."
1219  << "year.." << syear2 << "month..." << smon2 << "day..." << sday2 << "hour.." << shour2
1220  << "min.." << smin2 << "sec..." << second2 << endl;
1221 
1222  if (l1cinput->start_timeflag == 0) { // forcing to hms =0:0:0
1223  secs = 0.;
1224  }
1225 
1226  tfile_ini_sec = ymds2unix(syear, smon, sday, secs);
1227  tfile_end_sec =
1228  tfile_ini_sec + 49 * 60 * 3; // 49 minutes per half-orbit, 3x cause 1.5 orbit, forthy swaths (3 before)
1229 
1230  int16_t y_ini, mo_ini, d_ini, h_ini, mi_ini;
1231  double sec_ini;
1232  unix2ymdhms(tfile_ini_sec, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
1233  if(l1cinput->verbose)cout << "tfile_ini_sec.."
1234  << "YEAR.." << y_ini << "MONTH.." << mo_ini << "DAY.." << d_ini << "HOUR.." << h_ini << "MIN.."
1235  << mi_ini << "SEC.." << sec_ini << endl;
1236  unix2ymdhms(tfile_end_sec, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
1237  if(l1cinput->verbose)cout << "tfile_end_sec.."
1238  << "YEAR.." << y_ini << "MONTH.." << mo_ini << "DAY.." << d_ini << "HOUR.." << h_ini << "MIN.."
1239  << mi_ini << "SEC.." << sec_ini << endl;
1240 
1241  l1cfile->tfile_ini_sec = tfile_ini_sec;
1242  l1cfile->tfile_end_sec = tfile_end_sec;
1243 
1244  // open dimensions
1246  NcDim dimhp = nc_l1a->getDim("SC_hkt_pkts");
1247  NcDim dimtp = nc_l1a->getDim("max_SC_packet");
1248  NcDim dimor = nc_l1a->getDim("orb_records");
1249 
1250  number_hkpack = dimhp.getSize();
1251  number_scpack = dimtp.getSize();
1252  number_orecords = dimor.getSize();
1253 
1254  // total counters---
1255  number_hkpack_tot += number_hkpack;
1256  number_scpack_tot += number_scpack;
1257  nr += number_orecords;
1258 
1259  if(l1cinput->verbose)cout << "number_hkpack_tot.." << number_hkpack_tot << "number_scpack_tot.." << number_scpack_tot
1260  << "number of orbit records total.REAL." << nr << endl;
1261 
1262  // allocat mem
1263  hkpackets =
1264  allocate2d_uchar(number_hkpack, number_scpack); // NUMBER of ephem elements x max SIZE of package
1265  apids = (uint8_t*)calloc(number_hkpack, sizeof(uint8_t));
1266 
1267  orb_time = (double*)calloc(number_orecords, sizeof(double));
1268  orb_lat = (double*)calloc(number_orecords, sizeof(double));
1269  orb_lon = (double*)calloc(number_orecords, sizeof(double));
1270  orb_dir = (int*)calloc(number_orecords, sizeof(int));
1271 
1272  // open groups
1273  NcGroup telGrp = nc_l1a->getGroup("housekeeping_data");
1274  NcGroup navGrp = nc_l1a->getGroup("navigation_data");
1275 
1276  // open vars ids
1277  NcVar v1 = telGrp.getVar("SC_HKT_packets");
1278  v1.getVar(&hkpackets[0][0]);
1279 
1280  v1 = navGrp.getVar("orb_time");
1281  v1.getVar(&orb_time[0]);
1282  v1 = navGrp.getVar("orb_lat");
1283  v1.getVar(&orb_lat[0]);
1284  v1 = navGrp.getVar("orb_lon");
1285  v1.getVar(&orb_lon[0]);
1286 
1287  for (size_t hk = 0; hk < number_hkpack; hk++) {
1288  ubnum1 = (uint8_t)hkpackets[hk][0]; //-48;//48 is 0 ascii code
1289  ubnum2 = (uint8_t)hkpackets[hk][1];
1290 
1291  apids[hk] = (ubnum1 % 8) * 256 + ubnum2;
1292  if (apids[hk] == 128) {
1293  ix.push_back(hk); // packet index where ephemr
1294  }
1295  }
1296 
1297  if(l1cinput->verbose)cout << "#number of ephem elements...." << ix.size() << "for HKT file..." << ptstr << "#orecords..."
1298  << number_orecords << endl;
1299  n_ephem = ix.size();
1300 
1301  // Reverse the byte order from small endian to large endian, small endian last byte is stored first
1302  // convert to double and later swap endian
1303 
1304  // allocate mem
1305  orb_vel = (double*)calloc(n_ephem, sizeof(double));
1306  grn_vel = (double*)calloc(n_ephem, sizeof(double));
1307 
1308  sec = (double*)calloc(n_ephem, sizeof(double));
1309  year = (int16_t*)calloc(n_ephem, sizeof(int16_t));
1310  mon = (int16_t*)calloc(n_ephem, sizeof(int16_t));
1311  day = (int16_t*)calloc(n_ephem, sizeof(int16_t));
1312  hour = (int16_t*)calloc(n_ephem, sizeof(int16_t));
1313  min = (int16_t*)calloc(n_ephem, sizeof(int16_t));
1314 
1315  tai58_sec = (double*)calloc(n_ephem, sizeof(double));
1316 
1317  // process all packets----
1318  double tai58;
1319  c = 0;
1320  // TIME
1321 
1322  while (c < ix.size()) {
1323  double* tai_ptr = (double*)(hkpackets[ix[c]] + 16); // moving pointer to position 16
1324  swapc_bytes((char*)tai_ptr, 8,
1325  1); // swap 8 bytes once or 1 16-23 pos, lets say we have 16-32 indexes and we want
1326  // to swap 8 bytes at the time, some we need ntime=2
1327  tai58 = *tai_ptr;
1328  tai58unix = tai58_to_unix(tai58);
1329 
1330  unix2ymdhms(tai58unix, &oneyear, &onemon, &oneday, &onehour, &onemin, &onesec);
1331 
1332  year[c] = oneyear;
1333  mon[c] = onemon;
1334  day[c] = oneday;
1335  hour[c] = onehour;
1336  min[c] = onemin;
1337  sec[c] = onesec;
1338 
1339  tai58_sec[c] = tai58unix;
1340 
1341  c++;
1342  }
1343 
1344  // POS/VEL alloc mem-----
1345  orb_array2* posr = new orb_array2[n_ephem]();
1346  orb_array2* velr = new orb_array2[n_ephem]();
1347 
1348  double dp1, dp2, dp3;
1349 
1350  // computed orbital velocity from ephemeris
1351  c = 0;
1352  while (c < ix.size()) {
1353  double* posi = (double*)(hkpackets[ix[c]] + 120);
1354  double* veli = (double*)(hkpackets[ix[c]] + 144);
1355  double* ecmat = (double*)(hkpackets[ix[c]] + 176);
1356 
1357  swapc_bytes((char*)posi, 8, 3);
1358  swapc_bytes((char*)veli, 8, 3);
1359  swapc_bytes((char*)ecmat, 8, 9);
1360 
1361  // assuming ecmat index 1-3 first row, 4-6 second row and 7-9 third row
1362  // dot product
1363  // position
1364  dp1 = posi[0] * ecmat[0] + posi[1] * ecmat[1] + posi[2] * ecmat[2];
1365  dp2 = posi[0] * ecmat[3] + posi[1] * ecmat[4] + posi[2] * ecmat[5];
1366  dp3 = posi[0] * ecmat[6] + posi[1] * ecmat[7] + posi[2] * ecmat[8];
1367 
1368  posr[c][0] = dp1;
1369  posr[c][1] = dp2;
1370  posr[c][2] = dp3;
1371 
1372  // velocity
1373  dp1 = veli[0] * ecmat[0] + veli[1] * ecmat[1] + veli[2] * ecmat[2];
1374  dp2 = veli[0] * ecmat[3] + veli[1] * ecmat[4] + veli[2] * ecmat[5];
1375  dp3 = veli[0] * ecmat[6] + veli[1] * ecmat[7] + veli[2] * ecmat[8];
1376 
1377  velr[c][0] = dp1;
1378  velr[c][1] = dp2;
1379  velr[c][2] = dp3;
1380 
1381  velr[c][0] = velr[c][0] + posr[c][1] * omeg;
1382  velr[c][1] = velr[c][1] - posr[c][0] * omeg;
1383 
1384  // computing orbital velocity
1385  vxyz = sqrt(velr[c][0] * velr[c][0] + velr[c][1] * velr[c][1] +
1386  velr[c][2] * velr[c][2]); // units m/s
1387  orb_vel[c] = vxyz;
1388 
1389  pos_norm = sqrt(posr[c][0] * posr[c][0] + posr[c][1] * posr[c][1] + posr[c][2] * posr[c][2]);
1390  clatg2 = sqrt(posr[c][0] * posr[c][0] + posr[c][1] * posr[c][1]) / pos_norm;
1391  rl2 = Re * (1 - fe) / (sqrt(1 - (2 - fe) * fe * clatg2 * clatg2));
1392  grn_vel[c] = vxyz * rl2 / pos_norm;
1393 
1394  c++;
1395  }
1396 
1397  for (short n = 0; n < n_ephem; n++) {
1398  orb_time_tot[cc] = orb_time[n];
1399  orb_lon_tot[cc] = orb_lon[n];
1400  orb_lat_tot[cc] = orb_lat[n];
1401  orb_vel_tot[cc] = orb_vel[n];
1402  grn_vel_tot[cc] = grn_vel[n];
1403  velr_tot[cc][0] = velr[n][0];
1404  velr_tot[cc][1] = velr[n][1];
1405  velr_tot[cc][2] = velr[n][2];
1406  posr_tot[cc][0] = posr[n][0];
1407  posr_tot[cc][1] = posr[n][1];
1408  posr_tot[cc][2] = posr[n][2];
1409  year_tot[cc] = year[n];
1410  day_tot[cc] = day[n];
1411  hour_tot[cc] = hour[n];
1412  cc++;
1413  }
1414 
1415  if (hkpackets != nullptr)
1416  free (hkpackets);
1417  hkpackets = nullptr;
1418 
1419  if (apids != nullptr)
1420  free (apids);
1421  apids = nullptr;
1422 
1423  ix.clear();
1424 
1425  if (tai58_sec != nullptr)
1426  free (tai58_sec);
1427  tai58_sec = nullptr;
1428 
1429  if (sec != nullptr)
1430  free (sec);
1431  sec = nullptr;
1432 
1433  if (year != nullptr)
1434  free (year);
1435  year = nullptr;
1436 
1437  if (mon != nullptr)
1438  free (mon);
1439  mon = nullptr;
1440 
1441  if (day != nullptr)
1442  free (day);
1443  day = nullptr;
1444 
1445  if (hour != nullptr)
1446  free (hour);
1447  hour = nullptr;
1448 
1449  if (min != nullptr)
1450  free (min);
1451  min = nullptr;
1452 
1453  if (orb_lat != nullptr)
1454  free (orb_lat);
1455  orb_lat = nullptr;
1456 
1457  if (orb_lon != nullptr)
1458  free (orb_lon);
1459  orb_lon = nullptr;
1460 
1461  if (orb_time != nullptr)
1462  free (orb_time);
1463  orb_time = nullptr;
1464 
1465  if (orb_vel != nullptr)
1466  free (orb_vel);
1467  orb_vel = nullptr;
1468 
1469  if (grn_vel != nullptr)
1470  free (grn_vel);
1471  grn_vel = nullptr;
1472 
1473  if (orb_dir != nullptr)
1474  free (orb_dir);
1475  orb_dir = nullptr;
1476 
1477  if (posr != nullptr)
1478  delete[] (posr);
1479  posr = nullptr;
1480 
1481  if (velr != nullptr)
1482  delete[] (velr);
1483  velr = nullptr;
1484  if (senz != nullptr)
1485  free (senz);
1486  senz = nullptr;
1487 
1488  if (sena != nullptr)
1489  free (sena);
1490  sena = nullptr;
1491 
1492  nc_l1a->close();
1493 
1494  if(l1cinput->verbose)cout<<"total number or orbital records #"<<cc-1<<endl;
1495  number_orecords_tot=cc-1;
1496  l1cfile->norb_rec=number_orecords_tot;
1497  // determine ini/end time indexes for asc/desc passes
1498  for (size_t k = 0; k < number_orecords_tot - 1; k++) {
1499  if (orb_lat_tot[k + 1] > orb_lat_tot[k])
1500  orb_dir_tot[k] = 1; // ascending
1501  else
1502  orb_dir_tot[k] = 0; // descending
1503  }
1504 
1505  // find time limits for half-orbits
1506  int swt = 1;
1507  int tlimits_swt[2][2] = {{-1, -1}, {-1, -1}}; // 2 swath x time limits ini/end
1508  tlimits_swt[swt - 1][0] = 0; // only the record indexes
1509 
1510  for (size_t k = 0; k < number_orecords_tot - 1; k++) {
1511  if (orb_dir_tot[k] != orb_dir_tot[k + 1]) {
1512  if (swt == 1) {
1513  tlimits_swt[swt - 1][1] = k; // this time is a 'local' end cause may continue circling the
1514  // earth after ascending or descending
1515  swt++;
1516  tlimits_swt[swt - 1][0] = k + 1;
1517  } else if (swt == 2) { // another partial swath
1518  tlimits_swt[swt - 1][1] = k;
1519  tlimits_swt[swt - 2][1] = number_orecords_tot - 1;
1520  swt++;
1521  }
1522  }
1523  if (k == number_orecords_tot - 2 && swt == 2)
1524  tlimits_swt[swt - 1][1] = number_orecords_tot - 1;
1525  }
1526 
1527  nswath = swt;
1528 
1529  if(l1cinput->verbose)cout << "nswath...................." << nswath << "number_orecords_tot.." << number_orecords_tot << endl;
1530 
1531  // split half-orbits---- group lat/lon time and velocity vector
1532  for (size_t sw = 0; sw < 2; sw++) {
1533  if (nswath == 2) {
1534  if (sw == 0) {
1535  ix1 = tlimits_swt[sw][0];
1536  ix2 = tlimits_swt[sw][1];
1537  if(l1cinput->verbose)cout << "nswath#2......swat#1...ix1" << ix1 << "ix2.." << ix2 << endl;
1538  } else if (sw == 1) {
1539  ix3 = tlimits_swt[sw][0];
1540  // ix4 = tlimits_swt[sw][1];
1541  ix4 = number_orecords_tot - 1;
1542  if(l1cinput->verbose)cout << "nswath2.....swat#2...ix3" << ix3 << "ix4.." << ix4 << endl;
1543  }
1544  } else if (nswath == 3) { //>2 we have 4 limits for swath 1
1545  if (sw == 0) {
1546  ix1 = tlimits_swt[sw][0];
1547  ix2 = tlimits_swt[sw + 1][0] - 1;
1548  ix3 = tlimits_swt[sw + 1][1] + 1;
1549  ix4 = tlimits_swt[sw][1];
1550  if(l1cinput->verbose)cout << "nswath3....swat#1...ix1" << ix1 << "ix2.." << ix2 << "ix3.." << ix3 << "ix4..."
1551  << ix4 << endl;
1552  } else if (sw == 1) {
1553  ix5 = tlimits_swt[sw][0];
1554  //ix6 = tlimits_swt[sw][1];
1555 
1556  // This is a cludge to get the end granule to work.
1557  // Just use all of the HKT orbit records
1558  //ix5 = 2;
1559  ix6 = number_orecords_tot - 1;
1560 
1561  if(l1cinput->verbose)cout << "nswath3....swat#2...ix5" << ix5 << "ix6.." << ix6 << endl;
1562  }
1563  } else {
1564  if(l1cinput->verbose)cout << "number of swaths is less than 2..!! exit..." << endl;
1565  // exit(1);
1566  }
1567  }
1568 
1569  float tcross1 = -999., loncross1 = -999., tcross2 = -999., loncross2 = -999.;
1570 
1571 
1572  l1cfile->num_gridlines = 4800;//3860
1573  num_gridlines = l1cfile->num_gridlines;
1574  //--------------------------------------------------------------------------------------
1575  // nswath#2
1576  //----case type: half orbits are consecutive
1577  // swath1
1578  //--------------------------------------------------------------------------------------
1579  if (nswath == 2) {
1580  if (ix1 >= 0 && ix5 < 0) {
1581  norbs = ix2 - ix1 + 1;
1582  // compute solar zenith angle for checking day/night conditions for each swath ---
1583  //-------------------------------------------------------
1584  day_mode = sunz_swt(ix1, ix2, hour_tot, day_tot, year_tot, orb_lat_tot, orb_lon_tot);
1585 
1586  if(l1cinput->verbose)cout << "day_mode at nswath #2 SWATH1." << day_mode << endl;
1587 
1588  if (day_mode == 1) {
1589  // orbit direction
1590  l1cfile->orb_dir = orb_dir_tot[ix1];
1591 
1592  // ini and end UTC--this is swath orbital time, not interpolated time
1593  tswt = (double*)calloc(norbs, sizeof(double));
1594  latswt = (double*)calloc(norbs, sizeof(double));
1595  lonswt = (double*)calloc(norbs, sizeof(double));
1596 
1597  status = ect_swt(l1cfile, ix1, ix2, orb_time_tot, orb_lat_tot, orb_lon_tot, orb_vel_tot,
1598  grn_vel_tot, tswt, latswt, lonswt, &tcross1, &loncross1, &mov1, &mgv1);
1599 
1600  if(l1cinput->verbose)cout << "nswath==2 --tcross equat crossing in (s)..swath#1..." << tcross1 << "loncross1..."
1601  << loncross1 << endl;
1602 
1603  if (latswt != nullptr)
1604  free (latswt);
1605  latswt = nullptr;
1606  if (lonswt != nullptr)
1607  free (lonswt);
1608  lonswt = nullptr;
1609 
1610 
1611  if (status == 0) {
1612  tmgv1 = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
1613  tmgvf = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
1614 
1615  swtime_swt2(1, l1cinput, l1cfile, norbs, tswt, tcross1, mgv1, tmgv1);
1616  // ini/end times for swath
1617  num_gridlines = l1cfile->num_gridlines;
1618  if(l1cinput->verbose)cout << "number of across bins L1C grid...#.." << l1cfile->nbinx << "l1cfile->n_views..."
1619  << l1cfile->n_views << endl;
1620 
1621  lat_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
1622  lon_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
1623  alt = allocate2d_float(num_gridlines, l1cfile->nbinx);
1624 
1625  orb_to_latlon(ix1,ix4,num_gridlines, l1cfile->nbinx, orb_time_tot, posr_tot,
1626  velr_tot, mgv1, tmgv1, tmgvf, lat_gd, lon_gd, alt,FirsTerrain);
1627 
1628 
1629  if (tmgv1 != nullptr)
1630  free (tmgv1);
1631  tmgv1 = nullptr;
1632 
1633  l1cfile->num_gridlines = l1cfile->num_gridlines - 1;
1634 
1635  if (l1cinput->grantype == 1) {
1636  tswt_ini_sec = tfile_ini_sec + tmgvf[0];
1637  tswt_end_sec = tfile_ini_sec + tmgvf[num_gridlines - 2];
1638  create_time_swt(num_gridlines, tfile_ini_sec, tmgvf, tswt_ini_sec,tswt_end_sec,
1639  &tswt_ini, &tswt_ini_file,&tswt_mid, &tswt_end);
1640 
1641  l1cfile->tswt_ini = tswt_ini;
1642  l1cfile->tswt_mid = tswt_mid;
1643  l1cfile->tswt_end = tswt_end;
1644  l1cfile->tswt_ini_file = tswt_ini_file;
1645  create_SOCEA2(1, l1cinput, l1cfile, lat_gd, lon_gd, alt,
1646  tmgvf); // THIS IS SWATH PROCESSING
1647  } else if (l1cinput->grantype == 0) {
1648  //--------------------------------------------------------------
1649  // granule processing---------------------------------------------
1650  //-----------------------------------------------------------
1651  deltasec = tmgvf[num_gridlines - 2] - tmgvf[0] + 1;
1652  if(l1cinput->verbose)cout << "deltasec..swath." << deltasec << endl;
1653 
1654  if (tswt != nullptr)
1655  free (tswt);
1656  tswt = nullptr;
1657 
1658  // numgran = 144 * 2;//1 day
1659  numgran=12;
1660  l1cfile->numgran = numgran;
1661  gd_per_gran = round(num_gridlines / 10); // 10 granules per half orbit
1662  l1cfile->gd_per_gran = gd_per_gran;
1663  if(l1cinput->verbose)cout << "estimated # of granules to be processed..." << numgran << "gd_per_gran..."
1664  << gd_per_gran << "#gridlines.." << num_gridlines << endl;
1665 
1666  if(FirsTerrain) write_L1C_granule2(1, l1cfile, l1cinput, tmgvf, lat_gd, lon_gd, alt,orb_time_tot);
1667  } else {
1668  cout << "ERROR selecting grantype, must be 0: granules or 1: "
1669  "swath........................."
1670  << endl;
1671  exit(1);
1672  }
1673  //-----------------------------------------------------------------
1674  //-----------------------------------------------------------------
1675  if (lat_gd != nullptr)
1676  free (lat_gd);
1677  lat_gd = nullptr;
1678  if (lon_gd != nullptr)
1679  free (lon_gd);
1680  lon_gd = nullptr;
1681  if (tmgvf != nullptr)
1682  free (tmgvf);
1683  tmgvf = nullptr;
1684  if (alt != nullptr)
1685  free (alt);
1686  alt = nullptr;
1687 
1688  } // status =0
1689 
1690  } // end day_mode
1691  else {
1692  if(l1cinput->verbose)cout << "ERROR swath #1 day_mode = " << day_mode
1693  << "nightime...continue to swath2 (nswath#2).............." << endl;
1694  }
1695 
1696  // nswath=2
1697  //-------------------------------------------------------------------------------------------------------
1698  // swath#2
1699  //-------------------------------------------------------------------------------------------------------
1700  // if(l1cinput->swath_num==2){
1701  norbs = ix4 - ix3 + 1;
1702  // compute solar zenith angle for checking day/night conditions for each swath --
1703  day_mode = sunz_swt(ix3, ix4, hour_tot, day_tot, year_tot, orb_lat_tot, orb_lon_tot);
1704 
1705  if(l1cinput->verbose)cout << "day_mode at nswath #2 SWATH2" << day_mode << endl;
1706 
1707  if (day_mode == 1) {
1708  // orbit direction
1709  l1cfile->orb_dir = orb_dir_tot[ix3];
1710 
1711  tswt2 = (double*)calloc(norbs, sizeof(double));
1712  latswt2 = (double*)calloc(norbs, sizeof(double));
1713  lonswt2 = (double*)calloc(norbs, sizeof(double));
1714 
1715  status = ect_swt(l1cfile, ix3, ix4, orb_time_tot, orb_lat_tot, orb_lon_tot, orb_vel_tot,
1716  grn_vel_tot, tswt2, latswt2, lonswt2, &tcross2, &loncross2, &mov2, &mgv2);
1717 
1718  if(l1cinput->verbose)cout << "nswath==2 ---tcross equat crossing in (s)..swath#2..." << tcross2 << "loncross2..."
1719  << loncross2 << endl;
1720 
1721  if (latswt2 != nullptr)
1722  free (latswt2);
1723  latswt2 = nullptr;
1724  if (lonswt2 != nullptr)
1725  free (lonswt2);
1726  lonswt2 = nullptr;
1727 
1728 
1729 
1730  if (status == 0) {
1731  tmgv2 = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
1732  tmgvf2 = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
1733 
1734  swtime_swt2(2, l1cinput, l1cfile, norbs, tswt2, tcross2, mgv2, tmgv2);
1735  num_gridlines = l1cfile->num_gridlines;
1736 
1737  // granule processing ------------------------
1738  numgran = 144;
1739  l1cfile->numgran = numgran;
1740  gd_per_gran = round(num_gridlines / 10);
1741  l1cfile->gd_per_gran = gd_per_gran;
1742 
1743  if(l1cinput->verbose)cout << "# of granules to be processed..." << numgran << "gd_per_gran..." << gd_per_gran
1744  << endl;
1745 
1746  if (tswt2 != nullptr)
1747  free (tswt2);
1748  tswt2 = nullptr;
1749 
1750  if(l1cinput->verbose)cout << "number of across bins L1C grid...#.." << l1cfile->nbinx << endl;
1751  lat_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
1752  lon_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
1753  alt = allocate2d_float(num_gridlines, l1cfile->nbinx);
1754 
1755  orb_to_latlon(ix1,ix4,num_gridlines, l1cfile->nbinx, orb_time_tot, posr_tot,
1756  velr_tot, mgv2, tmgv2, tmgvf2, lat_gd, lon_gd, alt,FirsTerrain);
1757 
1758 
1759  if (tmgv2 != nullptr)
1760  free (tmgv2);
1761  tmgv2 = nullptr;
1762 
1763  l1cfile->num_gridlines = l1cfile->num_gridlines - 1;
1764 
1765  if (l1cinput->grantype == 1) {
1766  tswt_ini_sec = tfile_ini_sec + tmgvf2[0];
1767  tswt_end_sec = tfile_ini_sec + tmgvf2[num_gridlines - 2];
1768 
1769  create_time_swt(num_gridlines, tfile_ini_sec, tmgvf2, tswt_ini_sec, tswt_end_sec,
1770  &tswt_ini, &tswt_ini_file, &tswt_mid,&tswt_end);
1771 
1772  l1cfile->tswt_ini = tswt_ini;
1773  l1cfile->tswt_mid = tswt_mid;
1774  l1cfile->tswt_end = tswt_end;
1775  l1cfile->tswt_ini_file = tswt_ini_file;
1776  create_SOCEA2(2, l1cinput, l1cfile, lat_gd, lon_gd, alt,
1777  tmgvf2); // THIS IS SWATH PROCESSING
1778  } else if (l1cinput->grantype == 0) {
1779  //--------------------------------------------------------------
1780  // granule processing---------------------------------------------
1781  //---------------------------------------------------------------
1782  deltasec = tmgvf2[num_gridlines - 2] - tmgvf2[0] + 1;
1783  if(l1cinput->verbose)cout << "deltasec..swath." << deltasec << endl;
1784 
1785  numgran = 144 * 2;
1786  l1cfile->numgran = numgran;
1787  gd_per_gran = round(num_gridlines / 10); // 10 granules per half orbit
1788  l1cfile->gd_per_gran = gd_per_gran;
1789 
1790  if(l1cinput->verbose)cout << "estimated # of granules to be processed..." << numgran << "gd_per_gran..."
1791  << gd_per_gran << "#gridlines.." << num_gridlines << endl;
1792  if(FirsTerrain) write_L1C_granule2(2, l1cfile, l1cinput, tmgvf2, lat_gd, lon_gd, alt,orb_time_tot);
1793  } else {
1794  cout << "ERROR selecting grantype, must be 0: granules or 1: "
1795  "swath........................."
1796  << endl;
1797  exit(1);
1798  }
1799  //--------------------------------------------------------------------
1800  if (lat_gd != nullptr)
1801  free (lat_gd);
1802  lat_gd = nullptr;
1803  if (lon_gd != nullptr)
1804  free (lon_gd);
1805  lon_gd = nullptr;
1806  if (tmgvf2 != nullptr)
1807  free (tmgvf2);
1808  tmgvf2 = nullptr;
1809  if (alt != nullptr)
1810  free (alt);
1811  alt = nullptr;
1812  } // status 0
1813  else {
1814  if(l1cinput->verbose)cout << "ERROR swath #2 does not cross the equator..NO L1C grid for swath #2" << endl;
1815  }
1816 
1817  } else {
1818  if(l1cinput->verbose)cout << "day_mode = 0 nightime...no L1C grid produced--exit (swath#2 --- "
1819  "nswath#2).............."
1820  << day_mode << endl;
1821 // exit(1);
1822  }
1823 
1824  } // end index x1 and x5 condition
1825 
1826  } // end nswath=2
1827 
1828  //-------------HALF-ORBITS IN TWO NON-CONSECUTIVE PORTIONS ------
1829  // nswath =3
1830  // swath1
1831  //-------------------------------------------------------------
1832  c = 0;
1833  if (nswath == 3) {
1834  if (ix1 >= 0 && ix5 > 0) {
1835  norbs = ix2 - ix1 + 1;
1836 
1837  day_mode = sunz_swt(ix1, ix2, hour_tot, day_tot, year_tot, orb_lat_tot, orb_lon_tot);
1838 
1839  if(l1cinput->verbose)cout << "day_mode at nswath #3 SWATH1--segment#1: " << day_mode << endl;
1840 
1841  if (day_mode == 1) { // daylight
1842  l1cfile->orb_dir = orb_dir_tot[ix1];
1843 
1844  tswt = (double*)calloc(norbs, sizeof(double));
1845  latswt = (double*)calloc(norbs, sizeof(double));
1846  lonswt = (double*)calloc(norbs, sizeof(double));
1847 
1848  status1 = ect_swt(l1cfile, ix1, ix2, orb_time_tot, orb_lat_tot, orb_lon_tot, orb_vel_tot,
1849  grn_vel_tot, tswt, latswt, lonswt, &tcross1, &loncross1, &mov1, &mgv1);
1850  if(l1cinput->verbose)cout << "nswath==3 --tcross equat crossing in (s)..swath#1.(segment #1).." << tcross1
1851  << "loncross1..." << loncross1 << "orbit direction.0 is descending." << l1cfile->orb_dir
1852  << "mov1.." << mov1 << "mgv1.." << mgv1 << endl;
1853 
1854  c = 0;
1855  if (tcross1 < 0.) {
1856  if (tswt != nullptr)
1857  free (tswt);
1858  tswt = nullptr;
1859  if (latswt != nullptr)
1860  free (latswt);
1861  latswt = nullptr;
1862  if (lonswt != nullptr)
1863  free (lonswt);
1864  lonswt = nullptr;
1865 
1866  norbs = ix4 - ix3 + 1;
1867 
1868  tswt = (double*)calloc(norbs, sizeof(double));
1869  latswt = (double*)calloc(norbs, sizeof(double));
1870  lonswt = (double*)calloc(norbs, sizeof(double));
1871 
1872  tcross1 = -999, loncross1 = -999.;
1873 
1874  l1cfile->orb_dir = orb_dir_tot[ix3];
1875 
1876  day_mode = sunz_swt(ix3, ix4, hour_tot, day_tot, year_tot, orb_lat_tot, orb_lon_tot);
1877 
1878  if(l1cinput->verbose)cout << "day_mode at nswath #3 SWATH1--segment#2: " << day_mode << endl;
1879 
1880  status2 = ect_swt(l1cfile, ix3, ix4, orb_time_tot, orb_lat_tot, orb_lon_tot, orb_vel_tot,
1881  grn_vel_tot, tswt, latswt, lonswt, &tcross1, &loncross1, &mov1, &mgv1);
1882  if(l1cinput->verbose)cout << "nswath==3 -swath1---tcross equat crossing in (s)..swath#1.(segment #2).."
1883  << tcross1 << "loncross1..." << loncross1 << "orbit direction.0 is descending."
1884  << l1cfile->orb_dir << "mov1.." << mov1 << "mgv1.." << mgv1 << endl;
1885  } // end segment2
1886 
1887 
1888 
1889  if (status1 == 0 || status2 == 0 && day_mode == 1) {
1890  tmgv1 = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
1891 
1892  if(l1cinput->verbose)cout << "#gridlines..." << l1cfile->num_gridlines << "norbs.." << norbs << endl;
1893 
1894  for (int i = 0; i < l1cfile->num_gridlines; i++)
1895  tmgv1[i] = -1;
1896 
1897  swtime_swt2_segment(
1898  1, l1cinput, l1cfile, norbs, tswt, tcross1, mgv1,
1899  tmgv1); // number of gridlines may change here!!!!!!!!! not anymore 4000 and
1900  // assymetric around the equator so bina and binb not the same!!!
1901  num_gridlines = l1cfile->num_gridlines;
1902 
1903  double* tmgv1_segment = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
1904  tmgvf = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
1905 
1906  int ss = 0;
1907  for (int i = 0; i < num_gridlines; i++) {
1908  if (tmgv1[i] >= 0) {
1909  tmgv1_segment[ss] = tmgv1[i];
1910  ss++;
1911  }
1912  }
1913 
1914  if(l1cinput->verbose)cout << "#segment gridlines.ss counter.." << ss << "equal to num_gridlines..."
1915  << num_gridlines << endl;
1916 
1917  if (tswt != nullptr)
1918  free (tswt);
1919  tswt = nullptr;
1920  if (latswt != nullptr)
1921  free (latswt);
1922  latswt = nullptr;
1923  if (lonswt != nullptr)
1924  free (lonswt);
1925  lonswt = nullptr;
1926 
1927  if(l1cinput->verbose)cout << "number of across bins L1C grid...#.." << l1cfile->nbinx << "l1cfile->n_views..."
1928  << l1cfile->n_views << endl;
1929  lat_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
1930  lon_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
1931  alt = allocate2d_float(num_gridlines, l1cfile->nbinx);
1932 
1933  if(status1==0){
1934  ix_ini=ix1;
1935  ix_end=ix2;
1936  }
1937  else if(status2==0){
1938  ix_ini=ix3;
1939  ix_end=ix4;
1940  }
1941  orb_to_latlon(ix_ini,ix_end,num_gridlines, l1cfile->nbinx, orb_time_tot, posr_tot,
1942  velr_tot, mgv1, tmgv1_segment, tmgvf, lat_gd, lon_gd, alt,FirsTerrain);
1943 
1944 
1945  if (tmgv1 != nullptr)
1946  free (tmgv1);
1947  tmgv1 = nullptr;
1948  if (tmgv1_segment != nullptr)
1949  free (tmgv1_segment);
1950  tmgv1_segment = nullptr;
1951 
1952  l1cfile->num_gridlines = l1cfile->num_gridlines - 1;
1953 
1954  if (l1cinput->grantype == 1) {
1955  if(l1cinput->verbose)cout << "Processing swath--->" << endl;
1956  tswt_ini_sec = tfile_ini_sec + tmgvf[0];
1957  tswt_end_sec = tfile_ini_sec + tmgvf[num_gridlines - 2];
1958 
1959  create_time_swt(num_gridlines, tfile_ini_sec, tmgvf, tswt_ini_sec, tswt_end_sec,
1960  &tswt_ini, &tswt_ini_file, &tswt_mid,&tswt_end);
1961 
1962  l1cfile->tswt_ini = tswt_ini;
1963  l1cfile->tswt_mid = tswt_mid;
1964  l1cfile->tswt_end = tswt_end;
1965  l1cfile->tswt_ini_file = tswt_ini_file;
1966 
1967  create_SOCEA2(1, l1cinput, l1cfile, lat_gd, lon_gd, alt, tmgvf);
1968  } else if (l1cinput->grantype == 0) {
1969  if(l1cinput->verbose)cout << "Processing granules--->" << endl;
1970  numgran = 144 * 2;
1971  l1cfile->numgran = numgran;
1972  gd_per_gran = round(num_gridlines / 10);
1973  l1cfile->gd_per_gran = gd_per_gran;
1974 
1975  if(l1cinput->verbose)cout << "estimated # of granules to be processed..." << numgran << "gd_per_gran..."
1976  << gd_per_gran << "#gridlines.." << num_gridlines << endl;
1977  if(FirsTerrain) write_L1C_granule2(1, l1cfile, l1cinput, tmgvf, lat_gd, lon_gd, alt,orb_time_tot);
1978  } else {
1979  cout << "ERROR selecting grantype, must be 0: granules or 1: "
1980  "swath........................."
1981  << endl;
1982  exit(1);
1983  }
1984 
1985  //-----------------------------------------------------------------
1986  //-----------------------------------------------------------------
1987  if (lat_gd != nullptr)
1988  free (lat_gd);
1989  lat_gd = nullptr;
1990  if (lon_gd != nullptr)
1991  free (lon_gd);
1992  lon_gd = nullptr;
1993  if (tmgvf != nullptr)
1994  free (tmgvf);
1995  tmgvf = nullptr;
1996  if (alt != nullptr)
1997  free (alt);
1998  alt = nullptr;
1999  } // end tcross (all segments)
2000  else {
2001  if(l1cinput->verbose)
2002  {
2003  cout << "ERROR swath #1 does not cross the equator..." << endl;
2004  cout << "checking swath segment #2..." << endl;
2005  }
2006  }
2007  } // end day_mode
2008  else {
2009  if(l1cinput->verbose){
2010  cout << "day_mode==0 (nighttime)...." << day_mode << endl;
2011  cout << "checking swath segment #2..." << endl;
2012  }
2013  }
2014  } // end indexes x1 and x5 condition
2015 
2016  //-----------------------------------------------------------------------------------
2017  // nswath =3
2018  // swath 2
2019  //----------------------------------------------------------------------------------
2020 
2021  norbs = ix6 - ix5 + 1;
2022 
2023  day_mode = sunz_swt(ix5, ix6, hour_tot, day_tot, year_tot, orb_lat_tot, orb_lon_tot);
2024 
2025  if(l1cinput->verbose)cout << "day_mode at nswath #3 SWATH#2: " << day_mode << endl;
2026 
2027  if (day_mode == 1) {
2028  l1cfile->orb_dir = orb_dir_tot[ix5];
2029 
2030  tswt2 = (double*)calloc(norbs, sizeof(double));
2031  latswt2 = (double*)calloc(norbs, sizeof(double));
2032  lonswt2 = (double*)calloc(norbs, sizeof(double));
2033 
2034  status = ect_swt(l1cfile, ix5, ix6, orb_time_tot, orb_lat_tot, orb_lon_tot, orb_vel_tot,
2035  grn_vel_tot, tswt2, latswt2, lonswt2, &tcross2, &loncross2, &mov2, &mgv2);
2036  if(l1cinput->verbose)cout << "nswath==3 -swath2--tcross equat crossing in (s)..swath#2..." << tcross2 << "loncross2..."
2037  << loncross2 << "orbit direction.0 is descending." << l1cfile->orb_dir << "mov2.." << mov2
2038  << "mgv2.." << mgv2 << endl;
2039 
2040  if (latswt2 != nullptr)
2041  free (latswt2);
2042  latswt2 = nullptr;
2043  if (lonswt2 != nullptr)
2044  free (lonswt2);
2045  lonswt2 = nullptr;
2046 
2047 
2048 
2049  if (status == 0 && day_mode == 1) {//crossing and daytime
2050  tmgv2 = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
2051  tmgvf2 = (double*)calloc(l1cfile->num_gridlines, sizeof(double));
2052 
2053  swtime_swt2(2, l1cinput, l1cfile, norbs, tswt2, tcross2, mgv2, tmgv2);
2054  num_gridlines = l1cfile->num_gridlines;
2055 
2056  if (tswt2 != nullptr)
2057  free (tswt2);
2058  tswt2 = nullptr;
2059 
2060  if(l1cinput->verbose)cout << "number of across bins L1C grid...#.." << l1cfile->nbinx << "l1cfile->n_views..."
2061  << l1cfile->n_views << endl;
2062  lat_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
2063  lon_gd = allocate2d_float(num_gridlines, l1cfile->nbinx);
2064  alt = allocate2d_float(num_gridlines, l1cfile->nbinx);
2065 
2066  orb_to_latlon(ix5,ix6,num_gridlines, l1cfile->nbinx, orb_time_tot, posr_tot,
2067  velr_tot, mgv2, tmgv2, tmgvf2, lat_gd, lon_gd, alt,FirsTerrain);
2068 
2069 
2070  if (tmgv2 != nullptr)
2071  free (tmgv2);
2072  tmgv2 = nullptr;
2073 
2074  l1cfile->num_gridlines = l1cfile->num_gridlines - 1;
2075 
2076  if (l1cinput->grantype == 1) {
2077  tswt_ini_sec = tfile_ini_sec + tmgvf2[0];
2078  tswt_end_sec = tfile_ini_sec + tmgvf2[num_gridlines - 2];
2079 
2080  create_time_swt(num_gridlines, tfile_ini_sec, tmgvf2, tswt_ini_sec, tswt_end_sec,
2081  &tswt_ini, &tswt_ini_file, &tswt_mid,&tswt_end);
2082 
2083  l1cfile->tswt_ini = tswt_ini;
2084  l1cfile->tswt_mid = tswt_mid;
2085  l1cfile->tswt_end = tswt_end;
2086  l1cfile->tswt_ini_file = tswt_ini_file;
2087  create_SOCEA2(2, l1cinput, l1cfile, lat_gd, lon_gd, alt, tmgvf2);
2088  } else if (l1cinput->grantype == 0) {
2089  deltasec = tmgvf2[num_gridlines - 2] - tmgvf2[0] + 1;
2090  if(l1cinput->verbose)cout << "deltasec..swath." << deltasec << endl;
2091  numgran = 144 * 2;
2092 
2093  l1cfile->numgran = numgran;
2094  gd_per_gran = round(num_gridlines / 10); // 10 granules per half orbit
2095  l1cfile->gd_per_gran = gd_per_gran;
2096 
2097  if(l1cinput->verbose)cout << "estimated # of granules to be processed..." << numgran << "gd_per_gran..."
2098  << gd_per_gran << "#gridlines.." << num_gridlines << endl;
2099 
2100  if(FirsTerrain) write_L1C_granule2(2, l1cfile, l1cinput, tmgvf2, lat_gd, lon_gd, alt,orb_time_tot);
2101  } else {
2102  cout << "ERROR selecting grantype, must be 0: granules or 1: "
2103  "swath........................."
2104  << endl;
2105  exit(1);
2106  }
2107  if (lat_gd != nullptr)
2108  free (lat_gd);
2109  lat_gd = nullptr;
2110  if (lon_gd != nullptr)
2111  free (lon_gd);
2112  lon_gd = nullptr;
2113  if (tmgvf2 != nullptr)
2114  free (tmgvf2);
2115  tmgvf2 = nullptr;
2116  if (alt != nullptr)
2117  free (alt);
2118  alt = nullptr;
2119  } // status 0
2120  else {
2121  if(l1cinput->verbose) cout << "ERROR swath #2 does not cross the equator..NO L1C grid for swath #2" << endl;
2122  }
2123 
2124  } // end day_mode==1
2125  else {
2126  cout << "day_mode = 0 nightime...no L1C grid produced--exit (swath#2 ---nswath#3).............."
2127  << day_mode << endl;
2128  // exit(1);
2129  }
2130  } // end nswath =3
2131  } // end big loop files HKT batch processing
2132 
2133  free (year_tot);
2134  free (day_tot);
2135  free (hour_tot);
2136  free (orb_time_tot);
2137  free (orb_lon_tot);
2138  free (orb_lat_tot);
2139  free (orb_dir_tot);
2140  delete[] (posr_tot);
2141  delete[] (velr_tot);
2142  free (orb_vel_tot);
2143  free (grn_vel_tot);
2144 
2145  return 0;
2146 }
2147 
2148 
2149 int L1C::create_SOCEA2(int swtd, L1C_input* l1cinput, l1c_filehandle* l1cfile, float** lat_gd, float** lon_gd,
2150  float** altitude, double* time_nad) {
2151  string tswt_ini, tswt_end, tswt_ini_file, date_created;
2152  int asc_mode = -1;
2153  std::string extstr = ".nc";
2154  string dirstr, prodstr;
2155  const char* filename_lt;
2156  int Nwest = -1, Neast = -1, Ngring = -1, midix = -1, dp = -1;
2157  int p, ix = -1;
2158  float latemp = -1, lontemp1 = -1, lontemp2 = -1, dlat_gd = -1, dlon_gd = -1, dlat20 = -1, dlon20 = -1,
2159  lon360 = -1;
2160  int re = -1, rw = -1;
2161  int32_t iyear=0,iday=0,msec=0;
2162  string datezerotime,yearstr,monstr,daystr;
2163  int16_t syear, smon, sday;
2164  double tdate_ini,secs;
2165  int twodays=0;
2166 
2167  if(l1cinput->verbose)cout << "Creating SOCEA proj for the whole orbit......." << endl;
2168 
2169  if (l1cinput->sensor == 34) {
2170  //"SPEXONE";
2171  l1cfile->nbinx = 29;
2172  midix = 14;
2173  // NVIEWS=l1cfile->n_views;
2174  // NBANDS=400;
2175  } else if (l1cinput->sensor == 30 || l1cinput->sensor == 31) { // OCIS/OCIS
2176  //"OCI";
2177  l1cfile->nbinx = 519;
2178  midix = 259;
2179  // NVIEWS=2;
2180  // NBANDS=249;
2181  } else if (l1cinput->sensor == 35) {
2182  //"HARP2";
2183  l1cfile->nbinx = 457;
2184  midix = 228;
2185  // NVIEWS=l1cfile->n_views;
2186  // NBANDS=4;
2187  } else {
2188  cout << "sensor by default is OCI option 2....." << endl;
2189  //"OCI";
2190  l1cfile->nbinx = 519;
2191  midix = 259;
2192  // NVIEWS=2;
2193  // NBANDS=249;
2194  }
2195 
2196  for(int gd=0;gd<l1cfile->num_gridlines;gd++) {
2197  if(time_nad[gd]<0)
2198  twodays=1;
2199  }
2200  // time nadir in seconds
2201 
2202  tswt_ini_file = l1cfile->tswt_ini_file;
2203  prodstr = "PACE." + tswt_ini_file + ".L1C" + extstr;
2204 
2205  l1cfile->gridname = prodstr.c_str();
2206 
2207  filename_lt = prodstr.c_str();
2208  char* gridchar = strdup(filename_lt);
2209  string l1c_str = filename_lt;
2210 
2211  NcFile* nc_output;
2212  try {
2213  nc_output = new NcFile(filename_lt, NcFile::replace);
2214  } catch (NcException& e) {
2215  e.what();
2216  cerr << "l1cgen l1c_pflag= 5 : producing L1C grid: " + l1c_str << endl;
2217  exit(1);
2218  }
2219 
2220  bin_str binl1c;
2221  binl1c.sensor = l1cinput->sensor;
2222  binl1c.date_mid_grid=l1cfile->tswt_mid;
2223  meta_l1c_grid(gridchar, &binl1c, l1cfile->num_gridlines, nc_output);
2224  // gobal attrs--
2225  asc_mode = l1cfile->orb_dir;
2226  if (asc_mode == 1)
2227  dirstr = "Ascending";
2228  else if (asc_mode == 0)
2229  dirstr = "Descending";
2230  tswt_ini = l1cfile->tswt_ini;
2231  tswt_end = l1cfile->tswt_end;
2232 
2233  nc_output->putAtt("processing_version", l1cinput->pversion);
2234  nc_output->putAtt("history", l1cinput->history);
2235  nc_output->putAtt("product_name", prodstr);
2236  nc_output->putAtt("startDirection", dirstr);
2237  nc_output->putAtt("endDirection", dirstr);
2238  nc_output->putAtt("time_coverage_start", tswt_ini);
2239  nc_output->putAtt("time_coverage_end", tswt_end);
2240 
2241  isodate2ydmsec((char *) l1cfile->tswt_mid.c_str(), &iyear, &iday, &msec);
2242  double dist_es=esdist_(&iyear,&iday,&msec);
2243  nc_output->putAtt("sun_earth_distance", ncFloat, dist_es);
2244  if(l1cinput->verbose)cout << "sun_earth_distance -- mid gridline" <<dist_es<<endl;
2245 
2246  tdate_ini = isodate2unix(l1cinput->start_time);
2247  unix2ymds(tdate_ini, &syear, &smon, &sday, &secs);
2248  yearstr=std::to_string(syear);
2249  monstr=std::to_string(smon);
2250  if(monstr.size() == 1)
2251  monstr = "0" + monstr;
2252  daystr=std::to_string(sday);
2253  if(daystr.size() == 1)
2254  daystr = "0" + daystr;
2255  datezerotime="seconds since " + yearstr + "-" + monstr + "-" + daystr;
2256  // vars
2257  NcGroup ba_grp = nc_output->getGroup("bin_attributes");
2258  NcVar v1 = ba_grp.getVar("nadir_view_time");
2259 
2260  if(twodays) {
2261  for(int gd=0;gd<l1cfile->num_gridlines;gd++) {
2262  time_nad[gd]+=3600*24;
2263  }
2264  }
2265 
2266  v1.putVar(&time_nad[0]);
2267  v1.putAtt("units",datezerotime);
2268  // v1=ba_grp.getVar("view_time_offsets");
2269  // v1.putVar(&time_off[0][0][0]);
2270  NcGroup geo_grp = nc_output->getGroup("geolocation_data");
2271  v1 = geo_grp.getVar("latitude");
2272  v1.putVar(&lat_gd[0][0]);
2273  v1 = geo_grp.getVar("longitude");
2274  v1.putVar(&lon_gd[0][0]);
2275  v1 = geo_grp.getVar("height");
2276  v1.putVar(&altitude[0][0]);
2277 
2278  // GRING-------
2279  // determine the number of GCpoint indexes
2280  // default 6 for the swath sides + number of coordinates every 20 degrees latitude
2281 
2282  // ascending pass
2283  if (l1cfile->orb_dir == 1) {
2284  Nwest = round((lat_gd[l1cfile->num_gridlines - 1][0] - lat_gd[0][0]) / 20);
2285  Neast = round(
2286  (lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1] - lat_gd[0][l1cfile->nbinx - 1]) / 20);
2287  } else // descending
2288  {
2289  Neast = (round(lat_gd[0][0] - lat_gd[l1cfile->num_gridlines - 1][0]) / 20);
2290  Nwest = round(
2291  (lat_gd[0][l1cfile->nbinx - 1] - lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1]) / 20);
2292  }
2293 
2294  // first NGring estimate----
2295  Ngring = Nwest + Neast + 6;
2296 
2297  if (Ngring > 0) {
2298  float* latarr = (float*)calloc(Ngring, sizeof(float));
2299  float* lonarr = (float*)calloc(Ngring, sizeof(float));
2300  int* p_west = (int*)calloc(Ngring, sizeof(int));
2301  int* p_east = (int*)calloc(Ngring, sizeof(int));
2302 
2303  // corners--counterclockwise and ascending pass
2304  if (l1cfile->orb_dir == 1) {
2305  if (Ngring == 6) {
2306  latarr[0] = lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2307  latarr[1] = lat_gd[l1cfile->num_gridlines - 1][midix];
2308  latarr[2] = lat_gd[l1cfile->num_gridlines - 1][0];
2309  latarr[3] = lat_gd[0][0];
2310  latarr[4] = lat_gd[0][midix];
2311  latarr[5] = lat_gd[0][l1cfile->nbinx - 1];
2312 
2313  lonarr[0] = lon_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2314  lonarr[1] = lon_gd[l1cfile->num_gridlines - 1][midix];
2315  lonarr[2] = lon_gd[l1cfile->num_gridlines - 1][0];
2316  lonarr[3] = lon_gd[0][0];
2317  lonarr[4] = lon_gd[0][midix];
2318  lonarr[5] = lon_gd[0][l1cfile->nbinx - 1];
2319  } else {
2320  latarr[0] = lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2321  latarr[1] = lat_gd[l1cfile->num_gridlines - 1][midix];
2322  latarr[2] = lat_gd[l1cfile->num_gridlines - 1][0];
2323 
2324  lonarr[0] = lon_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2325  lonarr[1] = lon_gd[l1cfile->num_gridlines - 1][midix];
2326  lonarr[2] = lon_gd[l1cfile->num_gridlines - 1][0];
2327 
2328  latemp = latarr[2];
2329  latemp -= 20;
2330  p = 1;
2331  rw = 0;
2332  // west side
2333  while (latemp > lat_gd[0][0]) {
2334  latarr[2 + p] = latemp;
2335  if(l1cinput->verbose)cout << "p west--" << p << "point# = " << 3 + p << "lat20 = " << latemp << endl;
2336  p_west[rw] = 3 + p;
2337  latemp -= 20;
2338  p++;
2339  rw++;
2340  }
2341 
2342  p--;
2343 
2344  latarr[2 + p + 1] = lat_gd[0][0];
2345  latarr[2 + p + 2] = lat_gd[0][midix];
2346  latarr[2 + p + 3] = lat_gd[0][l1cfile->nbinx - 1];
2347 
2348  lonarr[2 + p + 1] = lon_gd[0][0];
2349  lonarr[2 + p + 2] = lon_gd[0][midix];
2350  lonarr[2 + p + 3] = lon_gd[0][l1cfile->nbinx - 1];
2351 
2352  latemp = latarr[2 + p + 3];
2353  latemp += 20;
2354  p++;
2355 
2356  int c = 1;
2357  re = 0;
2358  // east side
2359  while (latemp < lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1]) {
2360  latarr[5 + p] = latemp;
2361  if(l1cinput->verbose)cout << "p east--" << c << "point# = " << 6 + p << "lat20 = " << latemp << endl;
2362  p_east[re] = 6 + p;
2363  latemp += 20;
2364  p++;
2365  c++;
2366  re++;
2367  }
2368 
2369  p--;
2370 
2371  if(l1cinput->verbose){cout << "mid points west.." << rw << "mid points east.." << re << endl;
2372  cout << "# points in GRING = " << 6 + p << "# mid points west--" << rw
2373  << "# mid points east--" << re << endl;
2374  }
2375  if (Ngring == 6)
2376  dp = 6;
2377  else
2378  dp = rw + re + 6;
2379  // west
2380  for (int i = 0; i < rw; i++) {
2381  ix = p_west[i] - 1;
2382  if(l1cinput->verbose) cout<<"lat:"<<latarr[ix]<<"pwest--"<<p_west[i]<<endl;
2383  for (int row = 0; row < l1cfile->num_gridlines - 1; row++) {
2384  if (latarr[ix] > lat_gd[row][0] && latarr[ix] <= lat_gd[row + 1][0]) {
2385  if(l1cinput->verbose) cout<<"mid point west# = "<<i+1<<"found index between #row= "<<row+1<<"and row ="<<row+2<<endl;
2386 
2387  if (lon_gd[row][0] < 0.)
2388  lontemp1 = lon_gd[row][0] + 360;
2389  else
2390  lontemp1 = lon_gd[row][0];
2391  if (lon_gd[row + 1][0] < 0.)
2392  lontemp2 = lon_gd[row + 1][0] + 360;
2393  else
2394  lontemp2 = lon_gd[row + 1][0];
2395 
2396  dlat_gd = abs(lat_gd[row + 1][0] - lat_gd[row][0]);
2397  dlon_gd = abs(lontemp1 - lontemp2);
2398  dlat20 = abs(latarr[ix] - lat_gd[row + 1][0]);
2399  dlon20 = dlat20 * dlon_gd / dlat_gd;
2400 
2401  if (lontemp1 > lontemp2)
2402  lon360 = lontemp2 + dlon20;
2403  else
2404  lon360 = lontemp2 - dlon20;
2405  if (lon360 > 180)
2406  lon360 = lon360 - 360.;
2407  lonarr[ix] = lon360;
2408  if(l1cinput->verbose) cout<<"lon_gd row+1.."<<lon_gd[row+1][0]<<"lon_gd row.."<<lon_gd[row][0]<<"lon cring.."<<lonarr[ix]<<endl;
2409  break;
2410  }
2411  }
2412  }
2413  // east
2414  for (int i = 0; i < re; i++) {
2415  ix = p_east[i] - 1;
2416  if(l1cinput->verbose)cout << "lat:" << latarr[ix] << "peast--" << p_east[i] << endl;
2417  for (int row = 0; row < l1cfile->num_gridlines - 1; row++) {
2418  if (latarr[ix] > lat_gd[row][l1cfile->nbinx - 1] &&
2419  latarr[ix] <= lat_gd[row + 1][l1cfile->nbinx - 1]) {
2420  if(l1cinput->verbose) cout<<"mid point east# = "<<i+1<<"found index between #row= "<<row+1<<"and row ="<<row+2<<endl;
2421 
2422  if (lon_gd[row][l1cfile->nbinx - 1] < 0.)
2423  lontemp1 = lon_gd[row][l1cfile->nbinx - 1] + 360;
2424  else
2425  lontemp1 = lon_gd[row][l1cfile->nbinx - 1];
2426  if (lon_gd[row + 1][l1cfile->nbinx - 1] < 0.)
2427  lontemp2 = lon_gd[row + 1][l1cfile->nbinx - 1] + 360;
2428  else
2429  lontemp2 = lon_gd[row + 1][l1cfile->nbinx - 1];
2430 
2431  dlat_gd =
2432  abs(lat_gd[row + 1][l1cfile->nbinx - 1] - lat_gd[row][l1cfile->nbinx - 1]);
2433  dlon_gd = abs(lontemp1 - lontemp2);
2434  dlat20 = abs(latarr[ix] - lat_gd[row][l1cfile->nbinx - 1]);
2435  dlon20 = dlat20 * dlon_gd / dlat_gd;
2436  if (lontemp1 > lontemp2)
2437  lon360 = lontemp1 - dlon20;
2438  else
2439  lon360 = lontemp1 + dlon20;
2440  if (lon360 > 180)
2441  lon360 = lon360 - 360.;
2442  lonarr[ix] = lon360;
2443  if(l1cinput->verbose) cout<<"lon_gd row+1.."<<lon_gd[row+1][l1cfile->nbinx-1]<<"lon_gd row.."<<lon_gd[row][l1cfile->nbinx-1]<<"lon cring.."<<lonarr[ix]<<endl;
2444 
2445  break;
2446  }
2447  }
2448  }
2449  }
2450  } else // descending orb
2451  {
2452  if (Ngring == 6) {
2453  latarr[0] = lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2454  latarr[1] = lat_gd[l1cfile->num_gridlines - 1][midix];
2455  latarr[2] = lat_gd[l1cfile->num_gridlines - 1][0];
2456  latarr[3] = lat_gd[0][0];
2457  latarr[4] = lat_gd[0][midix];
2458  latarr[5] = lat_gd[0][l1cfile->nbinx - 1];
2459 
2460  lonarr[0] = lon_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2461  lonarr[1] = lon_gd[l1cfile->num_gridlines - 1][midix];
2462  lonarr[2] = lon_gd[l1cfile->num_gridlines - 1][0];
2463  lonarr[3] = lon_gd[0][0];
2464  lonarr[4] = lon_gd[0][midix];
2465  lonarr[5] = lon_gd[0][l1cfile->nbinx - 1];
2466  } else {
2467  latarr[0] = lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2468  latarr[1] = lat_gd[l1cfile->num_gridlines - 1][midix];
2469  latarr[2] = lat_gd[l1cfile->num_gridlines - 1][0];
2470 
2471  lonarr[0] = lon_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1];
2472  lonarr[1] = lon_gd[l1cfile->num_gridlines - 1][midix];
2473  lonarr[2] = lon_gd[l1cfile->num_gridlines - 1][0];
2474 
2475  latemp = latarr[2];
2476  latemp += 20;
2477  p = 1;
2478  int rw = 0;
2479  // west side
2480  while (latemp < lat_gd[0][0]) {
2481  latarr[2 + p] = latemp;
2482  if(l1cinput->verbose)cout << "p west--" << p << "point# = " << 3 + p << "lat20 = " << latemp << endl;
2483  p_west[rw] = 3 + p;
2484  latemp += 20;
2485  p++;
2486  rw++;
2487  }
2488 
2489  p--;
2490 
2491  latarr[2 + p + 1] = lat_gd[0][0];
2492  latarr[2 + p + 2] = lat_gd[0][midix];
2493  latarr[2 + p + 3] = lat_gd[0][l1cfile->nbinx - 1];
2494 
2495  lonarr[2 + p + 1] = lon_gd[0][0];
2496  lonarr[2 + p + 2] = lon_gd[0][midix];
2497  lonarr[2 + p + 3] = lon_gd[0][l1cfile->nbinx - 1];
2498  latemp = latarr[2 + p + 3];
2499  latemp -= 20;
2500  p++;
2501 
2502  int c = 1;
2503  int re = 0;
2504  // east side
2505  while (latemp > lat_gd[l1cfile->num_gridlines - 1][l1cfile->nbinx - 1]) {
2506  latarr[5 + p] = latemp;
2507  if(l1cinput->verbose)cout << "p east--" << c << "point# = " << 6 + p << "lat20 = " << latemp << endl;
2508  p_east[re] = 6 + p;
2509  latemp -= 20;
2510  p++;
2511  c++;
2512  re++;
2513  }
2514 
2515  p--;
2516 
2517  if(l1cinput->verbose){
2518  cout << "mid points west.." << rw << "mid points east.." << re << endl;
2519  cout << "# points in GRING = " << 6 + p << "# mid points west--" << rw
2520  << "# mid points east--" << re << endl;
2521  }
2522  if (Ngring == 6)
2523  dp = 6;
2524  else
2525  dp = rw + re + 6;
2526 
2527  // west side
2528  for (int i = 0; i < rw; i++) {
2529  ix = p_west[i] - 1;
2530  if(l1cinput->verbose) cout<<"lat:"<<latarr[ix]<<"pwest--"<<p_west[i]<<endl;
2531  for (int row = 0; row < l1cfile->num_gridlines - 1; row++) {
2532  if (latarr[ix] <= lat_gd[row][0] && latarr[ix] > lat_gd[row + 1][0]) {
2533  if(l1cinput->verbose) cout<<"mid point west# = "<<i+1<<"found index between #row= "<<row+1<<"and row ="<<row+2<<endl;
2534  if (lon_gd[row][0] < 0.)
2535  lontemp1 = lon_gd[row][0] + 360;
2536  else
2537  lontemp1 = lon_gd[row][0];
2538  if (lon_gd[row + 1][0] < 0.)
2539  lontemp2 = lon_gd[row + 1][0] + 360;
2540  else
2541  lontemp2 = lon_gd[row + 1][0];
2542 
2543  dlat_gd = abs(lat_gd[row][0] - lat_gd[row + 1][0]);
2544  dlon_gd = abs(lontemp1 - lontemp2);
2545  dlat20 = abs(latarr[ix] - lat_gd[row + 1][0]);
2546  dlon20 = dlat20 * dlon_gd / dlat_gd;
2547 
2548  if (lontemp1 > lontemp2)
2549  lon360 = lontemp2 + dlon20;
2550  else
2551  lon360 = lontemp2 - dlon20;
2552  if (lon360 > 180)
2553  lon360 = lon360 - 360.;
2554  lonarr[ix] = lon360;
2555  if(l1cinput->verbose) cout<<"lon_gd row+1.."<<lon_gd[row+1][0]<<"lon_gd row.."<<lon_gd[row][0]<<"lon cring.."<<lonarr[ix]<<endl;
2556  break;
2557  }
2558  }
2559  }
2560  // east side
2561  for (int i = 0; i < re; i++) {
2562  ix = p_east[i] - 1;
2563  if(l1cinput->verbose) cout<<"lat:"<<latarr[ix]<<"peast--"<<p_east[i]<<endl;
2564  for (int row = 0; row < l1cfile->num_gridlines - 1; row++) {
2565  if (latarr[ix] <= lat_gd[row][l1cfile->nbinx - 1] &&
2566  latarr[ix] > lat_gd[row + 1][l1cfile->nbinx - 1]) {
2567  if(l1cinput->verbose)cout<<"mid point east# = "<<i+1<<"found index between #row= "<<row+1<<"and row ="<<row+2<<endl;
2568  if (lon_gd[row][l1cfile->nbinx - 1] < 0.)
2569  lontemp1 = lon_gd[row][l1cfile->nbinx - 1] + 360;
2570  else
2571  lontemp1 = lon_gd[row][l1cfile->nbinx - 1];
2572  if (lon_gd[row + 1][l1cfile->nbinx - 1] < 0.)
2573  lontemp2 = lon_gd[row + 1][l1cfile->nbinx - 1] + 360;
2574  else
2575  lontemp2 = lon_gd[row + 1][l1cfile->nbinx - 1];
2576 
2577  dlat_gd =
2578  abs(lat_gd[row][l1cfile->nbinx - 1] - lat_gd[row + 1][l1cfile->nbinx - 1]);
2579  dlon_gd = abs(lontemp1 - lontemp2);
2580  dlat20 = abs(latarr[ix] - lat_gd[row][l1cfile->nbinx - 1]);
2581  dlon20 = dlat20 * dlon_gd / dlat_gd;
2582  if (lontemp1 > lontemp2)
2583  lon360 = lontemp1 - dlon20;
2584  else
2585  lon360 = lontemp1 + dlon20;
2586  if (lon360 > 180)
2587  lon360 = lon360 - 360.;
2588  lonarr[ix] = lon360;
2589  if(l1cinput->verbose) cout<<"lon_gd row+1.."<<lon_gd[row+1][l1cfile->nbinx-1]<<"lon_gd row.."<<lon_gd[row][l1cfile->nbinx-1]<<"lon cring.."<<lonarr[ix]<<endl;
2590 
2591  break;
2592  }
2593  }
2594  }
2595 
2596  } // else GRING>6
2597 
2598  } // end descending
2599 
2600  // sequence- GRING
2601 
2602  string onelat,onelon,onecoor,firstcoor;
2603  //GRing params
2604  string gring_polygon,gring_crs,gring_latmin,gring_latmax,gring_lonmin,gring_lonmax;
2605 
2606  for (int s = 0; s < dp; s++) {
2607  onelat=to_string(latarr[s]);
2608  onelon=to_string(lonarr[s]);
2609  if(s==0)
2610  {
2611  onecoor="POLYGON(("+onelat+" "+onelon+",";
2612  firstcoor=onelat+" "+onelon;
2613  }
2614  else onecoor=onelat+" "+onelon+",";
2615  gring_polygon+=onecoor;
2616  }
2617  //last coor
2618  gring_polygon+=firstcoor+"))";
2619  float latpmin=999,latpmax=-999,lonpmin=999,lonpmax=-999;
2620  for(int row=0;row<l1cfile->num_gridlines;row++)
2621  {
2622  for(int col=0;col<l1cfile->nbinx;col++)
2623  {
2624  if(lat_gd[row][col]<latpmin && lat_gd[row][col]!=BAD_FLT) latpmin=lat_gd[row][col];
2625  if(lat_gd[row][col]>latpmax && lat_gd[row][col]!=BAD_FLT) latpmax=lat_gd[row][col];
2626  if(lon_gd[row][col]<lonpmin && lon_gd[row][col]!=BAD_FLT) lonpmin=lon_gd[row][col];
2627  if(lon_gd[row][col]>lonpmax && lon_gd[row][col]!=BAD_FLT) lonpmax=lon_gd[row][col];
2628  }}
2629 
2630  gring_crs="EPSG:4326";
2631  nc_output->putAtt("geospatial_bounds", gring_polygon);
2632  nc_output->putAtt("geospatial_bounds_crs", gring_crs);
2633 
2634  //algo to find min/max coordinates inside the gring polygon
2635  nc_output->putAtt("geospatial_lat_min", ncFloat, latpmin);
2636  nc_output->putAtt("geospatial_lat_max", ncFloat, latpmax);
2637  nc_output->putAtt("geospatial_lon_min", ncFloat, lonpmin);
2638  nc_output->putAtt("geospatial_lon_max", ncFloat, lonpmax);
2639 
2640  free (latarr);
2641  free (lonarr);
2642  free (p_west);
2643  free (p_east);
2644 
2645  } else {
2646  cout << "ERROR EXTRACTING GRING coordinates!!-----" << endl;
2647  exit(1);
2648  }
2649 
2650  nc_output->close();
2651  return 0;
2652 }
2653 
2654 // fred/me implementation OK 3/22/2023
2655 int L1C::search_l1cgen(L1C_input* l1cinput, l1c_str* l1cstr, l1c_filehandle* l1cfile, short** gdindex) {
2656  int32_t num_gridlines, nbinx;
2657  float gnvm;
2658  float gnvec[3], gvec[3], bvec[3];
2659  int flag_out = -1;
2660  size_t pix;
2661  int32_t i;
2662  size_t num_pixels;
2663  int irow = -1, col = -1;
2664  float bmcm, bm = 100;
2665  double db;
2666  float c1, c2, c3;
2667  double fudge = 0.00001, dotprod, dot_firstline, dot_lastline;
2668  flag_out = -1;
2669 
2670  num_gridlines = l1cfile->num_gridlines;
2671  nbinx = l1cfile->nbinx;
2672 
2673  num_pixels = l1cfile->npix;
2674 
2675  db = (l1cinput->grid_resolution) / 6371 / 2; // Half of bin size in radians
2676 
2677  // big loop
2678  // dot product
2679  for (pix = 0; pix < num_pixels; pix++) {
2680  bvec[0] = cos(l1cstr->lonpix[pix] * M_PI / 180) * cos(l1cstr->latpix[pix] * M_PI / 180);
2681  bvec[1] = sin(l1cstr->lonpix[pix] * M_PI / 180) * cos(l1cstr->latpix[pix] * M_PI / 180);
2682  bvec[2] = sin(l1cstr->latpix[pix] * M_PI / 180);
2683 
2684  for (i = 0; i < num_gridlines; i++) {
2685  if (l1cstr->latpix[pix] > 90 || l1cstr->latpix[pix] < -90 || l1cstr->lonpix[pix] < -180 ||
2686  l1cstr->lonpix[pix] > 180) {
2687  /* cout << "ERROR in search_l1cgen --latitude longitude pixel out of the boundaries.."
2688  << "latpix>90 or <-90.." << l1cstr->latpix[pix] << "lonpix>180 or <-180.."
2689  << l1cstr->lonpix[pix] << endl;*/
2690  return 110;
2691  }
2692  // normal vectors for L1C rows
2693  if (l1cfile->lat_gd[i][nbinx - 1] > 90 || l1cfile->lat_gd[i][nbinx - 1] < -90 ||
2694  l1cfile->lon_gd[i][nbinx - 1] < -180 || l1cfile->lon_gd[i][nbinx - 1] > 180) {
2695  cout << "lat lon out of the boundaries.." << endl;
2696  exit(1);
2697  }
2698  if (l1cfile->lat_gd[i][0] > 90 || l1cfile->lat_gd[i][0] < -90 || l1cfile->lon_gd[i][0] < -180 ||
2699  l1cfile->lon_gd[i][0] > 180) {
2700  cout << "lat lon out of the boundaries.." << endl;
2701  exit(1);
2702  }
2703 
2704  gnvec[0] = sin(l1cfile->lon_gd[i][nbinx - 1] * M_PI / 180) *
2705  cos(l1cfile->lat_gd[i][nbinx - 1] * M_PI / 180) *
2706  sin(l1cfile->lat_gd[i][0] * M_PI / 180) -
2707  sin(l1cfile->lat_gd[i][nbinx - 1] * M_PI / 180) *
2708  sin(l1cfile->lon_gd[i][0] * M_PI / 180) * cos(l1cfile->lat_gd[i][0] * M_PI / 180);
2709  gnvec[1] = sin(l1cfile->lat_gd[i][nbinx - 1] * M_PI / 180) *
2710  cos(l1cfile->lon_gd[i][0] * M_PI / 180) * cos(l1cfile->lat_gd[i][0] * M_PI / 180) -
2711  cos(l1cfile->lon_gd[i][nbinx - 1] * M_PI / 180) *
2712  cos(l1cfile->lat_gd[i][nbinx - 1] * M_PI / 180) *
2713  sin(l1cfile->lat_gd[i][0] * M_PI / 180);
2714  gnvec[2] = cos(l1cfile->lon_gd[i][nbinx - 1] * M_PI / 180) *
2715  cos(l1cfile->lat_gd[i][nbinx - 1] * M_PI / 180) *
2716  sin(l1cfile->lon_gd[i][0] * M_PI / 180) * cos(l1cfile->lat_gd[i][0] * M_PI / 180) -
2717  sin(l1cfile->lon_gd[i][nbinx - 1] * M_PI / 180) *
2718  cos(l1cfile->lat_gd[i][nbinx - 1] * M_PI / 180) *
2719  cos(l1cfile->lon_gd[i][0] * M_PI / 180) * cos(l1cfile->lat_gd[i][0] * M_PI / 180);
2720 
2721  // vector norm
2722  gnvm = sqrt(gnvec[0] * gnvec[0] + gnvec[1] * gnvec[1] + gnvec[2] * gnvec[2]);
2723  if (isnan(gnvm) == 1) {
2724  cout << "NAN value for gnvm.." << endl;
2725  exit(1);
2726  }
2727  if (gnvm == 0) {
2728  cout << "ERROR gnvm == 0--- WE CANT NORMALIZE..." << endl;
2729  exit(1);
2730  }
2731  // normalization
2732  gnvec[0] = gnvec[0] / gnvm;
2733  gnvec[1] = gnvec[1] / gnvm;
2734  gnvec[2] = gnvec[2] / gnvm;
2735 
2736  // for each pixels
2737  // dot prod, orbital normaliz and transposed by pixel vector
2738  dotprod = gnvec[0] * bvec[0] + gnvec[1] * bvec[1] + gnvec[2] * bvec[2];
2739 
2740  if (i == 0) {
2741  dot_firstline = dotprod;
2742  }
2743  if (i == num_gridlines - 1) {
2744  dot_lastline = dotprod;
2745  }
2746 
2747  if (dotprod - fudge <= db && dotprod + fudge > -db) {
2748  gdindex[pix][0] = i + 1; // first found
2749  }
2750 
2751  } // end lines
2752  // for each pixels
2753  if (dot_firstline <= db && dot_lastline > -db) {
2754 
2755  // find column
2756  irow = gdindex[pix][0] - 1;
2757  if (irow < 0) {
2758  cout << "ERROR icol in search_l1c..."
2759  << "at pix#.." << pix + 1 << "and irow#.." << irow << endl;
2760  exit(1);
2761  }
2762 
2763  for (int j = 0; j < nbinx; j++) {
2764  gvec[0] =
2765  cos(l1cfile->lon_gd[irow][j] * M_PI / 180) * cos(l1cfile->lat_gd[irow][j] * M_PI / 180);
2766  gvec[1] =
2767  sin(l1cfile->lon_gd[irow][j] * M_PI / 180) * cos(l1cfile->lat_gd[irow][j] * M_PI / 180);
2768  gvec[2] = sin(l1cfile->lat_gd[irow][j] * M_PI / 180);
2769 
2770  c1 = bvec[0] - gvec[0];
2771  c2 = bvec[1] - gvec[1];
2772  c3 = bvec[2] - gvec[2];
2773 
2774  bmcm = sqrt(c1 * c1 + c2 * c2 + c3 * c3);
2775  if (bmcm < bm) {
2776  bm = bmcm;
2777  col = j + 1;
2778  }
2779  }
2780  if (col < 1) {
2781  cout << "ERROR col in search_l1c..."
2782  << "at pix#.." << pix + 1 << "and row#.." << irow + 1 << endl;
2783  exit(1);
2784  }
2785 
2786  gdindex[pix][1] = col;
2787  bm = 100;
2788  col = -1;
2789  } else {
2790  gdindex[pix][0] = -1; // actually these are row/col not irow/icol like in my searching algos!!!
2791  gdindex[pix][1] = -1;
2792  }
2793 
2794  } // end pixels
2795 
2796  flag_out = -1;
2797  for (pix = 0; pix < num_pixels; pix++) {
2798  if (gdindex[pix][0] > 0 && gdindex[pix][1] > 0) {
2799  flag_out = 0;
2800  } else {
2801  if(l1cinput->verbose) cout<<"THIS LINE WILL BE SKIPPED -- NO PIXELS BINNED.............."<<endl;
2802  }
2803  }
2804 
2805  if (flag_out == 0)
2806  return 0;
2807  else
2808  return 1;
2809 }
2810 
2811 
2812 // version after Fred proj
2813 int32_t L1C::openL1Cgrid3(l1c_str* l1cstr, l1c_filehandle* l1cfile, L1C_input* l1cinput) {
2814  int32_t num_gridlines, nbinx; // number of gridlines to be processed
2815  std::string str;
2816  const char* ptstr;
2817  int32_t NY = -1, NX = -1;
2818  std::string ifile_str;
2819  string gridname, azeast_name;
2820  string tswt_ini_file;
2821  std::string fname_out, senstr, monstr, daystr, yearstr, prodstr, gdstr, swtstr, extstr, granstr,
2822  timestr, azstr, missionstr, ofilestr;
2823  NcFile* nc_l1c;
2824 
2825  // char tmp[256];
2826  // getcwd(tmp, 256);
2827  // pathstr=tmp;
2828 
2829  if (l1cinput->projection == 0) { // socea
2830  gridname = l1cinput->l1c_grid;
2831  if(l1cinput->verbose)cout << "projection type is SOCEA!!!........for grid file..:" << gridname << endl;
2832  }
2833  if (l1cinput->projection == 1) { // fixed bearing projection ----
2834  gridname = "/accounts/mamontes/images/OCIS/sean/out/FB_L1Cgrid.nc";
2835  }
2836 
2837  if(l1cinput->verbose)cout << "gridname.." << gridname << endl;
2838  ptstr = gridname.c_str();
2839 
2840  if(l1cinput->verbose)cout << "Opening L1C grid ." << gridname << "for parallax-cloud corrections......." << endl;
2841  try {
2842  nc_l1c = new NcFile(ptstr, NcFile::read);
2843  } catch (NcException& e) {
2844  e.what();
2845  cerr << "l1cgen l1c_pflag= 8:: Failure write FULL L1C grid: " + gridname << endl;
2846  exit(1);
2847  }
2848 
2849  // global attributs
2850  NcGroupAtt i1 = nc_l1c->getAtt("instrument");
2851  i1.getValues(l1cfile->instrument);
2852  i1 = nc_l1c->getAtt("startDirection"); // NcGroupAtt is a global attr!!
2853  i1.getValues(l1cfile->start_dir);
2854  i1 = nc_l1c->getAtt("endDirection");
2855  i1.getValues(l1cfile->end_dir);
2856  i1 = nc_l1c->getAtt("time_coverage_start");
2857  i1.getValues(l1cfile->tswt_ini);
2858  i1 = nc_l1c->getAtt("time_coverage_end");
2859  i1.getValues(l1cfile->tswt_end);
2860  i1 = nc_l1c->getAtt("nadir_bin");
2861  i1.getValues(l1cfile->binstr);
2862 
2863  // dims
2864  NcDim ydim = nc_l1c->getDim("bins_along_track");
2865  NcDim xdim = nc_l1c->getDim("bins_across_track");
2866 
2867  //************************************************************
2868  // mem allocation for geolocation pointers---
2869  num_gridlines = ydim.getSize();
2870  nbinx = xdim.getSize();
2871 
2872  if(l1cinput->verbose)cout << "num_gridlines..." << num_gridlines << "nbinx....." << nbinx << endl;
2873 
2874  l1cfile->lat_gd = allocate2d_float(num_gridlines, nbinx);
2875  l1cfile->lon_gd = allocate2d_float(num_gridlines, nbinx);
2876  l1cfile->alt_gd = allocate2d_float(num_gridlines, nbinx);
2877  l1cfile->index_xy = allocate2d_short(num_gridlines, nbinx);
2878  l1cfile->lat_asort = allocate2d_float(num_gridlines, nbinx);
2879 
2880  NcGroup geo_grp = nc_l1c->getGroup("geolocation_data");
2881  NcVar v1 = geo_grp.getVar("latitude");
2882  v1.getVar(&l1cfile->lat_gd[0][0]);
2883  v1 = geo_grp.getVar("longitude");
2884  v1.getVar(&l1cfile->lon_gd[0][0]);
2885  v1 = geo_grp.getVar("altitude");
2886  v1.getVar(&l1cfile->alt_gd[0][0]);
2887 
2888  nc_l1c->close();
2889 
2890  NY = num_gridlines;
2891  NX = nbinx;
2892  l1cfile->num_gridlines = NY;
2893  l1cfile->nbinx = NX;
2894 
2895  //*********** SORTING LAT/LON ASCENDING AND TRACK INDEXES ****************************
2896  //***********************************************************************************
2897  vector<pair<float, int>> vp;
2898 
2899  // fill vector with lat_gd column--
2900  if(l1cinput->verbose)cout << "*********** SORTING LAT/LON GRIDPOINTS AND TRACKING INDEXES ********************" << endl;
2901 
2902  for (int col = 0; col < NX; col++) {
2903  for (int row = 0; row < NY; row++) {
2904  vp.push_back(make_pair(l1cfile->lat_gd[row][col] + 90., row));
2905  }
2906  // sort
2907  stable_sort(vp.begin(), vp.end());
2908 
2909  for (unsigned int i = 0; i < vp.size(); i++) {
2910  l1cfile->lat_asort[i][col] = vp[i].first;
2911  l1cfile->index_xy[i][col] = vp[i].second;
2912  }
2913  vp.clear();
2914  }
2915 
2916  return 0;
2917 }
2918 
2919 
2920 
2921 int32_t L1C::write_L1C_granule2(int swtd, l1c_filehandle* l1cfile, L1C_input* l1cinput, double* tmgv,
2922  float** lat_gd, float** lon_gd, float** alt_gd,double *orb_time_tot) {
2923  int32_t num_gridlines, nbinx, NY1 = -1, NY2 = -1;
2924  int32_t NX, NY;
2925  const char* filename_lt;
2926  float **lat_out = nullptr, **lon_out = nullptr, **alt_out = nullptr;
2927  double* time_nad_out = nullptr;
2928  int asc_mode = l1cfile->orb_dir;
2929  string senstr;
2930  double tg_ini, tg_end,tg_mid;
2931  int numgran;
2932  double tgridline=-1, tfile_ini_sec=-1, tfile_end_sec=-1;
2933  int16_t* granid = nullptr;
2934  int32_t gransize;
2935  short** gdindex = nullptr;
2936  double** gdtime = nullptr;
2937  int16_t y_ini, mo_ini, d_ini, h_ini, mi_ini, y_end, mo_end, d_end, h_end, mi_end, y_mid, mo_mid, d_mid, h_mid, mi_mid,y_zero, mo_zero, d_zero,h_zero,mi_zero;
2938  double sec_ini, sec_end,sec_mid,sec_zero;
2939  string tswt_ini, tswt_end, tswt_ini_file;
2940 
2941  int16_t syear, smon, sday, syear2, smon2, sday2;
2942  double secs, secs2;
2943  double tgran_ini, tgran_ini_sec, tgran_end, tgran_end_sec;
2944  string gfull;
2945  int16_t gtime;
2946  int Ngring = -1, dp = -1;
2947  int32_t iyear=0,iday=0,msec=0;
2948  int twodays=0,nadir_bin_index;
2949  Geodesic geod(Constants::WGS84_a(), Constants::WGS84_f());
2950  GeodesicLine line;
2951  size_t norb_rec=l1cfile->norb_rec;
2952  double tfile_ini_offset=-1;
2953 
2954  if (l1cinput->sensor == 34) {
2955  senstr = "SPEXONE";
2956  l1cfile->nbinx = 29;
2957  nadir_bin_index=14;
2958 
2959  } else if (l1cinput->sensor == 30 || l1cinput->sensor == 31) {
2960  senstr = "OCI"; // OCIS
2961  l1cfile->nbinx = 519;
2962  nadir_bin_index=259;
2963 
2964  } else if (l1cinput->sensor == 35) {
2965  senstr = "HARP2";
2966  l1cfile->nbinx = 457;
2967  nadir_bin_index=228;
2968  } else {
2969  if(l1cinput->verbose)cout << "sensor by default is OCI option 2....." << endl;
2970  senstr = "OCI";
2971  l1cfile->nbinx = 519;
2972  nadir_bin_index=259;
2973  }
2974 
2975 
2976  gransize = l1cfile->gransize;
2977 
2978  num_gridlines = l1cfile->num_gridlines;
2979  if(l1cinput->verbose)cout << "swath#....................................................................." << swtd
2980  << "asc_mode..." << asc_mode << endl;
2981 
2982  granid = (int16_t*)calloc(num_gridlines, sizeof(int16_t));
2983 
2984  numgran = l1cfile->numgran;//1 day of granules or 288x5 minutes/60 = 24 h
2985 
2986  nbinx = l1cfile->nbinx;
2987 
2988  gdtime = allocate2d_double(numgran, 2);
2989  gdindex = allocate2d_short(numgran, 2);
2990 
2991  tfile_ini_sec = l1cfile->tfile_ini_sec;
2992  tfile_ini_offset=tfile_ini_sec;
2993  tfile_end_sec = l1cfile->tfile_end_sec;
2994  double time_zero=tfile_ini_sec-24*3600;
2995  tg_ini = tfile_ini_sec; // always unix time
2996  tg_end = tg_ini + gransize * 60; // always in seconds
2997 
2998  // first granule----
2999  if (l1cinput->start_time[0] != '\0' && l1cinput->end_time[0] != '\0' && l1cinput->grantype == 0) {
3000  if(l1cinput->verbose)cout << "Processing L1C granules between ........................" << l1cinput->start_time
3001  << "and ........................." << l1cinput->end_time << endl;
3002  tgran_ini = isodate2unix(l1cinput->start_time);
3003  tgran_end = isodate2unix(l1cinput->end_time);
3004  unix2ymds(tgran_ini, &syear, &smon, &sday, &secs);
3005  unix2ymds(tgran_end, &syear2, &smon2, &sday2, &secs2);
3006  tgran_ini_sec = tgran_ini;
3007 
3008  tg_ini = tgran_ini_sec;
3009  tg_end = tg_ini + gransize * 60;
3010  tgran_end_sec = tgran_end;
3011 
3012  tfile_end_sec = tgran_end_sec + gransize * 60;
3013 
3014  unix2ymdhms(tgran_ini_sec, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
3015  if(l1cinput->verbose)cout << "tgran_ini_sec.."
3016  << "YEAR.." << y_ini << "MONTH.." << mo_ini << "DAY.." << d_ini << "HOUR.." << h_ini << "MIN.."
3017  << mi_ini << "SEC.." << sec_ini << endl;
3018  unix2ymdhms(tgran_end_sec, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
3019 
3020  if (tgran_end_sec > tfile_end_sec) {
3021  cout << "ERROR--tgran_end_sec>tfile_end_sec...wrong command line (gran_end_time).gran_end_time "
3022  "beyond swath time limits.."
3023  << endl;
3024  exit(1);
3025  }
3026 
3027  if(l1cinput->verbose)cout << "tgran_end_sec.."
3028  << "YEAR.." << y_ini << "MONTH.." << mo_ini << "DAY.." << d_ini << "HOUR.." << h_ini << "MIN.."
3029  << mi_ini << "SEC.." << sec_ini << endl;
3030  unix2ymdhms(tg_ini, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
3031  if(l1cinput->verbose)cout << "tg_ini.."
3032  << "YEAR.." << y_ini << "MONTH.." << mo_ini << "DAY.." << d_ini << "HOUR.." << h_ini << "MIN.."
3033  << mi_ini << "SEC.." << sec_ini << endl;
3034  unix2ymdhms(tg_end, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
3035  if(l1cinput->verbose)cout << "tg_end.."
3036  << "YEAR.." << y_ini << "MONTH.." << mo_ini << "DAY.." << d_ini << "HOUR.." << h_ini << "MIN.."
3037  << mi_ini << "SEC.." << sec_ini << endl;
3038  } else {
3039  if(l1cinput->verbose)cout << "ERROR Processing L1C granules between initial and final time --"
3040  << "start_time.." << l1cinput->start_time << "end_time.." << l1cinput->end_time << "grantype..."
3041  << l1cinput->grantype << endl;
3042  }
3043 
3044  for (int i = 0; i < num_gridlines; i++) {
3045  granid[i] = -1;
3046  }
3047 
3048  for (int i = 0; i < numgran; i++) {
3049  for (size_t j = 0; j < 2; j++) {
3050  gdtime[i][j] = -1;
3051  gdindex[i][j] = -1;
3052  }
3053  }
3054 
3055  int neg = 0, gmin = -1, c = 0;
3056  if(tmgv[0]<0) twodays=1;//the granule can start later even tough the grid contains gridlines from previous day
3057 
3058 
3059  for (int gran = 0; gran < numgran; gran++) {
3060  for (int i = 0; i < num_gridlines; i++) {
3061  tgridline =
3062  tfile_ini_sec + tmgv[i]; // seconds of the day are added to second since unix time reference
3063  if (tg_end > tgran_end_sec)
3064  tg_end = tgran_end_sec;
3065  if (tgridline >= tg_ini && tgridline < tg_end) {
3066 
3067  if (l1cinput->start_time[0] != '\0' || l1cinput->end_time[0] != '\0') {
3068  if (tg_ini >= tgran_ini_sec && tg_end <= tgran_end_sec + gransize * 60) {
3069  if (gmin < 0) { // first time of the granule or th_ini
3070  gdtime[c][0] = tg_ini;
3071  if (i == num_gridlines - 1)
3072  gdtime[c][1] = tgridline;
3073  else
3074  gdtime[c][1] = tg_end;
3075  gdindex[c][0] = i;
3076  gmin = 1;
3077  }
3078  gdindex[c][1] = i;
3079  if (i == num_gridlines - 1)
3080  gdtime[c][1] = tgridline;
3081  else
3082  gdtime[c][1] = tg_end;
3083  } // gran time
3084  } // command line
3085  else {
3086  if (gmin < 0) {
3087  gdtime[c][0] = tg_ini;
3088  if (i == num_gridlines - 1)
3089  gdtime[c][1] = tgridline;
3090  else
3091  gdtime[c][1] = tg_end;
3092  gdindex[c][0] = i;
3093  gmin = 1;
3094  }
3095  gdindex[c][1] = i;
3096  if (i == num_gridlines - 1)
3097  gdtime[c][1] = tgridline;
3098  else
3099  gdtime[c][1] = tg_end;
3100  }
3101  }
3102 
3103 
3104  if (tmgv[i] < 0 && gran == 0) { // previous day data
3105  neg++;
3106  }
3107  } // gridlines
3108  // new granule----
3109  tg_ini = tg_end;
3110  tg_end = tg_ini + gransize * 60;
3111  gmin = -1;
3112  c++;
3113 
3114  if (tg_ini > tgridline) {
3115  if(l1cinput->verbose)cout << "gridlines with negative time..." << neg << endl;
3116  break;
3117  }
3118  // granule selection constraints----------------
3119  if (l1cinput->start_time[0] != '\0' || l1cinput->end_time[0] != '\0') {
3120  if (tg_ini >= tg_end || tg_ini >= tgran_end_sec || tg_end > (tgran_end_sec + gransize * 60)) {
3121  if(l1cfile->verbose)cout<<"ERROR : GRAN # "<<c+1<<"tg_ini >= tg_end || tg_ini >= tgran_end_sec || tg_end > (tgran_end_sec + gransize * 60"<<endl;
3122  break;
3123  }
3124  }
3125 
3126  } // granules
3127 
3128  std::string timestr, missionstr, fname_out, fname_out_nopath, pathstr, monstr, daystr, yearstr, secstr,
3129  mistr, hstr, prodstr, gdstr, swtstr, extstr, ofilestr, dirstr1, dirstr2,datetimestr1, datetimestr2,datetimestr3,
3130  fdatetimestr1, date_created,datezerotime;
3131  std::string y_create, m_create, d_create, t_create;
3132 
3133  pathstr = "";
3134  missionstr = "PACE";
3135  extstr = ".nc";
3136  // write filenames to list----
3137  if(l1cinput->outlist[0]=='\0')
3138  strcpy(l1cinput->outlist, "l1c.tmp");
3139  string outxt(l1cinput->outlist);
3140 
3141  outxt = pathstr + outxt;
3142 
3143  std::ofstream outf;
3144 
3145  outf.open(outxt, std::ofstream::out);
3146 
3147  if (outf) {
3148  if(l1cinput->verbose)cout << "writing L1C granules to outfile..." << outxt << endl;
3149  } else {
3150  std::cerr << "output file.." << outxt << " could not be opened for writing!\n";
3151  return 1;
3152  }
3153 
3154  // Current date/time based on current system
3155  date_created = unix2isodate(now(), 'G');
3156 
3157  int16_t ngridlines, totlines = 0;
3158  // write files with granid>0 ---------------------
3159  for (int gran = 0; gran < numgran; gran++) {
3160  // twodays=0;
3161 
3162  NY1 = gdindex[gran][0];
3163  NY2 = gdindex[gran][1];
3164  /* for (int i = NY1; i < NY2 + 1; i++)
3165  {
3166  if(tmgv[i]<0) twodays=1;
3167  }
3168  */
3169  if(twodays)
3170  {
3171  tfile_ini_offset=time_zero;
3172  if(l1cinput->verbose==1) cout<<"twodays flag is ON!! : "<<twodays<<"tzero offset = "<< tfile_ini_offset<<endl;
3173  }
3174 
3175 
3176  if(l1cinput->verbose){
3177  tg_ini = gdtime[gran][0];
3178  tg_end = gdtime[gran][1];
3179  unix2ymdhms(tg_ini, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
3180  unix2ymdhms(tg_end, &y_end, &mo_end, &d_end, &h_end, &mi_end, &sec_end);
3181  cout<<"--------------------------"<<endl;
3182  cout<<"gran# "<<gran+1<<"day ini "<<d_ini<<"h ini "<<h_ini<<"mi ini "<<mi_ini<<"sec ini "<<sec_ini<<"tgini "<<tg_ini<<endl;
3183  cout<<"gran# "<<gran+1<<"day end "<<d_end<<"h end "<<h_end<<"mi end "<<mi_end<<"sec end "<<sec_end<<"tend "<<tg_end<<endl;
3184  tg_ini = tfile_ini_offset+orb_time_tot[0];
3185  tg_end = tfile_ini_sec+orb_time_tot[norb_rec-1];
3186  unix2ymdhms(tg_ini, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
3187  unix2ymdhms(tg_end, &y_end, &mo_end, &d_end, &h_end, &mi_end, &sec_end);
3188  cout<<"orbital time limits ----------------- "<<endl;
3189  cout<<"gran# "<<gran+1<<"day ini "<<d_ini<<"h ini "<<h_ini<<"mi ini "<<mi_ini<<"sec ini "<<sec_ini<<"tgini "<<tg_ini<<endl;
3190  cout<<"gran# "<<gran+1<<"day end "<<d_end<<"h end "<<h_end<<"mi end "<<mi_end<<"sec end "<<sec_end<<"tend "<<tg_end<<endl;
3191  cout<<"---------------------------"<<endl;
3192  }
3193 
3194  double torb_off=0;//5 more minutes after orb to avoid time truncation by granule
3195  //some orbital records from previous day but not in the crossing granule
3196  if(tfile_ini_offset+orb_time_tot[0]>tfile_ini_sec+orb_time_tot[norb_rec-1]) orb_time_tot[0]-=24*3600;
3197  if(gdtime[gran][1]>tfile_ini_sec+orb_time_tot[norb_rec-1]) torb_off=gdtime[gran][1]-tfile_ini_sec-orb_time_tot[norb_rec-1];
3198 
3199  if(torb_off>2*60) torb_off=0;
3200 
3201  if (gdtime[gran][0]>=(tfile_ini_offset+orb_time_tot[0]) && gdtime[gran][1]<=(tfile_ini_sec+orb_time_tot[norb_rec-1]+torb_off))
3202  {
3203  if(l1cinput->verbose)cout << "gran #..." << gran + 1 <<endl;
3204 
3205  tg_ini = gdtime[gran][0];
3206  tg_end = gdtime[gran][1];
3207  tg_mid=(tg_ini+tg_end)/2;
3208 
3209  unix2ymdhms(tg_ini, &y_ini, &mo_ini, &d_ini, &h_ini, &mi_ini, &sec_ini);
3210  unix2ymdhms(tg_end, &y_end, &mo_end, &d_end, &h_end, &mi_end, &sec_end);
3211  unix2ymdhms(tg_mid, &y_mid, &mo_mid, &d_mid, &h_mid, &mi_mid, &sec_mid);
3212 
3213  if (mi_end * 60 == 0)
3214  {
3215  gtime = ((60 * 60 + round(sec_end)) - mi_ini * 60 - round(sec_ini)) / 60;
3216  }
3217  else
3218  {
3219  gtime = ((mi_end * 60 + round(sec_end)) - mi_ini * 60 - round(sec_ini)) / 60;
3220  }
3221 
3222  if (l1cinput->verbose && gtime > gransize)
3223  {
3224  cout << "WARNING: rounding errors --gtime = " << gtime << " is greater than granule size = " << gransize << endl; }
3225  if ((gtime % gransize) == 0)
3226  gfull = "1";
3227  else
3228  gfull = "0";
3229  if(l1cinput->verbose)cout << "gtime.." << gtime << "gfull.." << gfull << endl;
3230 
3231  // # gridlines per granule
3232  // gridline index 0-num_gridlines
3233  ngridlines = gdindex[gran][1] - gdindex[gran][0] + 1;
3234  totlines += ngridlines;
3235  if(l1cinput->verbose)cout << "ngridlines..." << ngridlines << "tot gridlines..." << totlines << endl;
3236 
3237  if(twodays)//elapsed from previous day
3238  {
3239  h_zero=0;
3240  mi_zero=0;
3241  sec_zero=0;
3242  unix2ymdhms(time_zero, &y_zero, &mo_zero, &d_zero, &h_zero, &mi_zero, &sec_zero);
3243  }
3244  else //elapsed from the same day
3245  {
3246  d_zero=d_ini;
3247  mo_zero=mo_ini;
3248  y_zero=y_ini;
3249  }
3250 
3251  daystr = std::to_string(d_zero);
3252  if(daystr.size() == 1)
3253  daystr = "0" + daystr;
3254  monstr = std::to_string(mo_zero);
3255  if(monstr.size() == 1)
3256  monstr = "0" + monstr;
3257  yearstr = std::to_string(y_zero);
3258 
3259  datezerotime="seconds since " + yearstr + "-" + monstr + "-" + daystr;
3260  string datemid= unix2isodate(tg_mid, 'G');
3261  isodate2ydmsec((char *) datemid.substr(0,19).c_str(), &iyear, &iday, &msec);
3262  double dist_es=esdist_(&iyear,&iday,&msec);
3263 
3264  if(l1cinput->verbose)cout << "sun_earth_distance -- mid gridline" <<dist_es<<endl;
3265 
3266  // output file with list of granules
3267  string dateini= unix2isodate(tg_ini, 'G');
3268  string datend= unix2isodate(tg_end, 'G');
3269 
3270 
3271  yearstr=dateini.substr(0,4);
3272  monstr=dateini.substr(5,2);
3273  daystr=dateini.substr(8,2);
3274  hstr=dateini.substr(11,2);
3275  mistr=dateini.substr(14,2);
3276  secstr=dateini.substr(17,2);
3277 
3278  fdatetimestr1 = yearstr + monstr + daystr + "T" + hstr + mistr + secstr;
3279  fname_out = pathstr + "PACE." + fdatetimestr1 + ".L1C" + extstr;
3280  fname_out_nopath = "PACE." + fdatetimestr1 + ".L1C" + extstr;
3281 
3282  if(l1cinput->verbose)cout << "granule filename.." << fname_out << endl;
3283 
3284  outf << fname_out_nopath << "," << dateini.substr(0,19) << "," << datend.substr(0,19) << "," << gfull << "\n";
3285 
3286 
3287  l1cfile->gridname = fname_out.c_str();
3288  filename_lt = fname_out.c_str();
3289  char* gridchar = strdup(filename_lt);
3290  string l1c_str = filename_lt;
3291 
3292  NcFile* nc_output;
3293  try {
3294  nc_output = new NcFile(filename_lt, NcFile::replace);
3295  } catch (NcException& e) {
3296  e.what();
3297  cerr << "l1cgen l1c_pflag= 5 : producing L1C grid: " + l1c_str << endl;
3298  exit(1);
3299  }
3300 
3301  bin_str binl1c;
3302  binl1c.sensor = l1cinput->sensor;
3303  meta_l1c_grid(gridchar, &binl1c, ngridlines, nc_output);
3304 
3305  if(l1cinput->pversion[0])
3306  nc_output->putAtt("processing_version", l1cinput->pversion);
3307  if(l1cinput->doi[0]) {
3308  nc_output->putAtt("identifier_product_doi", l1cinput->doi);
3309  nc_output->putAtt("identifier_product_doi_authority", "http://dx.doi.org");
3310  }
3311  nc_output->putAtt("history", l1cinput->history);
3312  nc_output->putAtt("product_name", fname_out_nopath);
3313  nc_output->putAtt("time_coverage_start", dateini.substr(0,19) + "Z");
3314  nc_output->putAtt("time_coverage_end", datend.substr(0,19) + "Z");
3315  nc_output->putAtt("sun_earth_distance", ncFloat, dist_es);
3316 
3317  NY = ngridlines;
3318  NX = nbinx;
3319  lat_out = allocate2d_float(NY, NX);
3320  lon_out = allocate2d_float(NY, NX);
3321  alt_out = allocate2d_float(NY, NX);
3322  time_nad_out = (double*)calloc(NY, sizeof(double)); // time of the day in seconds
3323 
3324  int cc = 0,ep_shift;
3325 
3326  if(tmgv[NY1]<0) ep_shift = 3600*24;
3327  else ep_shift=0;
3328  if(l1cinput->verbose) cout<<"twodays flag : ep_shift :"<<ep_shift<<endl;
3329  for (int i = NY1; i < NY2 + 1; i++) {
3330  for (int j = 0; j < NX; j++) {
3331  lat_out[cc][j] = lat_gd[i][j];
3332  lon_out[cc][j] = lon_gd[i][j];
3333  alt_out[cc][j] = alt_gd[i][j];
3334  time_nad_out[cc] = tmgv[i]+ep_shift;
3335  }
3336  cc++;
3337  }
3338  if(lat_gd[NY1+1][nadir_bin_index]>lat_gd[NY1][nadir_bin_index]) {dirstr1="Ascending";l1cfile->orb_dir=1;} else {dirstr1= "Descending";l1cfile->orb_dir=0;}
3339  if(lat_gd[NY2][nadir_bin_index]>lat_gd[NY2-1][nadir_bin_index]) dirstr2="Ascending";else dirstr2= "Descending";
3340  nc_output->putAtt("startDirection", dirstr1);
3341  nc_output->putAtt("endDirection", dirstr2);
3342 
3343  // vars
3344  NcGroup ba_grp = nc_output->getGroup("bin_attributes");
3345  NcVar v1 = ba_grp.getVar("nadir_view_time");
3346  v1.putVar(&time_nad_out[0]);
3347  v1.putAtt("units",datezerotime);
3348  NcGroup geo_grp = nc_output->getGroup("geolocation_data");
3349  v1 = geo_grp.getVar("latitude");
3350  v1.putVar(&lat_out[0][0]);
3351  v1 = geo_grp.getVar("longitude");
3352  v1.putVar(&lon_out[0][0]);
3353  v1 = geo_grp.getVar("height");
3354  v1.putVar(&alt_out[0][0]);
3355 
3356  // GRING-------
3357  // determine the number of GCpoint indexes
3358  // default 6 for the swath sides + number of coordinates every 20 degrees latitude
3359 
3360  Ngring = 4;
3361  if(Ngring>0){
3362  float* latarr = (float*)calloc(Ngring, sizeof(float));
3363  float* lonarr = (float*)calloc(Ngring, sizeof(float));
3364  // ascending pass
3365  if (l1cfile->orb_dir == 1) {
3366  // corners--counterclockwise and ascending pass
3367  latarr[0] = lat_gd[NY2][l1cfile->nbinx - 1];
3368  latarr[1] = lat_gd[NY2][0];
3369  latarr[2] = lat_gd[NY1][0];
3370  latarr[3] = lat_gd[NY1][l1cfile->nbinx - 1];
3371 
3372  lonarr[0] = lon_gd[NY2][l1cfile->nbinx - 1];
3373  lonarr[1] = lon_gd[NY2][0];
3374  lonarr[2] = lon_gd[NY1][0];
3375  lonarr[3] = lon_gd[NY1][l1cfile->nbinx - 1];
3376 
3377  } else // descending orb
3378  {
3379  latarr[0] = lat_gd[NY2][l1cfile->nbinx - 1];
3380  latarr[1] = lat_gd[NY2][0];
3381  latarr[2] = lat_gd[NY1][0];
3382  latarr[3] = lat_gd[NY1][l1cfile->nbinx - 1];
3383 
3384  lonarr[0] = lon_gd[NY2][l1cfile->nbinx - 1];
3385  lonarr[1] = lon_gd[NY2][0];
3386  lonarr[2] = lon_gd[NY1][0];
3387  lonarr[3] = lon_gd[NY1][l1cfile->nbinx - 1];
3388 
3389  } // end descending
3390 
3391  string onelat,onelon,onecoor,firstcoor;
3392  //GRing params
3393  string gring_polygon,gring_crs,gring_latmin,gring_latmax,gring_lonmin,gring_lonmax;
3394  dp=4;
3395 
3396  for (int s = 0; s < dp; s++) {
3397  onelat=to_string(latarr[s]);
3398  onelon=to_string(lonarr[s]);
3399  if(s==0)
3400  {
3401  onecoor="POLYGON(("+onelat+" "+onelon+",";
3402  firstcoor=onelat+" "+onelon;
3403  }
3404  else onecoor=onelat+" "+onelon+",";
3405 
3406  gring_polygon+=onecoor;
3407  }
3408 
3409 
3410  //last coor
3411  gring_polygon+=firstcoor+"))";
3412 
3413  float latpmin=999,latpmax=-999,lonpmin=999,lonpmax=-999;
3414  for(int row=NY1;row<NY2+1;row++)
3415  {
3416  for(int col=0;col<nbinx;col++)
3417  {
3418 
3419  if(lat_gd[row][col]<latpmin && lat_gd[row][col]!= BAD_FLT) latpmin=lat_gd[row][col];
3420  if(lat_gd[row][col]>latpmax && lat_gd[row][col]!= BAD_FLT) latpmax=lat_gd[row][col];
3421  if(lon_gd[row][col]<lonpmin && lon_gd[row][col]!= BAD_FLT) lonpmin=lon_gd[row][col];
3422  if(lon_gd[row][col]>lonpmax && lon_gd[row][col]!= BAD_FLT) lonpmax=lon_gd[row][col];
3423  }}
3424 
3425  gring_crs="EPSG:4326";
3426  nc_output->putAtt("geospatial_bounds", gring_polygon);
3427  nc_output->putAtt("geospatial_bounds_crs", gring_crs);
3428 
3429  //algo to find min/max coordinates inside the gring polygon
3430  nc_output->putAtt("geospatial_lat_min", ncFloat, latpmin);
3431  nc_output->putAtt("geospatial_lat_max", ncFloat, latpmax);
3432  nc_output->putAtt("geospatial_lon_min", ncFloat, lonpmin);
3433  nc_output->putAtt("geospatial_lon_max", ncFloat, lonpmax);
3434 
3435  free (latarr);
3436  free (lonarr);
3437  } else {
3438  cout << "ERROR EXTRACTING GRING coordinates!!-----" << endl;
3439  exit(1);
3440  }
3441 
3442  nc_output->close();
3443 
3444  if (lat_out != nullptr)
3445  free (lat_out);
3446  if (lon_out != nullptr)
3447  free (lon_out);
3448  if (alt_out != nullptr)
3449  free (alt_out);
3450  if (time_nad_out != nullptr)
3451  free (time_nad_out);
3452 
3453  add_proc_group_l1c(l1cinput,l1cfile,fname_out.c_str());
3454  } // gdtime >0
3455  } // end granule loop
3456 
3457  if (gdtime != nullptr)
3458  free (gdtime);
3459  gdtime = nullptr;
3460 
3461  if (gdindex != nullptr)
3462  free (gdindex);
3463  gdindex = nullptr;
3464 
3465  if (granid != nullptr)
3466  free (granid);
3467  granid = nullptr;
3468 
3469  outf.close();
3470 
3471  return 0;
3472 }
3473 
3474 
3475 
3476 
3477 // this version includes info coming from Don open file routines--
3478 
3479 int32_t L1C::load_l1c_filehandle4(l1c_filehandle* l1cfile, L1C_input* l1cinput) {
3480  int nfiles;
3481  string filename;
3482  nfiles = l1cinput->files.size();
3483 
3484  for (int j = 0; j < nfiles; j++) {
3485  filename = l1cinput->files[j];
3486  l1cfile->ifiles.push_back(filename);
3487  }
3488 
3489  l1cfile->l1c_pflag = l1cinput->l1c_pflag;
3490  l1cfile->swath_num = l1cinput->swath_num;
3491 
3492  for (int j = 0; j < 10; j++) { // up to 10 granules
3493  l1cfile->selgran[j] = l1cinput->selgran[j];
3494  }
3495  // projection params
3496  l1cfile->gres = l1cinput->grid_resolution;
3497  l1cfile->verbose = l1cinput->verbose;
3498  l1cfile->proj_type = l1cinput->projection;
3499  l1cfile->cloud_corrected = l1cinput->cloud_correct;
3500  // multi attributes (view, pol, bands)
3501  l1cfile->overlap_vflag = l1cinput->overlap_vflag;
3502  l1cfile->overlap_pflag = l1cinput->overlap_pflag;
3503  l1cfile->overlap_bflag = l1cinput->overlap_bflag;
3504  // uncertainty params l1c merged products
3505  l1cfile->unc_meth = l1cinput->unc_meth;
3506  l1cfile->unc_thres_v = l1cinput->unc_thres_v;
3507  l1cfile->unc_thres_p = l1cinput->unc_thres_p;
3508  l1cfile->unc_thres_b = l1cinput->unc_thres_b;
3509 
3510  // calibration attributes
3511  // l1cfile->Fobar=file->Fobar;//nc
3512 
3513  cout << "ok transfering filehandle to l1c_filehandle info.." << endl;
3514  return 0;
3515 }
3516 
3517 }
l1c::l1c_filehandle::lat_asort
float ** lat_asort
Definition: l1c_filehandle.h:129
l1c::l1c_filehandle::start_dir
std::string start_dir
Definition: l1c_filehandle.h:96
value
int32 value
Definition: Granule.c:1235
l1c::L1C_input::overlap_bflag
bool overlap_bflag
Definition: l1c_input.h:80
l1c::l1c_filehandle::orbit_node_lon
float orbit_node_lon
Definition: l1c_filehandle.h:90
allocate3d.h
Utility functions for allocating and freeing three-dimensional arrays of various types.
NX
#define NX
Definition: main_biosmap.c:51
l1c::l1c_filehandle::nbinx
int32_t nbinx
Definition: l1c_filehandle.h:137
l1c::l1c_filehandle::unc_thres_p
float unc_thres_p
Definition: l1c_filehandle.h:165
j
int j
Definition: decode_rs.h:73
l1c::l1c_filehandle::progname
std::string progname
Definition: l1c_filehandle.h:34
day
int32_t day
Definition: atrem_corl1v2.h:308
status
int status
Definition: l1_czcs_hdf.c:32
l1c::l1c_filehandle::binstr
std::string binstr
Definition: l1c_filehandle.h:98
check_err
void check_err(const int stat, const int line, const char *file)
Definition: nc4utils.c:35
l1c::open_l1atol1c3
int32_t open_l1atol1c3(L1C_input *l1cinput, l1c_filehandle *l1cfile)
l1c::l1c_filehandle::proj_type
size_t proj_type
Definition: l1c_filehandle.h:147
min
These are used to scale the SD before writing it to the HDF4 file The default is and which means the product is not scaled at all Since the product is usually stored as a float inside of this is a way to write the float out as a integer l2prod min
Definition: HOWTO_Add_a_product.txt:76
l1c::L1C_input::verbose
bool verbose
Definition: l1c_input.h:42
l1c::l1c_filehandle::index_xy
short ** index_xy
Definition: l1c_filehandle.h:130
l1c::L1C_input::l2prod
char l2prod[2048]
Definition: l1c_input.h:41
l1c_input.h
l1c::l1c_filehandle::sensor
int32_t sensor
Definition: l1c_filehandle.h:138
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
bin_str::sensor
int sensor
Definition: l1c_latlongrid.h:28
l1c::L1C_input::overlap_vflag
bool overlap_vflag
Definition: l1c_input.h:78
l1c_str.h
l1c::L1C_input::unc_thres_p
float unc_thres_p
Definition: l1c_input.h:84
unix2ymds
void unix2ymds(double usec, int16_t *year, int16_t *mon, int16_t *day, double *secs)
Definition: unix2ymds.c:8
l1c::swtime_swt2_segment
int32_t swtime_swt2_segment(int swt, L1C_input *l1cinput, l1c_filehandle *l1cfile, int32_t norbs, double *tswt, double tcross, double mgv, double *tmgv)
l1c::l1c_str::latpix
float * latpix
Definition: l1c_str.h:99
l1c::l1c_filehandle::lat_gd
float ** lat_gd
Definition: l1c_filehandle.h:126
l1c::l1c_filehandle::orb_dir
size_t orb_dir
Definition: l1c_filehandle.h:89
l1c::l1c_filehandle::l1c_pflag
size_t l1c_pflag
Definition: l1c_filehandle.h:45
hico.proc_hico.valid_min
valid_min
Definition: proc_hico.py:246
msec
int32 * msec
Definition: l1_czcs_hdf.c:31
l1c::L1C_input::west
float west
Definition: l1c_input.h:49
esdist_
double esdist_(int32_t *year, int32_t *day, int32_t *msec)
l1c::L1C_input::selday
int32_t selday
Definition: l1c_input.h:65
syear
int syear
Definition: l1_czcs_hdf.c:15
l1c::swtime_swt2
int32_t swtime_swt2(int swt, L1C_input *l1cinput, l1c_filehandle *l1cfile, int32_t norbs, double *tswt, double tcross, double mgv, double *tmgv)
file_format
Definition: filetype.h:69
hico.proc_hico.valid_max
valid_max
Definition: proc_hico.py:247
l1c::l1c_filehandle::n_views
size_t n_views
Definition: l1c_filehandle.h:69
l1c::write_L1C_granule2
int32_t write_L1C_granule2(int swtd, l1c_filehandle *l1cfile, L1C_input *l1cinput, double *tmgv, float **lat_gd, float **lon_gd, float **alt_gd, double *orb_time_tot)
string
@ string
Definition: GEO_read_param_file.c:64
l1c::l1c_filehandle::gd_per_gran
int32_t gd_per_gran
Definition: l1c_filehandle.h:112
l1c::l1c_filehandle::eqt
double eqt
Definition: l1c_filehandle.h:91
l1c::l1c_str::lonpix
float * lonpix
Definition: l1c_str.h:100
l1c::L1C_input::east
float east
Definition: l1c_input.h:50
hico.proc_hico.senz
senz
Definition: proc_hico.py:239
M_PI
#define M_PI
Definition: dtranbrdf.cpp:19
swapc_bytes
int swapc_bytes(char *in, int nbyte, int ntime)
Definition: swapc_bytes.c:4
sday
int sday
Definition: l1_czcs_hdf.c:15
l1c::ect_swt
int32_t ect_swt(l1c_filehandle *l1cfile, int ix1, int ix2, double *tswt_tot, double *latswt_tot, double *lonswt_tot, double *ovel_tot, double *gvel_tot, double *tswt, double *latswt, double *lonswt, float *tcross, float *loncross, double *ovel, double *gvel)
l1c::L1C_input::cloud_height
int32_t cloud_height
Definition: l1c_input.h:72
l1c::L1C_input::ix_l1cprod
int32_t ix_l1cprod[3]
Definition: l1c_input.h:61
l1c::L1C_input::l1c_grid
char l1c_grid[FILENAME_MAX]
Definition: l1c_input.h:30
l1c::L1C_input::demcloud_flag
int demcloud_flag
Definition: l1c_input.h:76
l1c::L1C_input::terrain_correct
bool terrain_correct
Definition: l1c_input.h:73
l1c::l1c_filehandle::overlap_bflag
size_t overlap_bflag
Definition: l1c_filehandle.h:161
sunangs_
void sunangs_(int *, int *, float *, float *, float *, float *, float *)
Definition: sunangs.c:25
f
double precision function f(R1)
Definition: tmd.lp.f:1454
l1c::l1c_filehandle::npix
size_t npix
Definition: l1c_filehandle.h:77
l1c::L1C_input::north
float north
Definition: l1c_input.h:48
l1c::l1c_filehandle::numgran
int numgran
Definition: l1c_filehandle.h:105
l1c::l1c_filehandle
Definition: l1c_filehandle.h:25
l1c::l1c_filehandle::tswt_ini
std::string tswt_ini
Definition: l1c_filehandle.h:92
Re
float Re
Definition: l1c.cpp:57
l1c::cross_product_double
void cross_product_double(double vector_a[], double vector_b[], double temp[])
Definition: l1c.cpp:229
l1c::L1C_input::swath_num
int32_t swath_num
Definition: l1c_input.h:56
l1c::L1C_input::south
float south
Definition: l1c_input.h:47
allocate4d.h
Utility functions for allocating and freeing four-dimensional arrays of various types.
l1c::L1C_input::pversion
char pversion[256]
Definition: l1c_input.h:36
strdup
char * strdup(const char *)
getFormat
file_format getFormat(char *filename)
Definition: filetype.c:192
l1c::L1C_input::sensor
int32_t sensor
Definition: l1c_input.h:64
hico.proc_hico.sena
sena
Definition: proc_hico.py:241
l1c::L1C_input::l1c_anc
char l1c_anc[FILENAME_MAX]
Definition: l1c_input.h:31
l1c::l1c_filehandle::gridname
const char * gridname
Definition: l1c_filehandle.h:43
l1c::load_l1c_filehandle4
int32_t load_l1c_filehandle4(l1c_filehandle *l1cfile, L1C_input *l1cinput)
l1c::l1c_filehandle::cloud_corrected
size_t cloud_corrected
Definition: l1c_filehandle.h:117
l1c::L1C_input::cloud_type
int cloud_type
Definition: l1c_input.h:75
l1c::sunz_swt
int sunz_swt(int ix1_swt, int ix2_swt, int16_t *hour_swt, int16_t *day_swt, int16_t *year_swt, double *olat, double *olon)
Definition: l1c.cpp:890
l1c::add_proc_group_l1c
int32_t add_proc_group_l1c(L1C_input *l1cinput, l1c_filehandle *l1cfile, const char *filename)
l1aconvert_octs.NC_WRITE
string NC_WRITE
Definition: l1aconvert_octs.py:29
dp
float dp[MODELMAX]
Definition: atrem_corl1v2.h:248
l1c::l1c_filehandle::overlap_pflag
size_t overlap_pflag
Definition: l1c_filehandle.h:160
filename
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
allocate2d_uchar
unsigned char ** allocate2d_uchar(size_t h, size_t w)
Allocate a two-dimensional array of type unsigned char of a given size.
Definition: allocate2d.c:35
l1c::l1c_filehandle::tswt_end
std::string tswt_end
Definition: l1c_filehandle.h:93
l1c::l1c_filehandle::tfile_ini_sec
double tfile_ini_sec
Definition: l1c_filehandle.h:106
bin_str
Definition: l1c_latlongrid.h:15
fix_mac_rpath.line
line
Definition: fix_mac_rpath.py:56
hawkeye_methods.h
l1c
Definition: l1c.cpp:71
isodate2ydmsec
void isodate2ydmsec(char *date, int32_t *year, int32_t *day, int32_t *msec)
Definition: date2ydmsec.c:20
l1c::L1C_input::gransize
int gransize
Definition: l1c_input.h:53
l1c::l1c_filehandle::swath_num
size_t swath_num
Definition: l1c_filehandle.h:40
l1c::L1C_input::start_timeflag
int start_timeflag
Definition: l1c_input.h:57
l1c::L1C_input::doi
char doi[256]
Definition: l1c_input.h:37
bin_str::date_mid_grid
std::string date_mid_grid
Definition: l1c_latlongrid.h:93
l1c::L1C_input::bintype
int bintype
Definition: l1c_input.h:55
allocate2d.h
Utility functions for allocating and freeing two-dimensional arrays of various types.
l1c::L1C_input::files_l1c
std::vector< std::string > files_l1c
Definition: l1c_input.h:45
l1c::l1c_filehandle::norb_rec
size_t norb_rec
Definition: l1c_filehandle.h:150
l1c::L1C_input::files
std::vector< std::string > files
Definition: l1c_input.h:44
l1c::create_time_swt
int32_t create_time_swt(int num_gridlines, double tfile_ini_sec, double *tmgvf, double tswt_ini_sec, double tswt_end_sec, std::string *tswt_ini, std::string *tswt_ini_file, std::string *tswt_mid, std::string *tswt_end)
l1c::L1C_input::grid_resolution
float grid_resolution
Definition: l1c_input.h:63
l1c::l1c_filehandle::end_dir
std::string end_dir
Definition: l1c_filehandle.h:97
l1c::L1C_input::sort_method
int32_t sort_method
Definition: l1c_input.h:71
l1c::l1c_filehandle::verbose
bool verbose
Definition: l1c_filehandle.h:39
hico.proc_hico.units
units
Definition: proc_hico.py:245
dotprod
long dotprod(void *dp, signed short a[])
l1c::L1C_input::selmon
int32_t selmon
Definition: l1c_input.h:66
l1c::L1C_input::selyear
int32_t selyear
Definition: l1c_input.h:67
NY
#define NY
Definition: main_biosmap.c:52
l1c::l1c_filehandle::version
std::string version
Definition: l1c_filehandle.h:33
meta_l1c_grid
int meta_l1c_grid(char *gridname, bin_str *binl1c, int16_t num_gridlines, NcFile *nc_output)
Definition: l1c_latlongrid.cpp:913
l1c::L1C_input::grantype
int grantype
Definition: l1c_input.h:54
color_dtdb.bg
bg
Definition: color_dtdb.py:451
l1c_filehandle.h
BAD_FLT
#define BAD_FLT
Definition: jplaeriallib.h:19
re
#define re
Definition: l1_czcs.c:695
Point_t
bg::model::point< double, 2, bg::cs::geographic< bg::degree > > Point_t
Definition: l1c.cpp:67
geometry
subroutine geometry
Definition: atrem_app_refl_f90_cubeio.f:1329
l2_str.h
l1c::l1c_filehandle::tswt_ini_file
std::string tswt_ini_file
Definition: l1c_filehandle.h:95
l1c::l1c_filehandle::unc_thres_b
float unc_thres_b
Definition: l1c_filehandle.h:166
l1c::l1c_filehandle::unc_meth
size_t unc_meth
Definition: l1c_filehandle.h:163
l1c::l1c_filehandle::selgran
int16_t selgran[10]
Definition: l1c_filehandle.h:42
l1c::L1C_input::unc_thres_v
float unc_thres_v
Definition: l1c_input.h:83
l1c::l1c_str
Definition: l1c_str.h:19
tai58_to_unix
double tai58_to_unix(double tai58)
Definition: yds2tai.c:29
l1c::L1C_input::history
char history[200]
Definition: l1c_input.h:35
l1c::L1C_input::unc_thres_b
float unc_thres_b
Definition: l1c_input.h:85
l1c::l1c_filehandle::gransize
int gransize
Definition: l1c_filehandle.h:111
l1c::L1C_input::overlap_pflag
bool overlap_pflag
Definition: l1c_input.h:79
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
l1c::L1C_input
Definition: l1c_input.h:24
allocate2d_double
double ** allocate2d_double(size_t h, size_t w)
Allocate a two-dimensional array of type double of a given size.
Definition: allocate2d.c:145
l1c::l1c_filehandle::tswt_mid
std::string tswt_mid
Definition: l1c_filehandle.h:94
l1c::L1C_input::cloud_correct
int cloud_correct
Definition: l1c_input.h:74
l1c::l1c_filehandle::unc_thres_v
float unc_thres_v
Definition: l1c_filehandle.h:164
s
data_t s[NROOTS]
Definition: decode_rs.h:75
Polygon_t
bg::model::polygon< Point_t > Polygon_t
Definition: l1c.cpp:68
l1c::L1C_input::l1c_pflag
int32_t l1c_pflag
Definition: l1c_input.h:62
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)
l1c::L1C_input::start_time
char start_time[50]
Definition: l1c_input.h:32
l1c::L1C_input::unc_meth
int32_t unc_meth
Definition: l1c_input.h:82
c
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude resolving resolving GSFcd00179 Corrected handling of fill values for[Sensor|Solar][Zenith|Azimuth] resolving MODxl01751 Changed to validate LUT version against a value retrieved from the resolving MODxl02056 Changed to calculate Solar Diffuser angles without adjustment for estimated post launch changes in the MODIS orientation relative to incidentally resolving defects MODxl01766 Also resolves MODxl01947 Changed to ignore fill values in SCI_ABNORM and SCI_STATE rather than treating them as resolving MODxl01780 Changed to use spacecraft ancillary data to recognise when the mirror encoder data is being set by side A or side B and to change calculations accordingly This removes the need for seperate LUTs for Side A and Side B data it makes the new LUTs incompatible with older versions of the and vice versa Also resolves MODxl01685 A more robust GRing algorithm is being which will create a non default GRing anytime there s even a single geolocated pixel in a granule Removed obsolete messages from seed as required for compatibility with version of the SDP toolkit Corrected test output file names to end in per delivery and then split off a new MYD_PR03 pcf file for Aqua Added AssociatedPlatformInstrumentSensor to the inventory metadata in MOD01 mcf and MOD03 mcf Created new versions named MYD01 mcf and MYD03 where AssociatedPlatformShortName is rather than Terra The program itself has been changed to read the Satellite Instrument validate it against the input L1A and LUT and to use it determine the correct files to retrieve the ephemeris and attitude data from Changed to produce a LocalGranuleID starting with MYD03 if run on Aqua data Added the Scan Type file attribute to the Geolocation copied from the L1A and attitude_angels to radians rather than degrees The accumulation of Cumulated gflags was moved from GEO_validate_earth_location c to GEO_locate_one_scan c
Definition: HISTORY.txt:464
l1c::l1c_filehandle::lon_gd
float ** lon_gd
Definition: l1c_filehandle.h:127
l1c::L1C_input::outlist
char outlist[FILENAME_MAX]
Definition: l1c_input.h:29
l1c::openL1Cgrid3
int32_t openL1Cgrid3(l1c_str *l1cstr, l1c_filehandle *l1cfile, L1C_input *l1cinput)
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
l1c::cross_product_norm_double
double cross_product_norm_double(double vector_a[], double vector_b[])
Definition: l1c.cpp:235
l1c::L1C_input::projection
int32_t projection
Definition: l1c_input.h:70
l1c::L1C_input::selgran
int selgran[10]
Definition: l1c_input.h:60
abs
#define abs(a)
Definition: misc.h:90
l1c::l1c_filehandle::gres
float gres
Definition: l1c_filehandle.h:149
i
int i
Definition: decode_rs.h:71
l1c::create_SOCEA2
int32_t create_SOCEA2(int swtd, L1C_input *l1cinput, l1c_filehandle *l1cfile, float **lat_gd, float **lon_gd, float **altitude, double *tswt)
getanc_aquarius.gran
gran
Definition: getanc_aquarius.py:234
l1c::l1c_filehandle::num_gridlines
int32_t num_gridlines
Definition: l1c_filehandle.h:133
orb_to_latlon
int orb_to_latlon(size_t ix_swt_ini, size_t ix_swt_end, size_t num_gridlines, int nbinx, double *orb_time_tot, orb_array2 *p, orb_array2 *v, double mgv1, double *tmgv1, double *tmgvf, float **lat_gd, float **lon_gd, float **alt, int FirsTerrain)
Definition: hawkeye_methods.cpp:41
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
l1c::L1C_input::end_time
char end_time[50]
Definition: l1c_input.h:33
str
Definition: aerosol.c:136
verbose
int verbose
Definition: fmt_check.c:6
l1c::search_l1cgen
int search_l1cgen(L1C_input *l1cinput, l1c_str *l1cstr, l1c_filehandle *l1cfile, short **gdindex)
k
int k
Definition: decode_rs.h:73
l1c::l1c_filehandle::ifiles
std::vector< std::string > ifiles
Definition: l1c_filehandle.h:48
l1c::l1c_filehandle::overlap_vflag
size_t overlap_vflag
Definition: l1c_filehandle.h:159
p
float p[MODELMAX]
Definition: atrem_corl1.h:131
l1c::l1c_filehandle::tfile_end_sec
double tfile_end_sec
Definition: l1c_filehandle.h:107
str
#define str(s)
Definition: leveltwo2csv-main.c:19
Box_t
bg::model::box< Point_t > Box_t
Definition: l1c.cpp:69
l1c::l1c_filehandle::alt_gd
float ** alt_gd
Definition: l1c_filehandle.h:128
isodate2unix
double isodate2unix(const char *isodate)
Definition: unix2isodate.c:61
l1c::l1c_filehandle::instrument
std::string instrument
Definition: l1c_filehandle.h:99
orb_array2
double orb_array2[3]
Definition: hawkeye_methods.h:9


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