OB.DAAC Logo
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;
71  l2binmatch_read_options(list, argc, argv);
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 
83  printf("Opening %s\n", input->ifile[0]);
84  if (openL2(input->ifile[0], 0x0, &l2_str)) {
85  printf("-E- %s: Error reading L2 data %s.\n", argv[0], input->ifile[0]);
86  exit(EXIT_FAILURE);
87  }
88 
89  if ((meta_l2 = (meta_l2Type *) malloc(sizeof (meta_l2Type))) == NULL) {
90  printf("-E- %s: Error allocating memory for L2 metadata.\n", argv[0]);
91  exit(EXIT_FAILURE);
92  }
93 
94  if (readL2meta(meta_l2, 0)) {
95  printf("-E- %s: Error reading L2 metadata %s.\n", argv[0],
96  input->ifile[0]);
97  exit(EXIT_FAILURE);
98  }
99 
100  /* grab the target sensor name from the L2 meta-data */
101  l2sensorID = sensorName2SensorId(meta_l2->sensor_name);
102 
103  /*
104  * Check output file for matching sensorID
105  */
106  int existingSensorID;
107 
108  /* Save old HDF5 error handler */
109  H5E_auto_t old_func;
110  void *old_client_data;
111  H5Eget_auto(H5E_DEFAULT, &old_func, &old_client_data);
112 
113  /* Turn off error handling */
114  H5Eset_auto(H5E_DEFAULT, NULL, NULL);
115 
116  if (Hishdf(input->ofile[0]) == TRUE || H5Fis_hdf5(input->ofile[0]) == TRUE) {
117  ds_id = openDS(input->ofile[0]);
118  status = readAttr(ds_id, "sensorID", &existingSensorID);
119  if (status) {
120  printf("-E- %s: Problem reading output file: %s\n", argv[0],
121  input->ofile[0]);
122  exit(EXIT_FAILURE);
123  }
124  endDS(ds_id);
125  if (l2sensorID != existingSensorID) {
126  printf(
127  "-E- %s: Mixing L2 data from different sensors, %s and %s.\n",
128  argv[0], sensorId2SensorName(l2sensorID),
129  sensorId2SensorName(existingSensorID));
130  exit(EXIT_FAILURE);
131  }
132 
133  }
134  /* Restore previous HDF5 error handler */
135  H5Eset_auto(H5E_DEFAULT, old_func, old_client_data);
136 
137  // open the target L3 file
138  L3File* l3File = new L3File();
139  if (!l3File->open(input->tgtfile)) {
140  printf("-E- Could not open ifile=\"%s\".\n", input->tgtfile);
141  exit(EXIT_FAILURE);
142  }
143  l3sensorID = sensorName2SensorId(
144  l3File->getMetaData()->sensor_name);
145  if (l3sensorID == -1) {
146  l3sensorID = instrumentPlatform2SensorId(
147  l3File->getMetaData()->sensor,
148  l3File->getMetaData()->mission);
149  if (l3sensorID == -1) {
150  printf("-E- Unknown sensor name %s\n",
151  l3File->getMetaData()->sensor_name);
152  exit(EXIT_FAILURE);
153  }
154  }
155 
156  // Get L3 product list
157  nprod = l3File->getHdfBinObject()->nprod();
158  char *fullprodlist = (char *) malloc(l3File->getHdfBinObject()->query());
159 
160  l3File->getHdfBinObject()->query(fullprodlist);
161 
162  l3File->setActiveProductList(fullprodlist);
163 
164  //Setup flag mask
165  strcpy(buf, l2_str.flagnames);
166  setupflags(buf, input->flaguse, &flagusemask, &required_mask, &status);
167  if (status < 0) {
168  printf("-E- %s: Error reading L2 flags %s.\n", argv[0], input->ifile[0]);
169  exit(EXIT_FAILURE);
170  }
171 
172  // initialize the depth file e.g., ETOPO1
173  dem_init();
174 
175  nbands = rdsensorinfo(l2sensorID, 0, NULL, NULL);
176 
177  // Set up list of products to match with L3 file
178  int32_t *l2prodinx;
179  if ((l2prodinx = (int32_t *) malloc(nprod * sizeof (int32_t))) == NULL) {
180  printf("-E- %s: Error allocating memory for L2 product index.\n",
181  argv[0]);
182  exit(EXIT_FAILURE);
183  }
184 
185  const char *l3prod;
186  for (i = 0; i < nprod; i++) {
187  l3prod = l3File->getHdfBinObject()->getProdName(i);
188  for (j = 0; j < l2_str.nprod; j++) {
189  if (strcmp(l2_str.prodname[j], l3prod) == 0) {
190  l2prodinx[i] = j;
191  strcpy(outprod[i], l3prod);
192  break;
193  }
194  }
195  }
196  // Set up list of products to output as 1D arrays (i.e. NOT matched to a L3 product)
197  strcpy(vars1Dnames[0], "lon"); // lon and lat you get by default
198  strcpy(vars1Dnames[1], "lat");
199 
200  int32_t oneDprodinx[L1_MAXPROD];
201  k = 2; //start at 2 since lon and lat are defacto
202  for (j = 0; j < l2_str.nprod; j++) {
203  int skip = 0;
204  for (i = 0; i < nprod; i++) {
205  if (strcmp(l2_str.prodname[j], outprod[i]) == 0) {
206  skip = 1;
207  break;
208  }
209  }
210  if (skip == 0 && strcmp(l2_str.prodname[j], "pixnum") != 0
211  && strcmp(l2_str.prodname[j], "mside") != 0
212  && strcmp(l2_str.prodname[j], "detnum") != 0
213  && strcmp(l2_str.prodname[j], "l2_flags") != 0
214  && strcmp(l2_str.prodname[j], "latitude") != 0
215  && strcmp(l2_str.prodname[j], "longitude") != 0) {
216  oneDprodinx[k++] = j;
217  strcpy(vars1Dnames[nvars1d++], l2_str.prodname[j]);
218  }
219  }
220 
221  /* Set the end pixel if it was not set by command argument */
222  if (l1_input->epixl == -1 || l1_input->epixl > l2_str.nsamp)
223  l1_input->epixl = l2_str.nsamp;
224  if (l1_input->eline == -1 || l1_input->eline > l2_str.nrec)
225  l1_input->eline = l2_str.nrec;
226  if (l1_input->spixl < 1)
227  l1_input->spixl = 1;
228  if (l1_input->sline < 1)
229  l1_input->sline = 1;
230 
231  spix = MAX(l1_input->spixl - 1, 0);
232  epix = MIN(l1_input->epixl - 1, l2_str.nsamp - 1);
233  dpix = MAX(l1_input->dpixl, 1);
234  sscan = MAX(l1_input->sline - 1, 0);
235  escan = MIN(l1_input->eline - 1, l2_str.nrec - 1);
236  dscan = MAX(l1_input->dline, 1);
237 
238  if (sscan > escan || spix > epix) {
239  printf("-E- %s: scan and pixel limits make no sense.\n", argv[0]);
240  printf(" start scan = %d\n", sscan + 1);
241  printf(" end scan = %d\n", escan + 1);
242  printf(" start pixel = %d\n", spix + 1);
243  printf(" end pixel = %d\n", epix + 1);
244  exit(EXIT_FAILURE);
245  }
246 
247  for (iscan = sscan; iscan <= escan; iscan += dscan) {
248 
249  if ((iscan % 100) == 0)
250  printf("Processing scan %d of %d\n", iscan, l2_str.nrec);
251 
252  if (readL2(&(l2_str), 0, iscan, -1, NULL)) {
253  printf("%s: Error reading L2 data file %s at scan %d.\n", argv[0],
254  l2_str.filename, iscan);
255  continue;
256  }
257  for (ipix = spix; ipix <= epix; ipix += dpix) {
258  // Make sure pixels pass the flags defined as masks in flaguse
259  if ((l2_str.l2_flags[ipix] & flagusemask) == 0) {
260 
261  lat = l2_str.latitude[ipix];
262  lon = l2_str.longitude[ipix];
263 
264  L3Bin* l3Bin = l3File->getClosestBin(lat, lon);
265 
266  if (l3Bin) {
267  /* check:
268  * depth is greater than input depth
269  * minimum nscenes is met
270  * minimum nsamples is met
271  */
272 
273  if (input->vcal_depth > get_dem(lat, lon)
274  && l3Bin->getNscenes() >= input->vcal_min_nscene
275  && l3Bin->getNobs() >= input->vcal_min_nbin) {
276  pixrec = alloc_calrec(2, nbands, nprod, nvars1d);
277  pixrec->sensorID = l2sensorID;
278  pixrec->year = l2_str.year;
279  pixrec->day = l2_str.day;
280  pixrec->msec = l2_str.msec;
281  pixrec->iscan = iscan;
282  pixrec->mside = l2_str.mside[iscan];
283  pixrec->detnum = l2_str.detnum[iscan];
284  if (l2_str.pixnum)
285  pixrec->pixnum = l2_str.pixnum[ipix];
286  else
287  pixrec->pixnum = ipix;
288  pixrec->vars1D[0] = lon;
289  pixrec->vars1D[1] = lat;
290  // populate the "2D" matched pixel arrays
291  for (j = 0; j < nprod; j++) {
292  pixrec->data[j][0] = l2_str.l2_data[l2prodinx[j]][ipix];
293  pixrec->data[j][1] = l3Bin->getMean(j);
294  }
295  // populate the "1D" arrays
296  for (j = 2; j < nvars1d; j++) {
297  pixrec->vars1D[j] = l2_str.l2_data[oneDprodinx[j]][ipix];
298  }
299 
300  pixrecArray.push_back(pixrec);
301  }
302  }
303  }
304  }
305  }
306  npixs = pixrecArray.size();
307  printf("Number of pixels matched: %d\n", npixs);
308 
309  if (npixs > 2 * input->subsamp) {
310  j = 0;
311  npixsOut = npixs / input->subsamp;
312  ds_id = calfile_open(input->ofile[0], l2sensorID, npixsOut, 2, nprod, nvars1d,
313  outprod, vars1Dnames, &runCounter, BINMATCH);
314  for (i = 0; i < npixs; i += input->subsamp) {
315  if (j < npixsOut) {
316  pixrec = pixrecArray[i];
317  calfile_write(ds_id, pixrec, runCounter++, npixsOut, 2, nprod, nbands,
318  nvars1d, outprod, vars1Dnames, BINMATCH);
319  free_calrec(pixrec, nbands, nprod);
320  j++;
321  }
322  }
323 
324  if (calfile_close(ds_id) != 0) {
325  printf("Trouble closing file %s\n", input->ofile[0]);
326  exit(EXIT_FAILURE);
327  } else {
328  printf("Number of pixels output: %d\n", npixsOut);
329  printf("Done processing file %s\n", input->ofile[0]);
330  }
331  } // Free up some memory
332 
333  closeL2(&l2_str, 0);
334  freeL2meta(meta_l2);
335 
336  freeL2(&l2_str);
337 
338  exit(EXIT_SUCCESS);
339 
340 }
#define L1_MAXPROD
Definition: filehandle.h:20
float *** aots
int32_t openL2(const char *fname, char *plist, l2_prod *l2_str)
Definition: readL2scan.c:296
#define MAX(A, B)
Definition: swl0_utils.h:26
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:739
int j
Definition: decode_rs.h:73
int status
Definition: l1_czcs_hdf.c:32
list(APPEND LIBS ${PGSTK_LIBRARIES}) add_executable(atteph_info_modis atteph_info_modis.c) target_link_libraries(atteph_info_modis $
Definition: CMakeLists.txt:7
int instrumentPlatform2SensorId(const char *instrument, const char *platform)
Definition: sensorInfo.c:302
int l2binmatch_input(int argc, char **argv, clo_optionList_t *list)
idDS openDS(const char *filename)
Definition: wrapper.c:616
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:2081
int dem_init()
Definition: read_mask.c:43
virtual bool open(const char *fileName)
Definition: L3File.cpp:430
int sensorName2SensorId(const char *name)
Definition: sensorInfo.c:268
#define TRUE
Definition: rice.h:165
float * lat
int32_t freeL2meta(meta_l2Type *meta_l2)
Definition: readL2scan.c:2255
int32_t readL2(l2_prod *l2_str, int32_t ifile, int32_t recnum, int32_t iprod, unsigned char *scan_in_rowgroup)
Definition: readL2scan.c:1250
int main(int argc, char *argv[])
int32_t closeL2(l2_prod *l2_str, int32_t ifile)
Definition: readL2scan.c:1936
int readAttr(idDS ds_id, const char *nam, void *data)
virtual bool setActiveProductList(const char *prodStr)
Definition: L3File.cpp:591
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
void setupflags(char *flagdef, char *flaguse, uint32_t *flagusemask, uint32_t *required, int *status)
Definition: setupflags.c:5
int16_t fileID
virtual int32_t nprod()
Definition: hdf_bin.h:108
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:475
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:198
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:784
Definition: dfutils.h:28
float * lon
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:634
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:2026
int l2binmatch_init_options(clo_optionList_t *list)
int k
Definition: decode_rs.h:73