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
cloud_flag.c
Go to the documentation of this file.
1 #include "l1.h"
2 #include <nc4utils.h>
3 #include <hdf4utils.h>
4 
5 int get_sdps_cld_mask( l1str *l1rec, int32_t ip, char *cld_category ){
6 /*
7  get_sdps_cld_mask - get the cloud mask catergoy, and the cloud/no cloud
8  from the cloud mask file made for SDPS
9 
10  Returns: int the cloud value: 0 no cloud, 1 cloud
11  l1str *l1rec: supply the cloud mask file name, controls
12  int32_t ip: pixel #
13  char *cld_category: the cloud category: 0 - confident_cloudy,
14  1 - probablyclear, 2 - probablyclear, 3 - probablyclear
15  this will return the above category even for the new format, although it returns
16  only 0 for cloud, 3 for clear
17 
18  WDR, 29 Mar 2021
19  WDR, 9 Nov 2023, also be able to read the new format with SDS cloud_flag
20  in group geophysical_data
21  */
22  static size_t start[] = {0,0}, count[] = {0,0}, fil_np;
23  static int firstCall = 1;
24  static int ncid, d_id, dim_id;
25  static signed char* msk_lcl;
26  static int lastScan = -1;
27  static int new_cl_msk_fmt;
28  int grp_id;
29 
30  if (firstCall) {
31  char sdsname[HDF4_UTILS_MAX_NAME] = "";
32  char file [FILENAME_MAX] = "";
33 
34  if (!l1_input->cld_msk_file[0]) {
35  printf("-E- %s line %d: no cloud flag file provided\n",
36  __FILE__, __LINE__);
37  exit(1);
38  }
39 
40  strcpy(file, l1_input->cld_msk_file);
41  printf("Loading cloud mask flag file %s\n", file);
42 
43  /* Open the file */
44  if (nc_open(file, 0, &ncid) != NC_NOERR) {
45  fprintf(stderr,
46  "-E- %s %d: file: %s is not netcdf, not acceptable cloud maskfile\n",
47  __FILE__, __LINE__, file);
48  exit(1);
49  }
50  // get the # pixels
51  if (nc_inq_dimid(ncid, "pixels_per_line", &dim_id ) != NC_NOERR) {
52  fprintf(stderr,
53  "-E- %s %d: file: %s Could not read the pixels_per_line id\n",
54  __FILE__, __LINE__, file);
55  exit(1);
56  }
57 
58  if (nc_inq_dimlen(ncid, dim_id, &fil_np ) != NC_NOERR) {
59  fprintf(stderr,
60  "-E- %s %d: file: %s Could not read the pixels_per_line\n",
61  __FILE__, __LINE__, file);
62  exit(1);
63  }
64  // try reading the new group/sds name. if group fails,
65  // assume old format
66  if( nc_inq_grp_ncid( ncid, "geophysical_data", &grp_id ) != NC_NOERR) {
67  new_cl_msk_fmt = 0;
68  fprintf(stderr,
69  "-I- %s %d: file: %s Cloud mask read switch to old format\n",
70  __FILE__, __LINE__, file);
71  strcpy(sdsname, "CF_CATEGORY");
72  if( nc_inq_varid(ncid, sdsname, &d_id ) != NC_NOERR) {
73  fprintf(stderr,
74  "-E- %s %d: file: %s Could not find SDS CF_CATEGORY\n",
75  __FILE__, __LINE__, file);
76  exit(1);
77  }
78  } else {
79  // for new mask format
80  ncid = grp_id;
81  new_cl_msk_fmt = 1;
82  if(l1_input->cloud_mask_opt == 0) {
83  strcpy(sdsname, "cloud_flag");
84  } else if(l1_input->cloud_mask_opt == 1) {
85  strcpy(sdsname, "cloud_flag_dilated");
86  } else {
87  fprintf(stderr,
88  "-E- %s:%d - cloud_mask_opt=%d Illegal option\n",
89  __FILE__, __LINE__, l1_input->cloud_mask_opt);
90  exit(1);
91  }
92  if( nc_inq_varid(grp_id, sdsname, &d_id ) != NC_NOERR) {
93  fprintf(stderr,
94  "-E- %s %d: file: %s Could not find SDS CF_CATEGORY\n",
95  __FILE__, __LINE__, file);
96  exit(1);
97  }
98  }
99  // get the range to process and set it
100  start[1] = l1_input->spixl - 1; // spixl, epixl are 1-origin!
101  count[1] = l1_input->epixl - l1_input->spixl + 1;
102  count[0] = 1;
103  // start[0] will be the line to read
104  if( ( msk_lcl = (signed char *)
105  malloc( count[1] * sizeof( signed char *) ) ) == NULL ) {
106  fprintf(stderr,
107  "-E- %s %d: file: %s Could not allocate cloud mask storage\n",
108  __FILE__, __LINE__, file);
109  exit(1);
110  }
111  firstCall = 0;
112  // end setup
113  }
114 
115  if (lastScan != l1rec->iscan) {
116  start[0] = l1rec->iscan;
117  if( nc_get_vara_schar( ncid, d_id, start, count, msk_lcl ) != NC_NOERR ) {
118  fprintf(stderr,
119  "-E- %s %d: Could not read the cloud mask file line\n",
120  __FILE__, __LINE__);
121  exit(1);
122  }
123  lastScan = l1rec->iscan;
124  }
125 
126  if( new_cl_msk_fmt == 0 ) {
127  *cld_category = msk_lcl[ip];
128  } else {
129  if( msk_lcl[ip] == BAD_BYTE )
130  *cld_category = BAD_BYTE;
131  else if( msk_lcl[ip] == 0 )
132  *cld_category = 3;
133  else *cld_category = 0;
134  }
135 
136  if (*cld_category == BAD_BYTE )
137  return (BAD_BYTE);
138  else if (*cld_category < 2 )
139  return (1);
140  else
141  return (0);
142 }
143 
144 char get_cloudmask_meris(l1str *l1rec, int32_t ip) {
145  // Cloud Masking for MERIS and OLCI
146 
147  static int ib443, ib490, ib560, ib620, ib665, ib681, ib709, ib754, ib865, ib885, firstCall = 1;
148  int ipb;
149  float *rhos = l1rec->rhos, *cloud_albedo = l1rec->cloud_albedo;
150  float ftemp, cldtmp;
151  char flagcld;
152  static productInfo_t* rhos_product_info;
153 
154  if (firstCall == 1) {
155  ib443 = bindex_get(443);
156  ib490 = bindex_get(490);
157  ib560 = bindex_get(560);
158  ib620 = bindex_get(620);
159  ib665 = bindex_get(665);
160  ib681 = bindex_get(681);
161  ib709 = bindex_get(709);
162  ib754 = bindex_get(754);
163  ib865 = bindex_get(865);
164  ib885 = bindex_get(885);
165 
166  if (ib443 < 0 || ib490 < 0 || ib560 < 0 || ib620 < 0 || ib665 < 0 || ib681 < 0 || ib709 < 0 || ib754 < 0 || ib865 < 0 || ib885 < 0) {
167  printf("get_habs_cldmask: incompatible sensor wavelengths for this algorithm\n");
168  exit(EXIT_FAILURE);
169  }
170 
171  if (rhos_product_info == NULL) {
172  rhos_product_info = allocateProductInfo();
173  findProductInfo("rhos", l1rec->l1file->sensorID, rhos_product_info);
174  }
175 
176  firstCall = 0;
177  }
178  flagcld = 0;
179  ipb = l1rec->l1file->nbands*ip;
180 
181  if (rhos[ipb + ib443] >= 0.0 &&
182  rhos[ipb + ib560] >= 0.0 &&
183  rhos[ipb + ib620] >= 0.0 &&
184  rhos[ipb + ib665] >= 0.0 &&
185  rhos[ipb + ib681] >= 0.0 &&
186  rhos[ipb + ib709] >= 0.0 &&
187  rhos[ipb + ib754] >= 0.0) {
188  // turbidity signal in water
189  if ((rhos[ipb + ib443] <= rhos_product_info->validMax) &&
190  (rhos[ipb + ib665] <= rhos_product_info->validMax) &&
191  (rhos[ipb + ib620] <= rhos_product_info->validMax) &&
192  (rhos[ipb + ib681] <= rhos_product_info->validMax) &&
193  (rhos[ipb + ib754] <= rhos_product_info->validMax)) {
194 
195  ftemp = (rhos[ipb + ib620] + rhos[ipb + ib665] + rhos[ipb + ib681]) - 3 * rhos[ipb + ib443] - \
196  (rhos[ipb + ib754] - rhos[ipb + ib443]) / (754 - 443)*(620 + 665 + 681 - 3 * 443);
197  } else {
198  ftemp = 0;
199  }
200  //switch to rhos_865 where cldalb fails
201  if (cloud_albedo[ip] >= 0.0) {
202  cldtmp = cloud_albedo[ip];
203  } else {
204  cldtmp = rhos[ipb + ib865];
205  }
206 
207  //remove turbidity signal from cloud albedo
208  if (ftemp > 0) cldtmp = cldtmp - 3 * ftemp;
209 
210  if (cldtmp > 0.08) {
211  flagcld = 1;
212  }
213 
214  // to deal with scum look at relative of NIR and blue for lower albedos
215  if ((rhos[ipb + ib754] <= rhos_product_info->validMax) &&
216  (rhos[ipb + ib709] <= rhos_product_info->validMax) &&
217  (rhos[ipb + ib681] <= rhos_product_info->validMax) &&
218  (rhos[ipb + ib560] <= rhos_product_info->validMax)) {
219 
220  if ((rhos[ipb + ib443] <= rhos_product_info->validMax) &&
221  (rhos[ipb + ib490] <= rhos_product_info->validMax)) {
222  if ((rhos[ipb + ib754] + rhos[ipb + ib709]) > (rhos[ipb + ib443] + rhos[ipb + ib490]) && cldtmp < 0.1 && (rhos[ipb + ib560] > rhos[ipb + ib681])) flagcld = 0;
223  }
224  if (rhos[ipb + ib665] <= rhos_product_info->validMax) {
225  if (((rhos[ipb + ib754] + rhos[ipb + ib709]) - (rhos[ipb + ib665] + rhos[ipb + ib681])) > 0.01 && cldtmp < 0.15 && (rhos[ipb + ib560] > rhos[ipb + ib681])) flagcld = 0;
226  if ((((rhos[ipb + ib754] + rhos[ipb + ib709]) - (rhos[ipb + ib665] + rhos[ipb + ib681])) / cldtmp) > 0.1 && (rhos[ipb + ib560] > rhos[ipb + ib681])) flagcld = 0;
227  }
228  }
229  if (rhos[ipb + ib665] <= rhos_product_info->validMax) {
230  if ((rhos[ipb + ib665] > 0.1) && (cldtmp > 0.15)) {
231  flagcld = 1;
232  }
233  }
234  if ((rhos[ipb + ib865] >= rhos_product_info->validMin) &&
235  (rhos[ipb + ib865] <= rhos_product_info->validMax)) {
236  // flag high glint
237  if (rhos[ipb + ib865] - cloud_albedo[ip] > 0.25) {
238  flagcld = 1;
239  }
240  } else {
241  flagcld = 1;
242  }
243  }
244 
245  if (l1rec->iscan == l1rec->l1file->nscan) {
246  freeProductInfo(rhos_product_info);
247  }
248 
249  return (flagcld);
250 }
251 
252 char get_cloudmask_modis(l1str *l1rec, int32_t ip) {
253  // Cloud Masking for MODIS
254  static int firstCall = 1;
255  static int ib469, ib555, ib645, ib667, ib859, ib1240, ib2130;
256  int ipb;
257  float *rhos = l1rec->rhos, *cloud_albedo = l1rec->cloud_albedo;
258  float ftemp, ftemp2, ftemp3;
259  float cloudthr = 0.027;
260  char flagcld;
261 
262  if (firstCall == 1) {
263  ib469 = bindex_get(469);
264  ib555 = bindex_get(555);
265  ib645 = bindex_get(645);
266  ib667 = bindex_get(667);
267  ib859 = bindex_get(859);
268  ib1240 = bindex_get(1240);
269  ib2130 = bindex_get(2130);
270  if (ib469 < 0 || ib555 < 0 || ib645 < 0 || ib667 < 0 || ib859 < 0 || ib1240 < 0 || ib2130 < 0) {
271  printf("get_habs_cldmask: incompatible sensor wavelengths for this algorithm\n");
272  exit(EXIT_FAILURE);
273  }
274  firstCall = 0;
275  }
276 
277  ipb = l1rec->l1file->nbands*ip;
278  flagcld = 0;
279  ftemp = 0; //rhos[ipb+ib667];
280  // first correct for turbid water
281 
282  if (rhos[ipb + ib667] < 0.0) ftemp = 0.0;
283  ftemp2 = cloud_albedo[ip] - ftemp;
284 
285  if (ftemp2 > 0.027) flagcld = 1;
286 
287  // non-water check 1240 is bright relative to 859 and the combination is bright
288  // this may hit glint by accident, need to be checked.
289 
290  if (rhos[ipb + ib1240] / rhos[ipb + ib859] > 0.5 && (rhos[ipb + ib1240] + rhos[ipb + ib2130]) > 0.10) flagcld = 1;
291 
292 
293  // now try to correct for glint
294  // region check was thrown out {IF (region = "OM") cloudthr = 0.04} rjh 11/2/2015
295 
296  ftemp = rhos[ipb + ib645] - rhos[ipb + ib555] + (rhos[ipb + ib555] - rhos[ipb + ib859])*(645.0 - 555.0) / (859.0 - 555.0);
297  ftemp2 = cloud_albedo[ip] + ftemp;
298  if (ftemp2 < cloudthr) flagcld = 0;
299  if (rhos[ipb + ib859] / rhos[ipb + ib1240] > 4.0) flagcld = 0;
300 
301  // scum areas
302 
303  if ((rhos[ipb + ib859] - rhos[ipb + ib469]) > 0.01 && cloud_albedo[ip] < 0.30) flagcld = 0;
304  if ((rhos[ipb + ib859] - rhos[ipb + ib645]) > 0.01 && cloud_albedo[ip] < 0.15) flagcld = 0;
305  if (rhos[ipb + ib1240] < 0.2)
306  ftemp2 = ftemp2 - (rhos[ipb + ib859] - rhos[ipb + ib1240]) * fabs(rhos[ipb + ib859] - rhos[ipb + ib1240]) / cloudthr;
307 
308  ftemp3 = ftemp2;
309  if (ftemp2 < cloudthr * 2) {
310  if ((rhos[ipb + ib555] - rhos[ipb + ib1240]) > (rhos[ipb + ib469] - rhos[ipb + ib1240])) {
311  ftemp3 = ftemp2 - (rhos[ipb + ib555] - rhos[ipb + ib1240]);
312  } else {
313  ftemp3 = ftemp2 - (rhos[ipb + ib469] - rhos[ipb + ib1240]);
314  }
315  }
316 
317  if (ftemp3 < cloudthr) flagcld = 0;
318 
319  return (flagcld);
320 }
321 
322 char get_cldmask(l1str *l1rec, int32_t ip) {
323  //function for cloud mask by pixel
324  switch (l1rec->l1file->sensorID) {
325  case MERIS:
326  case OLCIS3A:
327  case OLCIS3B:
328  return (get_cloudmask_meris(l1rec, ip));
329  break;
330  case MODISA:
331  case MODIST:
332  return (get_cloudmask_modis(l1rec, ip));
333  break;
334  default:
335  printf("HABS cldmsk not supported for this sensor (%s).\n",
336  sensorId2SensorName(l1rec->l1file->sensorID));
337  exit(EXIT_FAILURE);
338  }
339  return (0);
340 }
get_cldmask
char get_cldmask(l1str *l1rec, int32_t ip)
Definition: cloud_flag.c:322
OLCIS3A
#define OLCIS3A
Definition: sensorDefs.h:32
freeProductInfo
void freeProductInfo(productInfo_t *info)
Definition: productInfo.cpp:145
BAD_BYTE
#define BAD_BYTE
Definition: genutils.h:25
NULL
#define NULL
Definition: decode_rs.h:63
l1rec
read l1rec
Definition: HOWTO_Add_a_sensor.txt:72
MERIS
#define MERIS
Definition: sensorDefs.h:22
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
bindex_get
int bindex_get(int32_t wave)
Definition: windex.c:45
allocateProductInfo
productInfo_t * allocateProductInfo()
Definition: productInfo.cpp:133
l1_input
l1_input_t * l1_input
Definition: l1_options.c:9
get_sdps_cld_mask
int get_sdps_cld_mask(l1str *l1rec, int32_t ip, char *cld_category)
Definition: cloud_flag.c:5
seadasutils.aquarius_utils.start
start
Definition: aquarius_utils.py:27
findProductInfo
int findProductInfo(const char *productName, int sensorId, productInfo_t *info)
Definition: productInfo.cpp:983
sensorId2SensorName
const char * sensorId2SensorName(int sensorId)
Definition: sensorInfo.c:273
get_cloudmask_modis
char get_cloudmask_modis(l1str *l1rec, int32_t ip)
Definition: cloud_flag.c:252
fabs
#define fabs(a)
Definition: misc.h:93
get_cloudmask_meris
char get_cloudmask_meris(l1str *l1rec, int32_t ip)
Definition: cloud_flag.c:144
HDF4_UTILS_MAX_NAME
#define HDF4_UTILS_MAX_NAME
Definition: hdf4utils.h:11
OLCIS3B
#define OLCIS3B
Definition: sensorDefs.h:41
strcpy
How many dimensions is the output array Default is Not sure if anything above will work correctly strcpy(l2prod->title, "no title yet")
MODISA
#define MODISA
Definition: sensorDefs.h:19
count
int count
Definition: decode_rs.h:79


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