NASA Logo
Ocean Color Science Software

ocssw V2022
main_l2binmatch.cpp
Go to the documentation of this file.
1 /* ========================================================================
2  * l2binmatch creates a file with L2 sensor validation data
3  *
4  * L2 files for a given time period from the sensor being validated
5  * L3 daily, 8-day, etc. bin file from the sensor used as validation truth
6  *
7  * ======================================================================== */
8 //todo implement product selection via the l2prod parameter
9 //todo implement bandshift option for handling different L2 and L3 sensors
10 #include <stdio.h>
11 #include <unistd.h>
12 #include <stdlib.h>
13 #include "l2binmatch_input.h"
14 #include <readL2scan.h>
15 #include <L3File.h>
16 #include "calfile_utils.h"
17 #include <setupflags.h>
18 
19 #include <vector>
20 
21 using namespace std;
22 using namespace l3;
23 
24 int16_t fileID = 0;
25 
26 float ***aots, ***nlws; //arrays to contain the values needed to avg over for inversion
27 instr *input; /* input parameters structure */
28 
29 int main(int argc, char* argv[]) {
30  long runCounter;
31  int32_t i, j, k;
32  int32_t iscan = 0; /* input scan number */
33  int32_t ipix; /* input pixel number */
34  int32_t npixs; /* number of matched pixels in the L2 file */
35  int32_t npixsOut; /* number of matched pixels to output */
36 
37  int32_t spix = 0; /* start pixel for subscene process */
38  int32_t epix = -1; /* end pixel for subscene process */
39  int32_t dpix;
40  int32_t sscan = 0; /* start scan for subscene process */
41  int32_t escan = -1; /* end scan for subscene process */
42  int32_t dscan;
43  int32_t nprod, nbands, status;
44 
45  float lon, lat;
46 
47  static l2_prod l2_str; /* input data for l2 calibrator */
48  meta_l2Type *meta_l2; /* input metadata for l2 calibrator */
49 
50  int32_t l2sensorID, l3sensorID;
51 
52  uint32 flagusemask;
53  uint32 required_mask;
54  char buf[5000];
55  idDS ds_id;
56 
57  // static instr input; /* input parameters structure */
58  calstr* pixrec; /* output detector-run structure */
59  std::vector<calstr*> pixrecArray; /* array of pixrec to output */
60  int nvars1d = 2; // initialize to 2, since lon and lat are always output
61  char vars1Dnames[L1_MAXPROD][32]; // floating point one dimensional data to include in output
62  char outprod[L1_MAXPROD][32];
63 
64  /* hold all of the command line options */
66 
67  if (argc == 1 || strcmp(argv[1], "-version") == 0
68  || strcmp(argv[1], "--version") == 0) {
69  want_verbose = 0;
73  exit(EXIT_FAILURE);
74  }
75  want_verbose = 1;
76 
77  /* Parse input parameters */
78  if (l2binmatch_input(argc, argv, list) != 0) {
79  printf("-E- %s: Error parsing input parameters.\n", argv[0]);
80  exit(EXIT_FAILURE);
81  }
82  printf("Opening %s\n", input->ifile[0]);
83  if (openL2(input->ifile[0], 0x0, &l2_str)) {
84  printf("-E- %s: Error reading L2 data %s.\n", argv[0], input->ifile[0]);
85  exit(EXIT_FAILURE);
86  }
87 
88  if ((meta_l2 = (meta_l2Type *) malloc(sizeof (meta_l2Type))) == NULL) {
89  printf("-E- %s: Error allocating memory for L2 metadata.\n", argv[0]);
90  exit(EXIT_FAILURE);
91  }
92 
93  if (readL2meta(meta_l2, 0)) {
94  printf("-E- %s: Error reading L2 metadata %s.\n", argv[0],
95  input->ifile[0]);
96  exit(EXIT_FAILURE);
97  }
98 
99  /* grab the target sensor name from the L2 meta-data */
100  l2sensorID = sensorName2SensorId(meta_l2->sensor_name);
101 
102  /*
103  * Check output file for matching sensorID
104  */
105  int existingSensorID;
106 
107  /* Save old HDF5 error handler */
108  H5E_auto_t old_func;
109  void *old_client_data;
110  H5Eget_auto(H5E_DEFAULT, &old_func, &old_client_data);
111 
112  /* Turn off error handling */
113  H5Eset_auto(H5E_DEFAULT, NULL, NULL);
114 
115  if (Hishdf(input->ofile[0]) == TRUE || H5Fis_hdf5(input->ofile[0]) == TRUE) {
116  ds_id = openDS(input->ofile[0]);
117  status = readAttr(ds_id, "sensorID", &existingSensorID);
118  if (status) {
119  printf("-E- %s: Problem reading output file: %s\n", argv[0],
120  input->ofile[0]);
121  exit(EXIT_FAILURE);
122  }
123  endDS(ds_id);
124  if (l2sensorID != existingSensorID) {
125  printf(
126  "-E- %s: Mixing L2 data from different sensors, %s and %s.\n",
127  argv[0], sensorId2SensorName(l2sensorID),
128  sensorId2SensorName(existingSensorID));
129  exit(EXIT_FAILURE);
130  }
131 
132  }
133  /* Restore previous HDF5 error handler */
134  H5Eset_auto(H5E_DEFAULT, old_func, old_client_data);
135 
136 
137 
138  // open the target L3 file
139  L3File* l3File = new L3File();
140  if (!l3File->open(input->tgtfile)) {
141  printf("-E- Could not open ifile=\"%s\".\n", input->tgtfile);
142  exit(EXIT_FAILURE);
143  }
144  l3sensorID = sensorName2SensorId(
145  l3File->getMetaData()->sensor_name);
146  if (l3sensorID == -1) {
147  l3sensorID = instrumentPlatform2SensorId(
148  l3File->getMetaData()->sensor,
149  l3File->getMetaData()->mission);
150  if (l3sensorID == -1) {
151  printf("-E- Unknown sensor name %s\n",
152  l3File->getMetaData()->sensor_name);
153  exit(EXIT_FAILURE);
154  }
155  }
156 
157  // Get L3 product list
158  nprod = l3File->getHdfBinObject()->nprod();
159  char *fullprodlist = (char *) malloc(l3File->getHdfBinObject()->query());
160 
161  l3File->getHdfBinObject()->query(fullprodlist);
162 
163  l3File->setActiveProductList(fullprodlist);
164 
165  //Setup flag mask
166  strcpy(buf, l2_str.flagnames);
167  setupflags(buf, input->flaguse, &flagusemask, &required_mask, &status,l2_str.l2_bits);
168  if (status < 0) {
169  printf("-E- %s: Error reading L2 flags %s.\n", argv[0], input->ifile[0]);
170  exit(EXIT_FAILURE);
171  }
172 
173  // initialize the depth file e.g., ETOPO1
174  dem_init();
175 
176  nbands = rdsensorinfo(l2sensorID, 0, NULL, NULL);
177 
178  // Set up list of products to match with L3 file
179  int32_t *l2prodinx;
180  if ((l2prodinx = (int32_t *) malloc(nprod * sizeof (int32_t))) == NULL) {
181  printf("-E- %s: Error allocating memory for L2 product index.\n",
182  argv[0]);
183  exit(EXIT_FAILURE);
184  }
185 
186  const char *l3prod;
187  for (i = 0; i < nprod; i++) {
188  l3prod = l3File->getHdfBinObject()->getProdName(i);
189  for (j = 0; j < l2_str.nprod; j++) {
190  if (strcmp(l2_str.prodname[j], l3prod) == 0) {
191  l2prodinx[i] = j;
192  strcpy(outprod[i], l3prod);
193  break;
194  }
195  }
196  }
197  // Set up list of products to output as 1D arrays (i.e. NOT matched to a L3 product)
198  strcpy(vars1Dnames[0], "lon"); // lon and lat you get by default
199  strcpy(vars1Dnames[1], "lat");
200 
201  int32_t oneDprodinx[L1_MAXPROD];
202  k = 2; //start at 2 since lon and lat are defacto
203  for (j = 0; j < l2_str.nprod; j++) {
204  int skip = 0;
205  for (i = 0; i < nprod; i++) {
206  if (strcmp(l2_str.prodname[j], outprod[i]) == 0) {
207  skip = 1;
208  break;
209  }
210  }
211  if (skip == 0 && strcmp(l2_str.prodname[j], "pixnum") != 0
212  && strcmp(l2_str.prodname[j], "mside") != 0
213  && strcmp(l2_str.prodname[j], "detnum") != 0
214  && strcmp(l2_str.prodname[j], "l2_flags") != 0
215  && strcmp(l2_str.prodname[j], "latitude") != 0
216  && strcmp(l2_str.prodname[j], "longitude") != 0) {
217  oneDprodinx[k++] = j;
218  strcpy(vars1Dnames[nvars1d++], l2_str.prodname[j]);
219  }
220  }
221 
222  /* Set the end pixel if it was not set by command argument */
223  if (l1_input->epixl == -1 || l1_input->epixl > l2_str.nsamp)
224  l1_input->epixl = l2_str.nsamp;
225  if (l1_input->eline == -1 || l1_input->eline > l2_str.nrec)
226  l1_input->eline = l2_str.nrec;
227  if (l1_input->spixl < 1)
228  l1_input->spixl = 1;
229  if (l1_input->sline < 1)
230  l1_input->sline = 1;
231 
232  spix = MAX(l1_input->spixl - 1, 0);
233  epix = MIN(l1_input->epixl - 1, l2_str.nsamp - 1);
234  dpix = MAX(l1_input->dpixl, 1);
235  sscan = MAX(l1_input->sline - 1, 0);
236  escan = MIN(l1_input->eline - 1, l2_str.nrec - 1);
237  dscan = MAX(l1_input->dline, 1);
238 
239  if (sscan > escan || spix > epix) {
240  printf("-E- %s: scan and pixel limits make no sense.\n", argv[0]);
241  printf(" start scan = %d\n", sscan + 1);
242  printf(" end scan = %d\n", escan + 1);
243  printf(" start pixel = %d\n", spix + 1);
244  printf(" end pixel = %d\n", epix + 1);
245  exit(EXIT_FAILURE);
246  }
247 
248  for (iscan = sscan; iscan <= escan; iscan += dscan) {
249 
250  if ((iscan % 100) == 0)
251  printf("Processing scan %d of %d\n", iscan, l2_str.nrec);
252 
253  if (readL2(&(l2_str), 0, iscan, -1, NULL)) {
254  printf("%s: Error reading L2 data file %s at scan %d.\n", argv[0],
255  l2_str.filename, iscan);
256  continue;
257  }
258  for (ipix = spix; ipix <= epix; ipix += dpix) {
259  // Make sure pixels pass the flags defined as masks in flaguse
260  if ((l2_str.l2_flags[ipix] & flagusemask) == 0) {
261 
262  lat = l2_str.latitude[ipix];
263  lon = l2_str.longitude[ipix];
264 
265  L3Bin* l3Bin = l3File->getClosestBin(lat, lon);
266 
267  if (l3Bin) {
268  /* check:
269  * depth is greater than input depth
270  * minimum nscenes is met
271  * minimum nsamples is met
272  */
273 
274  if (input->vcal_depth > get_dem(lat, lon)
275  && l3Bin->getNscenes() >= input->vcal_min_nscene
276  && l3Bin->getNobs() >= input->vcal_min_nbin) {
277  pixrec = alloc_calrec(2, nbands, nprod, nvars1d);
278  pixrec->sensorID = l2sensorID;
279  pixrec->year = l2_str.year;
280  pixrec->day = l2_str.day;
281  pixrec->msec = l2_str.msec;
282  pixrec->iscan = iscan;
283  pixrec->mside = l2_str.mside[iscan];
284  pixrec->detnum = l2_str.detnum[iscan];
285  if (l2_str.pixnum)
286  pixrec->pixnum = l2_str.pixnum[ipix];
287  else
288  pixrec->pixnum = ipix;
289  pixrec->vars1D[0] = lon;
290  pixrec->vars1D[1] = lat;
291  // populate the "2D" matched pixel arrays
292  for (j = 0; j < nprod; j++) {
293  pixrec->data[j][0] = l2_str.l2_data[l2prodinx[j]][ipix];
294  pixrec->data[j][1] = l3Bin->getMean(j);
295  }
296  // populate the "1D" arrays
297  for (j = 2; j < nvars1d; j++) {
298  pixrec->vars1D[j] = l2_str.l2_data[oneDprodinx[j]][ipix];
299  }
300 
301  pixrecArray.push_back(pixrec);
302  }
303  }
304  }
305  }
306  }
307  npixs = pixrecArray.size();
308  printf("Number of pixels matched: %d\n", npixs);
309 
310  if (npixs > 2 * input->subsamp) {
311  j = 0;
312  npixsOut = npixs / input->subsamp;
313  ds_id = calfile_open(input->ofile[0], l2sensorID, npixsOut, 2, nprod, nvars1d,
314  outprod, vars1Dnames, &runCounter, BINMATCH);
315  for (i = 0; i < npixs; i += input->subsamp) {
316  if (j < npixsOut) {
317  pixrec = pixrecArray[i];
318  calfile_write(ds_id, pixrec, runCounter++, npixsOut, 2, nprod, nbands,
319  nvars1d, outprod, vars1Dnames, BINMATCH);
320  free_calrec(pixrec, nbands, nprod);
321  j++;
322  }
323  }
324 
325  if (calfile_close(ds_id) != 0) {
326  printf("Trouble closing file %s\n", input->ofile[0]);
327  exit(EXIT_FAILURE);
328  } else {
329  printf("Number of pixels output: %d\n", npixsOut);
330  printf("Done processing file %s\n", input->ofile[0]);
331  }
332  } // Free up some memory
333 
334  closeL2(&l2_str, 0);
335  freeL2meta(meta_l2);
336 
337  freeL2(&l2_str);
338 
339  exit(EXIT_SUCCESS);
340 
341 }
#define L1_MAXPROD
Definition: filehandle.h:20
void setupflags(char *flagdef, char *flaguse, uint32_t *flagusemask, uint32_t *required, int *status, int32_t *BITS)
Definition: setupflags.c:7
float *** aots
#define MAX(A, B)
Definition: swl0_utils.h:25
int l2binmatch_read_options(clo_optionList_t *list, int argc, char *argv[])
#define MIN(x, y)
Definition: rice.h:169
#define EXIT_SUCCESS
Definition: GEO_basic.h:72
virtual L3Bin * getClosestBin(float lat, float lon)
Definition: L3File.cpp:763
int j
Definition: decode_rs.h:73
int status
Definition: l1_czcs_hdf.c:32
int instrumentPlatform2SensorId(const char *instrument, const char *platform)
Definition: sensorInfo.c:405
int l2binmatch_input(int argc, char **argv, clo_optionList_t *list)
idDS openDS(const char *filename)
Definition: wrapper.c:606
int calfile_write(idDS ds_id, calstr *calrec, int recnum, int ydim, int xdim, int nprods, int nbands, int nvars1d, char l2prods[L1_MAXPROD][32], char vars1Dnames[L1_MAXPROD][32], caltype ctype)
#define NULL
Definition: decode_rs.h:63
char mission[SM_ATTRSZ]
Definition: meta_l3b.h:21
int32_t readL2meta(meta_l2Type *meta_l2, int32_t ifile)
Definition: readL2scan.c:2151
int dem_init()
Definition: read_mask.c:43
virtual bool open(const char *fileName)
Definition: L3File.cpp:454
int sensorName2SensorId(const char *name)
Definition: sensorInfo.c:371
#define TRUE
Definition: rice.h:165
int32_t freeL2meta(meta_l2Type *meta_l2)
Definition: readL2scan.c:2325
int32_t readL2(l2_prod *l2_str, int32_t ifile, int32_t recnum, int32_t iprod, unsigned char *scan_in_rowgroup)
Definition: readL2scan.c:1336
int main(int argc, char *argv[])
int32_t closeL2(l2_prod *l2_str, int32_t ifile)
Definition: readL2scan.c:1995
int readAttr(idDS ds_id, const char *nam, void *data)
virtual bool setActiveProductList(const char *prodStr)
Definition: L3File.cpp:615
virtual int query()
Definition: bin_io.cpp:2764
char sensor[SM_ATTRSZ]
Definition: meta_l3b.h:23
clo_optionList_t * clo_createList()
Definition: clo.c:532
list(APPEND LIBS ${NETCDF_LIBRARIES}) find_package(GSL REQUIRED) include_directories($
Definition: CMakeLists.txt:8
int32_t openL2(const char *fname, const char *plist, l2_prod *l2_str)
Definition: readL2scan.c:316
int16_t fileID
virtual int32_t nprod()
Definition: hdf_bin.h:110
void clo_printUsage(clo_optionList_t *list)
Definition: clo.c:1988
void free_calrec(calstr *calrec, int nbands, int nprods)
l1_input_t * l1_input
Definition: l1_options.c:9
virtual meta_l3bType * getMetaData()
Definition: L3File.cpp:499
int want_verbose
@ BINMATCH
Definition: calfile_utils.h:21
float *** nlws
int32_t getNobs() const
Definition: L3File.cpp:194
float getMean(int32_t prodId=0) const
Definition: L3File.cpp:263
calstr * alloc_calrec(int ydim, int nbands, int nprods, int nvar1d)
int32 dpix
Definition: l1_czcs_hdf.c:22
int32_t nbands
const char * sensorId2SensorName(int sensorId)
Definition: sensorInfo.c:273
virtual const char * getProdName(int prodNum) const
Definition: bin_io.cpp:2831
int32 spix
Definition: l1_czcs_hdf.c:21
int32_t iscan
virtual Hdf::hdf_bin * getHdfBinObject() const
Definition: L3File.cpp:808
Definition: dfutils.h:28
int32 epix
Definition: l1_czcs_hdf.c:23
instr * input
int32_t rdsensorinfo(int32_t, int32_t, const char *, void **)
Definition: rdsensorinfo.c:69
Definition: L3File.cpp:10
int calfile_close(idDS ds_id)
float get_dem(float lat, float lon)
Definition: read_mask.c:112
int endDS(idDS ds_id)
Definition: wrapper.c:624
idDS calfile_open(char *ofile, int sensorID, int ydim, int xdim, int nprods, int nvars1d, char l2prods[L1_MAXPROD][32], char vars1Dnames[L1_MAXPROD][32], long *numExistingRecords, caltype ctype)
Definition: calfile_utils.c:48
char sensor_name[SM_ATTRSZ]
Definition: meta_l3b.h:19
int32_t getNscenes() const
Definition: L3File.cpp:202
int i
Definition: decode_rs.h:71
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
int32_t freeL2(l2_prod *l2_str)
Definition: readL2scan.c:2083
int l2binmatch_init_options(clo_optionList_t *list)
int k
Definition: decode_rs.h:73