OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
RsViirs.cpp
Go to the documentation of this file.
1 /*******************************************************************************
2  *
3  * NAME: RsViirs.cpp
4  *
5  * DESCRIPTION: Running the program on a granule will modify the data
6  * in-place ..
7  *
8  * Based on a technique described in:
9  * Remote Sensing 2016, 8(1), 79
10  * Article Titled "Improved VIIRS and MODIS SST Imagery"
11  * By Irina Gladkova, Alexander Ignatov, Fazlul Shahriar, Yury Kihai,
12  * Don Hillger and Boris Petrenko
13  *
14  *
15  * Created on: February 17, 2019
16  * Author: Sam Anderson
17  *
18  *******************************************************************************/
19 
20 
21 #include <netcdf>
22 #include <DDOptions.h>
23 #include "RsViirs.h"
24 
25 using namespace std;
26 using namespace netCDF;
27 using namespace netCDF::exceptions;
28 
29 // column break points
30 short SBP[NBREAKS+1] = {
31  0, 5, 87, 170, 358, 567, 720, 850, 997, 1120, 1275, 1600
32 };
33 
34 // relative row that the pixel comes from
36  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
37  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
38  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
39  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
40  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
41  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
42  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
43  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
44  {+8, +8, +8, 0, 0, 0, 0, 0, 0, 0, 0},
45  {+8, -1, -1, +7, 0, 0, 0, 0, 0, 0, 0},
46  {-2, +7, +7, -1, +6, 0, 0, 0, 0, 0, 0},
47  {+7, -2, -2, +6, -1, +5, 0, 0, 0, 0, 0},
48  {-3, +6, +6, -2, +5, -1, +4, 0, 0, 0, 0},
49  {+6, -3, -3, +5, -2, +4, -1, +3, 0, 0, 0},
50  {-4, +5, +5, -3, +4, -2, +3, -1, +2, 0, 0},
51  {+5, -4, -4, +4, -3, +3, -2, +2, -1, +1, 0},
52 };
53 
54 // relative row that the pixel comes from
56  {-5, +4, +4, -4, +3, -3, +2, -2, +1, -1, 0},
57  {+4, -5, -5, +3, -4, +2, -3, +1, -2, 0, 0},
58  {-6, +3, +3, -5, +2, -4, +1, -3, 0, 0, 0},
59  {+3, -6, -6, +2, -5, +1, -4, 0, 0, 0, 0},
60  {-7, +2, +2, -6, +1, -5, 0, 0, 0, 0, 0},
61  {+11, -7, -7, +1, -6, 0, 0, 0, 0, 0, 0},
62  {+1, +1, +1, -7, 0, 0, 0, 0, 0, 0, 0},
63  {-9, +9, -8, 0, 0, 0, 0, 0, 0, 0, 0},
64  {+9, -9, +8, 0, 0, 0, 0, 0, 0, 0, 0},
65  {-1, -1, -1, +7, 0, 0, 0, 0, 0, 0, 0},
66  {-11, +7, +7, -1, +6, 0, 0, 0, 0, 0, 0},
67  {+7, -2, -2, +6, -1, +5, 0, 0, 0, 0, 0},
68  {-3, +6, +6, -2, +5, -1, +4, 0, 0, 0, 0},
69  {+6, -3, -3, +5, -2, +4, -1, +3, 0, 0, 0},
70  {-4, +5, +5, -3, +4, -2, +3, -1, +2, 0, 0},
71  {+5, -4, -4, +4, -3, +3, -2, +2, -1, +1, 0},
72 };
73 
74 // relative row that the pixel comes from
76  {-5, +4, +4, -4, +3, -3, +2, -2, +1, -1, 0},
77  {+4, -5, -5, +3, -4, +2, -3, +1, -2, 0, 0},
78  {-6, +3, +3, -5, +2, -4, +1, -3, 0, 0, 0},
79  {+3, -6, -6, +2, -5, +1, -4, 0, 0, 0, 0},
80  {-7, +2, +2, -6, +1, -5, 0, 0, 0, 0, 0},
81  {+2, -7, -7, +1, -6, 0, 0, 0, 0, 0, 0},
82  {-8, +1, +1, -7, 0, 0, 0, 0, 0, 0, 0},
83  {-8, -8, -8, 0, 0, 0, 0, 0, 0, 0, 0},
84  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
85  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
86  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
87  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
88  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
89  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
90  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
91  {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
92 };
93 
94 
95 /**************************************************************************
96  * NAME: RsViirs()
97  *
98  * DESCRIPTION: Class Constructor
99  *
100  *************************************************************************/
101 
103 }
104 
105 RsViirs::RsViirs( int lines, int pixels ) {
106  lines_ = lines;
107  pixels_ = pixels;
108 }
109 
110 /**************************************************************************
111  * NAME: ~RsViirs()
112  *
113  * DESCRIPTION: Class Destructor
114  *
115  *************************************************************************/
116 
118 }
119 
120 /**************************************************************************
121  * NAME: process()
122  *
123  * DESCRIPTION: Processes data and produces output values for granule
124  *
125  *************************************************************************/
126 
128  int status = DD_SUCCESS;
129 
130  string filepath_l1b = get_option(INPUT_L1B);
131  string filepath_geo = get_option(INPUT_GEO);
132 
133  if (filepath_l1b.empty() || filepath_geo.empty()) {
134  return DD_FAIL;
135  }
136 
137  NcFile* nc_io;
138  try {
139  nc_io = new NcFile(filepath_l1b, NcFile::read);
140  } catch (NcException& e) {
141  e.what();
142  cerr << "RsInput:: Failure opening VIIRS L1B input file: "
143  + filepath_l1b << endl;
144  return DD_FAIL;
145  }
146  NcDim line_dim = nc_io->getDim("number_of_lines");
147  NcDim pixel_dim = nc_io->getDim("number_of_pixels");
148 
149  lines_ = line_dim.getSize();
150  pixels_ = pixel_dim.getSize();
151  delete nc_io;
152 
153  status = generate_sort_index();
154 
155  sio* sa = new sio;
156  memset( sa, 0, sizeof(sio));
157  uio* ua = new uio;
158  memset( ua, 0, sizeof(uio));
159  usio* usa = new usio;
160  memset( usa, 0, sizeof(usio));
161  fio* fa = new fio;
162  memset( fa, 0, sizeof(fio));
163 
164  try {
165  nc_io = new NcFile(filepath_geo, NcFile::write);
166  } catch (NcException& e) {
167  e.what();
168  cerr << "RsViirs:: Failure opening VIIRS Geolocation input file: "
169  + filepath_geo << endl;
170  return DD_FAIL;
171  }
172 
173  NcGroup nc_group = nc_io->getGroup("geolocation_data");
174  NcVar nc_var = nc_group.getVar("quality_flag");
175  nc_var.getVar(&sa->in[0][0]);
176  resort(sa);
177  nc_var.putVar(&sa->out[0][0]);
178  nc_var = nc_group.getVar("height");
179  nc_var.getVar(&sa->in[0][0]);
180  resort(sa);
181  nc_var.putVar(&sa->out[0][0]);
182  nc_var = nc_group.getVar("range");
183  nc_var.getVar(&sa->in[0][0]);
184  resort(sa);
185  nc_var.putVar(&sa->out[0][0]);
186  nc_var = nc_group.getVar("sensor_azimuth");
187  nc_var.getVar(&sa->in[0][0]);
188  resort(sa);
189  nc_var.putVar(&sa->out[0][0]);
190  nc_var = nc_group.getVar("sensor_zenith");
191  nc_var.getVar(&sa->in[0][0]);
192  resort(sa);
193  nc_var.putVar(&sa->out[0][0]);
194  nc_var = nc_group.getVar("solar_azimuth");
195  nc_var.getVar(&sa->in[0][0]);
196  resort(sa);
197  nc_var.putVar(&sa->out[0][0]);
198  nc_var = nc_group.getVar("solar_zenith");
199  nc_var.getVar(&sa->in[0][0]);
200  resort(sa);
201  nc_var.putVar(&sa->out[0][0]);
202  nc_var = nc_group.getVar("latitude");
203  nc_var.getVar(&fa->in[0][0]);
204  resort(fa);
205  nc_var.putVar(&fa->out[0][0]);
206  nc_var = nc_group.getVar("longitude");
207  nc_var.getVar(fa);
208  resort(fa);
209  nc_var.putVar(&fa->out[0][0]);
210  delete nc_io;
211 
212  try {
213  nc_io = new NcFile(filepath_l1b, NcFile::write);
214  } catch (NcException& e) {
215  e.what();
216  cerr << "RsViirs:: Failure opening VIIRS L1B input file: "
217  + filepath_l1b << endl;
218  return DD_FAIL;
219  }
220  nc_group = nc_io->getGroup("observation_data");
221  for (int ib = 0; ib < NMBANDS; ib++) {
222  if (ib==12) continue;
223  char str[20] = "";
224  sprintf(str, "M%02d", ib + 1);
225  nc_var = nc_group.getVar(string(str));
226  if (ib==12) {
227  nc_var.getVar(&ua->in[0][0]);
228  resort(ua);
229  fill_the_fills(ua);
230  nc_var.putVar(&ua->out[0][0]);
231  } else {
232  nc_var.getVar(&usa->in[0][0]);
233  resort(usa);
234  fill_the_fills(usa);
235  nc_var.putVar(&usa->out[0][0]);
236  }
237  sprintf(str, "M%02d_quality_flags", ib + 1);
238  nc_var = nc_group.getVar(string(str));
239  nc_var.getVar(&sa->in[0][0]);
240  resort(sa);
241  nc_var.putVar(&sa->out[0][0]);
242  }
243  delete nc_io;
244 
245  delete sa;
246  delete ua;
247  delete usa;
248  delete fa;
249 
250  return status;
251 }
252 
253 /**************************************************************************
254  * NAME: generate_sort_index()
255  *
256  * DESCRIPTION: Generate a image of latitude sorting indices.
257  *
258  *************************************************************************/
259 
261 
262  int ND = NDETECTORS;
263  for (int i = 0; i < NBREAKS; i++) {
264  for (int ix = SBP[i]; ix < SBP[i+1]; ix++) {
265  for (int iy = 0; iy < ND; iy++) {
266  si_[iy][ix] = iy + SORT_FIRST[iy][i];
267  si_[iy][pixels_-ix-1] = iy + SORT_FIRST[iy][i];
268  }
269  for (int iy = ND; iy < lines_ - ND; iy++) {
270  si_[iy][ix] = iy + SORT_MID[iy % ND][i];
271  si_[iy][pixels_-ix-1] = iy + SORT_MID[iy % ND][i];
272  }
273  for (int iy = lines_ - ND; iy < lines_; iy++) {
274  si_[iy][ix] = iy + SORT_LAST[iy % ND][i];
275  si_[iy][pixels_-ix-1] = iy + SORT_LAST[iy % ND][i];
276  }
277  }
278  }
279  return DD_SUCCESS;
280 }
281 
282 /**************************************************************************
283  * NAME: resort()
284  *
285  * DESCRIPTION: Returns the sorted image.
286  * Si is the image of sort indices.
287  *
288  *************************************************************************/
289 
291  for (int iy = 0; iy < lines_; iy++) {
292  for (int ix = 0; ix < pixels_; ix++) {
293  io->out[iy][ix] = io->in[si_[iy][ix]][ix];
294  }
295  }
296  return DD_SUCCESS;
297 }
298 
300  for (int iy = 0; iy < lines_; iy++) {
301  for (int ix = 0; ix < pixels_; ix++) {
302  io->out[iy][ix] = io->in[si_[iy][ix]][ix];
303  }
304  }
305  return DD_SUCCESS;
306 }
307 
309  for (int iy = 0; iy < lines_; iy++) {
310  for (int ix = 0; ix < pixels_; ix++) {
311  io->out[iy][ix] = io->in[si_[iy][ix]][ix];
312  }
313  }
314  return DD_SUCCESS;
315 }
316 
318  for (int iy = 0; iy < lines_; iy++) {
319  for (int ix = 0; ix < pixels_; ix++) {
320  io->out[iy][ix] = io->in[si_[iy][ix]][ix];
321  }
322  }
323  return DD_SUCCESS;
324 }
325 
326 /**************************************************************************
327  * NAME: fill_the_fills()
328  *
329  * DESCRIPTION: Replace fill values with average of above and below.
330  *
331  *************************************************************************/
332 
334  const unsigned short filltest = SOUB_UINT16_FILL;
335 
336  for (int iy = 1; iy < lines_ - 1; iy++) {
337  for (int ix = 0; ix < pixels_; ix++) {
338  if (io->out[iy][ix] >= filltest) {
339  if (io->out[iy - 1][ix] < filltest &&
340  io->out[iy + 1][ix] < filltest) {
341  io->out[iy][ix] = io->out[iy - 1][ix] / 2 +
342  io->out[iy + 1][ix] / 2;
343  } else if (io->out[iy - 1][ix] < filltest) {
344  io->out[iy][ix] = io->out[iy - 1][ix];
345  } else if (io->out[iy + 1][ix] < filltest) {
346  io->out[iy][ix] = io->out[iy + 1][ix];
347  }
348  }
349  }
350  }
351  for (int ix = 0; ix < pixels_; ix++) {
352  if (io->out[0][ix] >= filltest) {
353  if (io->out[1][ix] < filltest) {
354  io->out[0][ix] = io->out[1][ix];
355  }
356  }
357  }
358  for (int ix = 0; ix < pixels_; ix++) {
359  if (io->out[lines_-1][ix] >= filltest) {
360  if (io->out[lines_-2][ix] < filltest) {
361  io->out[lines_-1][ix] = io->out[lines_-2][ix];
362  }
363  }
364  }
365 
366  return DD_SUCCESS;
367 }
368 
370  const unsigned int filltest = 327680;
371 
372  for (int iy = 1; iy < lines_ - 1; iy++) {
373  for (int ix = 0; ix < pixels_; ix++) {
374  if (io->out[iy][ix] >= filltest) {
375  if (io->out[iy - 1][ix] < filltest &&
376  io->out[iy + 1][ix] < filltest) {
377  io->out[iy][ix] = io->out[iy - 1][ix] / 2 +
378  io->out[iy + 1][ix] / 2;
379  } else if (io->out[iy - 1][ix] < filltest) {
380  io->out[iy][ix] = io->out[iy - 1][ix];
381  } else if (io->out[iy + 1][ix] < filltest) {
382  io->out[iy][ix] = io->out[iy + 1][ix];
383  }
384  }
385  }
386  }
387  for (int ix = 0; ix < pixels_; ix++) {
388  if (io->out[0][ix] >= filltest) {
389  if (io->out[1][ix] < filltest) {
390  io->out[0][ix] = io->out[1][ix];
391  }
392  }
393  }
394  for (int ix = 0; ix < pixels_; ix++) {
395  if (io->out[lines_-1][ix] >= filltest) {
396  if (io->out[lines_-2][ix] < filltest) {
397  io->out[lines_-1][ix] = io->out[lines_-2][ix];
398  }
399  }
400  }
401 
402  return DD_SUCCESS;
403 }
404 
406  const unsigned short filltest = SOUB_UINT16_FILL;
407 
408  for (int iy = 1; iy < lines_ - 1; iy++) {
409  for (int ix = 0; ix < pixels_; ix++) {
410  if (io->out[iy][ix] < filltest) {
411  if (io->out[iy - 1][ix] > filltest &&
412  io->out[iy + 1][ix] > filltest) {
413  io->out[iy][ix] = io->out[iy - 1][ix] / 2 +
414  io->out[iy + 1][ix] / 2;
415  } else if (io->out[iy - 1][ix] > filltest) {
416  io->out[iy][ix] = io->out[iy - 1][ix];
417  } else if (io->out[iy + 1][ix] > filltest) {
418  io->out[iy][ix] = io->out[iy + 1][ix];
419  }
420  }
421  }
422  }
423  for (int ix = 0; ix < pixels_; ix++) {
424  if (io->out[0][ix] < filltest) {
425  if (io->out[1][ix] > filltest) {
426  io->out[0][ix] = io->out[1][ix];
427  }
428  }
429  }
430  for (int ix = 0; ix < pixels_; ix++) {
431  if (io->out[lines_-1][ix] < filltest) {
432  if (io->out[lines_-2][ix] > filltest) {
433  io->out[lines_-1][ix] = io->out[lines_-2][ix];
434  }
435  }
436  }
437 
438  return DD_SUCCESS;
439 }
440 
short SBP[NBREAKS+1]
Definition: RsViirs.cpp:30
RsViirs()
Definition: RsViirs.cpp:102
int resort(sio *io)
Definition: RsViirs.cpp:290
int status
Definition: l1_czcs_hdf.c:32
short SORT_LAST[NDETECTORS][NBREAKS]
Definition: RsViirs.cpp:75
short out[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:77
unsigned int in[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:84
string get_option(const string &name)
Definition: DDOptions.cpp:211
~RsViirs()
Definition: RsViirs.cpp:117
const unsigned short SOUB_UINT16_FILL
Definition: RsViirs.h:51
short SORT_FIRST[NDETECTORS][NBREAKS]
Definition: RsViirs.cpp:35
const string INPUT_L1B
Definition: DDOptions.cpp:40
unsigned short in[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:80
short SORT_MID[NDETECTORS][NBREAKS]
Definition: RsViirs.cpp:55
float in[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:72
const int NDETECTORS
Definition: RsViirs.h:65
float out[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:73
const string INPUT_GEO
Definition: DDOptions.cpp:41
int generate_sort_index()
Definition: RsViirs.cpp:260
Definition: RsViirs.h:79
short in[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:76
const int NMBANDS
Definition: RsViirs.h:62
int process()
Definition: RsViirs.cpp:127
const char * str
Definition: l1c_msi.cpp:35
const int NBREAKS
Definition: RsViirs.h:66
unsigned int out[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:85
Definition: RsViirs.h:71
Definition: RsViirs.h:75
Definition: RsViirs.h:83
int fill_the_fills(usio *io)
Definition: RsViirs.cpp:333
int i
Definition: decode_rs.h:71
unsigned short out[NSCANS *NDETECTORS][NPIXELS]
Definition: RsViirs.h:81