A sphere with three ocean waves in differing shades of blue. Next to the sphere, there is the bolded text "Ocean Color". Under that, there are two acronyms, separated by a vertical pipe: OB.DAAC (Ocean Biology Distributed Active Archive Center) and OBPG (Ocean Biology Processing Group). Lorem ipsum dolor sit amet
NASA Logo
Ocean Color Science Software

ocssw V2022
polcor.c
Go to the documentation of this file.
1 /* ========================================================================== */
2 /* module polcor.c - functions to read and apply polarization correction */
3 /* */
4 /* Written By: B. Franz, SAIC, November 2005. */
5 /* W. Robinson, SAIC, 29 May 2009 generalize for use with VIIRS */
6 /* ========================================================================== */
7 
8 #include "l12_proto.h"
9 #include "polcor_hawkeye.h"
10 #include "polcor_oci.h"
11 #include <xcal.h>
12 
13 #define NMIR 2 /* always 2 mirror sides */
14 #define NDET 40 /* big enough for largest detector array */
15 #define NANG 15 /* big enough for largest set of table AOI or scan angles */
16 
17 void polcor(l1str *l1rec, int32_t ip) {
18  static int firstCall = 1;
19 
20  typedef float m_array[NMIR][NANG][NDET];
21  static m_array *m12, *m13;
22  static float *detfac;
23  static float maxpix;
24  static int32 nang;
25  static float *ang = NULL, margin_s;
26  static int32_t nir_s;
27  static int32_t nir_l;
28 
29  static double *dm12;
30  static double *dm13;
31  static int *have_xcal;
32 
33  static float *pxwt; /* weights for each pixel in the line and */
34  static char *pxangix; /* index of low angle storage */
35 
36  float *polcor = l1rec->polcor;
37  float *delta_polcor=NULL;
38  float *dpol = l1rec->dpol;
39 
40  int32_t pixnum = l1rec->pixnum[ip];
41  int32_t detnum = l1rec->detnum;
42  int32_t mside = l1rec->mside;
43  int32_t nbands = l1rec->l1file->nbands;
44 
45  int32_t sensorID = l1rec->l1file->sensorID;
46 
47  float m1 [NANG], deriv_m1[NANG];
48  float alpha;
49  int32_t ib, ipb;
50  int32_t idet, iang;
51  float wt, zang;
52  float L_x, L_qp, L_up;
53 
54  int iagsm[] = {0, 640, 1376, 3152, 4928, 5664, 6304};
55  int iagpx[] = {0, 640, 1008, 1600, 2192, 2560, 3200};
56  int iseq, ix1, ix2, ipx, irng, step[] = {1, 2, 3, 3, 2, 1};
57  char ix, attrname[50];
58  double sind = 0.0003104;
59  float rad_2_deg = 180. / acos(-1);
60  int get_wt(float, float *, int, float *, char *);
61 
62  if(l1rec->uncertainty) {
63  delta_polcor=l1rec->uncertainty->derv_pol;
64  }
65 
66  if (mside < 0) mside = 0;
67  if (mside > 1) mside = 1;
68 
69  for (ib = 0; ib < nbands; ib++) {
70  l1rec->polcor[ip * nbands + ib] = 1.0;
71  }
72 
73  if (input->pol_opt == 0 || input->pol_opt > 99)
74  return;
75 
76  if (sensorID == OCI) {
77  polcor_oci(l1rec, ip);
78  return;
79  }
80 
81  if(sensorID == HAWKEYE) {
82  polcor_hawkeye(l1rec, ip);
83  return;
84  }
85 
86  if (firstCall) {
87 
88  char name [H4_MAX_NC_NAME] = "";
89  char sdsname[H4_MAX_NC_NAME] = "";
90  char file [FILENAME_MAX] = "";
91  int32 sd_id;
92  int32 sds_id;
93  int32 rank;
94  int32 nt;
95  int32 dims[H4_MAX_VAR_DIMS];
96  int32 nattrs;
97  int32 start[3] = {0, 0, 0};
98  int32 end [3] = {1, 1, 1};
99  int32 status;
100 
101  int32 ndets;
102  int32 nmside;
103  int32_t im, ia, id, ixb;
104 
105 
106  if ((detfac = (float *) calloc(l1rec->l1file->nbands, sizeof (float))) == NULL) {
107  printf("-E- : Error allocating memory to dray_for_i\n");
108  exit(FATAL_ERROR);
109  }
110 
111  if ((have_xcal = (int *) calloc(l1rec->l1file->nbands, sizeof (int))) == NULL) {
112  printf("-E- : Error allocating memory to dray_for_i\n");
113  exit(FATAL_ERROR);
114  }
115 
116  for (ixb = 0; ixb < l1_input->xcal_nwave; ixb++) {
117  if ((l1_input->xcal_opt[ixb] & XCALPOL) != 0) {
118  if ((ib = bindex_get(l1_input->xcal_wave[ixb])) < 0) {
119  printf("-E- %sline %d: xcal wavelength %f does not match sensor\n",
120  __FILE__, __LINE__, l1_input->xcal_wave[ixb]);
121  exit(1);
122  };
123  have_xcal[ib] = 1;
124  }
125  }
126 
127  if ((m12 = (m_array *) calloc(l1rec->l1file->nbands, sizeof (m_array))) == NULL) {
128  printf("-E- : Error allocating memory to dray_for_i\n");
129  exit(FATAL_ERROR);
130  }
131  if ((m13 = (m_array *) calloc(l1rec->l1file->nbands, sizeof (m_array))) == NULL) {
132  printf("-E- : Error allocating memory to dray_for_i\n");
133  exit(FATAL_ERROR);
134  }
135  /* load polfile for each sensor wavelength */
136  for (ib = 0; ib < nbands; ib++) {
137 
138  sprintf(file, "%s%s%d%s", input->polfile, "_", l1rec->l1file->iwave[ib], ".hdf");
139 
140  printf("Loading polarization file %s\n", file);
141 
142  /* Open the file */
143  sd_id = SDstart(file, DFACC_RDONLY);
144  if (sd_id == FAIL) {
145  fprintf(stderr, "-E- %s line %d: SDstart(%s, %d) failed.\n",
146  __FILE__, __LINE__, file, DFACC_RDONLY);
147  exit(1);
148  }
149 
150  // number of detectors in polcor file from global attr
151  status = SDreadattr(sd_id, SDfindattr(sd_id, "Number of Detectors"), &ndets);
152  if (status != 0) {
153  printf("-E- %s Line %d: Error reading global attribute %s from %s.\n",
154  __FILE__, __LINE__, "Number of Detectors", file);
155  exit(1);
156  }
157 
158  status = SDreadattr(sd_id, SDfindattr(sd_id, "Number of Mirror Sides"), &nmside);
159  if (status != 0) {
160  printf("-E- %s Line %d: Error reading global attribute %s from %s.\n",
161  __FILE__, __LINE__, "Number of Mirror Sides", file);
162  exit(1);
163  }
164 
165  if (ib == 0) {
166 
167  // decide on which angle information needs to be read
168  // TODO: rework this to either use sensor name via sensorId2InstrumentName function
169  // or better yet, redo the polcor files as netCDF and be consistent so this isn't required
170  if ((sensorID == VIIRSN) || (sensorID == VIIRSJ1) || (sensorID == VIIRSJ2)) {
171  strcpy(attrname, "Number of Scan angles");
172  strcpy(sdsname, "scanangle");
173  } else {
174  strcpy(attrname, "Number of AOIs");
175  strcpy(sdsname, "AOI");
176  }
177 
178  // MODIS - attrname = Number of AOIs
179  status = SDreadattr(sd_id, SDfindattr(sd_id, attrname), &nang);
180  if (status != 0) {
181  printf(
182  "-E- %s Line %d: Error reading global attribute %s from %s.\n",
183  __FILE__, __LINE__, attrname, file);
184  exit(1);
185  }
186 
187  if ((ang = (float *) malloc(nang * sizeof (float))) == NULL) {
188  printf("-E- %s Line %d: Error allocating memory for ang.\n",
189  __FILE__, __LINE__);
190  exit(1);
191  }
192  // MODIS - sdsname = AOI
193  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
194  status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs); // grabs AOI items
195  status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) ang);
196  if (status != 0) {
197  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
198  __FILE__, __LINE__, sdsname, file);
199  exit(1);
200  } else {
201  status = SDendaccess(sds_id);
202  }
203  }
204 
205 
206  strcpy(sdsname, "m12");
207  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
208  for (im = 0; im < nmside; im++) for (ia = 0; ia < nang; ia++) for (id = 0; id < ndets; id++) {
209  start[0] = im;
210  start[1] = ia;
211  start[2] = id;
212  // MODIS - read in m12 by [side][AOI][detector] save for each band ib
213  status = SDreaddata(sds_id, start, NULL, end, (VOIDP) & m12[ib][im][ia][id]);
214  if (status != 0) {
215  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
216  __FILE__, __LINE__, sdsname, file);
217  exit(1);
218  }
219  }
220 
221  status = SDendaccess(sds_id);
222 
223  strcpy(sdsname, "m13");
224  sds_id = SDselect(sd_id, SDnametoindex(sd_id, sdsname));
225  for (im = 0; im < nmside; im++) for (ia = 0; ia < nang; ia++) for (id = 0; id < ndets; id++) {
226  start[0] = im;
227  start[1] = ia;
228  start[2] = id;
229  status = SDreaddata(sds_id, start, NULL, end, (VOIDP) & m13[ib][im][ia][id]);
230  if (status != 0) {
231  printf("-E- %s Line %d: Error reading SDS %s from %s.\n",
232  __FILE__, __LINE__, sdsname, file);
233  exit(1);
234  }
235  }
236  status = SDendaccess(sds_id);
237  status = SDend(sd_id);
238 
239  /* we need to translate from actual detectors to replicated detectors or */
240  /* aggregated detectors, to handle the HMODIS resolution switching. Get */
241  /* the conversion factor by comparing sensor ndets to table ndets. */
242  if ((sensorID == VIIRSN) || (sensorID == VIIRSJ1) || (sensorID == VIIRSJ2))
243  detfac[ib] = 1.;
244  else
245  detfac[ib] = l1rec->l1file->ndets * 1.0 / ndets; /* 0.25 to 4 */
246  }
247 
248 
249  /* get the max # pixels */
250  if (sensorID == MODISA || sensorID == MODIST) {
251  switch (l1_input->resolution) {
252  case 250:
253  maxpix = 5416.0;
254  break;
255  case 500:
256  maxpix = 2708.0;
257  break;
258  default:
259  maxpix = 1354.0;
260  break;
261  }
262  } else if ((sensorID == VIIRSN) || (sensorID == VIIRSJ1) || (sensorID == VIIRSJ2)) {
263  margin_s = l1rec->margin_s;
264  if (l1rec->scn_fmt == 0)
265  maxpix = 3200.;
266  else
267  maxpix = 6304 + margin_s * 2;
268  } else {
269  printf("-E- %s line %d: Mirror geometry unknown for %s\n",
270  __FILE__, __LINE__, sensorId2SensorDir(sensorID));
271  exit(1);
272  }
273 
274  // precompute the weights and angle (AOI on mirror or scan angle)
275  // dimension index for each pixel in the scan
276 
277  if ((pxwt = (float *) malloc(maxpix * sizeof ( float))) == NULL) {
278  printf("-E- %s Line %d: Error allocating memory for weights.\n",
279  __FILE__, __LINE__);
280  exit(1);
281  }
282  if ((pxangix = (char *) malloc(maxpix * sizeof ( char))) == NULL) {
283  printf("-E- %s Line %d: Error allocating memory for angle index.\n",
284  __FILE__, __LINE__);
285  exit(1);
286  }
287 
288  // depending on instrument, get the angle and find the weight and index
289  // into pol table for all pixels
290 
291  if ((sensorID == VIIRSN) || (sensorID == VIIRSJ1) || (sensorID == VIIRSJ2)) {
292  if (l1rec->scn_fmt == 0) /* aggregated */ {
293  for (irng = 0; irng < 6; irng++) {
294  ix1 = irng;
295  ix2 = irng + 1;
296  for (ipx = iagpx[ix1]; ipx < iagpx[ix2]; ipx++) {
297  iseq = ipx - iagpx[ix1];
298  zang = iagsm[ix1] + ((2 * iseq + 1) * step[irng] - 1.) / 2.;
299  zang = (zang - 3151.5) * sind * rad_2_deg;
300 
301  // use zang in same way as below - I hate to repeat this though
302  get_wt(zang, ang, nang, &wt, &ix);
303  pxwt[ipx] = wt;
304  pxangix[ipx] = ix;
305  }
306  }
307  } else /* unaggregated */ {
308  for (ipx = 0; ipx < maxpix; ipx++) {
309  zang = ((float) ipx - 3151.5 - margin_s) * sind * rad_2_deg;
310  get_wt(zang, ang, nang, &wt, &ix);
311  pxwt[ipx] = wt;
312  pxangix[ipx] = ix;
313  }
314  }
315 
316  } else {
317 
318  for (ipx = 0; ipx < maxpix; ipx++) {
319  zang = 10.5 + ipx / maxpix * (65.5 - 10.5);
320  get_wt(zang, ang, nang, &wt, &ix);
321  pxwt[ipx] = wt;
322  pxangix[ipx] = ix;
323  }
324  }
325 
326  nir_s = bindex_get(input->aer_wave_short);
327  nir_l = bindex_get(input->aer_wave_long);
328 
329  firstCall = 0;
330  }
331 
332  /* apply sensitivites to radiances to get corrections */
333 
334  for (ib = 0; ib < nbands; ib++) {
335 
336  ipb = ip * nbands + ib;
337  idet = (int32_t) rint(detnum / detfac[ib]);
338  alpha = l1rec->alpha[ip] / RADEG;
339  L_x = l1rec->Lt[ipb] / l1rec->tg_sol[ipb] / l1rec->tg_sen[ipb];
340  if(l1rec->uncertainty) {
341  delta_polcor[ipb]= 1/ l1rec->tg_sol[ipb] / l1rec->tg_sen[ipb];
342  }
343  if (L_x > 0.0) {
344  L_qp = l1rec->L_q[ipb] * cos(2 * alpha) + l1rec->L_u[ipb] * sin(2 * alpha);
345  L_up = l1rec->L_u[ipb] * cos(2 * alpha) - l1rec->L_q[ipb] * sin(2 * alpha);
346  if (have_xcal[ib]) {
347  dm12 = get_xcal(l1rec, XM12, l1rec->l1file->iwave[ib]);
348  dm13 = get_xcal(l1rec, XM13, l1rec->l1file->iwave[ib]);
349  polcor[ipb] = 1.0 / (1.0 - dm12[ip] * L_qp / L_x - dm13[ip] * L_up / L_x);
350  } else {
351  for (iang = pxangix[pixnum]; iang <= (pxangix[pixnum] + 1); iang++) {
352  m1[iang] = 1.0 / (1.0
353  - m12[ib][mside][iang][idet] * L_qp / L_x
354  - m13[ib][mside][iang][idet] * L_up / L_x);
355 
356  if(l1rec->uncertainty) {
357  deriv_m1[iang]=-m1[iang]*m1[iang]*(m12[ib][mside][iang][idet] * L_qp + m13[ib][mside][iang][idet] * L_up)/(L_x*L_x);
358  deriv_m1[iang]*=delta_polcor[ipb];
359  }
360  }
361  polcor[ipb] = m1[(int) pxangix[pixnum]] * (1. - pxwt[pixnum]) + m1[(int) pxangix[pixnum] + 1] * pxwt[pixnum];
362  if(l1rec->uncertainty) {
363  delta_polcor[ipb]=deriv_m1[(int) pxangix[pixnum]] * (1. - pxwt[pixnum]) + deriv_m1[(int) pxangix[pixnum] + 1] * pxwt[pixnum];
364  }
365  }
366  dpol [ipb] = sqrt(pow(l1rec->L_q[ipb], 2.0) + pow(l1rec->L_u[ipb], 2.0)) / L_x;
367 
368  /* quick-fix to polcor of aerosol bands only */
369  if (input->pol_opt == 6 && ib != nir_l && ib != nir_s)
370  polcor[ipb] = 1.0;
371 
372  } else {
373  polcor[ipb] = 1.0;
374  dpol [ipb] = 0.0;
375  }
376  }
377 }
378 
379 int get_wt(float zang, float *ang, int nang, float *wt, char *ix)
380 /*******************************************************************
381 
382  get_wt
383 
384  purpose: derive a weight and index into angle array
385 
386  Returns type: int - 0 if all is OK
387 
388  Parameters: (in calling order)
389  Type Name I/O Description
390  ---- ---- --- -----------
391  float zang I angle to interpolate to
392  float * ang I array of angles in
393  interpolation table
394  int nang I # of angle tie points
395  float * wt O derived weight
396  char * ix O index of first angle to use in
397  weighting
398 
399  *******************************************************************/ {
400  int iang, ix2;
401 
402  if (zang <= ang[0])
403  *ix = 0;
404  else if (zang >= ang[nang - 1])
405  *ix = nang - 2;
406  else {
407  for (iang = 1; iang < nang; iang++) {
408  if (zang < ang[iang]) {
409  *ix = iang - 1;
410  break;
411  }
412  }
413  }
414  ix2 = *ix + 1;
415  *wt = (zang - ang[(int) *ix]) / (ang[ix2] - ang[(int) *ix]);
416  return 0;
417 }
sensorId2SensorDir
const char * sensorId2SensorDir(int sensorId)
Definition: sensorInfo.c:315
XCALPOL
#define XCALPOL
Definition: l1.h:87
status
int status
Definition: l1_czcs_hdf.c:32
polcor
void polcor(l1str *l1rec, int32_t ip)
Definition: polcor.c:17
NANG
#define NANG
Definition: polcor.c:15
OCI
#define OCI
Definition: sensorDefs.h:42
FAIL
#define FAIL
Definition: ObpgReadGrid.h:18
NULL
#define NULL
Definition: decode_rs.h:63
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
VIIRSN
#define VIIRSN
Definition: sensorDefs.h:23
polcor_oci.h
MODIST
#define MODIST
Definition: sensorDefs.h:18
file
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 file
Definition: HISTORY.txt:413
im
int32_t im
Definition: atrem_corl1.h:161
XM13
#define XM13
Definition: xcal.h:12
input
instr * input
Definition: main_l2binmatch.cpp:27
l12_proto.h
get_wt
int get_wt(float zang, float *ang, int nang, float *wt, char *ix)
Definition: polcor.c:379
bindex_get
int bindex_get(int32_t wave)
Definition: windex.c:45
polcor_oci
void polcor_oci(l1str *l1rec, int32_t currPix)
Definition: polcor_oci.cpp:314
mside
int32_t mside
Definition: l1_hmodis_hdf.c:752
HAWKEYE
#define HAWKEYE
Definition: sensorDefs.h:39
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
FATAL_ERROR
#define FATAL_ERROR
Definition: swl0_parms.h:5
polcor_hawkeye.h
polcor_hawkeye
void polcor_hawkeye(l1str *l1rec, int32_t ip)
Definition: polcor_hawkeye.cpp:24
RADEG
#define RADEG
Definition: czcs_ctl_pt.c:5
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
XM12
#define XM12
Definition: xcal.h:11
get_xcal
double * get_xcal(l1str *l1rec, int type, int wave)
Definition: xcal.c:33
nbands
int32_t nbands
Definition: atrem_corl1v3.h:190
rank
Extra metadata that will be written to the HDF4 file l2prod rank
Definition: HOWTO_Add_a_product.txt:80
NMIR
#define NMIR
Definition: polcor.c:13
VIIRSJ2
#define VIIRSJ2
Definition: sensorDefs.h:44
color_dtdb.alpha
alpha
Definition: color_dtdb.py:352
NDET
#define NDET
Definition: polcor.c:14
check_products.end
end
Definition: check_products.py:189
fix_mac_rpath.id
string id
Definition: fix_mac_rpath.py:28
rint
double rint(double)
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
sensorID
int32_t sensorID[MAXNFILES]
Definition: l2bin.cpp:91
MODISA
#define MODISA
Definition: sensorDefs.h:19
VIIRSJ1
#define VIIRSJ1
Definition: sensorDefs.h:37


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

How to Use Earthdata Search

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster