|
ocssw
1.0
|
00001 /* =================================================================== */ 00002 /* Algorithms for computing diffuse attenuation coefficient for MSl12. */ 00003 /* */ 00004 /* B. Franz, NASA Ocean Biology Processing Group, SAIC, March 2005. */ 00005 /* =================================================================== */ 00006 00007 #include <stdlib.h> 00008 #include <math.h> 00009 #include "l12_proto.h" 00010 00011 #define KD_MAX 6.4 00012 #define KD_MIN 0.016 00013 00014 static float kdbad = BAD_FLT; 00015 00016 /* ------------------------------------------------------------------- */ 00017 /* Kd490_KD2 - diffuse attenuation at 490nm (2-band polynomial). */ 00018 /* */ 00019 /* Inputs: */ 00020 /* l2rec - level-2 structure containing one complete scan after */ 00021 /* atmospheric correction. */ 00022 /* Outputs: */ 00023 /* k490 - diffuse attenuation coefficient, 1 value per pixel. */ 00024 /* */ 00025 /* Algorithm: P.J.Werdell, June 2009 */ 00026 /* ------------------------------------------------------------------- */ 00027 void Kd490_KD2(l2str *l2rec, float *Kd) 00028 { 00029 static int32_t *w = NULL; 00030 static float *a = NULL; 00031 static int ib1 = -1; 00032 static int ib2 = -1; 00033 00034 int32_t ip, ipb; 00035 float R; 00036 00037 if (w == NULL) { 00038 w = l2rec->input->kd2w; 00039 a = l2rec->input->kd2c; 00040 if (w[0] < 0 || w[1] < 0) { 00041 printf("Kd490_KD2: algorithm coefficients not provided for this sensor.\n"); 00042 exit(1); 00043 } 00044 ib1 = bindex_get(w[0]); 00045 ib2 = bindex_get(w[1]); 00046 if (ib1 < 0 || ib2 < 0) { 00047 printf("Kd490_KD2: incompatible sensor wavelengths for this algorithm\n"); 00048 exit(1); 00049 } 00050 } 00051 00052 for (ip=0; ip<l2rec->npix; ip++) { 00053 00054 ipb = NBANDS*ip; 00055 00056 if (l2rec->mask[ip] || 00057 l2rec->Rrs[ipb+ib1] <= 0.0 || l2rec->Rrs[ipb+ib2] <= 0.0) { 00058 Kd[ip] = kdbad; 00059 l2rec->flags[ip] |= PRODFAIL; 00060 } else { 00061 R = log10(l2rec->Rrs[ipb+ib1]/l2rec->Rrs[ipb+ib2]); 00062 if (isnan(R)) { 00063 Kd[ip] = kdbad; 00064 l2rec->flags[ip] |= PRODFAIL; 00065 } else { 00066 Kd[ip] = a[0]+pow(10.0,a[1] + R*(a[2]+R*(a[3]+R*(a[4]+R*a[5])))); 00067 if (Kd[ip] > KD_MAX) { 00068 Kd[ip] = KD_MAX; 00069 l2rec->flags[ip] |= PRODWARN; 00070 } 00071 } 00072 } 00073 } 00074 } 00075 00076 00077 /* ------------------------------------------------------------------- */ 00078 /* Kd490_mueller() - diffuse attenuation at 490nm (J. Mueller). */ 00079 /* */ 00080 /* Inputs: */ 00081 /* l2rec - level-2 structure containing one complete scan after */ 00082 /* atmospheric correction. */ 00083 /* Outputs: */ 00084 /* k490 - diffuse attenuation coefficient, 1 value per pixel. */ 00085 /* */ 00086 /* Algorithm provided by: J. Mueller */ 00087 /* OCTS coefficents: S. Bailey, 16 July 2001 */ 00088 /* ------------------------------------------------------------------- */ 00089 void Kd490_mueller(l2str *l2rec, float k490[]) 00090 { 00091 int32_t ip, ipb; 00092 00093 static float a = 0.15645; 00094 static float b = -1.5401; 00095 static float Kw = 0.016; 00096 00097 static float badval = 0.0; 00098 static float maxval = 6.4; 00099 00100 if (l2rec->sensorID == OCTS) { 00101 /* Coefs for 490/565 band combination */ 00102 a = 0.2166; 00103 b = -1.6355; 00104 } 00105 00106 for (ip=0; ip<l2rec->npix; ip++) { 00107 00108 ipb = NBANDS*ip; 00109 00110 if (l2rec->mask[ip]) { /* pixel was masked */ 00111 k490[ip] = kdbad; 00112 l2rec->flags[ip] |= PRODFAIL; 00113 } else if (l2rec->nLw[ipb+2] <= 0.0 ) { /* unknown, high attenuation */ 00114 k490[ip] = maxval; 00115 } else if (l2rec->nLw[ipb+4] <= 0.0 ) { /* unknown, low attenuation */ 00116 k490[ip] = Kw; 00117 } else { 00118 k490[ip] = a*pow(l2rec->nLw[ipb+2]/l2rec->nLw[ipb+4],b) + Kw; 00119 if (k490[ip] > maxval) { 00120 k490[ip] = maxval; 00121 } 00122 } 00123 } 00124 } 00125 00126 00127 /* ------------------------------------------------------------------- */ 00128 /* Kd490_obpg - diffuse attenuation at 490nm (Mueller & Werdell). */ 00129 /* */ 00130 /* Inputs: */ 00131 /* l2rec - level-2 structure containing one complete scan after */ 00132 /* atmospheric correction. */ 00133 /* Outputs: */ 00134 /* k490 - diffuse attenuation coefficient, 1 value per pixel. */ 00135 /* */ 00136 /* Algorithm provided by: J. Mueller */ 00137 /* New coefficents and updated form: P.J.Werdell, February 2005. */ 00138 /* ------------------------------------------------------------------- */ 00139 void Kd490_obpg(l2str *l2rec, float k490[]) 00140 { 00141 int32_t ip, ipb; 00142 00143 static float a; 00144 static float b; 00145 static int32_t ib490; 00146 static int32_t ib555; 00147 00148 static int firstCall = 1; 00149 00150 if (firstCall) { 00151 00152 /* select 490/555 or 490/565 fit, depending on proximity of */ 00153 /* sensor bands to fitted bands */ 00154 00155 if ((ib555 = bindex_get(551)) > 0) { 00156 a = 0.1853; 00157 b = -1.349; 00158 } else if ((ib555 = bindex_get(565)) > 0) { 00159 a = 0.1787; 00160 b = -1.122; 00161 } else { 00162 printf("Kd_obpg: incompatible sensor wavelengths (no 555 or 565).\n"); 00163 exit(1); 00164 } 00165 00166 if ((ib490 = bindex_get(490)) < 0) { 00167 printf("Kd_obpg: incompatible sensor wavelengths (no 490).\n"); 00168 exit(1); 00169 } 00170 00171 firstCall = 0; 00172 } 00173 00174 for (ip=0; ip<l2rec->npix; ip++) { 00175 00176 ipb = NBANDS*ip; 00177 00178 if (l2rec->mask[ip] || 00179 l2rec->nLw[ipb+ib490] <= 0.0 || l2rec->nLw[ipb+ib555] <= 0.0) { 00180 k490[ip] = kdbad; 00181 l2rec->flags[ip] |= CHLFAIL; 00182 l2rec->flags[ip] |= PRODFAIL; 00183 } else { 00184 k490[ip] = a*pow(l2rec->nLw[ipb+ib490]/l2rec->nLw[ipb+ib555],b); 00185 if (k490[ip] > KD_MAX) { 00186 k490[ip] = KD_MAX; 00187 l2rec->flags[ip] |= PRODWARN; 00188 } 00189 /* not until reprocessing 00190 if (k490[ip] < KD_MIN) { 00191 k490[ip] = KD_MIN; 00192 l2rec->flags[ip] |= PRODWARN; 00193 } 00194 */ 00195 } 00196 00197 } 00198 } 00199 00200 00201 /* ------------------------------------------------------------------- */ 00202 /* Kd490_morel() - diffuse attenuation at 490 using Morel (2007) */ 00203 /* */ 00204 /* Inputs: */ 00205 /* l2rec - level-2 structure containing one complete scan after */ 00206 /* atmospheric correction. */ 00207 /* band - waveband number (0 - nbands-1) at which Kd computed. */ 00208 /* */ 00209 /* Outputs: */ 00210 /* Kd - diffuse attenuation at 490nm, 1 value per pixel. */ 00211 /* */ 00212 /* Description: */ 00213 /* This produces the estimate of diffuse attenation at 490 */ 00214 /* using the satellite derived chlorophyll. */ 00215 /* */ 00216 /* Reference: */ 00217 /* */ 00218 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */ 00219 /* Consistency of products derived from various ocean color sensors: */ 00220 /* An examination before merging these products and extending their */ 00221 /* applications, Remote Sensing of Environment, to be submitted. */ 00222 /* */ 00223 /* New equation 8 (combined LOV and NOMAD derivation) */ 00224 /* */ 00225 /* Original Implementation: B. Franz, February 2007 */ 00226 /*---------------------------------------------------------------------*/ 00227 void Kd490_morel(l2str *l2rec, float *Kd) 00228 { 00229 static float Kw = 0.01660; 00230 static float X = 0.077746; 00231 static float e = 0.672846; 00232 00233 float chl; 00234 int32_t ip; 00235 00236 for (ip=0; ip<l2rec->npix; ip++) { 00237 00238 chl = l2rec->chl[ip]; 00239 00240 if (l2rec->mask[ip] || chl <= 0.0) { 00241 Kd[ip] = kdbad; 00242 l2rec->flags[ip] |= PRODFAIL; 00243 } else { 00244 Kd[ip] = Kw + X * pow(chl,e); 00245 if (Kd[ip] > KD_MAX) { 00246 Kd[ip] = KD_MAX; 00247 l2rec->flags[ip] |= PRODWARN; 00248 } else 00249 if (Kd[ip] < KD_MIN) { 00250 Kd[ip] = KD_MIN; 00251 l2rec->flags[ip] |= PRODWARN; 00252 } 00253 } 00254 } 00255 00256 return; 00257 } 00258 00259 00260 /* ------------------------------------------------------------------- */ 00261 /* Kd_PAR_morel() - spectrally integrated attenuation using Morel */ 00262 /* */ 00263 /* Inputs: */ 00264 /* l2rec - level-2 structure containing one complete scan after */ 00265 /* atmospheric correction. */ 00266 /* depth - layer depth 1=1/k490, 2=2/k490 */ 00267 /* */ 00268 /* Outputs: */ 00269 /* Kd - Kd(PAR), 1 value per pixel. */ 00270 /* */ 00271 /* Description: */ 00272 /* */ 00273 /* Reference: */ 00274 /* */ 00275 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */ 00276 /* Consistency of products derived from various ocean color sensors: */ 00277 /* An examination before merging these products and extending their */ 00278 /* applications, Remote Sensing of Environment, to be submitted. */ 00279 /* */ 00280 /* Original Implementation: B. Franz, November 2006 */ 00281 /*---------------------------------------------------------------------*/ 00282 void Kd_PAR_morel(l2str *l2rec, int depth, float *Kd) 00283 { 00284 static float *Kd490 = NULL; 00285 00286 int32_t ip; 00287 00288 if (Kd490 == NULL) { 00289 if ((Kd490 = malloc(l2rec->npix*sizeof(float))) == NULL) { 00290 printf("-E- %s: Error allocating space for %d records.\n",__FILE__,l2rec->npix); 00291 exit(1); 00292 } 00293 } 00294 00295 Kd490_morel(l2rec,Kd490); 00296 00297 for (ip=0; ip<l2rec->npix; ip++) { 00298 00299 if (l2rec->mask[ip] || Kd490[ip] <= 0.0) { 00300 Kd[ip] = kdbad; 00301 l2rec->flags[ip] |= PRODFAIL; 00302 } else { 00303 switch (depth) { 00304 case 1: 00305 Kd[ip] = 0.0864 + 0.884 * Kd490[ip] - 0.00137/Kd490[ip]; 00306 break; 00307 case 2: 00308 Kd[ip] = 0.0665 + 0.874 * Kd490[ip] - 0.00121/Kd490[ip]; 00309 break; 00310 default: 00311 printf("-E- %s: Invalid depth for Kd(PAR) (1 or 2).\n",__FILE__); 00312 exit(1); 00313 break; 00314 } 00315 } 00316 } 00317 00318 return; 00319 } 00320 00321 /* ------------------------------------------------------------------- */ 00322 /* Zhl_morel() - heated layer depth using Morel (2007) */ 00323 /* */ 00324 /* Inputs: */ 00325 /* l2rec - level-2 structure containing one complete scan after */ 00326 /* atmospheric correction. */ 00327 /* */ 00328 /* Outputs: */ 00329 /* Zhl - heated layer depth (m) */ 00330 /* */ 00331 /* Description: */ 00332 /* */ 00333 /* Reference: */ 00334 /* */ 00335 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */ 00336 /* Consistency of products derived from various ocean color sensors: */ 00337 /* An examination before merging these products and extending their */ 00338 /* applications, Remote Sensing of Environment, to be submitted. */ 00339 /* */ 00340 /* Original Implementation: B. Franz, November 2006 */ 00341 /*---------------------------------------------------------------------*/ 00342 void Zhl_morel(l2str *l2rec, float *Zhl) 00343 { 00344 static float *KdPAR = NULL; 00345 00346 int32_t ip; 00347 00348 /* need Kd(PAR) at 2nd optical depth */ 00349 if (KdPAR == NULL) { 00350 if ((KdPAR = malloc(l2rec->npix*sizeof(float))) == NULL) { 00351 printf("-E- %s: Error allocating space for %d records.\n",__FILE__,l2rec->npix); 00352 exit(1); 00353 } 00354 } 00355 Kd_PAR_morel(l2rec,2,KdPAR); 00356 00357 for (ip=0; ip<l2rec->npix; ip++) { 00358 00359 if (l2rec->mask[ip] || KdPAR[ip] <= 1e-5) { 00360 Zhl[ip] = -999; 00361 l2rec->flags[ip] |= PRODFAIL; 00362 } else 00363 Zhl[ip] = 2.0/KdPAR[ip]; 00364 } 00365 00366 return; 00367 } 00368 00369 00370 /* ------------------------------------------------------------------- */ 00371 /* Kd532() - spectral diffuse attenuation using Mueller/Austin&Petzold */ 00372 /* */ 00373 /* Inputs: */ 00374 /* l2rec - level-2 structure containing one complete scan after */ 00375 /* atmospheric correction. */ 00376 /* flag - return Kd in meters (1) or 1/meters */ 00377 /* Kd - diffuse attenuation at 490 nm. */ 00378 /* Outputs: */ 00379 /* Kd - diffuse attenuation at 532 nm. */ 00380 /* */ 00381 /* Description: */ 00382 /* This produces the estimate of diffuse attenation at 532 nm from */ 00383 /* the estimate of diffuse attenuation at 490 nm using the spectral */ 00384 /* K algorithm of Austin and Petzold. */ 00385 /* */ 00386 /* Reference: */ 00387 /* Austin, R. W and T. J. Petzold, "Spectral Dependence of the */ 00388 /* Diffuse Attenuation Coefficient of Light in Ocean Waters", */ 00389 /* SPIE Vol 489 Ocean Optics (1984) pp 168-178 */ 00390 /* */ 00391 /* Original Implementation: P. Martinolich, NRL/Stennis, May 2005 */ 00392 /*---------------------------------------------------------------------*/ 00393 00394 void Kd532(l2str *l2rec, int flag, float k532[]) 00395 { 00396 const float M532 = 0.68052; 00397 const float KW490 = 0.0224; 00398 const float KW532 = 0.05356; 00399 00400 static float maxval = 6.4; 00401 00402 int ip; 00403 00404 float temp; 00405 00406 for (ip=0; ip<l2rec->npix; ip++) { 00407 00408 if ( k532[ip] < 0.0 ) 00409 k532[ip] = kdbad; 00410 else if ( k532[ip] >= maxval ) 00411 k532[ip] = maxval; 00412 else if (l2rec->mask[ip]) 00413 k532[ip] = kdbad; 00414 else { 00415 temp = M532*(k532[ip]-KW490) + KW532; 00416 if ( flag > 0 ) 00417 k532[ip] = 1.0 / temp; 00418 else 00419 k532[ip] = temp; 00420 if ( k532[ip] < 0.0 ) 00421 k532[ip] = kdbad; 00422 if ( k532[ip] > maxval ) 00423 k532[ip] = maxval; 00424 } 00425 00426 } 00427 } 00428 00429 /* ------------------------------------------------------------------- */ 00430 /* Kd_lee() - spectral diffuse attenuation using Lee, et. (2005) */ 00431 /* */ 00432 /* Inputs: */ 00433 /* l2rec - level-2 structure containing one complete scan after */ 00434 /* atmospheric correction. */ 00435 /* band - waveband number (0 - nbands-1) at which Kd computed. */ 00436 /* */ 00437 /* Outputs: */ 00438 /* Kd - diffuse attenuation at specified band, 1 value per pixel. */ 00439 /* */ 00440 /* Description: */ 00441 /* This produces the estimate of diffuse attenation at the given band */ 00442 /* using the satellite derived total absorption and backscattering. */ 00443 /* */ 00444 /* Reference: */ 00445 /* Lee, Z.P., M. Darecki, K.L. Carder, C.O.Davis, */ 00446 /* D. Stramski, W.J. Rhea, "Diffuse Attenuation coefficient of */ 00447 /* downwelling irradiance: An evalution of remote sensing methods". */ 00448 /* */ 00449 /* Original Implementation: P. Martinolich, NRL/Stennis, March 2005 */ 00450 /*---------------------------------------------------------------------*/ 00451 void Kd_lee(l2str *l2rec, int band, float *Kd) 00452 { 00453 00454 const float m1 = 4.18; 00455 const float m2 = 0.52; 00456 const float m3 = -10.80; 00457 00458 float m0; 00459 int ip, ipb; 00460 00461 if (l2rec->input->iop_opt == IOPNONE) { 00462 printf("IOP-based Kd_lee product requires iop model selection (iop_opt). "); 00463 printf("Using default model.\n"); 00464 l2rec->input->iop_opt = IOPDEFAULT; 00465 get_iops(l2rec,l2rec->input->iop_opt); 00466 } 00467 00468 for (ip=0; ip<l2rec->npix; ip++) { 00469 00470 ipb = ip*NBANDS + band; 00471 00472 if (l2rec->mask[ip] || 00473 l2rec->a[ipb] <= 0.0 || l2rec->bb[ipb] <= 0.0 ) { 00474 Kd[ip] = kdbad; 00475 l2rec->flags[ip] |= PRODFAIL; 00476 } else { 00477 m0 = 1.0 + 0.005 * l2rec->solz[ip]; 00478 Kd[ip] = m0 * l2rec->a[ipb] 00479 + m1 * (1.0 - m2 * exp( m3 * l2rec->a[ipb])) * l2rec->bb[ipb]; 00480 if (Kd[ip] > KD_MAX) { 00481 Kd[ip] = KD_MAX; 00482 l2rec->flags[ip] |= PRODWARN; 00483 } else 00484 if (Kd[ip] < KD_MIN) { 00485 Kd[ip] = KD_MIN; 00486 l2rec->flags[ip] |= PRODWARN; 00487 } 00488 } 00489 } 00490 00491 return; 00492 } 00493 00494 00495 00496 /* ------------------------------------------------------------------- */ 00497 /* Kd_PAR_lee() - spectrally integrated attenuation using ZP Lee */ 00498 /* */ 00499 /* Inputs: */ 00500 /* l2rec - level-2 structure containing one complete scan after */ 00501 /* atmospheric correction. */ 00502 /* */ 00503 /* Outputs: */ 00504 /* Kd - Kd(PAR), 1 value per pixel. */ 00505 /* */ 00506 /* Description: */ 00507 /* */ 00508 /* Reference: */ 00509 /* */ 00510 /* Original Implementation: B. Franz, September 2007 */ 00511 /*---------------------------------------------------------------------*/ 00512 void Kd_PAR_lee(l2str *l2rec, float *Kd) 00513 { 00514 void Zphotic_lee(l2str *l2rec, l2prodstr *p, float Zp[]); 00515 00516 float *Zp = Kd; 00517 l2prodstr p; 00518 int32_t ip; 00519 00520 // get first optical depth from Lee 00521 p.cat_ix = CAT_Zphotic_lee; 00522 p.prod_ix = -3; 00523 Zphotic_lee(l2rec, &p, Zp); 00524 00525 // invert to get Kd(PAR) at 1st optical depth 00526 for (ip=0; ip<l2rec->npix; ip++) { 00527 if (Zp[ip] > 0.0) 00528 Kd[ip] = 1.0/Zp[ip]; 00529 else { 00530 Kd[ip] = kdbad; 00531 l2rec->flags[ip] |= PRODFAIL; 00532 } 00533 } 00534 } 00535 00536 /* ------------------------------------------------------------------- */ 00537 /* Kd_jamet() - spectral diffuse attenuation using Jamet et al, (2012) */ 00538 /* */ 00539 /* Inputs: */ 00540 /* l2rec - level-2 structure containing one complete scan after */ 00541 /* atmospheric correction. */ 00542 /* band - waveband number (0 - nbands-1) at which Kd computed. */ 00543 /* */ 00544 /* Outputs: */ 00545 /* Kd - diffuse attenuation at specified band, 1 value per pixel. */ 00546 /* */ 00547 /* Description: */ 00548 /* This derives the Kd using the full set of Rrs available and only */ 00549 /* for the SeaWiFS instrument (at this time). A neural network is */ 00550 /* used to perform this. */ 00551 /* */ 00552 /* Reference: */ 00553 /* Jamet, C., H. Loisel, and D. Dessailly (2012), Retrieval of the */ 00554 /* spectral diffuse attenuation coefficient Kd(lambda) in open and */ 00555 /* coastal ocean waters using a neural network inversion, J Geophys. */ 00556 /* Res., 117, C10023, doi: 10.1029/2012JC008076. */ 00557 /* */ 00558 /* Original Implementation: C. Jamet */ 00559 /*---------------------------------------------------------------------*/ 00560 void Kd_jamet(l2str *l2rec, int band, float *Kd) 00561 { 00562 static int firstCall = 1; 00563 int ip, ipb; 00564 static float alg_lambda[16], alg_inp[16]; 00565 int32 i, j, bad_prd; 00566 static float *b1, *w1, *w2, *mean, *sd, b2; 00567 static float *xnorm, *alayer; 00568 static int32 nrrs, n_input, n_norm, n_nodes; 00569 static int32 bnd_ptr[NBANDS]; 00570 float y; 00571 /* 00572 * Initialize the coefficients for the computation 00573 */ 00574 if( firstCall ) 00575 { 00576 FILE *fp; 00577 char *filedir; 00578 char filename[FILENAME_MAX]; 00579 char line [801]; 00580 int com_ln; 00581 int poub; 00582 00583 firstCall = 0; 00584 if( ( filedir = getenv("OCDATAROOT")) == NULL) 00585 { 00586 printf("-E- %s, %d: OCDATAROOT env variable undefined.\n", 00587 __FILE__, __LINE__ ); 00588 exit(1); 00589 } 00590 strcpy(filename, filedir); 00591 switch (l2rec->sensorID) 00592 { 00593 case SEAWIFS: 00594 strcat(filename, "/seawifs/iop/kd_jamet/seawifs_kd_jamet.dat"); 00595 break; 00596 case MERIS: 00597 strcat(filename, "/meris/iop/kd_jamet/meris_kd_jamet.dat"); 00598 break; 00599 default: 00600 printf( 00601 "-E- %s, %d: Kd_jamet can only be generated for the SEAWIFS instrument\n", 00602 __FILE__, __LINE__ ); 00603 exit(1); 00604 break; 00605 } 00606 printf("Loading Kd coefficient table %s\n", filename ); 00607 00608 if( ( fp = fopen(filename,"r")) == NULL ) 00609 { 00610 fprintf(stderr,"-E- %s, %d: unable to open %s for reading\n", 00611 __FILE__,__LINE__,filename); 00612 exit(1); 00613 } 00614 /* 00615 * read the # bands for the Rrs and their values 00616 */ 00617 com_ln = 1; 00618 while( com_ln == 1 ) 00619 { 00620 if( fgets(line,800,fp) == NULL ) { 00621 fprintf( stderr,"-E- %s, %d: Error reading %s\n", 00622 __FILE__, __LINE__, filename ); 00623 exit(1); 00624 } 00625 if (line[0] != ';' ) com_ln = 0; 00626 } 00627 if( sscanf( line, "%d", &nrrs ) != 1 ) { 00628 fprintf( stderr,"-E- %s, %d: Error reading %s\n", 00629 __FILE__, __LINE__, filename ); 00630 exit(1); 00631 } 00632 for( i = 0; i < nrrs; i++ ) 00633 { 00634 fgets(line,800,fp); 00635 sscanf( line, "%f", &alg_lambda[i] ); 00636 } 00637 n_input = nrrs + 1; 00638 n_norm = n_input + 1; 00639 /* 00640 * get the # nodes in the hidden layer 00641 */ 00642 com_ln = 1; 00643 while( com_ln == 1 ) 00644 { 00645 fgets(line,800,fp); 00646 if (line[0] != ';' ) com_ln = 0; 00647 } 00648 sscanf( line, "%d", &n_nodes ); 00649 00650 b1 = malloc( n_nodes * sizeof( float ) ); 00651 w1 = malloc( n_nodes * n_norm * sizeof( float ) ); 00652 w2 = malloc( n_nodes * sizeof( float ) ); 00653 mean = malloc( n_norm * sizeof( float ) ); 00654 sd = malloc( n_norm * sizeof( float ) ); 00655 xnorm = malloc( n_input *sizeof( float ) ); 00656 alayer = malloc( n_nodes * sizeof( float ) ); 00657 00658 if( ( b1 == NULL ) || ( w1 == NULL ) || ( w2 == NULL ) || 00659 ( mean == NULL ) || ( sd == NULL ) || 00660 ( xnorm == NULL ) || ( alayer == NULL ) ) 00661 { 00662 fprintf(stderr,"-E- %s, %d: allocation of Kd weights failed\n", 00663 __FILE__,__LINE__ ); 00664 exit(1); 00665 } 00666 /* 00667 * get the weights for the hidden and output layers 00668 */ 00669 com_ln = 1; 00670 while( com_ln == 1 ) 00671 { 00672 fgets(line,800,fp); 00673 if (line[0] != ';' ) com_ln = 0; 00674 } 00675 /* ( this line is a throw-away line) */ 00676 for(i=0; i<n_nodes; i++) 00677 { 00678 fgets(line,800,fp); 00679 sscanf(line,"%d %d %f",&poub,&poub,&b1[i]); 00680 } 00681 fgets(line,800,fp); 00682 sscanf(line,"%d %d %f",&poub,&poub,&b2); 00683 for(j=0; j<n_input; j++) 00684 for(i=0; i<n_nodes; i++) 00685 { 00686 fgets(line,800,fp); 00687 sscanf(line,"%d %d %f",&poub,&poub, ( w1 + i + n_nodes * j ) ); 00688 } 00689 for(j=0; j<n_nodes; j++) 00690 { 00691 fgets(line,800,fp); 00692 sscanf(line,"%d %d %f",&poub,&poub,&w2[j]); 00693 } 00694 /* 00695 * lastly, read the mean and standard deviation data 00696 */ 00697 com_ln = 1; 00698 while( com_ln == 1 ) 00699 { 00700 fgets(line,800,fp); 00701 if (line[0] != ';' ) com_ln = 0; 00702 } 00703 sscanf(line, "%f",&mean[0]); 00704 for(i=1; i< n_norm; i++) { 00705 fgets(line,800,fp); 00706 sscanf(line, "%f",&mean[i]); 00707 } 00708 com_ln = 1; 00709 while( com_ln == 1 ) 00710 { 00711 fgets(line,800,fp); 00712 if (line[0] != ';' ) com_ln = 0; 00713 } 00714 sscanf( line, "%f",&sd[0]); 00715 for(i=1; i< n_norm; i++) { 00716 fgets(line,800,fp); 00717 sscanf( line, "%f",&sd[i]); 00718 } 00719 fclose( fp ); 00720 /* 00721 * set up correct band locations 00722 */ 00723 printf( "Kd_jamet band assignment:\n" ); 00724 printf( "Alg wave inst wave inst band\n" ); 00725 for( i = 0; i < nrrs; i++ ) { 00726 bnd_ptr[i] = windex( alg_lambda[i], l2rec->fwave, l2rec->nbands ); 00727 printf( "%7.2f %7.2f %3d\n", alg_lambda[i], 00728 l2rec->fwave[bnd_ptr[i]], bnd_ptr[i] ); 00729 } 00730 printf( "\n" ); 00731 } 00732 /* 00733 * derive the Kd 00734 */ 00735 for (ip=0; ip<l2rec->npix; ip++) { 00736 ipb = ip * NBANDS; 00737 /* 00738 * Normalize the Rrs, lambda 00739 */ 00740 bad_prd = 0; 00741 if( l2rec->mask[ip] ) bad_prd = 1; 00742 /* 00743 for( i = 0; i < nrrs; i++ ) 00744 { 00745 alg_inp[i] = l2rec->Rrs[ ipb + bnd_ptr[i] ]; 00746 if( i <= 1 ) { 00747 if( alg_inp[i] <= -0.01 ) bad_prd = 1; 00748 } 00749 else { 00750 if( alg_inp[i] <= 0.0 ) bad_prd = 1; 00751 } 00752 } 00753 */ 00754 for( i = 0; i < nrrs; i++ ) 00755 { 00756 alg_inp[i] = l2rec->Rrs[ ipb + bnd_ptr[i] ]; 00757 if( alg_inp[i] < BAD_FLT + 1 ) bad_prd = 1; 00758 } 00759 alg_inp[nrrs] = l2rec->fwave[band]; 00760 00761 if( bad_prd == 1 ){ 00762 Kd[ip] = kdbad; 00763 l2rec->flags[ip] |= PRODFAIL; 00764 } 00765 else { 00766 for( i = 0; i < n_input; i++ ) 00767 xnorm[i] = ( ( 2./3. ) *( alg_inp[i] - mean[i] ) ) / sd[i]; 00768 /* 00769 * Apply the nn algorithm 00770 */ 00771 for( i = 0; i < n_nodes; i++ ){ 00772 alayer[i] = 0.0; 00773 for( j = 0; j < n_input; j++ ){ 00774 alayer[i] += ( xnorm[j] * *( w1 + i + n_nodes * j ) ); 00775 } 00776 alayer[i] = 1.715905 * (float)tanh( ( 2./3.) * 00777 ( double)( alayer[i] + b1[i] ) ); 00778 } 00779 00780 for( j = 0, y = 0.; j < n_nodes; j++ ) 00781 y += ( alayer[j] * w2[j] ); 00782 /* 00783 * De-normalize the log(Kd) and make Kd 00784 */ 00785 y = 1.5 * ( y + b2) * sd[n_input] + mean[n_input]; 00786 *( Kd + ip )= (float)pow( 10., (double)y ); 00787 } 00788 } 00789 return; 00790 } 00791 00792 /* ------------------------------------------------------------------- */ 00793 /* get_Kd() - l2_hdf_generic interface for Kd */ 00794 /* ------------------------------------------------------------------- */ 00795 void get_Kd(l2str *l2rec, l2prodstr *p, float prod[]) 00796 { 00797 switch (p->cat_ix) { 00798 case CAT_Kd_mueller: 00799 Kd490_mueller(l2rec,prod); 00800 break; 00801 case CAT_Kd_obpg: 00802 Kd490_obpg(l2rec,prod); 00803 break; 00804 case CAT_Kd_KD2: 00805 Kd490_KD2(l2rec,prod); 00806 break; 00807 case CAT_Kd_lee: 00808 Kd_lee(l2rec,p->prod_ix,prod); 00809 break; 00810 case CAT_Kd_532: 00811 Kd490_mueller(l2rec,prod); 00812 Kd532(l2rec,p->prod_ix,prod); 00813 break; 00814 case CAT_Kd_morel: 00815 Kd490_morel(l2rec,prod); 00816 break; 00817 case CAT_Kd_jamet: 00818 Kd_jamet(l2rec,p->prod_ix,prod); 00819 break; 00820 case CAT_KPAR_morel: 00821 Kd_PAR_morel(l2rec,p->prod_ix,prod); 00822 break; 00823 case CAT_KPAR_lee: 00824 Kd_PAR_lee(l2rec,prod); 00825 break; 00826 case CAT_Zhl_morel: 00827 Zhl_morel(l2rec,prod); 00828 break; 00829 default: 00830 printf("Error: %s : Unknown product specifier: %d\n",__FILE__,p->cat_ix); 00831 exit(1); 00832 break; 00833 } 00834 00835 return; 00836 } 00837 00838 00839 00840 00841 00842 /* =================================================================== */ 00843 /* defunct versions */ 00844 /* =================================================================== */ 00845 00846 00847 /* ------------------------------------------------------------------- */ 00848 /* Kd_morel() - spectral diffuse attenuation using Morel (2007) */ 00849 /* */ 00850 /* Inputs: */ 00851 /* l2rec - level-2 structure containing one complete scan after */ 00852 /* atmospheric correction. */ 00853 /* band - waveband number (0 - nbands-1) at which Kd computed. */ 00854 /* */ 00855 /* Outputs: */ 00856 /* Kd - diffuse attenuation at specified band, 1 value per pixel. */ 00857 /* */ 00858 /* Description: */ 00859 /* This produces the estimate of diffuse attenation at the given band */ 00860 /* using the satellite derived chlorophyll. */ 00861 /* */ 00862 /* Reference: */ 00863 /* */ 00864 /* Morel, A., Y. Huot, B. Gentili, P.J. Werdell, S.B. Hooker (2007). */ 00865 /* Consistency of products derived from various ocean color sensors: */ 00866 /* An examination before merging these products and extending their */ 00867 /* applications, Remote Sensing of Environment, to be submitted. */ 00868 /* */ 00869 /* Equation 8 */ 00870 /* */ 00871 /* Original Implementation: B. Franz, August 2006 */ 00872 /*---------------------------------------------------------------------*/ 00873 /* 00874 void Kd_morel(l2str *l2rec, int band, float *Kd) 00875 { 00876 typedef struct kd_morel_table { 00877 float wave; 00878 float Kw; 00879 float X; 00880 float e; 00881 } kdtabstr; 00882 00883 static kdtabstr *kdtab = NULL; 00884 static int ntab = 0; 00885 00886 float chl; 00887 float Kd1, Kd2; 00888 float wt; 00889 float wave; 00890 int32_t ip, itab, i1, i2; 00891 00892 if (kdtab == NULL) { 00893 FILE *fp = NULL; 00894 char *tmp_str; 00895 char file [FILENAME_MAX] = ""; 00896 if ((tmp_str = getenv("OCDATAROOT")) == NULL) { 00897 printf("OCDATAROOT environment variable is not defined.\n"); 00898 exit(1); 00899 } 00900 strcpy(file,tmp_str); strcat(file,"/common/kd_morel.dat"); 00901 if ( (fp = fopen(file,"r")) == NULL ) { 00902 printf("-E- %s: Error opening file %s.\n",__FILE__,file); 00903 exit(1); 00904 } 00905 printf("\nLoading Morel spectral Kd table\n %s\n",file); 00906 fscanf(fp,"%d",&ntab); 00907 if (ntab <= 0) { 00908 printf("-E- %s: Error reading %s.\n",__FILE__,file); 00909 exit(1); 00910 } 00911 if ((kdtab = malloc(ntab*sizeof(kdtabstr))) == NULL) { 00912 printf("-E- %s: Error allocating space for %d records.\n",__FILE__,ntab); 00913 exit(1); 00914 } 00915 for (itab=0; itab<ntab; itab++) { 00916 fscanf(fp,"%f %f %f %f", 00917 &kdtab[itab].wave, 00918 &kdtab[itab].Kw, 00919 &kdtab[itab].X, 00920 &kdtab[itab].e); 00921 printf("%f %f %f %f\n", 00922 kdtab[itab].wave, 00923 kdtab[itab].Kw, 00924 kdtab[itab].X, 00925 kdtab[itab].e); 00926 } 00927 fclose(fp); 00928 printf("\n"); 00929 } 00930 00931 if (band >= 0) 00932 wave = l2rec->fwave[band]; 00933 else 00934 wave = 490.; 00935 00936 for (itab=0; itab<ntab; itab++) 00937 if (wave <= kdtab[itab].wave) 00938 break; 00939 00940 if (wave == kdtab[itab].wave) { 00941 i1 = itab; 00942 i2 = itab; 00943 } else if (itab == ntab) { 00944 i1 = itab-2; 00945 i2 = itab-1; 00946 wt = (wave-kdtab[i1].wave)/(kdtab[i2].wave - kdtab[i1].wave); 00947 } else { 00948 i1 = itab; 00949 i2 = itab+1; 00950 wt = (wave-kdtab[i1].wave)/(kdtab[i2].wave - kdtab[i1].wave); 00951 } 00952 00953 for (ip=0; ip<l2rec->npix; ip++) { 00954 00955 chl = l2rec->chl[ip]; 00956 00957 if (l2rec->mask[ip] || chl <= 0.0) { 00958 Kd[ip] = kdbad; 00959 l2rec->flags[ip] |= PRODFAIL; 00960 } else { 00961 if (i1 == i2) { 00962 Kd[ip] = kdtab[i1].Kw + kdtab[i1].X * pow(chl,kdtab[i1].e); 00963 } else { 00964 Kd1 = kdtab[i1].Kw + kdtab[i1].X * pow(chl,kdtab[i1].e); 00965 Kd2 = kdtab[i2].Kw + kdtab[i2].X * pow(chl,kdtab[i2].e); 00966 Kd[ip] = Kd1 + (Kd2-Kd1)*wt; 00967 } 00968 if (Kd[ip] > KD_MAX) { 00969 Kd[ip] = KD_MAX; 00970 l2rec->flags[ip] |= PRODWARN; 00971 } else 00972 if (Kd[ip] < KD_MIN) { 00973 Kd[ip] = KD_MIN; 00974 l2rec->flags[ip] |= PRODWARN; 00975 } 00976 } 00977 } 00978 00979 return; 00980 } 00981 */
1.7.6.1