28 static int numFlagMPHPixels = -1;
29 static int numFlagHABSPixels = -1;
30 static uint8_t *flags_mph =
NULL;
31 static uint8_t *flags_habs =
NULL;
35 if((numFlagMPHPixels != numPix) || !flags_mph) {
36 numFlagMPHPixels = numPix;
39 flags_mph = (uint8_t*) malloc(numPix);
43 if((numFlagHABSPixels != numPix) || !flags_habs) {
44 numFlagHABSPixels = numPix;
47 flags_habs = (uint8_t*) malloc(numPix);
52 float rhos620,
float rhos665,
float rhos709,
float rhos865,
float *ci) {
54 float kd, kd_709, ss_560;
56 kd = (0.7 * ((((rhos620 + rhos665) / 2.) - rhos865) / (((rhos443 + rhos490) / 2.) - rhos865)));
59 kd_709 = (0.7 * ((rhos709 - rhos865) / (((rhos443 + rhos490) / 2.) - rhos865)));
62 if (kd_709 > kd) kd = kd_709;
65 ss_560 = rhos560 - rhos443 + (rhos443 - rhos620)*(560 - 443) / (620 - 443);
68 if (((((rhos620 + rhos665) / 2.) - rhos865) > 0) &&
69 ((((rhos443 + rhos490) / 2.) - rhos865) > 0)) {
71 ((rhos865 <= rhos490) ||
72 (rhos865 <= rhos665) ||
73 (rhos865 <= rhos709)) &&
82 int ib0, ib1, ib2, ib3, ib4, ib5, ib6, ib7;
83 float wav0, wav1, wav2, wav3, wav4, wav5, wav6, wav7, nonci,
fac;
86 static productInfo_t* ci_product_info;
88 switch (l2rec->l1rec->l1file->sensorID) {
103 switch (l2rec->l1rec->l1file->sensorID) {
131 if (l2rec->l1rec->l1file->sensorID !=
MERIS &&
132 l2rec->l1rec->l1file->sensorID !=
OLCIS3A &&
133 l2rec->l1rec->l1file->sensorID !=
OLCIS3B) {
134 printf(
"MCI not supported for this sensor (%s).\n",
144 printf(
"HABS_CI: Hmm, something's really messed up.\n");
148 if (ci_product_info ==
NULL) {
150 findProductInfo(
"CI_cyano", l2rec->l1rec->l1file->sensorID, ci_product_info);
157 if (ib1 < 0 || ib2 < 0 || ib3 < 0) {
158 printf(
"(M)CI_stumpf: incompatible sensor wavelengths for this algorithm\n");
162 for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
164 ipb = l2rec->l1rec->l1file->nbands*ip;
167 if ((l2rec->l1rec->height[ip] == 0 && l2rec->l1rec->elev[ip] < -10) ||
168 (l2rec->l1rec->flags[ip] &
mask) != 0 ||
169 l2rec->l1rec->rhos[ipb + ib1] <= 0.0 ||
170 l2rec->l1rec->rhos[ipb + ib2] <= 0.0 ||
171 l2rec->l1rec->rhos[ipb + ib3] <= 0.0) {
173 l2rec->l1rec->flags[ip] |=
PRODFAIL;
180 ci[ip] =
fac * ((l2rec->l1rec->rhos[ipb + ib3] - l2rec->l1rec->rhos[ipb + ib1])*(wav2 - wav1) / (wav3 - wav1)
181 - (l2rec->l1rec->rhos[ipb + ib2] - l2rec->l1rec->rhos[ipb + ib1]));
184 if (l2rec->l1rec->l1file->sensorID ==
MERIS || l2rec->l1rec->l1file->sensorID ==
OLCIS3A
185 || l2rec->l1rec->l1file->sensorID ==
OLCIS3B) {
188 l2rec->l1rec->rhos[ipb + ib6], l2rec->l1rec->rhos[ipb + ib0],
189 l2rec->l1rec->rhos[ipb + ib1], l2rec->l1rec->rhos[ipb + ib3],
190 l2rec->l1rec->rhos[ipb + ib7], &ci[ip]);
193 nonci = l2rec->l1rec->rhos[ipb + ib1]
194 - l2rec->l1rec->rhos[ipb + ib0]
195 + (l2rec->l1rec->rhos[ipb + ib0]
196 - l2rec->l1rec->rhos[ipb + ib2])
197 *(wav1 - wav0) / (wav2 - wav0);
201 l2rec->l1rec->flags[ip] |=
PRODFAIL;
205 if (nonci >= 0) ci[ip] = 0;
207 if (nonci < 0) ci[ip] = 0;
214 ci[ip] =
fac * (l2rec->l1rec->rhos[ipb + ib2] - l2rec->l1rec->rhos[ipb + ib1]
215 - (l2rec->l1rec->rhos[ipb + ib3] - l2rec->l1rec->rhos[ipb + ib1])*(wav2 - wav1) / (wav3 - wav1));
223 if (ci[ip] < ci_product_info->validMin) {
224 ci[ip] = ci_product_info->validMin;
230 for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
231 if (flags_habs[ip] != 0){
233 l2rec->l1rec->flags[ip] |=
PRODFAIL;
236 if (l2rec->l1rec->iscan == l2rec->l1rec->l1file->nscan) {
253 float wav6, wav7, wav8, wav9, wav10, wav14;
254 int ip, ipb, ib6, ib7, ib8, ib9, ib10, ib14;
255 float Rmax0, Rmax1, wavmax0, wavmax1, ndvi;
256 float sipf, sicf, bair, mph0, mph1;
257 float *rhos = l2rec->l1rec->rhos;
259 if ((l2rec->l1rec->l1file->sensorID !=
MERIS) &&
260 (l2rec->l1rec->l1file->sensorID !=
OLCIS3A) &&
261 (l2rec->l1rec->l1file->sensorID !=
OLCIS3B)) {
262 printf(
"MPH not supported for this sensor (%s).\n",
281 if (ib6 < 0 || ib7 < 0 || ib8 < 0 || ib9 < 0 || ib10 < 0 || ib14 < 0) {
282 printf(
"MPH_stumpf: incompatible sensor wavelengths for this algorithm\n");
286 for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
287 ipb = l2rec->l1rec->l1file->nbands*ip;
295 if (rhos[ipb + ib8] > rhos[ipb + ib9]) {
297 Rmax0 = rhos[ipb + ib8];
300 Rmax0 = rhos[ipb + ib9];
302 if (Rmax0 > rhos[ipb + ib10]) {
307 Rmax1 = rhos[ipb + ib10];
311 sipf = rhos[ipb + ib7] - rhos[ipb + ib6] - (rhos[ipb + ib8] - rhos[ipb + ib6]) * (664 - 619) / (681 - 619);
313 sicf = rhos[ipb + ib8] - rhos[ipb + ib7] - (rhos[ipb + ib9] - rhos[ipb + ib7]) * (681 - 664) / (709 - 664);
315 ndvi = (rhos[ipb + ib14] - rhos[ipb + ib7]) / (rhos[ipb + ib14] + rhos[ipb + ib7]);
317 bair = rhos[ipb + ib9] - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (709 - 664) / (885 - 664);
318 mph0 = Rmax0 - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (wavmax0 - 664) / (885 - 664);
319 mph1 = Rmax1 - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (wavmax1 - 664) / (885 - 664);
321 if (wavmax1 != wav10) {
323 if (sicf >= 0 || sipf <= 0 || bair <= 0.002) {
324 chl_mph[ip] = 5.24e9 * pow(mph0, 4) - 1.95e8 * pow(mph0, 3) + 2.46e6 * pow(mph0, 2) + 4.02e3 * mph0 + 1.97;
326 chl_mph[ip] = 22.44 * exp(35.79 * mph1);
329 if (mph1 >= 0.02 || ndvi >= 0.2) {
331 if (sicf < 0 && sipf > 0) {
332 chl_mph[ip] = 22.44 * exp(35.79 * mph1);
338 chl_mph[ip] = 5.24e9 * pow(mph0, 4) - 1.95e8 * pow(mph0, 3) + 2.46e6 * pow(mph0, 2) + 4.02e3 * mph0 + 1.97;
354 float wav6, wav7, wav8, wav9, wav10, wav14;
355 int ip, ipb, ib6, ib7, ib8, ib9, ib10, ib14;
356 float Rmax0, Rmax1, wavmax0, wavmax1, ndvi;
357 float sipf, sicf, bair, mph1, chl_mph;
358 float *rhos = l2rec->l1rec->rhos;
359 static float thresh = 350;
361 if ((l2rec->l1rec->l1file->sensorID !=
MERIS) &&
362 (l2rec->l1rec->l1file->sensorID !=
OLCIS3A) &&
363 (l2rec->l1rec->l1file->sensorID !=
OLCIS3B)) {
364 printf(
"MPH not supported for this sensor (%s).\n",
385 if (ib6 < 0 || ib7 < 0 || ib8 < 0 || ib9 < 0 || ib10 < 0 || ib14 < 0) {
386 printf(
"MPH_stumpf: incompatible sensor wavelengths for this algorithm\n");
390 for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
393 ipb = l2rec->l1rec->l1file->nbands*ip;
396 if (l2rec->l1rec->rhos[ipb + ib6] <= 0.0 ||
397 l2rec->l1rec->rhos[ipb + ib7] <= 0.0 ||
398 l2rec->l1rec->rhos[ipb + ib8] <= 0.0 ||
399 l2rec->l1rec->rhos[ipb + ib9] <= 0.0 ||
400 (l2rec->l1rec->flags[ip] &
LAND) != 0 ||
401 (l2rec->l1rec->flags[ip] &
NAVFAIL) != 0) {
402 l2rec->l1rec->flags[ip] |=
PRODFAIL;
405 if (rhos[ipb + ib8] > rhos[ipb + ib9]) {
407 Rmax0 = rhos[ipb + ib8];
410 Rmax0 = rhos[ipb + ib9];
412 if (Rmax0 > rhos[ipb + ib10]) {
417 Rmax1 = rhos[ipb + ib10];
421 sipf = rhos[ipb + ib7] - rhos[ipb + ib6] - (rhos[ipb + ib8] - rhos[ipb + ib6]) * (664 - 619) / (681 - 619);
423 sicf = rhos[ipb + ib8] - rhos[ipb + ib7] - (rhos[ipb + ib9] - rhos[ipb + ib7]) * (681 - 664) / (709 - 664);
425 ndvi = (rhos[ipb + ib14] - rhos[ipb + ib7]) / (rhos[ipb + ib14] + rhos[ipb + ib7]);
427 bair = rhos[ipb + ib9] - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (709 - 664) / (885 - 664);
429 mph1 = Rmax1 - rhos[ipb + ib7] - (rhos[ipb + ib14] - rhos[ipb + ib7]) * (wavmax1 - 664) / (885 - 664);
431 if (wavmax1 != wav10) {
433 if (sicf < 0 && sipf > 0 && bair > 0.002) {
435 chl_mph = 22.44 * exp(35.79 * mph1);
436 if (chl_mph > thresh)
440 if (mph1 >= 0.02 || ndvi >= 0.2) {
443 if (sicf < 0 && sipf > 0) {
445 chl_mph = 22.44 * exp(35.79 * mph1);
446 if (chl_mph > thresh)
464 int ib443, ib490, ib510, ib560, ib620, ib665, ib681, ib709, ib754, ib865, ib885;
466 float *rhos = l2rec->l1rec->rhos;
467 float mdsi, cv,
mean, sum, sdev;
472 if (l2rec->l1rec->l1file->sensorID !=
MERIS &&
473 l2rec->l1rec->l1file->sensorID !=
OLCIS3A &&
474 l2rec->l1rec->l1file->sensorID !=
OLCIS3B) {
475 printf(
"HABS flags not supported for this sensor (%s).\n",
492 if (ib443 < 0 || ib490 < 0 || ib510 < 0 || ib560 < 0 || ib620 < 0 || ib665 < 0 || ib681 < 0 || ib709 < 0 || ib754 < 0 || ib865 < 0 || ib885 < 0) {
493 printf(
"get_flags_habs: incompatible sensor wavelengths for this algorithm\n");
497 for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
499 ipb = l2rec->l1rec->l1file->nbands*ip;
503 if ((
input->evalmask | 2) == 2){
504 flags_habs[ip] = (uint8_t) l2rec->l1rec->cloud[ip];
511 if (rhos[ipb + ib443] >= 0.0 &&
512 rhos[ipb + ib490] >= 0.0 &&
513 rhos[ipb + ib510] >= 0.0 &&
514 rhos[ipb + ib560] >= 0.0 &&
515 rhos[ipb + ib620] >= 0.0 &&
516 rhos[ipb + ib665] >= 0.0 &&
517 rhos[ipb + ib681] >= 0.0 &&
518 rhos[ipb + ib709] >= 0.0 &&
519 rhos[ipb + ib754] >= 0.0 &&
520 rhos[ipb + ib885] >= 0.0) {
522 if (rhos[ipb + ib885] > rhos[ipb + ib620] &&
523 rhos[ipb + ib885] > rhos[ipb + ib709] &&
524 rhos[ipb + ib885] > rhos[ipb + ib754] &&
525 rhos[ipb + ib885] > 0.01) {
530 if ((rhos[ipb + ib620] > rhos[ipb + ib560]) &&
531 (rhos[ipb + ib560] > 0.15) &&
532 (rhos[ipb + ib885] > 0.15)){
537 mdsi = (rhos[ipb + ib865] - rhos[ipb + ib885]) / (rhos[ipb + ib865] + rhos[ipb + ib885]);
540 int b[7] = {ib443, ib490, ib510, ib560, ib620, ib665, ib681};
542 for(
i = 0;
i < n; ++
i) {
543 sum += rhos[ipb +
b[
i]];
548 for(
i = 0;
i < n; ++
i) {
549 sdev += pow((rhos[ipb +
b[
i]] -
mean),2);
551 cv = sqrt(sdev / (n - 1)) /
mean;
553 if ((mdsi > 0.01) && (rhos[ipb + ib885] > 0.15) && (cv < 0.1)) {
557 float mci = 1.0 * (rhos[ipb + ib709] - rhos[ipb + ib681]
558 + (rhos[ipb + ib681] - rhos[ipb + ib754])*(709 - 681) / (754 - 681));
559 float ci = 1.0 * ((rhos[ipb + ib709] - rhos[ipb + ib665])*(681 - 665) / (709 - 665)
560 - (rhos[ipb + ib681] - rhos[ipb + ib665]));
563 rhos[ipb + ib560], rhos[ipb + ib620],
564 rhos[ipb + ib665], rhos[ipb + ib709],
565 rhos[ipb + ib865], &ci);
567 if ((mci < 0) && (ci > 0)) {
571 if (rhos[ipb + ib620] < 0.0 ||
572 rhos[ipb + ib665] < 0.0 ||
573 rhos[ipb + ib681] < 0.0 ||
574 rhos[ipb + ib709] < 0.0) {
585 int ib469, ib555, ib645, ib667, ib859, ib1240, ib2130;
587 float *rhos = l2rec->l1rec->rhos, *Rrs = l2rec->Rrs, *cloud_albedo = l2rec->l1rec->cloud_albedo;
588 float ftemp, ftemp2, ftemp3;
589 float cloudthr = 0.027;
593 if (l2rec->l1rec->l1file->sensorID !=
MODISA && l2rec->l1rec->l1file->sensorID !=
MODIST) {
594 printf(
"HABS flags not supported for this sensor (%s).\n",
607 if (ib469 < 0 || ib555 < 0 || ib645 < 0 || ib667 < 0 || ib859 < 0 || ib1240 < 0 || ib2130 < 0) {
608 printf(
"get_flags_habs: incompatible sensor wavelengths for this algorithm\n");
612 for (ip = 0; ip < l2rec->l1rec->npix; ip++) {
614 ipb = l2rec->l1rec->l1file->nbands*ip;
618 if (l2rec->chl[ip] < 0)
619 ftemp = Rrs[ipb + ib667];
621 ftemp = Rrs[ipb + ib667]*(0.45 + l2rec->chl[ip]*0.005) / 4.3;
623 if (Rrs[ipb + ib667] < 0.0) ftemp = 0.0;
624 ftemp2 = cloud_albedo[ip] - ftemp;
626 if (ftemp2 > 0.027) flags_habs[ip] |=
HABS_CLOUD;
630 if (rhos[ipb + ib1240] / rhos[ipb + ib859] > 0.5 && (rhos[ipb + ib1240] + rhos[ipb + ib2130]) > 0.10) flags_habs[ip] |=
HABS_CLOUD;
635 ftemp = rhos[ipb + ib645] - rhos[ipb + ib555] + (rhos[ipb + ib555] - rhos[ipb + ib859])*(645.0 - 555.0) / (859.0 - 555.0);
636 ftemp2 = cloud_albedo[ip] + ftemp;
637 if (ftemp2 < cloudthr) flags_habs[ip] = 0;
638 if (rhos[ipb + ib859] / rhos[ipb + ib1240] > 4.0) flags_habs[ip] = 0;
642 if ((rhos[ipb + ib859] - rhos[ipb + ib469]) > 0.01 && cloud_albedo[ip] < 0.30) flags_habs[ip] = 0;
643 if ((rhos[ipb + ib859] - rhos[ipb + ib645]) > 0.01 && cloud_albedo[ip] < 0.15) flags_habs[ip] = 0;
644 if (rhos[ipb + ib1240] < 0.2)
645 ftemp2 = ftemp2 - (rhos[ipb + ib859] - rhos[ipb + ib1240]) *
fabs(rhos[ipb + ib859] - rhos[ipb + ib1240]) / cloudthr;
648 if (ftemp2 < cloudthr * 2) {
649 if ((rhos[ipb + ib555] - rhos[ipb + ib1240]) > (rhos[ipb + ib469] - rhos[ipb + ib1240])) {
650 ftemp3 = ftemp2 - (rhos[ipb + ib555] - rhos[ipb + ib1240]);
652 ftemp3 = ftemp2 - (rhos[ipb + ib469] - rhos[ipb + ib1240]);
656 if (ftemp3 < cloudthr) flags_habs[ip] = 0;
658 if (rhos[ipb + ib555] >= 0 && rhos[ipb + ib1240] >= 0.0 && rhos[ipb + ib1240] > rhos[ipb + ib555])
665 static int32_t currentLine = -1;
667 if (currentLine == l2rec->l1rec->iscan )
669 currentLine = l2rec->l1rec->iscan;
670 switch (l2rec->l1rec->l1file->sensorID) {
681 printf(
"HABS flags not supported for this sensor (%s).\n",