OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
make_L3_v1.1.c
Go to the documentation of this file.
1 /*
2  File: make_L3.c
3 
4  Compiler options:
5  HP:
6  cc -Ae +O2 +Ofastaccess +Onolimit make_L3.c -o make_L3 ./LIB/lib.a -lf -lm -I /usr/local/HDF4.0/include -L/usr/local/HDF4.0/lib -lmfhdf -ldf -ljpeg -lz
7  SGI:
8  cc -O2 -fullwarn -woff 1209,1210 -64 make_L3.c -o make_L3 ./LIB/lib.a -I./LIB -I$HDFINC -L$HDFLIB -lmfhdf -ldf -lz -ljpeg -lm -lmalloc
9  or cc -O2 -fullwarn -woff 1209,1210 -64 make_L3.c -o make_L3 ./LIB/lib.a -I./LIB -I$HDFINC -L/home/jack/libhdf64bit.4.1r3/ -lmfhdf -ldf -lz -lm -ljpeg
10 
11  Revision History:
12  Version 1.0 May 30, 2000, J. Descloitres
13  Version 1.0.1 May 31, 2000, J. Descloitres
14  Version 1.0.2 June 14, 2000, J. Descloitres
15  Version 1.1 August 8, 2001, J. Gales
16  Version 1.11 November 19, 2001, J. Gales
17 
18  References and Credits:
19  Jacques Descloitres
20  University of Maryland - Department of Geography
21  NASA / GSFC Code 922
22  Bldg 32 Room S31
23  Greenbelt, MD 20771
24  Phone: 301 614 5456
25  Email: jack@ltpmail.gsfc.nasa.gov
26  */
27 
28 #include <stdint.h>
29 #include <stdio.h>
30 #include <math.h>
31 #include <string.h>
32 #include <fcntl.h>
33 #include <unistd.h>
34 #include <sys/stat.h>
35 #include "path.h"
36 #include "proto.h"
37 #include "proj_proto.h"
38 #include "mfhdf.h"
39 
40 enum {
42 };
43 
44 enum {
46 };
47 
48 enum {
50 };
51 
52 enum {
54 };
55 
56 enum {
58 };
59 
60 #define DBL_MAX 1.7976931348623157E+308 /* max decimal value of a "double"*/
61 
62 #define DEFPROJECTION HAMMER /* Default output projection */
63 #define FAST FALSE
64 #define ATMCOR TRUE /* Do atmospheric correction */
65 #define UPDATE TRUE
66 #define TILTSCANS FALSE
67 #define MAXMEMDEF 500 /* Maximum memory to allocate (in Mbytes) */
68 #define STEP 1
69 #define DEFPIXSZ 10000. /* Pixel size (in m) for output image files (if compatible) */
70 #define UNDEF -999. /* Undefined value */
71 #define FILL_UINT8 255 /* Fill value */
72 #define FILL_INT16 -32767 /* Fill value */
73 #define MINREFL 0.00 /* Minimum reflectance allowed for min. blue compositing */
74 #define MAXSENZDEF 55. /* Maximum view angle */
75 #define UMASK 002 /* Set permission for created files */
76 #define BLUEBAND BAND2 /* Blue band for minimum blue compositing */
77 #define WATER_VAL -16384 /* Water (non-land) pixel value */
78 
79 #define WAVELENGTH { "412", "443", "490", "510", "555", "670", "765", "865"} /* Band wavelength */
80 #define IRRAD {171.07, 189.91, 194.33, 188.24, 185.93, 152.14, 123.55, 100.00} /* TOA solar irradiance */
81 #define TAURAY {0.3139, 0.2341, 0.1561, 0.1324, 0.0942, 0.0439, 0.0257, 0.0155} /* TOA molecular reflectance */
82 #define OZONE {.00103, .00400, .02536, .04200, .09338, .04685, .00837, .00485} /* Coef. for ozone absorption */
83 #define WVAPOR {0.0000, 0.0000, 0.0000, 0.0000, 0.0008, 0.0043, 0.0007, 0.0057} /* Coef. for water vapor abs. */
84 #define OXYGEN {0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000} /* Coef. for oxygen absorption */
85 #define UOZ 0.270 /* Ozone abundance (cm.atm) */
86 #define PSURF 1013. /* Sea surface pressure */
87 #define NLPSURF 2160 /* Number of rows in surface pressure file */
88 #define NPPSURF 4320 /* Number of columns in surface pressure file */
89 
90 #define TWO_THIRDS 0.666666666666666666667
91 #define FOUR_THIRDS 1.333333333333333333333
92 #define TWO_PI 6.28318530717958647692
93 #define RAD2DEG 57.29577951308232087684 /* 180.0 / M_PI */
94 #define NpixelsLAC 1285
95 #define OMF2 0.99330562 /* (1.0 - F) * (1.0 - F) */
96 #define R 6371007.181 /* Earth radius (m) */
97 #define MAXLENGTH 100
98 #define EPSILON 1e-6
99 
100 #define N_CONVOL 7
101 #define CONVOL_MAX 30.0
102 
103 #define CROSS_PRODUCT(x1, y1, x2, y2, x3, y3) ( ((x2) - (x1)) * ((y3) - (y1)) - ((y2) - (y1)) * ((x3) - (x1)) )
104 
105 #define VERSION "1.1.1"
106 
107 typedef struct {
108  int16 *reflTOA[NBANDS], *refl[NBANDS], *solz, *senz, *phi, *smoke, *ndvi, *evi, *gemi, *ifno;
109  char projection, projname[MAXLENGTH];
111  char *filelist, *cmd_line;
112  int Nl, Np, w, h, buflines, bufstep;
113  int l1, l2, p1, p2;
114  double lon1, lon2, lat1, lat2;
115  double xmin, xmax, ymin, ymax;
116  char domain;
117  double pixsz, lon_center;
118  int32 sd_id;
119  int startline, startday, Nfiles;
120  FILE* input_fp;
121  FILE* fileuse_fp;
122 } OUTPUT;
123 
124 typedef struct {
125  char ndvi, evi, gemi, angles, smoke, gzip, atmcor, TOA;
126  char band[NBANDS];
127  float maxsenz;
128  int maxmem;
129  char composite, compositename[MAXLENGTH], path[MAXLENGTH];
130  char daystr[14];
131  char dayend[14];
132 } PROCESSINFO;
133 
134 typedef struct {
135  float *lon, *lat, *mus, *muv, *tilt;
136  int16 *ndvi, *evi, *senz, *solz, *phi;
137  int Nscans, Npixels, startpix, subsampling;
138  float *rho[NBANDS], *rhoTOA[NBANDS];
139  int doy, doy_prev, year, orbno, step, old_iscan, corruption;
140  int32 sd_id;
141  double *lon1, *lat1, *lon2, *lat2;
142  double *lontop, *lonbot, *lattop, *latbot;
143  double *xtop, *xbot, *ytop, *ybot;
144  char newfile, intersect, title[100], level;
145  int32 *l2_flags;
146 } INPUT;
147 
148 int alloc_new_array(void **data, int size, int32 num_type, int maxmem);
149 void check_process_parameters(PROCESSINFO *process, char **wl);
150 int get_geolocation_scan(int iscan, INPUT *input, PROCESSINFO process);
151 int get_lonlat(double xproj, double yproj, double *lon, double *lat, int projtype, double lon_center);
152 int get_xyproj(double lon, double lat, double *xproj, double *yproj, int projtype, double lon_center);
153 void init_output_file(OUTPUT *proj, INPUT input, PROCESSINFO process);
154 void init_output_projection(OUTPUT *proj, PROCESSINFO process);
155 void memcpy_block(int offset2, int offset1, int Nrows, int Ncols, OUTPUT *proj);
156 void parse_command_line(int argc, char **argv, OUTPUT *proj, PROCESSINFO *process);
157 void process_pixel(int ipix, int iline_out, int ipix_out, char **wl, INPUT input, OUTPUT *proj, PROCESSINFO process);
158 int project_scan(int iscan, char **wl, INPUT *input, OUTPUT *proj, PROCESSINFO process);
159 int read_file_block(char *filename, char **wl, int offset, int startline, int buflines, OUTPUT *proj);
160 int read_l2(int iscan, INPUT *input, char **wl);
161 void resample_geo_scan(int Npixels, float *inlon, float *inlat, double *outlon, double *outlat);
162 void set_buffer_parameters(OUTPUT *proj, PROCESSINFO *process);
163 void set_proj_parameters(OUTPUT *proj);
164 int update_file(char *filename, char **wl, int offset, int startline, int buflines, OUTPUT *proj, char verbose, char gzip);
165 void update_pixel(INPUT input, int16 ndvi, int ipix, int idx_out, int ifno, OUTPUT *proj);
166 void write_cmd_line(int32 sd_id, char *cmd_line);
167 int16 EVI(float blue, float red, float nir);
168 int16 NDVI(float red, float nir);
169 int16 GEMI(float red, float nir);
170 int16 SMOKE(float **rho, int ipix);
171 unsigned char read_mask(double lon, double lat, int16 close);
172 
173 float64 sf_refl = 0.0001, sf_angle = 0.01, sf_smoke = 0.0001, sf_vi = 0.0001;
174 
175 int main(int argc, char **argv) {
176 
177  char prefix[MAXLENGTH];
178  char filename[MAXLENGTH], inputfile[MAXLENGTH];
179  char string[100];
180  int badorb;
181  int ipix, iscan, j;
182  char *ptr;
183 
184  OUTPUT proj;
185  INPUT input;
186  PROCESSINFO process;
187 
188  int16 *psurf = NULL;
189  int iline_psurf, ipix_psurf, idx_psurf;
190  float psurf0 = PSURF;
191 
192  char *wl[NBANDS] = WAVELENGTH;
193  float taur0[NBANDS] = TAURAY;
194  float aoz[NBANDS] = OZONE;
195  float awv[NBANDS] = WVAPOR;
196  float ao2[NBANDS] = OXYGEN;
197  float Es[NBANDS] = IRRAD;
198  float taur[NBANDS], rhoray[NBANDS];
199  double trup[NBANDS], trdown[NBANDS];
200  double coefcalib[NBANDS];
201  float m, Ttotrayu, Ttotrayd, rho_pix, dsol, uoz = 0.270;
202  int ib;
203  float sphalb[5000];
204 
205  int32 *scan_year;
206  int32 *scan_day;
207  int32 *scan_msec;
208  int32 start[2];
209  int32 edges[2];
210  int32 sds_id;
211  int hr, min, sec;
212  char scantime_str[14];
213  int i;
214  int k;
215  int l;
216  float32 convol;
217  float32 convol_diff;
218  unsigned char noisy_scan;
219  float32 *center_lat;
220 
221  char* tmpstr;
222 
223  FILE *fd;
224 
226 
227 
228  printf("%s %s (%s %s)\n", "LANDBIN", VERSION, __DATE__, __TIME__);
229 
230  sphalb[0] = 0;
231  for (j = 1; j < 5000; j++) sphalb[j] = csalbr(j / 10000.);
232 
233  input.lon = input.lat = input.mus = input.muv = NULL;
234  input.senz = input.solz = input.phi = NULL;
235  input.lon1 = input.lat1 = input.lon2 = input.lat2 = NULL;
236  input.lontop = input.lonbot = input.lattop = input.latbot = NULL;
237  input.xtop = input.xbot = input.ytop = input.ybot = NULL;
238  input.l2_flags = NULL;
239 
240  for (ib = 0; ib < NBANDS; ib++)
241  input.rho[ib] = NULL;
242  for (ib = 0; ib < NBANDS; ib++)
243  input.rhoTOA[ib] = NULL;
244  input.ndvi = NULL;
245  input.evi = NULL;
246  input.step = STEP;
247 
248  proj.pixsz = DEFPIXSZ;
249  proj.lon_center = 0;
250  proj.domain = GLOBAL;
251  proj.startday = 0;
252  proj.sd_id = -1;
253  proj.bufstep = -1;
254  proj.lon1 = proj.lon2 = proj.lat1 = proj.lat2 = UNDEF;
255  proj.smoke = proj.ndvi = proj.evi = proj.gemi = proj.solz = proj.senz = proj.phi = NULL;
256  proj.filelist = NULL;
257  proj.Nfiles = -1;
258  proj.startline = 0;
259  proj.filename[0] = 0;
260  proj.input_fp = NULL;
261  proj.fileuse_fp = NULL;
262 
263  for (ib = 0; ib < NBANDS; ib++)
264  process.band[ib] = FALSE;
265  process.angles = FALSE;
266  process.ndvi = FALSE;
267  process.evi = FALSE;
268  process.gemi = FALSE;
269  process.maxmem = MAXMEMDEF;
270  process.maxsenz = MAXSENZDEF;
271  process.composite = MINSENZ;
272  process.smoke = FALSE;
273  process.gzip = FALSE;
274  process.TOA = FALSE;
275  process.atmcor = ATMCOR;
276 
277 
278  strcpy(process.daystr, "0000000000000");
279  strcpy(process.dayend, "9999999999999");
280 
281  if (argc == 1) exit(0);
282 
283  parse_command_line(argc, argv, &proj, &process);
284 
285  check_process_parameters(&process, wl);
286 
287  set_proj_parameters(&proj);
288 
289  printf("Initialization...(120901) \n");
290 
291  set_buffer_parameters(&proj, &process);
292 
293  init_output_projection(&proj, process);
294 
295  umask(UMASK);
296 
297 
298  while (1) {
299 
300 
301  if (proj.input_fp != NULL) {
302  tmpstr = fgets(inputfile, MAXLENGTH, proj.input_fp);
303 
304  if (tmpstr == NULL) goto bail;
305 
306  inputfile[strlen(inputfile) - 1] = 0;
307 
308  } else {
309 
310  i = scanf("%s", inputfile);
311 
312  if (i == EOF) goto bail;
313  }
314 
315 
316  /* open input L1 HDF file */
317  printf("Opening %s ...\n", inputfile);
318  if ((input.sd_id = SDstart(inputfile, DFACC_RDONLY)) == -1) {
319  fprintf(stderr, "Cannot open file %s\n", inputfile);
320  continue;
321  }
322  input.newfile = TRUE;
323 
324  /* get global attributes */
325  if (get_attributes(input.sd_id, prefix, &input.Nscans, &input.Npixels, &input.startpix, &input.subsampling, &input.year, &input.doy, &input.orbno, input.title) == -1) {
326  SDend(input.sd_id);
327  continue;
328  }
329 
330  if (strcmp(input.title, "SeaWiFS Level-1A Data") == 0) {
331  input.level = L1A;
332  if (ATMCOR && !process.atmcor) {
333  printf("This file is Level-1A data. Atmospheric correction reactivated.\n");
334  process.atmcor = TRUE;
335  }
336  } else if (strcmp(input.title, "SeaWiFS Level-2 Data") == 0) {
337  input.level = L2;
338  if (process.atmcor) {
339  printf("This file is Level 2 data. Atmospheric correction inhibited.\n");
340  process.atmcor = FALSE;
341  }
342  } else {
343  fprintf(stderr, "Product type not supported (%s)\n", input.title);
344  SDend(input.sd_id);
345  continue;
346  }
347 
348  if (process.atmcor) {
349  printf("Reading terrain elevation...\n");
350  alloc_new_array((void *) &psurf, NLPSURF * NPPSURF, DFNT_INT16, process.maxmem);
351  sprintf(process.path, "%s", ANC_DIR);
352  sprintf(filename, "%s/tbase.img", process.path);
353  make_psurf(psurf, filename, NLPSURF, NPPSURF, psurf0);
354  }
355 
356  /* Truncate granule name */
357  *(prefix + 14) = 0;
358 
359  if (proj.sd_id < 0) {
360  init_output_file(&proj, input, process);
361  input.doy_prev = input.doy;
362  }
363 
364  printf("\n");
365  printf("Year = %d\n", input.year);
366  printf("Day of Year = %d\n", input.doy);
367  printf("Orbit Number = %d\n", input.orbno);
368  printf("Number of scans = %d\n", input.Nscans);
369  printf("Number of pixels = %d\n", input.Npixels);
370 
371  if ((fd = fopen(BADORBLIST, "r")) == NULL)
372  fprintf(stderr, "Cannot read file %s\n", BADORBLIST);
373  else {
374  while ((ptr = fgets(string, 100, fd)) != NULL)
375  if (sscanf(string, "%d", &badorb) == 1 &&
376  input.orbno == badorb)
377  break;
378  fclose(fd);
379  if (ptr) {
380  printf("Orbit %d is listed in the bad orbit list. Skipped.\n", input.orbno);
381  SDend(input.sd_id);
382  continue;
383  }
384  }
385 
386  if (input.doy != proj.startday)
387  fprintf(stderr, "\nWARNING !!!!!: Compositing different days.\n\n");
388 
389  if (UPDATE &&
390  input.doy != input.doy_prev) {
391 
392  if (update_file(proj.filename, wl, 0, proj.startline, proj.buflines, &proj, TRUE, process.gzip) == -1)
393  exit(-1);
394  printf("\n");
395 
396  } /* if (input.doy != input.doy_prev) */
397 
398  input.doy_prev = input.doy;
399 
400  /* Precompute constants */
401  dsol = .9856 * (input.doy - 4) * M_PI / 180.;
402  dsol = 1. / pow(1. - .01673 * cos(dsol), 2);
403  for (ib = 0; ib < NBANDS; ib++)
404  coefcalib[ib] = M_PI / Es[ib] / dsol;
405 
406  for (ib = 0; ib < NBANDS; ib++)
407  if (process.band[ib])
408  input.rho[ib] = (float *) realloc(input.rho[ib], input.Npixels * sizeof (float));
409  if (input.level == L2) {
410  for (ib = 0; ib < NBANDS; ib++)
411  if (process.band[ib] && process.TOA)
412  input.rhoTOA[ib] = (float *) realloc(input.rhoTOA[ib], input.Npixels * sizeof (float));
413  if (process.ndvi)
414  input.ndvi = (int16 *) realloc(input.ndvi, input.Npixels * sizeof (int16));
415  if (process.evi)
416  input.evi = (int16 *) realloc(input.evi, input.Npixels * sizeof (int16));
417  }
418 
419  input.lon = (float *) realloc(input.lon, input.Npixels * sizeof (float));
420  input.lat = (float *) realloc(input.lat, input.Npixels * sizeof (float));
421  input.mus = (float *) realloc(input.mus, input.Npixels * sizeof (float));
422  input.muv = (float *) realloc(input.muv, input.Npixels * sizeof (float));
423  input.senz = (int16 *) realloc(input.senz, input.Npixels * sizeof (int16));
424  input.solz = (int16 *) realloc(input.solz, input.Npixels * sizeof (int16));
425  input.phi = (int16 *) realloc(input.phi, input.Npixels * sizeof (int16));
426  input.tilt = (float *) malloc(input.Nscans * sizeof (float));
427 
428  input.lon1 = (double *) realloc(input.lon1, (input.Npixels + 1) * sizeof (double));
429  input.lat1 = (double *) realloc(input.lat1, (input.Npixels + 1) * sizeof (double));
430  input.lon2 = (double *) realloc(input.lon2, (input.Npixels + 1) * sizeof (double));
431  input.lat2 = (double *) realloc(input.lat2, (input.Npixels + 1) * sizeof (double));
432  input.lontop = (double *) realloc(input.lontop, (input.Npixels + 1) * sizeof (double));
433  input.lattop = (double *) realloc(input.lattop, (input.Npixels + 1) * sizeof (double));
434  input.lonbot = (double *) realloc(input.lonbot, (input.Npixels + 1) * sizeof (double));
435  input.latbot = (double *) realloc(input.latbot, (input.Npixels + 1) * sizeof (double));
436  input.xtop = (double *) realloc(input.xtop, (input.Npixels + 1) * sizeof (double));
437  input.ytop = (double *) realloc(input.ytop, (input.Npixels + 1) * sizeof (double));
438  input.xbot = (double *) realloc(input.xbot, (input.Npixels + 1) * sizeof (double));
439  input.ybot = (double *) realloc(input.ybot, (input.Npixels + 1) * sizeof (double));
440 
441  input.l2_flags = (int32 *) realloc(input.l2_flags, input.Npixels * sizeof (int32));
442 
443  if (input.level == L2)
444  printf("Reading file %s ...\n", inputfile);
445  else {
446  printf("Reading SDS \"l1a_data\" in file %s ...\n", inputfile);
447  if (process.atmcor)
448  printf("Calibration -> Geolocation -> Computation of surface reflectance (atmospheric correction) -> Projection ...\n");
449  else
450  printf("Calibration -> Geolocation -> Computation of top-of-the-atmosphere reflectance (no atmospheric correction) -> Projection ...\n");
451  }
452 
453 
454  input.intersect = FALSE;
455  input.old_iscan = -10;
456  if (read_sds(input.sd_id, "tilt", input.tilt) == -1) continue;
457 
458  if (proj.filelist)
459  proj.filelist = (char *) realloc(proj.filelist, (strlen(proj.filelist) + strlen(inputfile) + 8) * sizeof (char));
460  else {
461  proj.filelist = (char *) malloc((strlen(inputfile) + 6) * sizeof (char)); /* +1 for initial "\n", +3 for "0: ", */
462  strcpy(proj.filelist, "\n"); /* +1 for "\n", +1 for NULL */
463  }
464  proj.Nfiles++;
465  sprintf(proj.filelist + strlen(proj.filelist), "%d: %s\n", proj.Nfiles, inputfile);
466 
467 
468  /* Read scan year, day, and, msec */
469  /* ------------------------------ */
470  scan_year = (int32 *) calloc(input.Nscans, sizeof (int32));
471  scan_day = (int32 *) calloc(input.Nscans, sizeof (int32));
472  scan_msec = (int32 *) calloc(input.Nscans, sizeof (int32));
473  center_lat = (float32 *) calloc(input.Nscans, sizeof (float32));
474 
475  start[0] = 0;
476  edges[0] = input.Nscans;
477 
478  sds_id = SDselect(input.sd_id, SDnametoindex(input.sd_id, "year"));
479  if (SDreaddata(sds_id, start, NULL, edges, scan_year) < 0) {
480  fprintf(stderr, "can't read SDS %s\n", "year");
481  return -1;
482  }
483  SDendaccess(sds_id);
484 
485  sds_id = SDselect(input.sd_id, SDnametoindex(input.sd_id, "day"));
486  if (SDreaddata(sds_id, start, NULL, edges, scan_day) < 0) {
487  fprintf(stderr, "can't read SDS %s\n", "day");
488  return -1;
489  }
490  SDendaccess(sds_id);
491 
492  sds_id = SDselect(input.sd_id, SDnametoindex(input.sd_id, "msec"));
493  if (SDreaddata(sds_id, start, NULL, edges, scan_msec) < 0) {
494  fprintf(stderr, "can't read SDS %s\n", "msec");
495  return -1;
496  }
497  SDendaccess(sds_id);
498 
499 
500  start[1] = input.Npixels / 2;
501  edges[1] = 1;
502 
503  sds_id = SDselect(input.sd_id, SDnametoindex(input.sd_id, "latitude"));
504  if (SDreaddata(sds_id, start, NULL, edges, center_lat) < 0) {
505  fprintf(stderr, "can't read SDS %s\n", "latitude");
506  return -1;
507  }
508  SDendaccess(sds_id);
509 
510 
511  for (iscan = 0; iscan < input.Nscans; iscan += input.step) {
512 
513  if (iscan % 500 == 0) printf("Processing scan no. %d...\n", iscan);
514 
515  noisy_scan = 0;
516 
517  /* Checking that scan is after daystr and before dayend */
518  /* ---------------------------------------------------- */
519  hr = (int) (scan_msec[iscan] / 1000.0 / 3600.0);
520  min = (int) (((scan_msec[iscan] / 1000.0 / 3600.0) - hr) * 60.0);
521  sec = (int) (((((scan_msec[iscan] / 1000.0 / 3600.0) - hr) * 60.0) - min) * 60.0);
522 
523  sprintf(scantime_str, "%4d%3d%2d%2d%2d", scan_year[iscan], scan_day[iscan],
524  hr, min, sec);
525  for (i = 0; i < 13; i++) if (scantime_str[i] == ' ') scantime_str[i] = '0';
526 
527  if (strcmp(scantime_str, process.daystr) < 0) {
528  /* printf("%s %s\n", process.daystr, scantime_str);*/
529  continue;
530  }
531 
532  if (strcmp(scantime_str, process.dayend) > 0) {
533  /* printf("%s %s\n", process.dayend, scantime_str);*/
534  break;
535  }
536 
537 
538  if (!TILTSCANS &&
539  (fabs(input.tilt[iscan] - input.tilt[iscan - 1 + 2 * (iscan == 0)]) > 0.1 ||
540  fabs(input.tilt[iscan] - input.tilt[iscan + 1 - 2 * (iscan == input.Nscans - 1)]) > 0.1))
541  continue;
542 
543 
544  /* Skip "latitude" jumps */
545  /* --------------------- */
546  if ((fabs(center_lat[iscan] - center_lat[iscan - 1 + 2 * (iscan == 0)]) > 0.1 ||
547  fabs(center_lat[iscan] - center_lat[iscan + 1 - 2 * (iscan == input.Nscans - 1)]) > 0.1))
548  continue;
549 
550 
551  /* Read and calibrate data */
552  if ((input.level == L1A && (input.corruption = compute_l1b(input.sd_id, iscan, input.rho, CALIBTABLE, input.newfile)) == -1) ||
553  (input.level == L2 && read_l2(iscan, &input, wl) == -1))
554  continue;
555 
556  /* Get navigation info */
557  if (input.old_iscan != iscan - input.step) {
558  if (get_geolocation_scan(iscan, &input, process) == -1)
559  continue;
560  resample_geo_scan(input.Npixels, input.lon, input.lat, input.lon2, input.lat2);
561  }
562 
563  /* Assumption: bounding box can be directly derived from 1st and last pixels */
564  /* To be checked */
565  /* Must be refined anyway to not discard scans that intersect
566  the region although no pixel center falls in the region */
567  if (proj.lon1 != UNDEF &&
568  (input.lat[input.Npixels - 1] > proj.lat1 ||
569  input.lat[0] < proj.lat2 ||
570  (input.lon[0] < input.lon[input.Npixels - 1] &&
571  (input.lon[0] > proj.lon2 ||
572  input.lon[input.Npixels - 1] < proj.lon1))))
573  continue;
574 
575 
576  /* L1A processing (if needed) */
577  /* -------------------------- */
578  if (input.level == L1A &&
579  input.corruption != 1) {
580 
581  /* Computation of reflectance & Atmospheric correction */
582 
583  for (ipix = 0; ipix < input.Npixels; ipix += input.step) {
584  if (input.solz[ipix] == FILL_INT16) continue;
585  if (process.atmcor) {
586  iline_psurf = (int) ((NLPSURF - 1) / 2. - input.lat[ipix] * (NLPSURF - 1) / M_PI + 0.5);
587  ipix_psurf = (int) ((NPPSURF - 1) / 2. + input.lon[ipix] * (NPPSURF - 1) / TWO_PI + 0.5);
588  idx_psurf = iline_psurf * NPPSURF + ipix_psurf;
589  for (ib = 0; ib < NBANDS; ib++)
590  if (process.band[ib]) taur[ib] = taur0[ib] * psurf[idx_psurf] / psurf0;
591  chand(input.phi[ipix] * sf_angle, input.muv[ipix], input.mus[ipix],
592  taur, rhoray, trup, trdown, NBANDS, process.band);
593  m = 1 / input.mus[ipix] + 1 / input.muv[ipix];
594  }
595  for (ib = 0; ib < NBANDS; ib++) {
596  if (process.band[ib]) {
597  if ((rho_pix = input.rho[ib][ipix]) > 0) {
598  rho_pix *= coefcalib[ib] / input.mus[ipix];
599  if (process.atmcor) {
600  Ttotrayu = ((TWO_THIRDS + input.muv[ipix]) + (TWO_THIRDS - input.muv[ipix]) *
601  trup[ib]) / (FOUR_THIRDS + taur[ib]);
602  Ttotrayd = ((TWO_THIRDS + input.mus[ipix]) + (TWO_THIRDS - input.mus[ipix]) *
603  trdown[ib]) / (FOUR_THIRDS + taur[ib]);
604  rho_pix = (rho_pix / expf(-m * uoz * aoz[ib]) - rhoray[ib]) / Ttotrayd / Ttotrayu;
605  if (ao2[ib] > 0) rho_pix /= expf(-m * ao2[ib]);
606  if (awv[ib] > 0) rho_pix /= expf(-m * awv[ib]);
607  rho_pix /= (1 + sphalb[(int) (taur[ib] * 10000 + 0.5)] * rho_pix);
608  }
609  input.rho[ib][ipix] = rho_pix;
610  } else
611  input.rho[ib][ipix] = UNDEF;
612  }
613  }
614  } /* end atmospheric correction for (ipix=0 ... */
615 
616  } /* if L1A ... */
617 
618 
619  /* Find and flag "noisy" scans */
620  /* --------------------------- */
621  for (ib = 0; ib < NBANDS; ib++) {
622  if (input.rho[ib] != 0x0) {
623  convol_diff = 0.0;
624  for (ipix = 0; ipix < input.Npixels; ipix += input.step) {
625  convol = 0.0;
626  for (k = -N_CONVOL / 2; k <= N_CONVOL / 2; k++) {
627  l = ipix + k;
628  if (l < 0) l = 0;
629  if (l >= input.Npixels) l = input.Npixels - 1;
630 
631  convol += input.rho[ib][l];
632  }
633  convol_diff += fabs((convol / N_CONVOL) - input.rho[ib][ipix]);
634  } /* ipix loop */
635 
636  if (convol_diff >= CONVOL_MAX) {
637  noisy_scan = 1;
638  break;
639  }
640 
641  } /* if (input.rho[ib] ... */
642  } /* ib loop */
643 
644  if (noisy_scan) continue;
645 
646 
647  /* Project scan and update output projection */
648  if (project_scan(iscan, wl, &input, &proj, process) == -1)
649  continue;
650 
651  input.old_iscan = iscan;
652  input.newfile = FALSE;
653 
654  } /* for (iscan=0 ... */
655 
656 
657  /* Close HDF input file */
658  SDend(input.sd_id);
659 
660 
661  if (!input.intersect) {
662  printf("Note: file %s did not contribute at all to the output projection.\n", inputfile);
663  } else {
664  printf("Note: file %s actually contributed to the output projection.\n", inputfile);
665  if (proj.fileuse_fp != NULL) {
666  fprintf(proj.fileuse_fp, "%s\n", inputfile);
667  }
668  }
669 
670  printf("\n");
671 
672  } /* while (scanf("%s", inputfile) != EOF) */
673 
674 bail:
675 
676  /* Close input filelist */
677  /* -------------------- */
678  if (proj.input_fp != NULL) fclose(proj.input_fp);
679 
680  /* Close fileuse file */
681  /* ------------------ */
682  if (proj.fileuse_fp != NULL) fclose(proj.fileuse_fp);
683 
684 
685  if (update_file(proj.filename, wl, 0, proj.startline, proj.buflines, &proj, TRUE, process.gzip) == -1)
686  exit(-1);
687 
688 
689  free(input.lat);
690  free(input.lon);
691  free(input.senz);
692  free(input.solz);
693  free(input.phi);
694  free(input.mus);
695  free(input.muv);
696  free(input.tilt);
697 
698  free(input.lat1);
699  free(input.lon1);
700  free(input.lat2);
701  free(input.lon2);
702  free(input.lattop);
703  free(input.lontop);
704  free(input.latbot);
705  free(input.lonbot);
706  free(input.xtop);
707  free(input.ytop);
708  free(input.xbot);
709  free(input.ybot);
710  free(input.l2_flags);
711 
712  for (ib = 0; ib < NBANDS; ib++) {
713  if (input.rho[ib]) free(input.rho[ib]);
714  if (input.rhoTOA[ib]) free(input.rhoTOA[ib]);
715  if (proj.refl[ib]) free(proj.refl[ib]);
716  if (proj.reflTOA[ib]) free(proj.reflTOA[ib]);
717  }
718  if (proj.smoke) free(proj.smoke);
719  if (proj.ndvi) free(proj.ndvi);
720  if (proj.evi) free(proj.evi);
721  if (proj.gemi) free(proj.gemi);
722  if (process.atmcor) free(psurf);
723  if (proj.solz) {
724  free(proj.solz);
725  free(proj.senz);
726  free(proj.phi);
727  }
728  free(proj.ifno);
729  /*read_mask((double) 0, (double) 0, (int16) 1); (JMG) */
730 
731  free(scan_year);
732  free(scan_day);
733  free(scan_msec);
734  free(center_lat);
735 
736  exit(0);
737 
738 } /* end */
739 
740 
741 
742 
743 
744 
745 
746 #define IFOV 0.0015835 /* Instantaneous FOV (rad) */
747 
748 int get_geolocation_scan(int iscan, INPUT *input, PROCESSINFO process) {
749  double rmtq[3], geovec[3], up[3], ea[3], no[3];
750  double sv, sn, se, sunv, sunn, sune;
751  double scan, a, b, c, r, q, Qx, Qz, tmp, uxy, temp, upxy;
752  int j, ipix;
753  float32 rm[9], pos[3], sun[3], coef[6];
754  static double *cosa = NULL, *sina = NULL;
755  float sena, senz, sola, solz, phi;
756 
757  if (input->level == L2)
758  return 0;
759 
760  /* Precompute the cosines and sines for all scan angles */
761  if (input->newfile) {
762  cosa = (double *) realloc(cosa, input->Npixels * sizeof (double));
763  sina = (double *) realloc(sina, input->Npixels * sizeof (double));
764  for (ipix = 0; ipix < input->Npixels; ipix++) {
765  scan = (double) ((input->startpix - 1 + input->subsampling * ipix - (NpixelsLAC - 1) / 2.) * IFOV);
766  cosa[ipix] = cos(scan);
767  sina[ipix] = sin(scan);
768  }
769  }
770 
771  /* get SDS related to navigation and time */
772  if (get_navig_sds_line(input->sd_id, iscan, rm, pos, sun, coef, input->newfile) == -1) {
773  fprintf(stderr, "Unable to get geolocation for scan %d\n", iscan);
774  return -1;
775  }
776 
777  for (ipix = 0; ipix < input->Npixels; ipix++) {
778 
779  /* Compute the sensor-to-surface vectors for all scan angles */
780  a = coef[0] * cosa[ipix] * cosa[ipix] + coef[1] * cosa[ipix] * sina[ipix] + coef[2] * sina[ipix] * sina[ipix];
781  b = coef[3] * cosa[ipix] + coef[4] * sina[ipix];
782  c = coef[5];
783 
784  /* Begin the solution of the quadratic equation */
785  r = b * b - 4.0 * a * c;
786 
787  /* Check for scan angles past the edge of the Earth */
788  if (r < 0) {
789 
790  input->lat[ipix] = UNDEF;
791  input->lon[ipix] = UNDEF;
792  input->senz[ipix] = FILL_INT16;
793  input->solz[ipix] = FILL_INT16;
794  input->phi[ipix] = FILL_INT16;
795 
796  } else {
797 
798  /* Solve for the magnitude of the sensor-to-pixel vector and compute its components */
799  q = (-b - sqrt(r)) / (2 * a);
800  Qx = q * cosa[ipix];
801  Qz = q * sina[ipix];
802 
803  /* Transform the vector from the sensor frame to the geocentric frame */
804  for (j = 0; j < 3; j++) {
805  rmtq[j] = Qx * rm[3 * j] + Qz * rm[3 * j + 2];
806  geovec[j] = rmtq[j] + pos[j];
807  }
808  /* Compute the geodetic latitude and longitude for the scan line */
809  tmp = sqrt(geovec[0] * geovec[0] + geovec[1] * geovec[1]) * OMF2;
810  input->lat[ipix] = atan2(geovec[2], tmp);
811  input->lon[ipix] = atan2(geovec[1], geovec[0]);
812 
813  /* Compute the local vertical, East, and North unit vectors */
814  uxy = geovec[0] * geovec[0] + geovec[1] * geovec[1];
815  temp = sqrt(geovec[2] * geovec[2] + OMF2 * OMF2 * uxy);
816  up[0] = OMF2 * geovec[0] / temp;
817  up[1] = OMF2 * geovec[1] / temp;
818  up[2] = geovec[2] / temp;
819  upxy = sqrt(up[0] * up[0] + up[1] * up[1]);
820  ea[0] = -up[1] / upxy;
821  ea[1] = up[0] / upxy;
822  ea[2] = 0;
823  no[0] = -up[2] * ea[1];
824  no[1] = up[2] * ea[0];
825  no[2] = up[0] * ea[1] - up[1] * ea[0];
826 
827  /* Decompose the spacecraft vector into vertical (up), North (no), and East (ea) components */
828  sv = sn = se = 0;
829  for (j = 0; j < 3; j++) {
830  sv += -rmtq[j] * up[j];
831  sn += -rmtq[j] * no[j];
832  se += -rmtq[j] * ea[j];
833  }
834  /* Compute the spacecraft zenith and azimuth angles */
835  senz = RAD2DEG * atan2(sqrt(sn * sn + se * se), sv);
836  input->senz[ipix] = senz / sf_angle + 0.5;
837  if (senz > process.maxsenz) {
838  input->solz[ipix] = FILL_INT16;
839  continue;
840  }
841  sena = RAD2DEG * atan2(se, sn);
842  if (sena < 0) sena += 360;
843  input->muv[ipix] = sv / sqrt(sn * sn + se * se + sv * sv);
844 
845  /* Decompose the sun vector into vertical (up), North (no), and East (ea) components */
846  sunv = sunn = sune = 0;
847  for (j = 0; j < 3; j++) {
848  sunv += sun[j] * up[j];
849  sunn += sun[j] * no[j];
850  sune += sun[j] * ea[j];
851  }
852  /* Compute the solar zenith and azimuth angles */
853  solz = RAD2DEG * atan2(sqrt(sunn * sunn + sune * sune), sunv);
854  if (solz >= 90)
855  input->solz[ipix] = FILL_INT16;
856  else
857  input->solz[ipix] = solz / sf_angle + 0.5;
858  sola = RAD2DEG * atan2(sune, sunn);
859  if (sola < 0) sola += 360;
860  phi = sola - sena;
861  if (phi > +180) phi -= 360;
862  if (phi < -180) phi += 360;
863  input->phi[ipix] = floor(phi / sf_angle + 0.5);
864  input->mus[ipix] = sunv;
865 
866  } /* if (r < 0) ... */
867 
868  } /* for (ipix=0 ... */
869 
870  return 0;
871 
872 }
873 
874 void resample_geo_scan(int Npixels, float *inlon, float *inlat, double *outlon, double *outlat) {
875  int ipix;
876 
877  outlon[0] = -0.5 * inlon[1] + 1.5 * inlon[0];
878  outlat[0] = -0.5 * inlat[1] + 1.5 * inlat[0];
879  if (inlon[1] - inlon[0] >= M_PI) outlon[0] += 1.5 * TWO_PI;
880  if (inlon[0] - inlon[1] >= M_PI) outlon[0] -= 0.5 * TWO_PI;
881  if (outlon[0] > M_PI) outlon[0] -= TWO_PI;
882  if (outlon[0] < -M_PI) outlon[0] += TWO_PI;
883 
884  for (ipix = 0; ipix < Npixels - 1; ipix++) {
885  outlon[ipix + 1] = 0.5 * inlon[ipix] + 0.5 * inlon[ipix + 1];
886  outlat[ipix + 1] = 0.5 * inlat[ipix] + 0.5 * inlat[ipix + 1];
887  if (inlon[ipix] - inlon[ipix + 1] >= M_PI) outlon[ipix + 1] += 0.5 * TWO_PI;
888  if (inlon[ipix + 1] - inlon[ipix] >= M_PI) outlon[ipix + 1] -= 0.5 * TWO_PI;
889  if (outlon[ipix + 1] > M_PI) outlon[ipix + 1] -= TWO_PI;
890  if (outlon[ipix + 1] < -M_PI) outlon[ipix + 1] += TWO_PI;
891  }
892 
893  outlon[Npixels] = 1.5 * inlon[Npixels - 1] - 0.5 * inlon[Npixels - 2];
894  outlat[Npixels] = 1.5 * inlat[Npixels - 1] - 0.5 * inlat[Npixels - 2];
895  if (inlon[Npixels - 1] - inlon[Npixels - 2] >= M_PI) outlon[Npixels] += -0.5 * TWO_PI;
896  if (inlon[Npixels - 2] - inlon[Npixels - 1] >= M_PI) outlon[Npixels] -= -0.5 * TWO_PI;
897  if (outlon[Npixels] > M_PI) outlon[Npixels] -= TWO_PI;
898  if (outlon[Npixels] < -M_PI) outlon[Npixels] += TWO_PI;
899 
900 }
901 
902 int read_file_block(char *filename, char **wl, int offset, int startline, int buflines, OUTPUT *proj) {
903  int32 sd_id, sds_id;
904  int ib;
905  char SDSname[H4_MAX_NC_NAME];
906 
907  printf("Reading %s (lines %d to %d)...\n", filename, startline, startline + buflines - 1);
908  if ((sd_id = SDstart(filename, DFACC_RDWR)) == -1) {
909  fprintf(stderr, "Cannot open file %s\n", filename);
910  return -1;
911  }
912  for (ib = 0; ib < NBANDS; ib++) {
913  if (proj->refl[ib]) {
914  sprintf(SDSname, "refl_%s", wl[ib]);
915  if ((sds_id = open_sds(sd_id, SDSname)) == -1 ||
916  read_sds_block(sds_id, startline, buflines, proj->refl[ib] + offset) == -1)
917  return -1;
918  }
919  if (proj->reflTOA[ib]) {
920  sprintf(SDSname, "reflTOA_%s", wl[ib]);
921  if ((sds_id = open_sds(sd_id, SDSname)) == -1 ||
922  read_sds_block(sds_id, startline, buflines, proj->reflTOA[ib] + offset) == -1)
923  return -1;
924  }
925  }
926  sprintf(SDSname, "input_file");
927  if ((sds_id = open_sds(sd_id, SDSname)) == -1 ||
928  read_sds_block(sds_id, startline, buflines, proj->ifno + offset) == -1)
929  return -1;
930  if (proj->senz) {
931  if (sprintf(SDSname, "senz") == -1 ||
932  (sds_id = open_sds(sd_id, SDSname)) == -1 ||
933  read_sds_block(sds_id, startline, buflines, proj->senz + offset) == -1)
934  return -1;
935  }
936  if (proj->solz) {
937  if (sprintf(SDSname, "solz") == -1 ||
938  (sds_id = open_sds(sd_id, SDSname)) == -1 ||
939  read_sds_block(sds_id, startline, buflines, proj->solz + offset) == -1 ||
940  sprintf(SDSname, "phi") == -1 ||
941  (sds_id = open_sds(sd_id, SDSname)) == -1 ||
942  read_sds_block(sds_id, startline, buflines, proj->phi + offset) == -1)
943  return -1;
944  }
945  if (proj->smoke) {
946  sprintf(SDSname, "smoke_index");
947  if ((sds_id = open_sds(sd_id, SDSname)) == -1 ||
948  read_sds_block(sds_id, startline, buflines, proj->smoke + offset) == -1)
949  return -1;
950  }
951  if (proj->ndvi) {
952  sprintf(SDSname, "NDVI");
953  if ((sds_id = open_sds(sd_id, SDSname)) == -1 ||
954  read_sds_block(sds_id, startline, buflines, proj->ndvi + offset) == -1)
955  return -1;
956  }
957  if (proj->evi) {
958  sprintf(SDSname, "EVI");
959  if ((sds_id = open_sds(sd_id, SDSname)) == -1 ||
960  read_sds_block(sds_id, startline, buflines, proj->evi + offset) == -1)
961  return -1;
962  }
963  if (proj->gemi) {
964  sprintf(SDSname, "GEMI");
965  if ((sds_id = open_sds(sd_id, SDSname)) == -1 ||
966  read_sds_block(sds_id, startline, buflines, proj->gemi + offset) == -1)
967  return -1;
968  }
969  SDend(sd_id);
970 
971  return 0;
972 
973 }
974 
975 int update_file(char *filename, char **wl, int offset, int startline, int buflines, OUTPUT *proj, char verbose, char gzip) {
976  int32 sd_id;
977  int ib;
978  char SDSname[H4_MAX_NC_NAME];
979  int16 fill_int16 = FILL_INT16;
980 
981  if (startline + buflines > proj->h) buflines = proj->h - startline;
982  if (startline >= 0)
983  printf("Updating %s (lines %d to %d)...\n", filename, startline, startline + buflines - 1);
984  if ((sd_id = SDstart(filename, DFACC_RDWR)) == -1) {
985  fprintf(stderr, "Cannot open file %s\n", filename);
986  return -1;
987  }
988  for (ib = 0; ib < NBANDS; ib++) {
989  if (proj->refl[ib]) {
990  sprintf(SDSname, "refl_%s", wl[ib]);
991  write_sds(sd_id, SDSname, proj->refl[ib] + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_refl, verbose, gzip);
992  }
993  }
994  for (ib = 0; ib < NBANDS; ib++) {
995  if (proj->reflTOA[ib]) {
996  sprintf(SDSname, "reflTOA_%s", wl[ib]);
997  write_sds(sd_id, SDSname, proj->reflTOA[ib] + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_refl, verbose, gzip);
998  }
999  }
1000  sprintf(SDSname, "input_file");
1001  write_sds(sd_id, SDSname, proj->ifno + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, 0, verbose, gzip);
1002  if (proj->senz) {
1003  sprintf(SDSname, "senz");
1004  write_sds(sd_id, SDSname, proj->senz + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_angle, verbose, gzip);
1005  }
1006  if (proj->solz) {
1007  sprintf(SDSname, "solz");
1008  write_sds(sd_id, SDSname, proj->solz + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_angle, verbose, gzip);
1009  sprintf(SDSname, "phi");
1010  write_sds(sd_id, SDSname, proj->phi + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_angle, verbose, gzip);
1011  }
1012  if (proj->smoke) {
1013  sprintf(SDSname, "smoke_index");
1014  write_sds(sd_id, SDSname, proj->smoke + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_smoke, verbose, gzip);
1015  }
1016  if (proj->ndvi) {
1017  sprintf(SDSname, "NDVI");
1018  write_sds(sd_id, SDSname, proj->ndvi + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_vi, verbose, gzip);
1019  }
1020  if (proj->evi) {
1021  sprintf(SDSname, "EVI");
1022  write_sds(sd_id, SDSname, proj->evi + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_vi, verbose, gzip);
1023  }
1024  if (proj->gemi) {
1025  sprintf(SDSname, "GEMI");
1026  write_sds(sd_id, SDSname, proj->gemi + offset, startline, buflines, proj->h, proj->w, DFNT_INT16, &fill_int16, sf_vi, verbose, gzip);
1027  }
1028 
1029  if (SDsetattr(sd_id, "File list", DFNT_CHAR8, (int) strlen(proj->filelist) + 1, proj->filelist) == -1)
1030  fprintf(stderr, "Cannot write attribute in file %s\n", filename);
1031 
1032  SDend(sd_id);
1033  return 0;
1034 }
1035 
1036 void memcpy_block(int offset2, int offset1, int Nrows, int Ncols, OUTPUT *proj) {
1037  int ib, k, irow;
1038  if (offset1 < 0 || offset2 < 0) return;
1039  for (k = 0; k < Nrows; k++) {
1040  if (offset2 < offset1) irow = k;
1041  else irow = Nrows - 1 - k;
1042  for (ib = 0; ib < NBANDS; ib++) {
1043  if (proj->refl[ib])
1044  memcpy(proj->refl[ib] + (offset2 + irow) * Ncols, proj->refl[ib] + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1045  if (proj->reflTOA[ib])
1046  memcpy(proj->reflTOA[ib] + (offset2 + irow) * Ncols, proj->reflTOA[ib] + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1047  }
1048  memcpy(proj->ifno + (offset2 + irow) * Ncols, proj->ifno + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1049  if (proj->senz)
1050  memcpy(proj->senz + (offset2 + irow) * Ncols, proj->senz + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1051  if (proj->solz) {
1052  memcpy(proj->solz + (offset2 + irow) * Ncols, proj->solz + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1053  memcpy(proj->phi + (offset2 + irow) * Ncols, proj->phi + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1054  }
1055  if (proj->smoke)
1056  memcpy(proj->smoke + (offset2 + irow) * Ncols, proj->smoke + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1057  if (proj->ndvi)
1058  memcpy(proj->ndvi + (offset2 + irow) * Ncols, proj->ndvi + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1059  if (proj->evi)
1060  memcpy(proj->evi + (offset2 + irow) * Ncols, proj->evi + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1061  if (proj->gemi)
1062  memcpy(proj->gemi + (offset2 + irow) * Ncols, proj->gemi + (offset1 + irow) * Ncols, Ncols * sizeof (int16));
1063  }
1064 }
1065 
1066 void process_pixel(int ipix, int iline_out, int ipix_out, char **wl, INPUT input, OUTPUT *proj, PROCESSINFO process) {
1067  int idx_out, old_startline, offset;
1068  int16 ndvi;
1069  char better;
1070 
1071  iline_out -= proj->l1;
1072  ipix_out -= proj->p1;
1073 
1074  if (iline_out < proj->startline ||
1075  iline_out > proj->startline + proj->buflines - 1) {
1076  old_startline = proj->startline;
1077  if (proj->startline < 0)
1078  proj->startline = iline_out - iline_out % proj->bufstep;
1079  while (iline_out < proj->startline)
1080  proj->startline -= proj->bufstep;
1081  while (iline_out > proj->startline + proj->buflines - 1)
1082  proj->startline += proj->bufstep;
1083  if (proj->startline < 0)
1084  proj->startline = 0;
1085  if (proj->startline + proj->buflines > proj->h)
1086  proj->startline = proj->h - proj->buflines;
1087 
1088  if (proj->startline > old_startline &&
1089  proj->startline < old_startline + proj->buflines - 1) {
1090  offset = proj->startline - old_startline;
1091  if (update_file(proj->filename, wl, 0, old_startline, offset, proj, FALSE, process.gzip) == -1)
1092  exit(-1);
1093  memcpy_block(0, offset, proj->buflines - offset, proj->w, proj);
1094  if (read_file_block(proj->filename, wl, (proj->buflines - offset) * proj->w, proj->startline + proj->buflines - offset,
1095  MIN(offset, proj->h - proj->startline - proj->buflines + offset), proj) == -1)
1096  exit(-1);
1097  } else if (old_startline > proj->startline &&
1098  old_startline < proj->startline + proj->buflines - 1) {
1099  offset = old_startline - proj->startline;
1100  if (update_file(proj->filename, wl, (proj->buflines - offset) * proj->w, old_startline + proj->buflines - offset,
1101  offset, proj, FALSE, process.gzip) == -1)
1102  exit(-1);
1103  memcpy_block(offset, 0, proj->buflines - offset, proj->w, proj);
1104  if (read_file_block(proj->filename, wl, 0, proj->startline,
1105  MIN(offset, proj->h - proj->startline), proj) == -1)
1106  exit(-1);
1107  } else {
1108  if (update_file(proj->filename, wl, 0, old_startline, proj->buflines, proj, FALSE, process.gzip) == -1)
1109  exit(-1);
1110  if (read_file_block(proj->filename, wl, 0, proj->startline,
1111  MIN(proj->buflines, proj->h - proj->startline), proj) == -1)
1112  exit(-1);
1113  }
1114  }
1115 
1116  idx_out = (iline_out - proj->startline) * proj->w + ipix_out;
1117  better = FALSE;
1118 
1119  if (proj->ifno[idx_out] == FILL_INT16) /* First observation for this pixel */
1120 
1121  better = TRUE;
1122 
1123  else {
1124 
1125  switch (process.composite) {
1126 
1127  case MINBLUE:
1128  if (input.rho[BLUEBAND][ipix] >= MINREFL &&
1129  (input.rho[BLUEBAND][ipix] < proj->refl[BLUEBAND][idx_out] * sf_refl ||
1130  proj->refl[BLUEBAND][idx_out] == FILL_INT16))
1131  better = TRUE;
1132  break;
1133 
1134  case MAXNDVI:
1135  if (input.level == L1A)
1136  ndvi = NDVI(input.rho[BAND6][ipix], input.rho[BAND8][ipix]);
1137  else
1138  ndvi = input.ndvi[ipix];
1139  if (ndvi != FILL_INT16 &&
1140  (ndvi > proj->ndvi[idx_out] ||
1141  proj->ndvi[idx_out] == FILL_INT16))
1142  better = TRUE;
1143  break;
1144 
1145  case MINSENZ:
1146  if (input.senz[ipix] < proj->senz[idx_out])
1147  better = TRUE;
1148  break;
1149 
1150  case LASTIN:
1151  better = TRUE;
1152  break;
1153 
1154  }
1155  }
1156 
1157  if (better) {
1158  if (proj->ndvi &&
1159  (proj->ifno[idx_out] == FILL_INT16 ||
1160  process.composite != MAXNDVI)) {
1161  if (input.level == L1A)
1162  ndvi = NDVI(input.rho[BAND6][ipix], input.rho[BAND8][ipix]);
1163  else
1164  ndvi = input.ndvi[ipix];
1165  }
1166  update_pixel(input, ndvi, ipix, idx_out, proj->Nfiles, proj);
1167  }
1168 
1169 }
1170 
1171 int get_xyproj(double lon, double lat, double *xproj, double *yproj, int projtype, double lon_center) {
1172 
1173  lon -= lon_center;
1174  if (lon < -M_PI) lon += TWO_PI;
1175  if (lon > M_PI) lon -= TWO_PI;
1176 
1177  switch (projtype) {
1178  case PLATECARREE: *xproj = lon * RAD2DEG; /* PLATECARREE was not initialized with RAD2DEG factor */
1179  *yproj = lat * RAD2DEG;
1180  break;
1181  case SINUSOIDAL: *xproj = lon * cos(lat) * RAD2DEG; /* SINUSOIDAL was not initialized with RAD2DEG factor */
1182  *yproj = lat * RAD2DEG;
1183  break;
1184  case MOLLWEIDE: if (molwfor(lon, lat, xproj, yproj) == -1) return -1;
1185  break;
1186  case ROBINSON: if (robfor(lon, lat, xproj, yproj) == -1) return -1;
1187  break;
1188  case HAMMER: if (hamfor(lon, lat, xproj, yproj) == -1) return -1;
1189  break;
1190  default: return -1;
1191  }
1192  return 0;
1193 }
1194 
1195 int get_lonlat(double xproj, double yproj, double *lon, double *lat, int projtype, double lon_center) {
1196 
1197  switch (projtype) {
1198  case PLATECARREE: *lon = xproj / RAD2DEG;
1199  *lat = yproj / RAD2DEG;
1200  break;
1201  case SINUSOIDAL: *lat = yproj / RAD2DEG;
1202  if (fabs(*lat) >= M_PI_2) return -1;
1203  *lon = xproj / cos(*lat) / RAD2DEG;
1204  if (fabs(*lon) > M_PI) return -1;
1205  break;
1206  case MOLLWEIDE: if (molwinv(xproj, yproj, lon, lat) == -1) return -1;
1207  break;
1208  case ROBINSON: if (robinv(xproj, yproj, lon, lat) == -1) return -1;
1209  break;
1210  case HAMMER: if (haminv(xproj, yproj, lon, lat) == -1) return -1;
1211  break;
1212  default: return -1;
1213  }
1214 
1215  *lon += lon_center;
1216  if (*lon < -M_PI) *lon += TWO_PI;
1217  if (*lon > M_PI) *lon -= TWO_PI;
1218 
1219  return 0;
1220 
1221 }
1222 
1223 void update_pixel(INPUT input, int16 ndvi, int ipix, int idx_out, int ifno, OUTPUT *proj) {
1224  int ib;
1225 
1226  for (ib = 0; ib < NBANDS; ib++) {
1227  if (input.rho[ib]) {
1228  if (input.rho[ib][ipix] != UNDEF &&
1229  input.rho[ib][ipix] / sf_refl <= SHRT_MAX)
1230  proj->refl[ib][idx_out] = (int16) floor(input.rho[ib][ipix] / sf_refl + 0.5);
1231  else
1232  proj->refl[ib][idx_out] = FILL_INT16;
1233  }
1234 
1235  if (input.rhoTOA[ib]) {
1236  proj->reflTOA[ib][idx_out] = (int16) floor(input.rhoTOA[ib][ipix] / sf_refl + 0.5);
1237  }
1238  }
1239 
1240 
1241  if (proj->senz)
1242  proj->senz[idx_out] = input.senz[ipix];
1243  if (proj->solz)
1244  proj->solz[idx_out] = input.solz[ipix];
1245  if (proj->phi)
1246  proj->phi[idx_out] = input.phi[ipix];
1247 
1248  /*
1249  mask_lonlat = read_mask((double) (input.lon[ipix]*RAD2DEG),
1250  (double) (input.lat[ipix]*RAD2DEG), (int16) 0);
1251 
1252 
1253  mask_lonlat = (input.l2_flags[ipix] & 0x2) / 0x2;
1254 
1255  if (mask_lonlat == 0 && proj->ifno[idx_out] == FILL_INT16) proj->ifno[idx_out] = WATER_VAL;
1256  if (mask_lonlat == 1) proj->ifno[idx_out] = ifno;
1257  */
1258 
1259  proj->ifno[idx_out] = ifno;
1260 
1261  if (proj->smoke)
1262  proj->smoke[idx_out] = SMOKE(input.rho, ipix);
1263 
1264  if (proj->ndvi)
1265  proj->ndvi[idx_out] = ndvi;
1266 
1267  if (proj->evi) {
1268  if (input.level == L1A)
1269  proj->evi[idx_out] = EVI(input.rho[BAND2][ipix], input.rho[BAND6][ipix], input.rho[BAND8][ipix]);
1270  else
1271  proj->evi[idx_out] = input.evi[ipix];
1272  }
1273 
1274  if (proj->gemi)
1275  proj->gemi[idx_out] = GEMI(input.rho[BAND6][ipix], input.rho[BAND8][ipix]);
1276 
1277 }
1278 
1279 int16 SMOKE(float **rho, int ipix) {
1280  if (rho[BAND1][ipix] < 0 ||
1281  rho[BAND6][ipix] > 0.4)
1282  return FILL_INT16;
1283  else if (rho[BAND2][ipix] < rho[BAND4][ipix] ||
1284  rho[BAND5][ipix] < rho[BAND6][ipix])
1285  return 0;
1286  else
1287  return (int16) (rho[BAND1][ipix] * (rho[BAND2][ipix] - rho[BAND4][ipix]) * (rho[BAND5][ipix] - rho[BAND6][ipix]) * 10000);
1288 }
1289 
1290 int16 NDVI(float red, float nir) {
1291  float ndvi;
1292 
1293  if (red != UNDEF &&
1294  nir != UNDEF &&
1295  red + nir != 0) {
1296  ndvi = (nir - red) / (nir + red);
1297  if (ndvi < -2.) ndvi = -2.;
1298  if (ndvi > 2.) ndvi = 2.;
1299  return (int16) floor(ndvi / sf_vi + 0.5);
1300  } else
1301  return FILL_INT16;
1302 }
1303 
1304 int16 EVI(float blue, float red, float nir) {
1305  float evi;
1306  static float L = 1, c1 = 6, c2 = 7.5;
1307  double val;
1308 
1309  if (red == UNDEF ||
1310  nir == UNDEF)
1311  return FILL_INT16;
1312  else {
1313  if (blue != UNDEF && /* Most cases - EVI formula */
1314  (blue <= red || red <= nir)) {
1315  if ((val = L + nir + c1 * red - c2 * blue) == 0)
1316  return FILL_INT16;
1317  else {
1318  evi = 2.5 * (nir - red) / val;
1319  if (evi < -2.) evi = -2.;
1320  if (evi > 2.) evi = 2.;
1321  }
1322  } else { /* Backup - SAVI formula */
1323  if ((val = 0.5 + nir + red) == 0)
1324  return FILL_INT16;
1325  else {
1326  evi = 1.5 * (nir - red) / val;
1327  if (evi < -2.) evi = -2.;
1328  if (evi > 2.) evi = 2.;
1329  }
1330  }
1331  return (int16) floor(evi / sf_vi + 0.5);
1332  }
1333 
1334 }
1335 
1336 #if 0
1337 
1338 int16 EVI(float blue, float red, float nir) {
1339  float evi;
1340  static float L = 1, c1 = 6, c2 = 7.5;
1341  double val;
1342 
1343  if (blue != UNDEF &&
1344  red != UNDEF &&
1345  nir != UNDEF &&
1346  (val = L + nir + c1 * red - c2 * blue) != 0) {
1347  evi = 2.5 * (nir - red) / val;
1348  if (evi < -2.) evi = -2.;
1349  if (evi > 2.) evi = 2.;
1350  return (int16) floor(evi / sf_vi + 0.5);
1351  } else
1352  return FILL_INT16;
1353 }
1354 #endif
1355 
1356 int16 GEMI(float red, float nir) {
1357  float gemi;
1358  double x;
1359  if (red != UNDEF &&
1360  nir != UNDEF &&
1361  red + nir + 0.5 != 0 &&
1362  nir != 1.) {
1363  x = (2 * (nir * nir - red * red) + 1.5 * nir + 0.5 * red) / (red + nir + 0.5);
1364  gemi = x * (1 - 0.25 * x) - (red - 0.125) / (1 - red);
1365  if (gemi < -2.) gemi = -2.;
1366  if (gemi > 2.) gemi = 2.;
1367  return (int16) floor(gemi / sf_vi + 0.5);
1368  } else
1369  return FILL_INT16;
1370 }
1371 
1372 int alloc_new_array(void **data, int size, int32 num_type, int maxmem) {
1373  int sizeoftype;
1374 
1375  if ((*data = (void *) malloc(size * (sizeoftype = DFKNTsize(num_type)))) == NULL) { /* Allocate a new array */
1376  printf(" Sorry, can't allocate %.1f Mbytes of memory.\n", size * sizeoftype / 1e6);
1377  exit(-1);
1378  }
1379  printf(" (%.1f Mbytes allocated: %d items of %d bytes)\n", size * sizeoftype / 1e6, size, sizeoftype);
1380  return 0;
1381 }
1382 
1383 void write_cmd_line(int32 sd_id, char *cmd_line) {
1384  char *old_cmd_line, *new_cmd_line;
1385  int32 attr_index, num_type, count;
1386  char name[H4_MAX_NC_NAME] = "Command line";
1387 
1388  new_cmd_line = strdup(cmd_line); /* Default value is current command line */
1389 
1390  if ((attr_index = SDfindattr(sd_id, name)) != -1 &&
1391  SDattrinfo(sd_id, attr_index, name, &num_type, &count) != -1) {
1392 
1393  old_cmd_line = (char *) malloc((count + 1) * sizeof (char));
1394  if (SDreadattr(sd_id, attr_index, old_cmd_line) != -1) {
1395  new_cmd_line = (char *) realloc(new_cmd_line, (count + strlen(cmd_line) + 1 + 1) * sizeof (char));
1396  sprintf(new_cmd_line, "%s%s", old_cmd_line, cmd_line);
1397  }
1398  free(old_cmd_line);
1399 
1400  }
1401 
1402  if (SDsetattr(sd_id, name, DFNT_CHAR8, (int) strlen(new_cmd_line) + 1, new_cmd_line) == -1)
1403  printf("Unable to write command line in file attributes.\n");
1404 
1405  free(new_cmd_line);
1406 
1407 }
1408 
1410 
1411  printf("Pixel size for output projection: %g m\n", proj->pixsz);
1412  proj->pixsz /= R / RAD2DEG;
1413 
1414  switch (proj->projection) {
1415  case PLATECARREE: sprintf(proj->projname, "Plate Carree");
1416  proj->Nl = ceil(180. / proj->pixsz);
1417  proj->Np = ceil(360. / proj->pixsz);
1418  break;
1419  case SINUSOIDAL: sprintf(proj->projname, "Sinusoidal");
1420  proj->Nl = ceil(180. / proj->pixsz);
1421  proj->Np = ceil(360. / proj->pixsz);
1422  break;
1423  case MOLLWEIDE: sprintf(proj->projname, "Mollweide");
1424  proj->Nl = ceil(2 * sqrt(2.) / M_PI * 180. / proj->pixsz);
1425  proj->Np = ceil(2 * sqrt(2.) / M_PI * 360. / proj->pixsz);
1426  molwforint(RAD2DEG, 0., 0., 0.);
1427  molwinvint(RAD2DEG, 0., 0., 0.);
1428  break;
1429  case ROBINSON: sprintf(proj->projname, "Robinson");
1430  proj->Nl = ceil(180. / proj->pixsz);
1431  proj->Np = ceil(0.9858 * 360. / proj->pixsz);
1432  robforint(RAD2DEG, 0., 0., 0.);
1433  robinvint(RAD2DEG, 0., 0., 0.);
1434  break;
1435  case HAMMER: sprintf(proj->projname, "Hammer-Aitoff");
1436  proj->Nl = ceil(2 * RAD2DEG * sqrt(2.) / proj->pixsz);
1437  proj->Np = ceil(4 * RAD2DEG * sqrt(2.) / proj->pixsz);
1438  hamforint(RAD2DEG, 0., 0., 0.);
1439  haminvint(RAD2DEG, 0., 0., 0.);
1440  break;
1441  default: fprintf(stderr, "Wrong selection for output projection\n");
1442  exit(-1);
1443  }
1444 
1445  printf("Output projection: %s\n", proj->projname);
1446 
1447  if (proj->domain == USERBOX &&
1448  (proj->l1 < 0 ||
1449  proj->p1 < 0 ||
1450  proj->l2 > proj->Nl - 1 ||
1451  proj->p2 > proj->Np - 1)) {
1452  if (proj->l1 < 0) proj->l1 = 0;
1453  if (proj->p1 < 0) proj->p1 = 0;
1454  if (proj->l2 > proj->Nl - 1) proj->l2 = proj->Nl - 1;
1455  if (proj->p2 > proj->Np - 1) proj->p2 = proj->Np - 1;
1456  printf("Requested subwindow cannot exceed the projection dimensions.\n");
1457  printf("Output subwindow has been adjusted: (column %d, row %d) to (column %d, row %d)\n",
1458  proj->p1, proj->l1, proj->p2, proj->l2);
1459  }
1460 
1461  if (proj->domain == GLOBAL) { /* No bounding box has been requested */
1462  printf("No bounding box requested. Output will be a global projection.\n");
1463  proj->l1 = 0;
1464  proj->p1 = 0;
1465  proj->l2 = proj->Nl - 1;
1466  proj->p2 = proj->Np - 1;
1467  }
1468 
1469  proj->h = proj->l2 - proj->l1 + 1;
1470  proj->w = proj->p2 - proj->p1 + 1;
1471  printf("Output projection size: %dx%d\n", proj->w, proj->h);
1472 
1473  proj->xmin = proj->pixsz * (proj->p1 - (proj->Np - 1) / 2. - 0.5);
1474  proj->xmax = proj->pixsz * (proj->p2 - (proj->Np - 1) / 2. + 0.5);
1475  proj->ymin = proj->pixsz * ((proj->Nl - 1) / 2. - proj->l2 - 0.5);
1476  proj->ymax = proj->pixsz * ((proj->Nl - 1) / 2. - proj->l1 + 0.5);
1477 
1478 }
1479 
1480 void parse_command_line(int argc, char **argv, OUTPUT *proj, PROCESSINFO *process) {
1481  char *arg, *ptr;
1482  int j, ib, projection;
1483 
1484  proj->cmd_line = (char *) malloc((strlen(argv[0]) + 1 + 1) * sizeof (char));
1485  sprintf(proj->cmd_line, "\n%s", argv[0]);
1486 
1487  for (j = 1; j < argc; j++) {
1488 
1489  arg = argv[j];
1490 
1491  proj->cmd_line = (char *) realloc(proj->cmd_line, (strlen(proj->cmd_line) + strlen(arg) + 1 + 1) * sizeof (char));
1492  strcat(proj->cmd_line, " ");
1493  strcat(proj->cmd_line, arg);
1494 
1495  if (strstr(arg, "-geobox=") == arg) {
1496  if (sscanf(arg + strlen("-geobox="), "%lg,%lg,%lg,%lg",
1497  &proj->lon1, &proj->lat1, &proj->lon2, &proj->lat2) == 4) {
1498  printf("Geographic area selected: (lon,lat)=(%g,%g)->(%g,%g)\n",
1499  proj->lon1, proj->lat1, proj->lon2, proj->lat2);
1500  if (proj->lon1 >= proj->lon2 ||
1501  proj->lat1 <= proj->lat2 ||
1502  proj->lon1 < -180 ||
1503  proj->lon2 > 180 ||
1504  proj->lat1 > 90 ||
1505  proj->lat2 < -90) {
1506  printf("Bad specification for geographic area. Request ignored.\n");
1507  proj->lon1 = proj->lon2 = proj->lat1 = proj->lat2 = UNDEF;
1508  } else {
1509  proj->lon1 /= RAD2DEG;
1510  proj->lon2 /= RAD2DEG;
1511  proj->lat1 /= RAD2DEG;
1512  proj->lat2 /= RAD2DEG;
1513  }
1514  printf("\n");
1515  }
1516  continue;
1517  }
1518 
1519  if (strstr(arg, "-box=") == arg) {
1520  if (sscanf(arg + strlen("-box="), "%d,%d,%d,%d",
1521  &proj->p1, &proj->l1, &proj->p2, &proj->l2) == 4) {
1522  printf("Output subwindow selected: (column %d, row %d) to (column %d, row %d)\n",
1523  proj->p1, proj->l1, proj->p2, proj->l2);
1524  if (proj->l1 >= proj->l2 ||
1525  proj->p1 >= proj->p2)
1526  printf("Bad specification for output subwindow. Request ignored.\n");
1527  else
1528  proj->domain = USERBOX;
1529  printf("\n");
1530  }
1531  continue;
1532  }
1533 
1534  if (strstr(arg, "-bands=") == arg) {
1535  arg += strlen("-bands=");
1536  while ((ptr = strchr(arg, 44))) *ptr = 32;
1537  ptr = arg;
1538  while (sscanf(ptr, "%d", &ib) == 1) {
1539  if (ib >= 1 && ib <= NBANDS) process->band[ib - 1] = TRUE;
1540  if ((ptr = strchr(ptr, 32))) ptr++;
1541  else ptr = strchr(arg, 0);
1542  }
1543  continue;
1544  }
1545 
1546  if (strstr(arg, "-center=") == arg) {
1547  arg += strlen("-center=");
1548  if (sscanf(arg, "%lg", &proj->lon_center) == 1)
1549  printf("Center longitude for output projection: %g deg.\n", proj->lon_center);
1550  proj->lon_center /= RAD2DEG;
1551  continue;
1552  }
1553 
1554  if (strstr(arg, "-pixsz=") == arg) {
1555  arg += strlen("-pixsz=");
1556  sscanf(arg, "%lg", &proj->pixsz);
1557  continue;
1558  }
1559 
1560  if (strstr(arg, "-maxsenz=") == arg) {
1561  arg += strlen("-maxsenz=");
1562  sscanf(arg, "%g", &process->maxsenz);
1563  printf("Maximum sensor zenith angle allowed: %g deg.\n", process->maxsenz);
1564  continue;
1565  }
1566 
1567  if (strstr(arg, "-maxmem=") == arg) {
1568  arg += strlen("-maxmem=");
1569  sscanf(arg, "%d", &process->maxmem);
1570  printf("Maximum memory allocation allowed: %d Mbytes\n", process->maxmem);
1571  continue;
1572  }
1573 
1574  if (strstr(arg, "-daystr=") == arg) {
1575  arg += strlen("-daystr=");
1576  strncpy(&process->daystr[0], arg, 13);
1577  continue;
1578  }
1579 
1580  if (strstr(arg, "-dayend=") == arg) {
1581  arg += strlen("-dayend=");
1582  strncpy(&process->dayend[0], arg, 13);
1583  continue;
1584  }
1585 
1586  if (strcmp(arg, "-smoke") == 0) {
1587  printf("Smoke index requested\n");
1588  process->smoke = TRUE;
1589  continue;
1590  }
1591 
1592  if (strcmp(arg, "-ndvi") == 0) {
1593  printf("NDVI requested\n");
1594  process->ndvi = TRUE;
1595  continue;
1596  }
1597 
1598  if (strcmp(arg, "-evi") == 0) {
1599  printf("EVI requested\n");
1600  process->evi = TRUE;
1601  continue;
1602  }
1603 
1604  if (strcmp(arg, "-gemi") == 0) {
1605  printf("GEMI requested\n");
1606  process->gemi = TRUE;
1607  continue;
1608  }
1609 
1610  if (strcmp(arg, "-angles") == 0) {
1611  printf("View and sun angles requested\n");
1612  process->angles = TRUE;
1613  continue;
1614  }
1615 
1616  if (strcmp(arg, "-gzip") == 0) {
1617  printf("Compression of output SDSs requested\n");
1618  process->gzip = TRUE;
1619  continue;
1620  }
1621 
1622  if (strcmp(arg, "-TOA") == 0) {
1623  printf("TOA reflectance requested in output\n");
1624  process->TOA = TRUE;
1625  continue;
1626  }
1627 
1628  if (strstr(arg, "-method=") == arg) {
1629  arg += strlen("-method=");
1630  if (strcmp(arg, "minblue") == 0)
1631  process->composite = MINBLUE;
1632  else if (strcmp(arg, "maxndvi") == 0)
1633  process->composite = MAXNDVI;
1634  else if (strcmp(arg, "minsenz") == 0)
1635  process->composite = MINSENZ;
1636  else if (strcmp(arg, "lastin") == 0)
1637  process->composite = LASTIN;
1638  else
1639  printf("Compositing method \"%s\" not supported\n", arg);
1640  continue;
1641  }
1642 
1643  if (strstr(arg, "-proj=") == arg) {
1644  arg += strlen("-proj=");
1645  projection = DEFPROJECTION;
1646  if (sscanf(arg, "%d", &projection) == 1) {
1647  if (projection < 0 ||
1648  projection > NUMBEROFPROJECTIONS - 1) { /* char is >= 0 anyway */
1649  printf("Invalid projection requested (%d). Request ignored.\n", projection);
1650  projection = DEFPROJECTION;
1651  }
1652  proj->projection = projection;
1653  }
1654  continue;
1655  }
1656 
1657  if (strstr(arg, "-bufstep=") == arg) {
1658  arg += strlen("-bufstep=");
1659  if (sscanf(arg, "%d", &proj->bufstep) == 1)
1660  printf("Buffer step requested: %d lines\n", proj->bufstep);
1661  continue;
1662  }
1663 
1664  if (strstr(arg, "-of=") == arg) {
1665  arg += strlen("-of=");
1666  strcpy(proj->filename, arg);
1667  printf("Output filename requested: %s\n", proj->filename);
1668  continue;
1669  }
1670 
1671 
1672  if (strstr(arg, "-input=") == arg) {
1673  arg += strlen("-input=");
1674 
1675  proj->input_fp = fopen(arg, "r");
1676  if (proj->input_fp == NULL) {
1677  printf("Input listing file: \"%s\" not found.\n", arg);
1678  exit(1);
1679  }
1680  printf("Input filelist filename requested: %s\n", arg);
1681  continue;
1682  }
1683 
1684  if (strstr(arg, "-fileuse=") == arg) {
1685  arg += strlen("-fileuse=");
1686  proj->fileuse_fp = fopen(arg, "w");
1687  printf("File Use file requested: %s\n", arg);
1688  continue;
1689  }
1690 
1691  }
1692 
1693 }
1694 
1696  int nbytes, ib;
1697 
1698  nbytes = 0;
1699 
1700  for (ib = 0; ib < NBANDS; ib++)
1701  if (process->band[ib]) {
1702  nbytes += DFKNTsize(DFNT_INT16); /* Each individual band */
1703  if (process->TOA) nbytes += DFKNTsize(DFNT_INT16); /* TOA reflectance */
1704  }
1705 
1706  nbytes += DFKNTsize(DFNT_INT16); /* Input file SDS */
1707 
1708  if (process->angles ||
1709  process->composite == MINSENZ)
1710  nbytes += DFKNTsize(DFNT_INT16); /* View zenith angle */
1711 
1712  if (process->angles)
1713  nbytes += 2 * DFKNTsize(DFNT_INT16); /* Solar zenith angle and relative azimuth angle */
1714 
1715  if (process->smoke)
1716  nbytes += DFKNTsize(DFNT_INT16); /* Smoke index */
1717 
1718  if (process->ndvi ||
1719  process->composite == MAXNDVI)
1720  nbytes += DFKNTsize(DFNT_INT16); /* NDVI */
1721 
1722  if (process->evi)
1723  nbytes += DFKNTsize(DFNT_INT16); /* EVI */
1724 
1725  if (process->gemi)
1726  nbytes += DFKNTsize(DFNT_INT16); /* GEMI */
1727 
1728  if (process->atmcor)
1729  proj->buflines = MIN(proj->h,
1730  (int) ((1e6 * process->maxmem - NLPSURF * NPPSURF * DFKNTsize(DFNT_INT16)) / proj->w / nbytes));
1731  else
1732  proj->buflines = MIN(proj->h,
1733  (int) (1e6 * process->maxmem / proj->w / nbytes));
1734 
1735  if (proj->buflines < proj->h) {
1736 
1737  printf("Output projection has to be bufferized: %d lines out of %d\n", proj->buflines, proj->h);
1738 
1739  /* bufstep is the shift (in # of lines) to apply to the buffer when pixel is out of the current buffer */
1740  /* If bufstep is too small the buffer will be updated too often */
1741  /* If bufstep is too big the buffer will be flushed before all lines are */
1742  /* processed, and the buffer will go back and forth in the output projection */
1743 
1744  if (proj->bufstep == -1) {
1745  proj->bufstep = proj->buflines / 5 + 1;
1746  printf("No buffer step requested. %d will be used.\n", proj->bufstep);
1747  } else {
1748  if (proj->bufstep > proj->h - proj->buflines) {
1749  proj->bufstep = proj->h - proj->buflines;
1750  printf("Buffer step adjusted to %d lines\n", proj->bufstep);
1751  }
1752  }
1753  if (process->gzip) {
1754  printf("Compression cannot be performed on bufferized output SDS. Request ignored.\n");
1755  process->gzip = FALSE;
1756  }
1757 
1758  }
1759 
1760 }
1761 
1762 void check_process_parameters(PROCESSINFO *process, char **wl) {
1763  int ib;
1764 
1765  printf("Requested bands: ");
1766  for (ib = 0; ib < NBANDS; ib++)
1767  if (process->band[ib])
1768  break;
1769  if (ib == NBANDS) /* No band was specified - Process all bands */
1770  for (ib = 0; ib < NBANDS; ib++)
1771  process->band[ib] = TRUE;
1772  for (ib = 0; ib < NBANDS; ib++)
1773  if (process->band[ib])
1774  printf("%s ", wl[ib]);
1775  printf("\n");
1776 
1777  switch (process->composite) {
1778  case MINBLUE: sprintf(process->compositename, "minimum band %d (%s nm)", BLUEBAND + 1, wl[BLUEBAND]);
1779  break;
1780  case MAXNDVI: strcpy(process->compositename, "maximum NDVI");
1781  break;
1782  case MINSENZ: strcpy(process->compositename, "minimum view angle");
1783  break;
1784  case LASTIN: strcpy(process->compositename, "last in on top");
1785  break;
1786  default: fprintf(stderr, "Invalid composite criterion (%d).\n", process->composite);
1787  exit(-1);
1788  }
1789  printf("Compositing criterion: %s\n", process->compositename);
1790 
1791  if (!process->angles &&
1792  process->composite == MINSENZ)
1793  printf("View angle will be computed.\n");
1794 
1795  printf("\n");
1796 
1797  if (process->composite == MINBLUE &&
1798  !process->band[BLUEBAND]) {
1799  fprintf(stderr, "Cannot make composite without band %d (%s nm)\n", BLUEBAND + 1, wl[BLUEBAND]);
1800  exit(-1);
1801  }
1802 
1803  if (process->smoke) {
1804  if (!process->band[BAND1] ||
1805  !process->band[BAND2] ||
1806  !process->band[BAND4] ||
1807  !process->band[BAND5] ||
1808  !process->band[BAND6]) {
1809  fprintf(stderr, "Cannot compute smoke index without bands %s %s %s %s %s\n",
1810  wl[BAND1], wl[BAND2], wl[BAND4], wl[BAND5], wl[BAND6]);
1811  exit(-1);
1812  }
1813  }
1814 
1815  if (process->evi) {
1816  if (!process->band[BAND2] ||
1817  !process->band[BAND6] ||
1818  !process->band[BAND8]) {
1819  fprintf(stderr, "Cannot compute EVI without bands %s, %s and %s\n", wl[BAND2], wl[BAND6], wl[BAND8]);
1820  exit(-1);
1821  }
1822  }
1823 
1824  if (process->ndvi ||
1825  process->composite == MAXNDVI) {
1826  if (!process->band[BAND6] ||
1827  !process->band[BAND8]) {
1828  fprintf(stderr, "Cannot compute NDVI without bands %s and %s\n", wl[BAND6], wl[BAND8]);
1829  exit(-1);
1830  }
1831  }
1832 
1833  if (process->gemi) {
1834  if (!process->band[BAND6] ||
1835  !process->band[BAND8]) {
1836  fprintf(stderr, "Cannot compute GEMI without bands %s and %s\n", wl[BAND6], wl[BAND8]);
1837  exit(-1);
1838  }
1839  }
1840 
1841  if (!process->band[BAND8]) {
1842  fprintf(stderr, "Cannot estimate cloudiness without band %s\n", wl[BAND8]);
1843  exit(-1);
1844  }
1845 
1846 }
1847 
1849  int ib, j;
1850 
1851  for (ib = 0; ib < NBANDS; ib++) {
1852  proj->refl[ib] = NULL;
1853  if (process.band[ib]) {
1854  alloc_new_array((void *) &proj->refl[ib], proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1855  for (j = 0; j < proj->buflines * proj->w; j++)
1856  proj->refl[ib][j] = FILL_INT16;
1857  }
1858  }
1859 
1860  for (ib = 0; ib < NBANDS; ib++) {
1861  proj->reflTOA[ib] = NULL;
1862  if (process.TOA &&
1863  process.band[ib]) {
1864  alloc_new_array((void *) &proj->reflTOA[ib], proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1865  for (j = 0; j < proj->buflines * proj->w; j++)
1866  proj->reflTOA[ib][j] = FILL_INT16;
1867  }
1868  }
1869 
1870  alloc_new_array((void *) &proj->ifno, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1871  for (j = 0; j < proj->buflines * proj->w; j++)
1872  proj->ifno[j] = FILL_INT16;
1873 
1874  if (process.angles ||
1875  process.composite == MINSENZ) {
1876  alloc_new_array((void *) &proj->senz, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1877  for (j = 0; j < proj->buflines * proj->w; j++)
1878  proj->senz[j] = FILL_INT16;
1879  }
1880 
1881  if (process.angles) {
1882  alloc_new_array((void *) &proj->solz, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1883  alloc_new_array((void *) &proj->phi, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1884  for (j = 0; j < proj->buflines * proj->w; j++)
1885  proj->solz[j] = FILL_INT16;
1886  for (j = 0; j < proj->buflines * proj->w; j++)
1887  proj->phi[j] = FILL_INT16;
1888  }
1889 
1890  if (process.smoke) {
1891  alloc_new_array((void *) &proj->smoke, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1892  for (j = 0; j < proj->buflines * proj->w; j++)
1893  proj->smoke[j] = FILL_INT16;
1894  }
1895 
1896  if (process.ndvi ||
1897  process.composite == MAXNDVI) {
1898  alloc_new_array((void *) &proj->ndvi, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1899  for (j = 0; j < proj->buflines * proj->w; j++)
1900  proj->ndvi[j] = FILL_INT16;
1901  }
1902 
1903  if (process.evi) {
1904  alloc_new_array((void *) &proj->evi, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1905  for (j = 0; j < proj->buflines * proj->w; j++)
1906  proj->evi[j] = FILL_INT16;
1907  }
1908 
1909  if (process.gemi) {
1910  alloc_new_array((void *) &proj->gemi, proj->buflines * proj->w, DFNT_INT16, process.maxmem);
1911  for (j = 0; j < proj->buflines * proj->w; j++)
1912  proj->gemi[j] = FILL_INT16;
1913  }
1914 
1915 }
1916 
1917 void init_output_file(OUTPUT *proj, INPUT input, PROCESSINFO process) {
1918  char path[MAXLENGTH];
1919  float64 attr;
1920 
1921  proj->startday = input.doy;
1922 
1923  /* Open output L3 HDF file */
1924  if (proj->filename[0] == 0) {
1925  sprintf(path, "%s", OUT_DIR);
1926  sprintf(proj->filename, "%s/%4.4d%3.3d.hdf", path, input.year, proj->startday);
1927  }
1928  printf("Creating %s ...\n", proj->filename);
1929  if ((proj->sd_id = SDstart(proj->filename, DFACC_CREATE)) == -1) {
1930  fprintf(stderr, "Cannot create file %s\n", proj->filename);
1931  exit(-1);
1932  }
1933 
1934  if (SDsetattr(proj->sd_id, "Projection", DFNT_CHAR8, (int) strlen(proj->projname) + 1, proj->projname) == -1)
1935  fprintf(stderr, "Cannot write attribute \"%s\" in file %s\n", "Projection", proj->filename);
1936 
1937  attr = proj->pixsz * R / RAD2DEG;
1938  if (SDsetattr(proj->sd_id, "Pixel size (meters)", DFNT_FLOAT64, 1, &attr) == -1)
1939  fprintf(stderr, "Cannot write attribute \"%s\" in file %s\n", "Pixel size", proj->filename);
1940 
1941  if (SDsetattr(proj->sd_id, "Maximum sensor zenith angle (degrees)", DFNT_FLOAT32, 1, &process.maxsenz) == -1)
1942  fprintf(stderr, "Cannot write attribute \"%s\" in file %s\n", "Maximum sensor zenith angle", proj->filename);
1943 
1944  if (SDsetattr(proj->sd_id, "Compositing criterion", DFNT_CHAR8,
1945  (int) strlen(process.compositename) + 1, process.compositename) == -1)
1946  fprintf(stderr, "Cannot write attribute \"%s\" in file %s\n", "Compositing criterion", proj->filename);
1947 
1948  if ((proj->p1 != 0) || (proj->l1 != 0)) {
1949  if (SDsetattr(proj->sd_id, "Subset Scan Range", DFNT_INT32, 2, &proj->l1) == -1)
1950  fprintf(stderr, "Cannot write attribute \"%s\" in file %s\n", "Subset Scan Range", proj->filename);
1951  if (SDsetattr(proj->sd_id, "Subset Pixel Range", DFNT_INT32, 2, &proj->p1) == -1)
1952  fprintf(stderr, "Cannot write attribute \"%s\" in file %s\n", "Subset Pixel Range", proj->filename);
1953  }
1954 
1955  write_cmd_line(proj->sd_id, proj->cmd_line);
1956 
1957 
1958  SDend(proj->sd_id);
1959 
1960 }
1961 
1962 int project_scan(int iscan, char **wl, INPUT *input, OUTPUT *proj, PROCESSINFO process) {
1963  int ipix, k, iline1, iline2, ipix1, ipix2;
1964  int iline_out, ipix_out;
1965  double xproj, yproj;
1966  double xarr[6], yarr[6];
1967  double xarrmin, xarrmax, yarrmin, yarrmax;
1968 
1969  if (FAST) {
1970  for (ipix = 0; ipix < input->Npixels; ipix += input->step) {
1971  if (input->solz[ipix] == FILL_INT16) continue;
1972  if (get_xyproj((double) input->lon[ipix], (double) input->lat[ipix], &xproj, &yproj,
1973  proj->projection, proj->lon_center) == -1) continue;
1974  iline_out = (proj->Nl - 1) / 2. - yproj / proj->pixsz + 0.5;
1975  ipix_out = (proj->Np - 1) / 2. + xproj / proj->pixsz + 0.5;
1976  if (iline_out >= proj->l1 &&
1977  iline_out <= proj->l2 &&
1978  ipix_out >= proj->p1 &&
1979  ipix_out <= proj->p2) {
1980  process_pixel(ipix, iline_out, ipix_out, wl, *input, proj, process);
1981  input->intersect = TRUE;
1982  }
1983  }
1984  }
1985 
1986  if (get_geolocation_scan(iscan + input->step - 2 * input->step *
1987  (iscan + input->step > input->Nscans - 1),
1988  input, process) == -1)
1989  return -1;
1990 
1991  if (FAST) return -1;
1992 
1993  memcpy(input->lon1, input->lon2, (input->Npixels + 1) * sizeof (double));
1994  memcpy(input->lat1, input->lat2, (input->Npixels + 1) * sizeof (double));
1995  resample_geo_scan(input->Npixels, input->lon, input->lat, input->lon2, input->lat2);
1996 
1997  if (input->old_iscan != iscan - input->step) {
1998  for (ipix = 0; ipix < input->Npixels + 1; ipix += input->step) {
1999  input->lontop[ipix] = -1.5 * input->lon1[ipix] + 0.5 * input->lon2[ipix];
2000  if (input->lon1[ipix] - input->lon2[ipix] > M_PI) input->lontop[ipix] += 0.5 * TWO_PI;
2001  if (input->lon2[ipix] - input->lon1[ipix] > M_PI) input->lontop[ipix] -= 1.5 * TWO_PI;
2002  input->lattop[ipix] = -1.5 * input->lat1[ipix] + 0.5 * input->lat2[ipix];
2003  if (get_xyproj(input->lontop[ipix], input->lattop[ipix], &input->xtop[ipix],
2004  &input->ytop[ipix], proj->projection, proj->lon_center) == -1) continue;
2005  }
2006  } else {
2007  memcpy(input->lontop, input->lonbot, (input->Npixels + 1) * sizeof (double));
2008  memcpy(input->lattop, input->latbot, (input->Npixels + 1) * sizeof (double));
2009  memcpy(input->xtop, input->xbot, (input->Npixels + 1) * sizeof (double));
2010  memcpy(input->ytop, input->ybot, (input->Npixels + 1) * sizeof (double));
2011  }
2012 
2013  if (iscan + input->step > input->Nscans - 1) {
2014  for (ipix = 0; ipix < input->Npixels + 1; ipix += input->step) {
2015  input->lonbot[ipix] = 1.5 * input->lon1[ipix] - 0.5 * input->lon2[ipix];
2016  if (input->lon1[ipix] - input->lon2[ipix] > M_PI) input->lonbot[ipix] -= 0.5 * TWO_PI;
2017  if (input->lon2[ipix] - input->lon1[ipix] > M_PI) input->lonbot[ipix] += 1.5 * TWO_PI;
2018  input->latbot[ipix] = 1.5 * input->lat1[ipix] - 0.5 * input->lat2[ipix];
2019  if (get_xyproj(input->lonbot[ipix], input->latbot[ipix], &input->xbot[ipix],
2020  &input->ybot[ipix], proj->projection, proj->lon_center) == -1) continue;
2021  }
2022  } else {
2023  for (ipix = 0; ipix < input->Npixels + 1; ipix += input->step) {
2024  input->lonbot[ipix] = 0.5 * input->lon1[ipix] + 0.5 * input->lon2[ipix];
2025  if (fabs(input->lon1[ipix] - input->lon2[ipix]) > M_PI)
2026  input->lonbot[ipix] += 0.5 * TWO_PI;
2027  input->latbot[ipix] = 0.5 * input->lat1[ipix] + 0.5 * input->lat2[ipix];
2028  if (get_xyproj(input->lonbot[ipix], input->latbot[ipix], &input->xbot[ipix],
2029  &input->ybot[ipix], proj->projection, proj->lon_center) == -1) continue;
2030  }
2031  }
2032 
2033  if (fabs(input->tilt[iscan] - input->tilt[iscan - 1 + 2 * (iscan == 0)]) > 0.1)
2034  for (ipix = 0; ipix < input->Npixels + 1; ipix += input->step)
2035  get_xyproj(input->lon1[ipix], input->lat1[ipix], &input->xtop[ipix],
2036  &input->ytop[ipix], proj->projection, proj->lon_center);
2037 
2038  if (fabs(input->tilt[iscan] - input->tilt[iscan + 1 - 2 * (iscan == input->Nscans - 1)]) > 0.1)
2039  for (ipix = 0; ipix < input->Npixels + 1; ipix += input->step)
2040  get_xyproj(input->lon1[ipix], input->lat1[ipix], &input->xbot[ipix],
2041  &input->ybot[ipix], proj->projection, proj->lon_center);
2042 
2043 
2044  /* Update output projection */
2045 
2046  for (ipix = 0; ipix < input->Npixels; ipix += input->step) {
2047  if (input->solz[ipix] == FILL_INT16) continue;
2048  if (proj->lon1 != UNDEF &&
2049  ((input->latbot[ipix] > proj->lat1 &&
2050  input->latbot[ipix + input->step] > proj->lat1) ||
2051  (input->lattop[ipix] < proj->lat2 &&
2052  input->lattop[ipix + input->step] < proj->lat2) ||
2053  (input->lontop[ipix] < proj->lon1 &&
2054  input->lontop[ipix + input->step] < proj->lon1 &&
2055  input->lonbot[ipix] < proj->lon1 &&
2056  input->lonbot[ipix + input->step] < proj->lon1) ||
2057  (input->lontop[ipix] > proj->lon2 &&
2058  input->lontop[ipix + input->step] > proj->lon2 &&
2059  input->lonbot[ipix] > proj->lon2 &&
2060  input->lonbot[ipix + input->step] > proj->lon2)))
2061  continue;
2062  xarr[0] = input->xtop[ipix];
2063  yarr[0] = input->ytop[ipix];
2064  xarr[1] = input->xtop[ipix + input->step];
2065  yarr[1] = input->ytop[ipix + input->step];
2066  xarr[2] = input->xbot[ipix + input->step];
2067  yarr[2] = input->ybot[ipix + input->step];
2068  xarr[3] = input->xbot[ipix];
2069  yarr[3] = input->ybot[ipix];
2070  xarr[4] = xarr[0];
2071  yarr[4] = yarr[0];
2072  xarr[5] = xarr[1];
2073  yarr[5] = yarr[1];
2074 
2075  xarrmin = DBL_MAX;
2076  xarrmax = -DBL_MAX;
2077  yarrmin = DBL_MAX;
2078  yarrmax = -DBL_MAX;
2079  for (k = 0; k < 4; k++) {
2080  if (xarr[k] < xarrmin) xarrmin = xarr[k];
2081  if (xarr[k] > xarrmax) xarrmax = xarr[k];
2082  if (yarr[k] < yarrmin) yarrmin = yarr[k];
2083  if (yarr[k] > yarrmax) yarrmax = yarr[k];
2084  }
2085 
2086  if (xarrmin > proj->xmax ||
2087  xarrmax < proj->xmin ||
2088  yarrmin > proj->ymax ||
2089  yarrmax < proj->ymin)
2090  continue;
2091 
2092  for (k = 0; k < 4; k++)
2093  if (CROSS_PRODUCT(xarr[k + 1], yarr[k + 1],
2094  xarr[k], yarr[k],
2095  xarr[k + 2], yarr[k + 2]) < 0 ||
2096  (xarr[k] * xarr[k + 1] < 0 &&
2097  fabs(xarr[k] - xarr[k + 1]) > 20 * proj->pixsz))
2098  break;
2099  if (k < 4) continue;
2100 
2101  iline1 = (proj->Nl - 1) / 2. - yarrmax / proj->pixsz + 0.5;
2102  iline2 = (proj->Nl - 1) / 2. - yarrmin / proj->pixsz + 0.5;
2103  ipix1 = (proj->Np - 1) / 2. + xarrmin / proj->pixsz + 0.5;
2104  ipix2 = (proj->Np - 1) / 2. + xarrmax / proj->pixsz + 0.5;
2105 
2106  for (iline_out = iline1; iline_out <= iline2; iline_out++)
2107  for (ipix_out = ipix1; ipix_out <= ipix2; ipix_out++) {
2108  if (iline_out < proj->l1 ||
2109  iline_out > proj->l2 ||
2110  ipix_out < proj->p1 ||
2111  ipix_out > proj->p2)
2112  continue;
2113  xproj = (ipix_out - (proj->Np - 1) / 2.) * proj->pixsz;
2114  yproj = ((proj->Nl - 1) / 2. - iline_out) * proj->pixsz;
2115  for (k = 0; k < 4; k++) {
2116  if (CROSS_PRODUCT(xarr[k + 1], yarr[k + 1],
2117  xproj, yproj,
2118  xarr[k], yarr[k]) > EPSILON ||
2119  CROSS_PRODUCT(xarr[k + 1], yarr[k + 1],
2120  xproj, yproj,
2121  xarr[k + 2], yarr[k + 2]) < -EPSILON) break;
2122  }
2123  if (k >= 4) {
2124  process_pixel(ipix, iline_out, ipix_out, wl, *input, proj, process);
2125  input->intersect = TRUE;
2126  }
2127  }
2128  }
2129 
2130  return 0;
2131 
2132 }
2133 
2134 typedef struct {
2135  char sds_name[H4_MAX_NC_NAME];
2136  int32 sds_id, num_type;
2137  int32 start[H4_MAX_VAR_DIMS], edges[H4_MAX_VAR_DIMS];
2138  float32 slope, intercept;
2139  void *data;
2140 } SDS;
2141 
2142 #define L2FIELDS (2 * NBANDS + 9)
2143 
2144 int read_l2(int iscan, INPUT *input, char **wl) {
2145  static int32 sds_index, rank, num_type, nattrs;
2146  static int32 start[H4_MAX_VAR_DIMS], edges[H4_MAX_VAR_DIMS], dim_sizes[H4_MAX_VAR_DIMS];
2147  static int16 *sola = NULL, *sena = NULL;
2148  float phi;
2149  int ib, j, k;
2150  static char first_time = TRUE;
2151  static SDS L2sds[L2FIELDS];
2152 
2153  if (first_time) {
2154 
2155  for (ib = 0; ib < NBANDS; ib++) {
2156  k = ib;
2157  sprintf(L2sds[k].sds_name, "rhos_%s", wl[ib]);
2158  L2sds[k].data = input->rho[ib];
2159  ;
2160  L2sds[k].num_type = DFNT_FLOAT32;
2161  }
2162  for (ib = 0; ib < NBANDS; ib++) {
2163  k = ib + NBANDS;
2164  sprintf(L2sds[k].sds_name, "rhot_%s", wl[ib]);
2165  L2sds[k].data = input->rhoTOA[ib];
2166  L2sds[k].num_type = DFNT_FLOAT32;
2167  }
2168 
2169  k++;
2170  strcpy(L2sds[k].sds_name, "l2_flags");
2171  L2sds[k].data = input->l2_flags;
2172  L2sds[k].num_type = DFNT_INT32;
2173 
2174  k++;
2175  strcpy(L2sds[k].sds_name, "longitude");
2176  L2sds[k].data = input->lon;
2177  L2sds[k].num_type = DFNT_FLOAT32;
2178  k++;
2179  strcpy(L2sds[k].sds_name, "latitude");
2180  L2sds[k].data = input->lat;
2181  L2sds[k].num_type = DFNT_FLOAT32;
2182  k++;
2183  strcpy(L2sds[k].sds_name, "ndvi");
2184  L2sds[k].data = input->ndvi;
2185  L2sds[k].num_type = DFNT_INT16;
2186  k++;
2187  strcpy(L2sds[k].sds_name, "evi");
2188  L2sds[k].data = input->evi;
2189  L2sds[k].num_type = DFNT_INT16;
2190  k++;
2191  strcpy(L2sds[k].sds_name, "solz");
2192  L2sds[k].data = input->solz;
2193  L2sds[k].num_type = DFNT_INT16;
2194  k++;
2195  strcpy(L2sds[k].sds_name, "senz");
2196  L2sds[k].data = input->senz;
2197  L2sds[k].num_type = DFNT_INT16;
2198  k++;
2199  strcpy(L2sds[k].sds_name, "sola");
2200  L2sds[k].data = NULL;
2201  L2sds[k].num_type = DFNT_INT16;
2202  k++;
2203  strcpy(L2sds[k].sds_name, "sena");
2204  L2sds[k].data = NULL;
2205  L2sds[k].num_type = DFNT_INT16;
2206 
2207  first_time = FALSE;
2208 
2209  }
2210 
2211  if (input->newfile) {
2212 
2213  if (input->phi) {
2214  sena = (int16 *) realloc(sena, input->Npixels * sizeof (int16));
2215  sola = (int16 *) realloc(sola, input->Npixels * sizeof (int16));
2216  L2sds[L2FIELDS - 2].data = sola;
2217  L2sds[L2FIELDS - 1].data = sena;
2218  }
2219  for (k = 0; k < L2FIELDS; k++) {
2220  if (L2sds[k].data) {
2221  if ((sds_index = SDnametoindex(input->sd_id, L2sds[k].sds_name)) == -1) {
2222  fprintf(stderr, "Cannot locate SDS %s %d\n", L2sds[k].sds_name, k);
2223  return -1;
2224  }
2225  if ((L2sds[k].sds_id = SDselect(input->sd_id, sds_index)) == -1) {
2226  fprintf(stderr, "Cannot select SDS %s\n", L2sds[k].sds_name);
2227  return -1;
2228  }
2229  if (SDgetinfo(L2sds[k].sds_id, L2sds[k].sds_name, &rank, dim_sizes, &num_type, &nattrs) < 0) {
2230  fprintf(stderr, "can't get info for SDS %s\n", L2sds[k].sds_name);
2231  return -1;
2232  }
2233  if (strcmp(L2sds[k].sds_name, "latitude") == 0 ||
2234  strcmp(L2sds[k].sds_name, "longitude") == 0) {
2235  L2sds[k].slope = 1 / RAD2DEG;
2236  L2sds[k].intercept = 0;
2237  } else if (strcmp(L2sds[k].sds_name, "l2_flags") == 0) {
2238  L2sds[k].slope = 0;
2239  L2sds[k].intercept = 0;
2240  } else if (L2sds[k].num_type == DFNT_FLOAT32) {
2241  L2sds[k].slope = 1.;
2242  L2sds[k].intercept = 0;
2243  } else if (read_attr(L2sds[k].sds_id, "slope", &L2sds[k].slope) < 0 ||
2244  read_attr(L2sds[k].sds_id, "intercept", &L2sds[k].intercept) < 0) {
2245  fprintf(stderr, "can't read scaling parameters of SDS %s\n", L2sds[k].sds_name);
2246  return -1;
2247  }
2248  if (num_type != L2sds[k].num_type) {
2249  fprintf(stderr, "SDS %s has the wrong data type (%d)\n", L2sds[k].sds_name, num_type);
2250  return -1;
2251  }
2252  if (dim_sizes[0] != input->Nscans ||
2253  dim_sizes[1] != input->Npixels) {
2254  fprintf(stderr, "SDS %s has the wrong dimensions (%dx%d)\n", L2sds[k].sds_name, input->Npixels, input->Nscans);
2255  return -1;
2256  }
2257  }
2258  }
2259 
2260  }
2261 
2262  start[0] = iscan;
2263  start[1] = 0;
2264  edges[0] = 1;
2265  edges[1] = input->Npixels;
2266  for (k = 0; k < L2FIELDS; k++) {
2267  if (L2sds[k].data) {
2268  if (SDreaddata(L2sds[k].sds_id, start, NULL, edges, L2sds[k].data) < 0) {
2269  fprintf(stderr, "can't read SDS %s\n", L2sds[k].sds_name);
2270  return -1;
2271  }
2272  if (k < NBANDS) /* Surface Reflectance: int16 -> float32 */
2273  for (j = input->Npixels - 1; j >= 0; j--)
2274  ((float32 *) L2sds[k].data)[j] = L2sds[k].slope * ((float32 *) L2sds[k].data)[j] + L2sds[k].intercept;
2275  else if (k > 2 * NBANDS && k < L2FIELDS - 6) /* TOA refl, longitude and latitude: float32 -> float32 */
2276  for (j = 0; j < input->Npixels; j++)
2277  ((float32 *) L2sds[k].data)[j] = L2sds[k].slope * ((float32 *) L2sds[k].data)[j] + L2sds[k].intercept;
2278  }
2279  }
2280 
2281  if (input->ndvi) /* NDVI, EVI, angles: int16 -> int16 */
2282  for (j = 0; j < input->Npixels; j++)
2283  input->ndvi[j] = floor((L2sds[L2FIELDS - 6].slope * input->ndvi[j] + L2sds[L2FIELDS - 6].intercept) / sf_vi + 0.5);
2284  if (input->evi)
2285  for (j = 0; j < input->Npixels; j++)
2286  input->evi[j] = floor((L2sds[L2FIELDS - 5].slope * input->evi[j] + L2sds[L2FIELDS - 5].intercept) / sf_vi + 0.5);
2287  if (input->solz)
2288  for (j = 0; j < input->Npixels; j++)
2289  input->solz[j] = (L2sds[L2FIELDS - 4].slope * input->solz[j] + L2sds[L2FIELDS - 4].intercept) / sf_angle + 0.5;
2290  if (input->senz)
2291  for (j = 0; j < input->Npixels; j++)
2292  input->senz[j] = (L2sds[L2FIELDS - 3].slope * input->senz[j] + L2sds[L2FIELDS - 3].intercept) / sf_angle + 0.5;
2293  if (input->phi)
2294  for (j = 0; j < input->Npixels; j++) {
2295  phi = L2sds[L2FIELDS - 2].slope * sola[j] + L2sds[L2FIELDS - 2].intercept
2296  - L2sds[L2FIELDS - 1].slope * sena[j] - L2sds[L2FIELDS - 1].intercept;
2297  if (phi > +180) phi -= 360;
2298  if (phi < -180) phi += 360;
2299  input->phi[j] = floor(phi / sf_angle + 0.5);
2300  }
2301  return 0;
2302 
2303 }
int hamfor(double lon, double lat, double *x, double *y)
Definition: proj_hamfor.c:60
double xmin
Definition: make_L3_v1.1.c:115
int robfor(double lon, double lat, double *x, double *y)
Definition: proj_robfor.c:108
double * lon2
Definition: make_L3_v1.1.c:141
@ MAXNDVI
Definition: make_L3_v1.1.c:45
@ Nfiles
Definition: hybrid.c:20
double ymin
Definition: make_L3_v1.1.c:115
#define STEP
Definition: make_L3_v1.1.c:68
data_t q
Definition: decode_rs.h:74
integer, parameter int16
Definition: cubeio.f90:3
int16 GEMI(float red, float nir)
void check_process_parameters(PROCESSINFO *process, char **wl)
#define MIN(x, y)
Definition: rice.h:169
int r
Definition: decode_rs.h:73
#define CONVOL_MAX
Definition: make_L3_v1.1.c:101
void red(int iz1, int iz2, int jz1, int jz2, int jm1, int jm2, int jmf, int ic1, int jc1, int jcf, int kc, float ***c, float **s)
int robinv(double x, double y, double *lon, double *lat)
Definition: proj_robinv.c:111
#define BADORBLIST
Definition: path.h:4
#define ATMCOR
Definition: make_L3_v1.1.c:64
#define OMF2
Definition: make_L3_v1.1.c:95
@ PLATECARREE
Definition: make_L3_v1.1.c:41
int j
Definition: decode_rs.h:73
void make_psurf(signed short *psurf, char *file, int Nlines, int Npixels, float psurf0)
Definition: make_psurf.c:6
int16 SMOKE(float **rho, int ipix)
#define L(lambda, T)
Definition: PreprocessP.h:185
float32 intercept
@ BAND5
Definition: make_L3_v1.1.c:53
int16 * senz
Definition: make_L3_v1.1.c:108
#define CROSS_PRODUCT(x1, y1, x2, y2, x3, y3)
Definition: make_L3_v1.1.c:103
@ MOLLWEIDE
Definition: make_L3_v1.1.c:41
void parse_command_line(int argc, char **argv, OUTPUT *proj, PROCESSINFO *process)
#define FAST
Definition: make_L3_v1.1.c:63
int read_sds_block(int32 sds_id, int iline, int buflines, void *dataP)
Definition: read_write.c:19
#define DEFPIXSZ
Definition: make_L3_v1.1.c:69
#define VERSION
Definition: make_L3_v1.1.c:105
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
int16 * evi
Definition: make_L3_v1.1.c:108
float * tilt
Definition: make_L3_v1.1.c:135
int year
Definition: make_L3_v1.1.c:139
#define FALSE
Definition: rice.h:164
@ BAND1
Definition: make_L3_v1.1.c:53
int main(int argc, char **argv)
Definition: make_L3_v1.1.c:175
#define L2FIELDS
#define NULL
Definition: decode_rs.h:63
int Nfiles
Definition: make_L3_v1.1.c:119
MOD_PR01 Production producing one five minute granule of output data in each run It can be configured to produce as many as three five minute granules per run Each execution with one construction record and one date file for each dataset In normal these are created by which splits them out of the hour datasets For LANCE they are created by which merges all session MODIS L0 datasets overlapping the requested time and extracts from the merged data those packets which fall within that time period Each scan of data is stored in the L1A granule that covers the start time of that scan
Definition: MOD_PR01_pr.txt:19
#define RAD2DEG
Definition: make_L3_v1.1.c:93
void resample_geo_scan(int Npixels, float *inlon, float *inlat, double *outlon, double *outlat)
Definition: make_L3_v1.1.c:874
char * cmd_line
Definition: make_L3_v1.1.c:111
int molwforint(double r, double center_long, double false_east, double false_north)
Definition: proj_molwfor.c:37
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 band
#define FILL_INT16
Definition: make_L3_v1.1.c:72
int robinvint(double r, double center_long, double false_east, double false_north)
Definition: proj_robinv.c:46
char daystr[14]
Definition: make_L3_v1.1.c:130
#define TAURAY
Definition: make_L3_v1.1.c:81
#define MAXLENGTH
Definition: make_L3_v1.1.c:97
#define EPSILON
Definition: make_L3_v1.1.c:98
#define MINREFL
Definition: make_L3_v1.1.c:73
int32 sds_id
float64 sf_refl
Definition: make_L3_v1.1.c:173
float h[MODELMAX]
Definition: atrem_corl1.h:131
float32 * pos
Definition: l1_czcs_hdf.c:35
@ BAND3
Definition: make_L3_v1.1.c:53
#define TRUE
Definition: rice.h:165
float * lat
char sds_name[H4_MAX_NC_NAME]
HDF4 data type of the output SDS Default is DFNT_FLOAT32 Common types used DFNT_INT32
int write_sds(int32 sd_id, char *SDSname, void *data, int startline, int buflines, int Nl, int Np, int32 numtype, void *fillvalue, float64 scale_factor, char verbose, char gzip)
Definition: read_write.c:75
double pixsz
Definition: make_L3_v1.1.c:117
char domain
Definition: make_L3_v1.1.c:116
@ GLOBAL
Definition: make_L3_v1.1.c:49
int update_file(char *filename, char **wl, int offset, int startline, int buflines, OUTPUT *proj, char verbose, char gzip)
Definition: make_L3_v1.1.c:975
int16 * ndvi
Definition: make_L3_v1.1.c:108
int get_attributes(int32 sd_id, char *prefix, int *Nlines, int *Npixels, int *startpix, int *subsampling, int *year, int *doy, int *orbno, char *title)
Definition: get_attributes.c:4
FILE * input_fp
Definition: make_L3_v1.1.c:120
#define BLUEBAND
Definition: make_L3_v1.1.c:76
int get_navig_sds_line(int32 sd_id, int iline, float32 *sen_mat, float32 *orb_vec, float32 *sun_ref, float32 *scan_ell, char newfile)
Definition: get_navig_sds.c:4
int alloc_new_array(void **data, int size, int32 num_type, int maxmem)
int16 * solz
Definition: make_L3_v1.1.c:108
#define TILTSCANS
Definition: make_L3_v1.1.c:66
double lon2
Definition: make_L3_v1.1.c:114
void set_buffer_parameters(OUTPUT *proj, PROCESSINFO *process)
#define WVAPOR
Definition: make_L3_v1.1.c:83
instr * input
char path[MAXLENGTH]
Definition: make_L3_v1.1.c:129
character(len=1000) if
Definition: names.f90:13
int16 * gemi
Definition: make_L3_v1.1.c:108
int startline
Definition: make_L3_v1.1.c:119
void init_output_file(OUTPUT *proj, INPUT input, PROCESSINFO process)
int setlinebuf(FILE *stream)
int32 sd_id
Definition: make_L3_v1.1.c:118
int read_l2(int iscan, INPUT *input, char **wl)
char dayend[14]
Definition: make_L3_v1.1.c:131
int robforint(double r, double center_long, double false_east, double false_north)
Definition: proj_robfor.c:43
float64 sf_vi
Definition: make_L3_v1.1.c:173
@ NUMBEROFPROJECTIONS
Definition: make_L3_v1.1.c:41
data_t tmp
Definition: decode_rs.h:74
void set_proj_parameters(OUTPUT *proj)
int subsampling
Definition: make_L3_v1.1.c:137
double ymax
Definition: make_L3_v1.1.c:115
#define UPDATE
Definition: make_L3_v1.1.c:65
#define N_CONVOL
Definition: make_L3_v1.1.c:100
#define NpixelsLAC
Definition: make_L3_v1.1.c:94
float32 slope
void update_pixel(INPUT input, int16 ndvi, int ipix, int idx_out, int ifno, OUTPUT *proj)
string path
Definition: color_dtdb.py:221
int16 * solz
Definition: make_L3_v1.1.c:136
#define TWO_PI
Definition: make_L3_v1.1.c:92
char * strdup(const char *)
#define IRRAD
Definition: make_L3_v1.1.c:80
int32 * l2_flags
Definition: make_L3_v1.1.c:145
char band[NBANDS]
Definition: make_L3_v1.1.c:126
FILE * fileuse_fp
Definition: make_L3_v1.1.c:121
void process_pixel(int ipix, int iline_out, int ipix_out, char **wl, INPUT input, OUTPUT *proj, PROCESSINFO process)
float32 slope[]
Definition: l2lists.h:30
@ LASTIN
Definition: make_L3_v1.1.c:45
void memcpy_block(int offset2, int offset1, int Nrows, int Ncols, OUTPUT *proj)
int read_attr(int32 sd_id, char *attrname, void *attr)
Definition: read_write.c:61
#define WAVELENGTH
Definition: make_L3_v1.1.c:79
void write_cmd_line(int32 sd_id, char *cmd_line)
@ USERBOX
Definition: make_L3_v1.1.c:49
int get_xyproj(double lon, double lat, double *xproj, double *yproj, int projtype, double lon_center)
char projname[MAXLENGTH]
Definition: make_L3_v1.1.c:109
#define M_PI
Definition: pml_iop.h:15
@ BAND4
Definition: make_L3_v1.1.c:53
@ BAND6
Definition: make_L3_v1.1.c:53
#define UNDEF
Definition: make_L3_v1.1.c:70
@ L1A
Definition: make_L3_v1.1.c:57
int32 open_sds(int32 sd_id, char *sds_name)
Definition: read_write.c:6
@ BAND2
Definition: make_L3_v1.1.c:53
char filename[FILENAME_MAX]
Definition: atrem_corl1.h:122
float32 intercept[]
Definition: l2lists.h:44
integer, parameter double
int haminv(double x, double y, double *lon, double *lat)
Definition: proj_haminv.c:59
@ L2
Definition: make_L3_v1.1.c:57
HDF4 data type of the output SDS Default is DFNT_FLOAT32 Common types used DFNT_INT16
#define CALIBTABLE
Definition: path.h:3
int32 sd_id
Definition: make_L3_v1.1.c:140
no change in intended resolving MODur00064 Corrected handling of bad ephemeris attitude data
Definition: HISTORY.txt:356
int get_geolocation_scan(int iscan, INPUT *input, PROCESSINFO process)
Definition: make_L3_v1.1.c:748
void * data
Definition: hybrid.c:31
#define FOUR_THIRDS
Definition: make_L3_v1.1.c:91
char * filelist
Definition: make_L3_v1.1.c:111
float64 sf_angle
Definition: make_L3_v1.1.c:173
double lon1
Definition: make_L3_v1.1.c:114
int16 * phi
Definition: make_L3_v1.1.c:108
data_t b[NROOTS+1]
Definition: decode_rs.h:77
int16 EVI(float blue, float red, float nir)
@ INPUT
Definition: hybrid.c:24
level
Definition: mapgen.py:186
int16 * refl[NBANDS]
Definition: make_L3_v1.1.c:108
#define MAXSENZDEF
Definition: make_L3_v1.1.c:74
subroutine csalbr(xtau, xalb)
Definition: 6sm1.f:3618
int molwinvint(double r, double center_long, double false_east, double false_north)
Definition: proj_molwinv.c:35
#define UMASK
Definition: make_L3_v1.1.c:75
#define NLPSURF
Definition: make_L3_v1.1.c:87
char filename[MAXLENGTH]
Definition: make_L3_v1.1.c:110
#define OUT_DIR
Definition: path.h:1
#define fabs(a)
Definition: misc.h:93
int32_t iscan
Extra metadata that will be written to the HDF4 file l2prod rank
int get_lonlat(double xproj, double yproj, double *lon, double *lat, int projtype, double lon_center)
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)
@ HAMMER
Definition: make_L3_v1.1.c:41
unsigned char read_mask(double lon, double lat, int16 close)
Definition: read_mask.c:14
int molwfor(double lon, double lat, double *x, double *y)
Definition: proj_molwfor.c:53
float * lon
char compositename[MAXLENGTH]
Definition: make_L3_v1.1.c:129
@ BAND8
Definition: make_L3_v1.1.c:53
#define DBL_MAX
Definition: make_L3_v1.1.c:60
#define TWO_THIRDS
Definition: make_L3_v1.1.c:90
for(i=0;i< NROOTS;i++) s[i]
Definition: decode_rs.h:85
l2prod offset
HDF4 data type of the output SDS Default is DFNT_FLOAT32 Common types used DFNT_FLOAT32
#define R
Definition: make_L3_v1.1.c:96
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
int16 * reflTOA[NBANDS]
Definition: make_L3_v1.1.c:108
#define MAXMEMDEF
Definition: make_L3_v1.1.c:67
double xmax
Definition: make_L3_v1.1.c:115
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] angles
Definition: HISTORY.txt:366
#define DEFPROJECTION
Definition: make_L3_v1.1.c:62
#define ANC_DIR
Definition: path.h:2
int molwinv(double x, double y, double *lon, double *lat)
Definition: proj_molwinv.c:52
int compute_l1b(int32 sd_id, int iline, float *l1b_data[], char *calibtable, char newfile)
Definition: compute_l1b.c:9
#define NPPSURF
Definition: make_L3_v1.1.c:88
int startday
Definition: make_L3_v1.1.c:119
double lat1
Definition: make_L3_v1.1.c:114
int i
Definition: decode_rs.h:71
int buflines
Definition: make_L3_v1.1.c:112
Definition: hybrid.c:27
msiBandIdx val
Definition: l1c_msi.cpp:34
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
float64 sf_smoke
Definition: make_L3_v1.1.c:173
int16 * ifno
Definition: make_L3_v1.1.c:108
PGE01 indicating that PGE02 PGE01 V6 for and PGE01 V2 for MOD03 were used to produce the granule By convention adopted in all MODIS Terra PGE02 code versions are The fourth digit of the PGE02 version denotes the LUT version used to produce the granule The source of the metadata environment variable ProcessingCenter was changed from a QA LUT value to the Process Configuration A sign used in error in the second order term was changed to a
Definition: HISTORY.txt:424
#define OZONE
Definition: make_L3_v1.1.c:82
@ BAND7
Definition: make_L3_v1.1.c:53
@ ROBINSON
Definition: make_L3_v1.1.c:41
int32 num_type
Definition: hybrid.c:29
#define PSURF
Definition: make_L3_v1.1.c:86
int verbose
Definition: fmt_check.c:6
@ NBANDS
Definition: make_L3_v1.1.c:53
int k
Definition: decode_rs.h:73
int haminvint(double r, double center_long, double false_east, double false_north)
Definition: proj_haminv.c:40
double lat2
Definition: make_L3_v1.1.c:114
int32 read_sds(l1info_struct l1info, char *arr_name, int32 *exp_ntyp, void *array)
Definition: read_sds.c:5
subroutine chand(xphi, xmuv, xmus, xtau, xrray)
Definition: 6sm1.f:3548
@ MINSENZ
Definition: make_L3_v1.1.c:45
int project_scan(int iscan, char **wl, INPUT *input, OUTPUT *proj, PROCESSINFO process)
int16 * smoke
Definition: make_L3_v1.1.c:108
int16 NDVI(float red, float nir)
@ SINUSOIDAL
Definition: make_L3_v1.1.c:41
int read_file_block(char *filename, char **wl, int offset, int startline, int buflines, OUTPUT *proj)
Definition: make_L3_v1.1.c:902
@ MINBLUE
Definition: make_L3_v1.1.c:45
double * ytop
Definition: make_L3_v1.1.c:143
double * lontop
Definition: make_L3_v1.1.c:142
int bufstep
Definition: make_L3_v1.1.c:112
int hamforint(double r, double center_long, double false_east, double false_north)
Definition: proj_hamfor.c:41
#define OXYGEN
Definition: make_L3_v1.1.c:84
double lon_center
Definition: make_L3_v1.1.c:117
#define IFOV
Definition: make_L3_v1.1.c:746
void init_output_projection(OUTPUT *proj, PROCESSINFO process)
char projection
Definition: make_L3_v1.1.c:109
int count
Definition: decode_rs.h:79