OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
prismbil2oci.cpp
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <string.h>
5 #include <unistd.h>
6 
7 #include <sstream>
8 #include <iomanip>
9 #include <stdint.h>
10 #include <iostream>
11 #include <fstream>
12 
13 #include "cdl_utils.h"
14 #include "netcdf.h"
15 #include "nc4utils.h"
16 #include "nc_gridutils.h"
17 #include "prism.h"
18 #include "rdprism.h"
19 #include <gsl/gsl_errno.h>
20 #include <gsl/gsl_spline.h>
21 #include <gsl/gsl_sort_double.h>
22 #include <timeutils.h>
23 
24 #define SENSOR "ocip"
25 #define CDLFILE "OCIP_Level-1B_Data_Structure.cdl"
26 #define NBANDS 285 // number of OCI bands
27 
28 extern "C" void handle_error(int status);
29 
30 using namespace std;
31 // Modification history:
32 // Programmer Organization Date Ver Description of change
33 // ---------- ------------ ---- --- ---------------------
34 // Joel Gales FutureTech 12/18/14 0.01 Original development
35 // Joel Gales FutureTech 01/12/15 0.10 Add support for writing
36 // interpolated ORCA Lt values
37 // and navigation
38 // Rick Healy SAIC 06/06/16 1.00 Integrated L2gen L1 reader for Aviris
39 // Rick Healy SAIC 07/15/16 1.00 Integrated L2gen L1 reader for Prism
40 
41 
42 #define VERSION "1.00"
43 
44 size_t interpLt(size_t startJ, size_t npixels, size_t k, size_t nbMax,
45  float ociaWavelength, size_t iwave, prism4ocia_t* data,
46  gsl_interp_accel* acc, float* Lt) {
47  size_t j, n;
48  size_t i0 = startJ, i1 = startJ + 5;
49  float minWave, maxWave;
50 
51  double* x = (double*) (calloc(nbMax, sizeof (double)));
52  double* y = (double*) (calloc(nbMax, sizeof (double)));
53  size_t iter = 0;
54  while (i0 >= 0 && i1 < nbMax) {
55  n = 0;
56  minWave = 99999;
57  maxWave = -1;
58  for (j = i0; j < i1; j++) {
59  if (data->Lt[j * npixels + k] > 0) {
60  x[n] = data->wave[j];
61  y[n] = data->Lt[j * npixels + k];
62  n++;
63  if (minWave > data->wave[j]) minWave = data->wave[j];
64  if (maxWave < data->wave[j]) maxWave = data->wave[j];
65  }
66  }
67 
68  if (n > 2 && minWave <= ociaWavelength && maxWave >= ociaWavelength) break;
69  if (minWave > ociaWavelength || n < 3)
70  i0--;
71  if (maxWave < ociaWavelength || n < 3)
72  i1++;
73  iter++;
74  }
75  if (n > 2) {
76  // Allocate space for spline
77  gsl_spline* spline = gsl_spline_alloc(gsl_interp_cspline, n);
78  gsl_sort2(x, 1, y, 1, n);
79  if (ociaWavelength > x[0]) {
80  // Sort wavelengths and corresponding Lt values
81  gsl_sort2(x, 1, y, 1, n);
82  // Initiate spline
83  gsl_spline_init(spline, x, y, n);
84  // Generate interpolated values
85  double yi = gsl_spline_eval(spline, ociaWavelength, acc);
86  Lt[iwave * npixels + k] = yi;
87  } else {
88  Lt[iwave * npixels + k] = -9999;
89  }
90  iwave++;
91  gsl_spline_free(spline);
92  }
93  free(x);
94  free(y);
95  return iwave;
96 }
97 
98 int main(int argc, char* argv[]) {
99  int status;
100  size_t sline = 0, eline = 0;
101  char filename[FILENAME_MAX];
102  char *filedir;
103  double scantime_first, scantime_last;
104  static int first_good_scantime = 0;
105  char isodatetime_first[30], isodatetime_last[30];
106 
107  cout << "prismbil2oci " << VERSION << " ("
108  << __DATE__ << " " << __TIME__ << ")" << endl;
109 
110  if (argc < 3) {
111  cout << endl <<
112  "prismbil2oci input_PRISM_file output_ORCA_file <sline (default=0)> <eline (default=last row)>"
113  << endl;
114  return 0;
115  }
116 
117  if (argc >= 4) {
118  sscanf(argv[3], "%ld", &sline);
119  if (sline < 0) sline = 0;
120  }
121 
122  if (argc >= 5) {
123  sscanf(argv[4], "%ld", &eline);
124  if (eline < 0) eline = 0;
125  }
126  // prism4orca_t *&data;
127 
128  if ((filedir = getenv("OCDATAROOT")) == NULL) {
129  printf("-E- %s: OCDATAROOT env variable undefined.\n", __FILE__);
130  return (-1);
131  }
132  strcpy(filename, filedir);
133 
134  strcat(filename, "/");
135  strcat(filename, SENSOR);
136  strcat(filename, "/");
137 
138  strcat(filename, CDLFILE);
139 
140  prism4ocia_t **data, *temp;
141  float fwhm[NBANDS];
142  data = (prism4ocia_t **) malloc(sizeof (prism4ocia_t *));
143  temp = open_prism(argv[1], data);
144  size_t nscans, npixels = temp->npix;
145 
146  if (eline <= 0) eline = temp->nscan;
147 
148  nscans = eline - sline;
149 
150  int ncid, grpID;
151  //char** dim_names, uint32_t* dim_size, size_t n_dims
152  const char *dimNames[] = {"number_of_lines", "pixels_per_line", "number_of_bands"};
153  size_t dimSize[] = {nscans, npixels, NBANDS};
154  size_t n_dims = 3;
155  char time_coverage_start[64] = "yyyy-mm-ddThh:mm:ssZ";
156 
157  // Create ORCA output file
158  static ncdfFile ociafile;
159  ociafile.cdlCreateDim(argv[2], filename, dimNames, dimSize, n_dims, nscans);
160 
161  // nc_def_dim(orcafile.ncid, "number_of_lines", nscans, &dimid[0]);
162  // nc_def_dim(orcafile.ncid, "pixels_per_line", nscans, &dimid[1]);
163  //
164  // Copy prism nav fields to ORCA
165  //
166  ncid = ociafile.getNcid();
167 
168  sprintf(time_coverage_start, "%4.4d-%2.2d-%2.2dT%2.2d:%2.2d:%2.2dZ", temp->year, temp->month, temp->day, temp->hour, temp->min, (int) temp->sec);
169  status = nc_put_att_text(ncid, NC_GLOBAL, "time_coverage_start", strlen(time_coverage_start) + 1, time_coverage_start);
170  // status = nc_put_att_text(ncid, NC_GLOBAL, "time_coverage_end", strlen(time_coverage_start)+1, time_coverage_start);
171  if (status != NC_NOERR) {
172  printf("-E- %s %d: %s for %s\n",
173  __FILE__, __LINE__, nc_strerror(status), "time_coverage_start");
174  exit(1);
175  }
176 
177  sprintf(time_coverage_start, "%4.4d-%2.2d-%2.2dT%2.2d:%2.2d:%2.2dZ", temp->year, temp->month, temp->day, temp->hour, temp->min, (int) temp->sec);
178  status = nc_put_att_text(ncid, NC_GLOBAL, "time_coverage_end", strlen(time_coverage_start) + 1, time_coverage_start);
179 
180 
181 
182  size_t nbands_av = temp->numBands;
183  size_t nbands_ocia = NBANDS;
184 
185  size_t n = 0, iwave = 0;
186  istringstream istr;
187  string str;
188 
189 
190  // Allocate space for prism Lt (short) values
191 
192  gsl_interp_accel *acc = gsl_interp_accel_alloc();
193 
194  float *Lt = (float *) calloc(npixels*nbands_ocia, sizeof (float));
195  float *ociaWavelength = (float *) calloc(nbands_ocia, sizeof (float));
196  double* utc = (double*) (calloc(nscans, sizeof (double)));
197  size_t j;
198 
199  const char *navfieldsOCIP[] = {"sena", "senz", "sola", "solz", "lat", "lon"};
200 
201  int varID;
202 
203  size_t startOCIP[2] = {0, 0};
204  size_t countOCIP[2] = {1, npixels};
205  size_t i;
206 
207 
208  for (size_t iscan = sline; iscan < eline; iscan++) {
209 
210  grpID = ociafile.getGid("observation_data");
211 
212  // Read prism Lt's for this scan
213  if (read_prism(temp, iscan)) {
214 
215  cout << "Can't read file " << argv[1] << endl;
216  exit(1);
217 
218  }
219 
220  if ((iscan % 100) == 0) {
221  cout << "Processing scan: " << iscan << " out of " << nscans << endl;
222  printf("year=%d doy=%d msec=%d month=%d day=%d hour=%d min=%d sec=%f\n ", temp->year, temp->doy, temp->msec, temp->month, temp->day, temp->hour, temp->min, temp->sec);
223  }
224  // For each pixel in scan ...
225  for (size_t k = 0; k < npixels; k++) {
226 
227  // Get number of good prism wavelengths for each pixel
228  n = 0;
229  j = 0;
230  for (j = 0; j < nbands_av; j++) {
231  if (temp->Lt[j * npixels + k] > 0) n++;
232  }
233 
234  // Bail if no good prism Lt values
235  if (n == 0) continue;
236 
237 
238  // Generate x & y arrays for spline
239  // x - prism wavelengths
240  // y - prism Lt values
241  iwave = 0;
242  for (j = 0; j < nbands_av; j++) {
243  ociaWavelength[iwave] = temp->wave[j];
244  Lt[iwave * npixels + k] = temp->Lt[k * nbands_av + j];
245  // Lt[iwave*npixels + k] = temp->Lt[j*npixels+k];
246  fwhm[iwave] = temp->fwhm[j];
247  iwave++;
248  }
249  // ociaWavelength[iwave] = 936.436901;
250  // fwhm[iwave] = 32.4702967;
251  // Lt[iwave*npixels + k] = (11.9269624*temp->Lt[222*npixels+k] +
252  // 10.78153 *temp->Lt[223*npixels+k] +
253  // 9.76180425*temp->Lt[224*npixels+k])/fwhm[iwave];
254  // iwave++;
255  //
256  // ociaWavelength[iwave] = 1071.049452;
257  // fwhm[iwave] = 20.3179033;
258  // Lt[iwave*npixels + k] = (10.1448143*temp->Lt[271*npixels+k] +
259  // 10.173089 *temp->Lt[272*npixels+k])/fwhm[iwave];
260  // iwave++;
261  //
262  // Lt[iwave*npixels + k] = temp->Lt[273*npixels+k];
263  // fwhm[iwave] = temp->fwhm[273];
264  // ociaWavelength[iwave] = temp->wave[273];
265  // iwave++;
266  //
267  // ociaWavelength[iwave] = 1085.272670;
268  // fwhm[iwave] = 53.6097915;
269  // Lt[iwave*npixels + k] = (10.7227049*temp->Lt[274*npixels+k] +
270  // 10.7221941*temp->Lt[275*npixels+k] +
271  // 10.7218208*temp->Lt[276*npixels+k] +
272  // 10.721585 *temp->Lt[277*npixels+k] +
273  // 10.7214867*temp->Lt[278*npixels+k])/fwhm[iwave];
274  // iwave++;
275  // //iwave = interpLt(133, npixels, k, nbands_av, orcaWavelength[iwave], iwave, temp, acc,Lt);
276  //
277 
278  } // pixel loop
279 
280  if (NBANDS != iwave) {
281  printf("Oops, something's wrong iwave (%d) != nbands (%d) \n", (int) iwave, NBANDS);
282  exit(-1);
283  }
284  // Write Lt's
285  for (size_t iwave = 0; iwave < nbands_ocia; iwave++) {
286  int varID;
287 
288  // printf("Lambda(%d) = %d\n",iwave+1,(int)(orcaWavelength[iwave]+.5));
289  status = nc_inq_varid(grpID, "Lt", &varID);
290  check_err(status, __LINE__, __FILE__);
291 
292  size_t start[3] = {iwave, iscan - sline, 0};
293  size_t count[3] = {1, 1, npixels};
294  status = nc_put_vara_float(grpID, varID, start, count,
295  &Lt[iwave * npixels]);
296  check_err(status, __LINE__, __FILE__);
297  }
298 
299  utc[iscan] = temp->scantime;
300  //if (utc[iscan]>0) printf("utc=%f\n",utc[iscan]);
301 
302  grpID = ociafile.getGid("navigation_data");
303 
304  startOCIP[0] = iscan - sline;
305 
306  for (i = 0; i < 6; i++) {
307  // cout << "Copying " << navfieldsPRISM[i] << " to " <<
308  // navfieldsOCIP[i] << endl;
309 
310  status = nc_inq_varid(grpID, navfieldsOCIP[i], &varID);
311  check_err(status, __LINE__, __FILE__);
312 
313  switch (i) {
314  case 0:
315  status = nc_put_vara_float(grpID, varID, startOCIP, countOCIP, &temp->sena[0]);
316  check_err(status, __LINE__, __FILE__);
317  break;
318  case 1:
319  status = nc_put_vara_float(grpID, varID, startOCIP, countOCIP, &temp->senz[0]);
320  check_err(status, __LINE__, __FILE__);
321  break;
322  case 2:
323  status = nc_put_vara_float(grpID, varID, startOCIP, countOCIP, &temp->sola[0]);
324  check_err(status, __LINE__, __FILE__);
325  break;
326  case 3:
327  status = nc_put_vara_float(grpID, varID, startOCIP, countOCIP, &temp->solz[0]);
328  check_err(status, __LINE__, __FILE__);
329  break;
330  case 4:
331  status = nc_put_vara_double(grpID, varID, startOCIP, countOCIP, &temp->lat[0]);
332  check_err(status, __LINE__, __FILE__);
333  break;
334  case 5:
335  status = nc_put_vara_double(grpID, varID, startOCIP, countOCIP, &temp->lon[0]);
336  check_err(status, __LINE__, __FILE__);
337  break;
338  default:
339  break;
340  }
341  }
342  if (temp->scantime > 0) {
343  scantime_last = temp->scantime;
344  if (!first_good_scantime) {
345  first_good_scantime = 1;
346  scantime_first = temp->scantime;
347  }
348  }
349 
350  } // scan loop
351 
352  if (scantime_first < scantime_last) {
353  strncpy(isodatetime_first, unix2isodate(scantime_first, 'G'), 30);
354  strncpy(isodatetime_last, unix2isodate(scantime_last, 'G'), 30);
355  } else {
356  strncpy(isodatetime_last, unix2isodate(scantime_first, 'G'), 30);
357  strncpy(isodatetime_first, unix2isodate(scantime_last, 'G'), 30);
358  }
359  status = nc_put_att_text(ncid, NC_GLOBAL, "time_coverage_start", strlen(isodatetime_first) + 1, isodatetime_first);
360  if (status != NC_NOERR) {
361  printf("-E- %s %d: %s for %s\n",
362  __FILE__, __LINE__, nc_strerror(status), "time_coverage_start");
363  exit(1);
364  }
365  status = nc_put_att_text(ncid, NC_GLOBAL, "time_coverage_end", strlen(isodatetime_last) + 1, isodatetime_last);
366  if (status != NC_NOERR) {
367  printf("-E- %s %d: %s for %s\n",
368  __FILE__, __LINE__, nc_strerror(status), "time_coverage_end");
369  exit(1);
370  }
371 
372  gsl_interp_accel_free(acc);
373 
374  size_t startw[2] = {0, 0};
375  size_t countw[2] = {nbands_ocia, 0};
376 
377  grpID = ociafile.getGid("sensor_band_parameters");
378 
379  status = nc_inq_varid(grpID, "wavelength", &varID);
380 
381  status = nc_put_vara_float(grpID, varID, startw, countw, ociaWavelength);
382  check_err(status, __LINE__, __FILE__);
383 
384  status = nc_inq_varid(grpID, "fwhm", &varID);
385 
386  status = nc_put_vara_float(grpID, varID, startw, countw, fwhm);
387  check_err(status, __LINE__, __FILE__);
388 
389  size_t starts[2] = {0, 0};
390  size_t counts[2] = {nscans, 0};
391 
392  const char *nscanfieldsPRISM[] = {"utc", "alt"};
393  const char *nscanfieldsOCIP[] = {"scan_start_time", "altitude"};
394  short* alt = (short*) (calloc(nscans, sizeof (short)));
395 
396  for (i = 0; i < nscans; i++) {
397  //utc[i] = temp->scantime;
398  alt[i] = temp->alt;
399  }
400  grpID = ociafile.getGid("scan_line_attributes");
401  countOCIP[1] = 0;
402  for (i = 0; i < 2; i++) {
403  cout << "Copying " << nscanfieldsPRISM[i] << " to " <<
404  nscanfieldsOCIP[i] << endl;
405 
406  status = nc_inq_varid(grpID, nscanfieldsOCIP[i], &varID);
407  check_err(status, __LINE__, __FILE__);
408 
409  switch (i) {
410  case 0:
411  status = nc_put_vara_double(grpID, varID, starts, counts, utc);
412  check_err(status, __LINE__, __FILE__);
413  break;
414  case 1:
415  status = nc_put_vara_short(grpID, varID, starts, counts, alt);
416  check_err(status, __LINE__, __FILE__);
417  break;
418  default:
419  break;
420  }
421  }
422 
423  // outfile.close();
424 
425  status = close_prism(temp);
426 
427  return 0;
428 }
#define NBANDS
int j
Definition: decode_rs.h:73
int status
Definition: l1_czcs_hdf.c:32
int getNcid()
Definition: cdl_utils.cpp:746
void check_err(const int stat, const int line, const char *file)
Definition: nc4utils.c:35
void handle_error(int status)
Definition: nccmp.c:219
#define NULL
Definition: decode_rs.h:63
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
uint8 * counts
Definition: l1_czcs_hdf.c:30
#define VERSION
int close_prism(prism4ocia_t *data)
Definition: read_prism.c:509
float fwhm[NBANDS]
Definition: atrem_corl1.h:144
int main(int argc, char *argv[])
size_t interpLt(size_t startJ, size_t npixels, size_t k, size_t nbMax, float ociaWavelength, size_t iwave, prism4ocia_t *data, gsl_interp_accel *acc, float *Lt)
int getGid(const char *grpName)
Definition: cdl_utils.cpp:736
prism4ocia_t * open_prism(char *filename, prism4ocia_t **data)
Definition: read_prism.c:74
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
int read_prism(prism4ocia_t *data, int32_t recnum)
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
SENSOR
Definition: DDProcess.h:81
const char * str
Definition: l1c_msi.cpp:35
int32_t iscan
char * unix2isodate(double dtime, char zone)
Definition: unix2isodate.c:10
PARAM_TYPE_NONE Default value No parameter is buried in the product name name_prefix is case insensitive string compared to the product name PARAM_TYPE_VIS_WAVE The visible wavelength bands from the sensor are buried in the product name The product name is compared by appending and name_suffix ie aph_412_giop where prod_ix will be set to PARAM_TYPE_IR_WAVE same search method as PARAM_TYPE_VIS_WAVE except only wavelength above are looped through but prod_ix is still based ie aph_2_giop for the second and prod_ix set to PARAM_TYPE_INT name_prefix is compared with the beginning of the product name If name_suffix is not empty the it must match the end of the product name The characters right after the prefix are read as an integer and prod_ix is set to that number strncpy(l2prod->name_prefix, "myprod", UNITLEN)
int cdlCreateDim(char *l1_filename, char *cdl_filename, const char **dim_names, size_t *dim_size, size_t n_dims, size_t numScans)
Definition: cdl_utils.cpp:294
#define CDLFILE
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 k
Definition: decode_rs.h:73
int count
Definition: decode_rs.h:79