|
ocssw
1.0
|
00001 /* 00002 * get_qaa.c 00003 * 00004 * MSl12 wrapper for Quasi-Analytic Algorithm 00005 * 00006 * Naval Research Laboratory 00007 * Stennis Space Center, MS 00008 */ 00009 00010 #include <stdlib.h> 00011 #include <math.h> 00012 #include "qaa.h" 00013 #include "l12_proto.h" 00014 #include "l2prod.h" 00015 00016 static int QaaRecNum = -1; 00017 00018 static int nbands; /* number of bands computed for line */ 00019 00020 /* one scan line of data */ 00021 00022 static float *atot; /* total absorption coefficient */ 00023 static float *adg; /* detrital absorption coefficient */ 00024 static float *aph; /* phytoplankton absorption coefficient */ 00025 static float *bb; /* backscatter coefficient */ 00026 static unsigned char *flags; /* per-pixel flags */ 00027 00028 /* pixel data set during initialization */ 00029 00030 static float *bbw; /* pure-water backscattering */ 00031 static float *aw; /* pure-water total absorption */ 00032 static float *fwave; /* wavelength (nm) */ 00033 00034 /* pixel data */ 00035 00036 static float *Rrs; /* above-water remote sensing reflectance */ 00037 static float *rrs; /* below-water remote sensing reflectance */ 00038 static float *u; 00039 static float *a_qaa; /* total absorption coefficient */ 00040 static float *bb_qaa; /* backscatter coefficient */ 00041 static float *bbp_qaa; /* backscatter coefficient */ 00042 static float *adg_qaa; /* detrital absorption coefficient */ 00043 static float *aph_qaa; /* phytoplankton absorption coefficient */ 00044 00045 static int ib410, ib440, ib490, ib555, ib670; 00046 static int do_decomp = 1; 00047 00048 /* have we run for this scan line? */ 00049 00050 static int qaa_ran(int recnum) 00051 { 00052 if ( recnum == QaaRecNum ) 00053 return 1; 00054 else 00055 return 0; 00056 } 00057 00058 /* allocate private arrays for a single scan line (412 ... 555) */ 00059 00060 static void qaa_alloc(int npix, int nbands) 00061 { 00062 atot = (float*) calloc(npix*nbands,sizeof(float)); 00063 aph = (float*) calloc(npix*nbands,sizeof(float)); 00064 adg = (float*) calloc(npix*nbands,sizeof(float)); 00065 bb = (float*) calloc(npix*nbands,sizeof(float)); 00066 flags = (unsigned char*) calloc(npix,sizeof(unsigned char)); 00067 } 00068 00069 /* allocate private arrays for a single pixel (may include a 640nm) */ 00070 00071 static void qaa_pixel_alloc(int nbands) 00072 { 00073 00074 fwave = (float*) calloc(nbands,sizeof(float)); 00075 bbw = (float*) calloc(nbands,sizeof(float)); 00076 aw = (float*) calloc(nbands,sizeof(float)); 00077 Rrs = (float*) calloc(nbands,sizeof(float)); 00078 rrs = (float*) calloc(nbands,sizeof(float)); 00079 u = (float*) calloc(nbands,sizeof(float)); 00080 a_qaa = (float*) calloc(nbands,sizeof(float)); 00081 bb_qaa = (float*) calloc(nbands,sizeof(float)); 00082 bbp_qaa = (float*) calloc(nbands,sizeof(float)); 00083 adg_qaa = (float*) calloc(nbands,sizeof(float)); 00084 aph_qaa = (float*) calloc(nbands,sizeof(float)); 00085 00086 } 00087 00088 /* NOTE: nbands and l2rec->nbands are not always equal */ 00089 00090 static void run_qaa(l2str *l2rec) 00091 { 00092 static int firstCall = 1; 00093 00094 int ip,ib,ipb; 00095 int i; 00096 unsigned char flags_qaa; 00097 00098 if (firstCall) { 00099 00100 int *w; 00101 00102 firstCall = 0; 00103 00104 /* limit to visible wavelengths */ 00105 nbands = rdsensorinfo(l2rec->sensorID,l2rec->input->evalmask,"NbandsVIS", NULL); 00106 printf("QAA v6 processing for %d bands\n", nbands ); 00107 00108 qaa_alloc(l2rec->npix,nbands); 00109 qaa_pixel_alloc(nbands); 00110 00111 for (ib=0; ib<nbands; ib++) 00112 fwave[ib] = l2rec->fwave[ib]; 00113 get_aw_bbw(l2rec,fwave,nbands,aw,bbw); 00114 00115 w = l2rec->input->qaa_wave; 00116 if (w[0] < 0 || w[1] < 0 || w[2] < 0 || w[3] < 0) { 00117 printf("qaa: algorithm coefficients not provided for this sensor.\n"); 00118 exit(1); 00119 } 00120 00121 // qaaf_v6 needs these bands at the very minimum 00122 00123 ib440 = bindex_get(w[1]); 00124 ib490 = bindex_get(w[2]); 00125 ib555 = bindex_get(w[3]); 00126 ib670 = bindex_get(w[4]); 00127 if ( ib440 < 0 || ib490 < 0 || ib555 < 0 || ib670 < 0) { 00128 printf("get_qaa: missing minimum required wavelengths " 00129 "(need 440,490,555,670).\n"); 00130 printf("get_qaa: qaa_wave[1] =%d, qaa_wave[2] =%d, " 00131 "qaa_wave[3] = %d, qaa_wave[4] = %d.\n", 00132 w[1], w[2], w[3], w[4]); 00133 exit(1); 00134 } 00135 00136 // if we want compute the aph/adg (call qaaf_decomp) we need this additional band 00137 00138 ib410 = bindex_get(w[0]); 00139 if ( ib410 < 0 ) { 00140 printf("get_qaa: incompatible sensor wavelengths for aph/adg (need 410).\n"); 00141 do_decomp = 0; 00142 } 00143 printf("QAA v6 bands: (%d) %d nm, (%d) %d nm, (%d) %d nm, (%d) %d nm, (%d) %d nm\n", 00144 ib410, w[0], ib440, w[1], ib490, w[2], ib555, w[3], ib670, w[4] ); 00145 printf("QAA v6 wav :"); 00146 for ( ib = 0; ib < nbands; ib++ ) 00147 printf(" %10.6f", fwave[ib]); 00148 printf("\n"); 00149 00150 printf("QAA v6 aw :"); 00151 for ( ib = 0; ib < nbands; ib++ ) 00152 printf(" %10.6f", aw[ib]); 00153 printf("\n"); 00154 00155 printf("QAA v6 bbw :"); 00156 for ( ib = 0; ib < nbands; ib++ ) 00157 printf(" %10.6f", bbw[ib]); 00158 printf("\n"); 00159 00160 qaa_init( ib410, ib440, ib490, ib555, ib670 ); 00161 qaa_set_param( QAA_S_PARAM, l2rec->input->qaa_adg_s); 00162 } 00163 00164 for (ip=0; ip<l2rec->npix; ip++) { 00165 00166 /* clear static globals */ 00167 for (ib=0; ib<nbands; ib++) { 00168 ipb = ip*nbands+ib; 00169 bb [ipb] = -0.1; 00170 atot[ipb] = -0.1; 00171 adg [ipb] = -0.1; 00172 aph [ipb] = -0.1; 00173 } 00174 flags[ip] = 0; 00175 00176 if ( !l2rec->mask[ip] ) { 00177 00178 flags_qaa = 0; 00179 for ( i = 0; i < nbands; i++ ) 00180 Rrs[i] = l2rec->Rrs[ip*NBANDS+i]; 00181 00182 // Version 6 00183 qaaf_v6( nbands, fwave, Rrs, aw, bbw, rrs, u, a_qaa, bb_qaa, &flags_qaa ); 00184 if ( do_decomp ) 00185 qaaf_decomp( nbands, fwave, rrs, a_qaa, aw, adg_qaa, aph_qaa, 00186 &flags_qaa ); 00187 00188 /* store results for this pixel in static globals */ 00189 00190 for (ib=0; ib<nbands; ib++) { 00191 00192 ipb = ip*nbands+ib; 00193 00194 if ( finite(bb_qaa[ib]) ) 00195 bb[ipb] = bb_qaa[ib]; 00196 else { 00197 bb[ipb] = -0.1; 00198 l2rec->flags[ip] |= PRODFAIL; 00199 } 00200 00201 if ( finite(a_qaa[ib]) ) 00202 atot[ipb] = a_qaa[ib]; 00203 else { 00204 atot[ipb] = -0.1; 00205 l2rec->flags[ip] |= PRODFAIL; 00206 } 00207 00208 if ( finite(adg_qaa[ib]) ) 00209 adg[ipb] = adg_qaa[ib]; 00210 else { 00211 adg[ipb] = -0.1; 00212 l2rec->flags[ip] |= PRODFAIL; 00213 } 00214 00215 if ( finite(aph_qaa[ib]) ) 00216 aph[ipb] = aph_qaa[ib]; 00217 else { 00218 aph[ipb] = -0.1; 00219 l2rec->flags[ip] |= PRODFAIL; 00220 } 00221 00222 /* ZP Lee, 17 August 2007 */ 00223 if (atot[ipb] > 0.0) atot[ipb] = MAX(atot[ipb],aw [ib]*1.05); 00224 if (bb [ipb] > 0.0) bb [ipb] = MAX(bb [ipb],bbw[ib]*1.05); 00225 } 00226 flags[ip] = flags_qaa; 00227 00228 } 00229 } 00230 00231 QaaRecNum = l2rec->iscan; 00232 00233 return; 00234 } 00235 00236 00237 /* interface to l2_hdf_generic() to return QAA flags */ 00238 00239 unsigned char *get_flags_qaa(l2str *l2rec) 00240 { 00241 if ( !qaa_ran(l2rec->iscan) ) 00242 run_qaa(l2rec); 00243 00244 return (flags); 00245 } 00246 00247 void get_qaa(l2str *l2rec, l2prodstr *p, float prod[]) 00248 { 00249 int prodID = p->cat_ix; 00250 int ib = p->prod_ix; 00251 int ip, ipb; 00252 00253 if ( !qaa_ran(l2rec->iscan) ) 00254 run_qaa(l2rec); 00255 00256 for (ip=0; ip<l2rec->npix; ip++) { 00257 00258 ipb = ip*nbands+ib; 00259 00260 switch (prodID) { 00261 00262 case CAT_a_qaa : 00263 if ( atot[ipb] > 0.0 ) 00264 prod[ip] = atot[ipb]; 00265 else 00266 prod[ip] = p->badData; 00267 break; 00268 00269 case CAT_adg_qaa : 00270 if ( adg[ipb] > 0.0 ) 00271 prod[ip] = adg[ipb]; 00272 else 00273 prod[ip] = p->badData; 00274 break; 00275 00276 case CAT_aph_qaa : 00277 if ( aph[ipb] > 0.0 ) 00278 prod[ip] = aph[ipb]; 00279 else 00280 prod[ip] = p->badData; 00281 break; 00282 00283 case CAT_bb_qaa : 00284 if ( bb[ipb] > 0.0 ) 00285 prod[ip] = bb[ipb]; 00286 else 00287 prod[ip] = p->badData; 00288 break; 00289 00290 case CAT_bbp_qaa : 00291 if ( (bb[ipb]-bbw[ib]) > 0.0 ) 00292 prod[ip] = bb[ipb] - bbw[ib]; 00293 else 00294 prod[ip] = p->badData; 00295 break; 00296 00297 case CAT_b_qaa : 00298 if ( bb[ipb] > 0.0 ) 00299 prod[ip] = bb[ipb] * 53.56857 + 0.00765; 00300 else 00301 prod[ip] = p->badData; 00302 break; 00303 00304 case CAT_c_qaa : 00305 if ( bb[ipb] > 0.0 && atot[ipb] > 0.0 ) 00306 prod[ip] = bb[ipb] * 53.56857 + 0.00765 + atot[ipb]; 00307 else 00308 prod[ip] = p->badData; 00309 break; 00310 00311 /* 00312 case CAT_chl_qaa : 00313 prod[ip] = (float) chl[ip]; 00314 break; 00315 */ 00316 00317 default: 00318 printf("-E- %s line %d : erroneous product ID %d passed to get_qaa().\n", 00319 __FILE__,__LINE__,prodID); 00320 exit(1); 00321 } 00322 } 00323 00324 return; 00325 } 00326 00327 00328 /* interface to convl12() to return QAA iops */ 00329 00330 void iops_qaa(l2str *l2rec) 00331 { 00332 int ib, ip; 00333 00334 if ( !qaa_ran(l2rec->iscan) ) 00335 run_qaa(l2rec); 00336 00337 for (ip=0; ip<l2rec->npix; ip++) { 00338 for (ib=0; ib<nbands; ib++) { 00339 l2rec->a [ip*NBANDS+ib] = atot[ip*nbands+ib]; 00340 l2rec->bb[ip*NBANDS+ib] = bb [ip*nbands+ib]; 00341 } 00342 } 00343 00344 return; 00345 } 00346 00347 00348 /* ------------------------------------------------------------------------------- */ 00349 /* interface to giop() */ 00350 /* ------------------------------------------------------------------------------- */ 00351 int get_bbp_qaa(l2str *l2rec, int ip, float tab_wave[], float tab_bbp[], int tab_nwave) 00352 { 00353 int ipb, ib; 00354 00355 if ( !qaa_ran(l2rec->iscan) ) 00356 run_qaa(l2rec); 00357 00358 for (ib=0; ib<nbands; ib++) { 00359 ipb = ip*nbands+ib; 00360 bbp_qaa[ib] = bb[ipb]-bbw[ib]; 00361 if (bbp_qaa[ib] < 0) 00362 return(0); 00363 } 00364 00365 ipb = ip*nbands; 00366 for (ib=0; ib<tab_nwave; ib++) { 00367 tab_bbp[ib] = linterp(fwave,bbp_qaa,nbands,tab_wave[ib]); 00368 } 00369 00370 return(1); 00371 }
1.7.6.1