|
ocssw
1.0
|
00001 #include "l12_proto.h" 00002 00003 #define NTAB 100 00004 #define MAXLINE 1024 00005 00006 static float badval = BAD_FLT; 00007 static int32_t LastRecNum = -1; 00008 static float adg_s = 0.018; 00009 static float *chl; 00010 static float *idx; 00011 00012 00013 /* ---------------------------------------------------------------------------- */ 00014 /* run_cdom_morel - computes A. Morel's CDOM index and associated prodsucts. */ 00015 /* */ 00016 /* Reference: */ 00017 /* */ 00018 /* A. Morel, G. Gentili, A simple band ratio technique to quantify the colored */ 00019 /* disolved and detrital organic material from ocean color remotely sensed data.*/ 00020 /* Rem. Sens. Env., 2009. */ 00021 /* */ 00022 /* Implementation: */ 00023 /* */ 00024 /* B. Franz, NASA/OBPG, April 2009. */ 00025 /* */ 00026 /* ---------------------------------------------------------------------------- */ 00027 void run_cdom_morel(l2str *l2rec) 00028 { 00029 static int firstCall = 1; 00030 static float idxtab [NTAB][NTAB]; // CDOM index 00031 static float chltab [NTAB][NTAB]; // Chlorophyll 00032 static float xtab [NTAB]; // R412:R443 ratio 00033 static float ytab [NTAB]; // R490:R555 ratio 00034 static int nx; 00035 static int ny; 00036 static int ib412 = -1; 00037 static int ib443 = -1; 00038 static int ib490 = -1; 00039 static int ib555 = -1; 00040 00041 float Q0[NBANDS]; 00042 float R412; 00043 float R443; 00044 float R490; 00045 float R555; 00046 float xrat, yrat; 00047 float t,u,w[4],wt; 00048 int i,j,ip,ipb; 00049 00050 if (firstCall) { 00051 00052 char filename[FILENAME_MAX]; 00053 char *delim = ";"; 00054 FILE *fp; 00055 char *tmp_str, *path; 00056 char line [MAXLINE]; 00057 char buff [MAXLINE]; 00058 00059 if ((path = getenv("OCDATAROOT")) == NULL) { 00060 printf("OCDATAROOT environment variable is not defined.\n"); 00061 exit(1); 00062 } 00063 00064 // number of visible bands and specific band indices 00065 00066 ib412 = bindex_get(412); 00067 ib443 = bindex_get(443); 00068 ib490 = bindex_get(490); 00069 ib555 = bindex_get(551); 00070 ib555 = bindex_get(545); 00071 if (ib555 < 0) ib555 = bindex_get(550); 00072 if (ib555 < 0) ib555 = bindex_get(555); 00073 if (ib555 < 0) ib555 = bindex_get(560); 00074 if (ib412 < 0 || ib443 < 0 || ib490 < 0 || ib555 < 0) { 00075 printf("cdom_morel: incompatible sensor wavelengths for this algorithm\n"); 00076 exit(1); 00077 } 00078 00079 // read look-up table for band ratio to chl 00080 00081 strcpy( filename, path); 00082 strcat( filename, "/common/morel_chl_R490_R555.dat"); 00083 fp = fopen( filename,"r"); 00084 if (!fp) { 00085 printf("-E- %s line %d: error opening (%s) file", __FILE__,__LINE__,filename); 00086 exit(1); 00087 } 00088 00089 printf("\nLoading Morel CDOM table from file %s\n", filename); 00090 00091 j=0; 00092 while (fgets( line, MAXLINE, fp ) != NULL) { 00093 strcpy(buff,line); 00094 tmp_str = strtok(buff,delim); 00095 if (strcmp(tmp_str,"R412/R443") == 0) { 00096 i=0; 00097 while ((tmp_str = strtok(NULL,delim)) != NULL) { 00098 xtab[i] = atof(tmp_str); 00099 i++; 00100 } 00101 nx = i; 00102 } else if (strcmp(tmp_str,"R490/R555") == 0) { 00103 ; 00104 } else { 00105 ytab[j] = atof(tmp_str); 00106 i = 0; 00107 while ((tmp_str = strtok(NULL,delim)) != NULL) { 00108 if (strcmp(tmp_str,"NA") != 0) 00109 chltab[i][j] = atof(tmp_str); 00110 else 00111 chltab[i][j] = -999; 00112 i++; 00113 } 00114 if (i != nx) { 00115 printf("-E- %s line %d: error reading (%s) file", __FILE__,__LINE__,filename); 00116 exit(1); 00117 } 00118 j++; 00119 } 00120 } 00121 ny = j; 00122 00123 printf("Read %d x %d entries.\n",nx,ny); 00124 00125 // read look-up table for band ratio to CDOM index 00126 00127 strcpy( filename, path); 00128 strcat( filename, "/common/morel_cdm_index.dat"); 00129 fp = fopen( filename,"r"); 00130 if (!fp) { 00131 printf("-E- %s line %d: error opening (%s) file", __FILE__,__LINE__,filename); 00132 exit(1); 00133 } 00134 00135 printf("\nLoading Morel CDOM table from file %s\n", filename); 00136 00137 j=0; 00138 while (fgets( line, MAXLINE, fp ) != NULL) { 00139 strcpy(buff,line); 00140 tmp_str = strtok(buff,delim); 00141 if (strcmp(tmp_str,"R412/R443") == 0) { 00142 i=0; 00143 while ((tmp_str = strtok(NULL,delim)) != NULL) { 00144 xtab[i] = atof(tmp_str); 00145 i++; 00146 } 00147 nx = i; 00148 } else if (strcmp(tmp_str,"R490/R555") == 0) { 00149 ; 00150 } else { 00151 ytab[j] = atof(tmp_str); 00152 i = 0; 00153 while ((tmp_str = strtok(NULL,delim)) != NULL) { 00154 if (strcmp(tmp_str,"NA") != 0) 00155 idxtab[i][j] = atof(tmp_str); 00156 else 00157 idxtab[i][j] = -999; 00158 i++; 00159 } 00160 if (i != nx) { 00161 printf("-E- %s line %d: error reading (%s) filename", __FILE__,__LINE__,filename); 00162 exit(1); 00163 } 00164 j++; 00165 } 00166 } 00167 ny = j; 00168 00169 printf("Read %d x %d entries.\n",nx,ny); 00170 00171 // allocate space for static arrays 00172 00173 if ((chl = calloc(l2rec->npix,sizeof(float))) == NULL) { 00174 printf("-E- %s line %d : error allocating memory for Morel CDOM.\n", 00175 __FILE__,__LINE__); 00176 exit(1); 00177 } 00178 00179 if ((idx = calloc(l2rec->npix,sizeof(float))) == NULL) { 00180 printf("-E- %s line %d : error allocating memory for Morel CDOM.\n", 00181 __FILE__,__LINE__); 00182 exit(1); 00183 } 00184 00185 firstCall = 0; 00186 } 00187 00188 // table look-up and interpolation for CDOM index and chl 00189 00190 for (ip=0; ip<l2rec->npix; ip++) { 00191 ipb = ip*NBANDS; 00192 idx[ip] = badval; 00193 chl[ip] = badval; 00194 if (l2rec->chl[ip] > 0.0) { 00195 // get Rrs and convert to irradiance reflectance 00196 qint_morel(l2rec->fwave,l2rec->nbands,0.0,l2rec->chl[ip],Q0); 00197 R412 = l2rec->Rrs[ipb+ib412] * Q0[ib412]; 00198 R443 = l2rec->Rrs[ipb+ib443] * Q0[ib443]; 00199 R490 = l2rec->Rrs[ipb+ib490] * Q0[ib490]; 00200 R555 = conv_rrs_to_555(l2rec->Rrs[ipb+ib555],l2rec->fwave[ib555]) * Q0[ib555]; 00201 if (R443 > 0.0 && R555 > 0.0) { 00202 // compute ratios and locate bounding table entries 00203 xrat = MAX(MIN(R412/R443,xtab[nx-1]),xtab[0]); 00204 yrat = MAX(MIN(R490/R555,ytab[ny-1]),ytab[0]); 00205 for (i=0; i<nx-2; i++) 00206 if (xtab[i] > xrat) 00207 break; 00208 for (j=0; j<ny-2; j++) 00209 if (ytab[j] > yrat) 00210 break; 00211 // bilinearly interpolate, weigh missing table values to zero 00212 t = (xrat-xtab[i])/(xtab[i+1]-xtab[i]); 00213 u = (yrat-ytab[j])/(ytab[j+1]-ytab[j]); 00214 00215 w[0] = (idxtab[i ][j ] > -1) * (1-t)*(1-u); 00216 w[1] = (idxtab[i ][j+1] > -1) * t*(1-u); 00217 w[2] = (idxtab[i+1][j+1] > -1) * t*u; 00218 w[3] = (idxtab[i+1][j ] > -1) * (1-t)*u; 00219 00220 wt = w[0] + w[1] + w[2] + w[3]; 00221 00222 if (wt > 0.0) { 00223 idx[ip] = (idxtab[i ][j ] * w[0] 00224 + idxtab[i ][j+1] * w[1] 00225 + idxtab[i+1][j+1] * w[2] 00226 + idxtab[i+1][j ] * w[3]) / wt; 00227 chl[ip] = (chltab[i ][j ] * w[0] 00228 + chltab[i ][j+1] * w[1] 00229 + chltab[i+1][j+1] * w[2] 00230 + chltab[i+1][j ] * w[3]) / wt; 00231 } else { 00232 l2rec->flags[ip] |= PRODFAIL; // missing table values 00233 } 00234 } else { 00235 l2rec->flags[ip] |= PRODFAIL; // bad input Rrs 00236 } 00237 } else { 00238 l2rec->flags[ip] |= PRODFAIL; // bad input chl 00239 } 00240 } 00241 00242 LastRecNum = l2rec->iscan; 00243 } 00244 00245 00246 /* ------------------------------------------------------------------- */ 00247 /* test if this line has been processed */ 00248 /* ------------------------------------------------------------------- */ 00249 int cdom_morel_ran(int recnum) 00250 { 00251 if ( recnum == LastRecNum ) 00252 return 1; 00253 else 00254 return 0; 00255 } 00256 00257 00258 /* ------------------------------------------------------------------- */ 00259 /* compute absorption (adg) for given chl and CDOM index */ 00260 /* ------------------------------------------------------------------- */ 00261 float adg_morel(float chl, float idx, float wave) 00262 { 00263 if (chl > 0.0 && idx > badval+1) 00264 return (idx*0.065*pow(chl,0.63)*exp(-adg_s*(wave-400))); 00265 else 00266 return (badval); 00267 } 00268 00269 00270 /* ------------------------------------------------------------------- */ 00271 /* compute percent CDOM for given chl and CDOM index */ 00272 /* ------------------------------------------------------------------- */ 00273 float pcdom_morel(float chl, float idx) 00274 { 00275 float adg; 00276 float aph; 00277 00278 if (chl > 0.0 && idx > badval+1) { 00279 aph = 0.03782*pow(chl,0.635); 00280 adg = adg_morel(chl,idx,440.0); 00281 return (100*adg/(adg+aph)); 00282 } else 00283 return (badval); 00284 } 00285 00286 00287 /* ------------------------------------------------------------------- */ 00288 /* compute cdom-corrected chlorophyll */ 00289 /* ------------------------------------------------------------------- */ 00290 float chl_cdomcorr_morel(float chl, float idx) 00291 { 00292 float A[] = {-73.65,-35.92,15.3,14.8}; 00293 float X, Y; 00294 00295 if (chl > 0.0 && idx > 0.0) { 00296 idx = MIN(MAX(idx,0.1),10.0); 00297 X = log10(idx); 00298 Y = X*(A[0]+X*(A[1]+X*(A[2]+X*A[3]))); 00299 return (chl*(1+Y/100)); 00300 } else 00301 return (badval); 00302 } 00303 00304 00305 /* ------------------------------------------------------------------- */ 00306 /* interface to convl12() */ 00307 /* ------------------------------------------------------------------- */ 00308 void get_cdom_morel(l2str *l2rec, l2prodstr *p, float prod[]) 00309 { 00310 int prodID = p->cat_ix; 00311 int ib = p->prod_ix; 00312 00313 int32_t ip, ipb; 00314 00315 if (!cdom_morel_ran(l2rec->iscan)) 00316 run_cdom_morel(l2rec); 00317 00318 for (ip=0; ip<l2rec->npix; ip++) { 00319 00320 ipb = ip*NBANDS+ib; 00321 00322 switch (prodID) { 00323 00324 case CAT_iCDOM_morel : 00325 prod[ip] = (float) idx[ip]; 00326 break; 00327 00328 case CAT_chl_morel : 00329 prod[ip] = (float) chl[ip]; 00330 break; 00331 00332 case CAT_adg_morel : 00333 prod[ip] = adg_morel(chl[ip],idx[ip],l2rec->fwave[ib]); 00334 break; 00335 00336 case CAT_pCDOM_morel : 00337 prod[ip] = pcdom_morel(chl[ip],idx[ip]); 00338 break; 00339 00340 case CAT_chl_cdomcorr_morel : 00341 prod[ip] = chl_cdomcorr_morel(l2rec->chl[ip],idx[ip]); 00342 break; 00343 00344 default: 00345 printf("-E- %s line %d : erroneous product ID %d passed to CDOM morel.\n", 00346 __FILE__,__LINE__,prodID); 00347 exit(1); 00348 } 00349 } 00350 00351 return; 00352 } 00353 00354
1.7.6.1