ocssw V2020
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[MAXPROD][32]; // floating point one dimensional data to include in output
62  char outprod[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  elev_init(input->elevfile, NULL);
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[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 (input->epixl == -1 || input->epixl > l2_str.nsamp)
223  input->epixl = l2_str.nsamp;
224  if (input->eline == -1 || input->eline > l2_str.nrec)
225  input->eline = l2_str.nrec;
226  if (input->spixl < 1)
227  input->spixl = 1;
228  if (input->sline < 1)
229  input->sline = 1;
230 
231  spix = MAX(input->spixl - 1, 0);
232  epix = MIN(input->epixl - 1, l2_str.nsamp - 1);
233  dpix = MAX(input->dpixl, 1);
234  sscan = MAX(input->sline - 1, 0);
235  escan = MIN(input->eline - 1, l2_str.nrec - 1);
236  dscan = MAX(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_elev(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]
293  * l2_str.nsamp + 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]
299  * l2_str.nsamp + ipix];
300  }
301 
302  pixrecArray.push_back(pixrec);
303  }
304  }
305  }
306  }
307  }
308  npixs = pixrecArray.size();
309  printf("Number of pixels matched: %d\n", npixs);
310 
311  if (npixs > 2 * input->subsamp) {
312  j = 0;
313  npixsOut = npixs / input->subsamp;
314  ds_id = calfile_open(input->ofile[0], l2sensorID, npixsOut, 2, nprod, nvars1d,
315  outprod, vars1Dnames, &runCounter, BINMATCH);
316  for (i = 0; i < npixs; i += input->subsamp) {
317  if (j < npixsOut) {
318  pixrec = pixrecArray[i];
319  calfile_write(ds_id, pixrec, runCounter++, npixsOut, 2, nprod, nbands,
320  nvars1d, outprod, vars1Dnames, BINMATCH);
321  free_calrec(pixrec, nbands, nprod);
322  j++;
323  }
324  }
325 
326  if (calfile_close(ds_id) != 0) {
327  printf("Trouble closing file %s\n", input->ofile[0]);
328  exit(EXIT_FAILURE);
329  } else {
330  printf("Number of pixels output: %d\n", npixsOut);
331  printf("Done processing file %s\n", input->ofile[0]);
332  }
333  } // Free up some memory
334 
335  closeL2(&l2_str, 0);
336  freeL2meta(meta_l2);
337 
338  freeL2(&l2_str);
339 
340  exit(EXIT_SUCCESS);
341 
342 }
float *** aots
#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:31
int instrumentPlatform2SensorId(const char *instrument, const char *platform)
Definition: sensorInfo.c:279
int l2binmatch_input(int argc, char **argv, clo_optionList_t *list)
float get_elev(float lat, float lon)
Definition: elev.c:108
idDS openDS(const char *filename)
Definition: wrapper.c:612
#define NULL
Definition: decode_rs.h:63
char mission[SM_ATTRSZ]
Definition: meta_l3b.h:21
idDS calfile_open(char *ofile, int sensorID, int ydim, int xdim, int nprods, int nvars1d, char l2prods[MAXPROD][32], char vars1Dnames[MAXPROD][32], long *numExistingRecords, caltype ctype)
Definition: calfile_utils.c:48
int32 readL2meta(meta_l2Type *meta_l2, int32 ifile)
Definition: readL2scan.c:2016
virtual bool open(const char *fileName)
Definition: L3File.cpp:430
int sensorName2SensorId(const char *name)
Definition: sensorInfo.c:245
list(APPEND LIBS ${NETCDF_LIBRARIES}) list(APPEND LIBS $
Definition: CMakeLists.txt:9
#define TRUE
Definition: rice.h:165
float * lat
int main(int argc, char *argv[])
virtual bool setActiveProductList(const char *prodStr)
Definition: L3File.cpp:591
virtual int query()
Definition: bin_io.cpp:2763
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
int32_t rdsensorinfo(int32_t sensorID, int32_t evalmask, const char *pname, void **pval)
Definition: rdsensorinfo.c:69
int16_t fileID
virtual int32_t nprod()
Definition: hdf_bin.h:108
void clo_printUsage(clo_optionList_t *list)
Definition: clo.c:1955
void free_calrec(calstr *calrec, int nbands, int nprods)
int32 readL2(l2_prod *l2_str, int32 ifile, int32 recnum, int32 iprod, unsigned char *scan_in_rowgroup)
Definition: readL2scan.c:1197
virtual meta_l3bType * getMetaData()
Definition: L3File.cpp:475
int want_verbose
#define MAXPROD
Definition: l2prod.h:4
@ 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:21
int32_t nbands
const char * sensorId2SensorName(int sensorId)
Definition: sensorInfo.c:175
virtual const char * getProdName(int prodNum) const
Definition: bin_io.cpp:2830
int32 spix
Definition: l1_czcs_hdf.c:20
int32_t iscan
virtual Hdf::hdf_bin * getHdfBinObject() const
Definition: L3File.cpp:784
Definition: dfutils.h:29
float * lon
int readAttr(idDS ds_id, const char *nam, VOIDP data)
Definition: wrapper.c:71
int32 epix
Definition: l1_czcs_hdf.c:22
int32 openL2(const char *fname, char *plist, l2_prod *l2_str)
Definition: readL2scan.c:276
instr * input
int32 freeL2meta(meta_l2Type *meta_l2)
Definition: readL2scan.c:2186
Definition: L3File.cpp:10
int calfile_write(idDS ds_id, calstr *calrec, int recnum, int ydim, int xdim, int nprods, int nbands, int nvars1d, char l2prods[MAXPROD][32], char vars1Dnames[MAXPROD][32], caltype ctype)
int calfile_close(idDS ds_id)
int32 freeL2(l2_prod *l2_str)
Definition: readL2scan.c:1963
int endDS(idDS ds_id)
Definition: wrapper.c:630
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")
int l2binmatch_init_options(clo_optionList_t *list)
int k
Definition: decode_rs.h:73
int32 closeL2(l2_prod *l2_str, int32 ifile)
Definition: readL2scan.c:1873
void elev_init(char *elevGlobalFilename, char *elevAuxFilename)
Definition: elev.c:66