|
ocssw
1.0
|
00001 #include "viirs_sim_sdr.h" 00002 #include <math.h> 00003 00004 int viirs_decal( ctl_struc *ctl, in_rec_struc *in_rec ) 00005 /*----------------------------------------------------------------------------- 00006 Program: viirs_decal.c 00007 00008 Description: convert VIIRS radiances for a scan into dn values 00009 The dn can be either dark subtracted (dn) or not (DN) based on the 00010 cal_dark_sub value (if 1 [default] make it dark sutract, 0 no subtract) 00011 00012 Returns type: int - 0 if good 00013 00014 Arguments: 00015 Type Name I/O Description 00016 ---- ---- --- ----------- 00017 ctl_struc * ctl I input controls 00018 in_rec_struc * in_rec I/O controls for input record reading 00019 00020 Modification history: 00021 00022 W. Robinson, SAIC 16 Sep 2010 Original development 00023 00024 ----------------------------------------------------------------------------*/ 00025 { 00026 static int first_time = 0, rvs_note = 0; 00027 static float dn_max = 4095.; /* max representable dn in 12 bits */ 00028 static float dn_trans_pct = 0.83; /* general % of max dn when gain 00029 transition occurs, NICST memo 00030 NICST_MEMO_10_012_*.doc */ 00031 int ilin, idet, ismp, ipix, ibnd, igain, iham, dn_sat; 00032 int npix, nbnd; 00033 float c0, c1, c2, a0, a1, a2, dark, lt, aoi, rvs_val, dn, big_dn; 00034 float min_rvs[MAX_BND], max_rvs[MAX_BND]; 00035 static vir_gain_struc gain; 00036 static vir_rvs_struc rvs; 00037 /* 00038 * set up the gain and rvs strucs from the source files 00039 * only do the vis NIR (7 bands) 00040 */ 00041 iham = in_rec->ham_side; 00042 nbnd = in_rec->nbnd; 00043 if( first_time == 0 ) 00044 { 00045 if( vset_cal_gain( ctl->count_decal_gain_file, &gain ) != 0 ) 00046 return 1; 00047 /* 00048 * check gain struct so it satisfys the needs 00049 */ 00050 if( gain.nham != 2 ) 00051 { 00052 printf( "%s, %d: # of HAM sides in gain struct (%d) != 2\n", __FILE__, 00053 __LINE__, gain.nham ); 00054 printf( "\tgain file: %s\n", ctl->count_decal_gain_file ); 00055 return 1; 00056 } 00057 if( gain.ndet != NDET ) 00058 { 00059 printf( 00060 "%s, %d: # of detectors in gain struct (%d) != the VIIRS # (%d)\n", 00061 __FILE__, __LINE__, gain.ndet, NDET ); 00062 printf( "\tgain file: %s\n", ctl->count_decal_gain_file ); 00063 return 1; 00064 } 00065 if( gain.ngain != 2 ) 00066 { 00067 printf( "%s, %d: # of gain states in gain struct (%d) != 2\n", __FILE__, 00068 __LINE__, gain.nham ); 00069 printf( "\tgain file: %s\n", ctl->count_decal_gain_file ); 00070 return 1; 00071 } 00072 if( gain.nbnd < N_VNIR_BND ) 00073 { 00074 printf( 00075 "%s, %d: # of bands in gain struct (%d) not >= # VIS / NIR bands (%d)\n", 00076 __FILE__, __LINE__, gain.nbnd, N_VNIR_BND ); 00077 printf( "\tgain file: %s\n", ctl->count_decal_gain_file ); 00078 return 1; 00079 } 00080 00081 if( vset_cal_rvs( ctl->count_decal_rvs_file, &rvs ) != 0 ) 00082 return 1; 00083 /* 00084 * check rvs struct so it satisfys the needs 00085 */ 00086 if( rvs.nham != 2 ) 00087 { 00088 printf( "%s, %d: # of HAM sides in RVS struct (%d) != 2\n", __FILE__, 00089 __LINE__, rvs.nham ); 00090 printf( "\trvs file: %s\n", ctl->count_decal_rvs_file ); 00091 return 1; 00092 } 00093 if( rvs.ndet != NDET ) 00094 { 00095 printf( 00096 "%s, %d: # of detectors in RVS struct (%d) != the VIIRS # (%d)\n", 00097 __FILE__, __LINE__, rvs.ndet, NDET ); 00098 printf( "\trvs file: %s\n", ctl->count_decal_rvs_file ); 00099 return 1; 00100 } 00101 if( rvs.nbnd < N_VNIR_BND ) 00102 { 00103 printf( "%s, %d: # of bands in RVS struct (%d) not >= # needed # (%d)\n", 00104 __FILE__, __LINE__, rvs.nbnd, nbnd ); 00105 printf( "\trvs file: %s\n", ctl->count_decal_rvs_file ); 00106 return 1; 00107 } 00108 00109 for( ibnd = 0; ibnd < N_VNIR_BND; ibnd++ ) 00110 { 00111 min_rvs[ibnd] = 999.; 00112 max_rvs[ibnd] = -999.; 00113 } 00114 00115 first_time = 1; 00116 } 00117 /* 00118 * loop through the pixels, lines, and bands to get the lt 00119 */ 00120 npix = in_rec->npix; 00121 for( ipix = 0; ipix < npix; ipix++ ) 00122 { 00123 /* get sample */ 00124 ismp = ipix - in_rec->margin[1]; 00125 00126 /* get the aoi from the sample */ 00127 vir_xf_scan( (float)ismp, VIR_SCAN_UASMP, VIR_SCAN_AOI, &aoi ); 00128 if( ( rvs_note == 0 ) && ( ( ipix % 200 ) == 0 ) ) 00129 { 00130 float scn_ang; 00131 vir_xf_scan( (float)ismp, VIR_SCAN_UASMP, VIR_SCAN_ANG, &scn_ang ); 00132 printf( "%s, %d: info, sample: %d, scan angle: %f, aoi: %f\n", 00133 __FILE__, __LINE__, ipix, scn_ang, aoi ); 00134 } 00135 00136 for( ilin = 0; ilin < in_rec->ndet_scan; ilin++ ) 00137 { 00138 /* get the detector from the line */ 00139 idet = ilin - in_rec->margin[0]; 00140 if( idet < 0 ) idet = 0; 00141 if( idet >= NDET ) idet = NDET - 1; 00142 00143 for( ibnd = 0; ibnd < N_VNIR_BND; ibnd++ ) 00144 { 00145 /* 00146 * if the lt is not flagged, get it and de-cal 00147 */ 00148 if( *( in_rec->bnd_q[ibnd] + ipix + npix * ilin ) == 0 ) 00149 { 00150 /* 00151 * Assume a high gain range and check later 00152 */ 00153 lt = *( in_rec->bnd_lt[ibnd] + ipix + npix * ilin ); 00154 igain = 1; 00155 dn_sat = 0; 00156 dark = *( gain.dark + ibnd + gain.nbnd * ( igain + gain.ngain * 00157 ( idet + gain.ndet * iham ) ) ); 00158 /* 00159 * get the gain, rvs coeffs 00160 */ 00161 c0 = *( gain.c0 + ibnd + gain.nbnd * ( igain + gain.ngain * 00162 ( idet + gain.ndet * iham ) ) ); 00163 c1 = *( gain.c1 + ibnd + gain.nbnd * ( igain + gain.ngain * 00164 ( idet + gain.ndet * iham ) ) ); 00165 c2 = *( gain.c2 + ibnd + gain.nbnd * ( igain + gain.ngain * 00166 ( idet + gain.ndet * iham ) ) ); 00167 00168 a0 = *( rvs.a0 + ibnd + rvs.nbnd * ( idet + rvs.ndet * iham ) ); 00169 a1 = *( rvs.a1 + ibnd + rvs.nbnd * ( idet + rvs.ndet * iham ) ); 00170 a2 = *( rvs.a2 + ibnd + rvs.nbnd * ( idet + rvs.ndet * iham ) ); 00171 00172 rvs_val = a0 + a1 * aoi + a2 * aoi * aoi; 00173 if( rvs_note == 0 ) 00174 { 00175 /* 00176 * record the min and max rvs for each band 00177 */ 00178 if( rvs_val < min_rvs[ibnd] ) 00179 min_rvs[ibnd] = rvs_val; 00180 if( rvs_val > max_rvs[ibnd] ) 00181 max_rvs[ibnd] = rvs_val; 00182 } 00183 /* 00184 * apply the gain, rvs to de cal 00185 */ 00186 if( c2 == 0 ) /* just linear fit */ 00187 { 00188 dn = ( lt * rvs_val - c0 ) / c1; 00189 } 00190 else /* a quadratic fit */ 00191 { 00192 dn = ( -c1 + sqrt( c1 * c1 - 4. * c2 * ( c0 - lt * rvs_val ) ) ) 00193 / ( 2. * c2 ); 00194 } 00195 /* 00196 * if a dual gain band, see if dn is above the transition 00197 */ 00198 if( gain.gain_ct[ibnd] == 2 ) 00199 { 00200 if( dn > ( dn_max * dn_trans_pct ) ) 00201 { 00202 igain = 0; 00203 dark = *( gain.dark + ibnd + gain.nbnd * ( igain + gain.ngain * 00204 ( idet + gain.ndet * iham ) ) ); 00205 /* 00206 * get the low gain coefficients 00207 */ 00208 c0 = *( gain.c0 + ibnd + gain.nbnd * ( igain + gain.ngain * 00209 ( idet + gain.ndet * iham ) ) ); 00210 c1 = *( gain.c1 + ibnd + gain.nbnd * ( igain + gain.ngain * 00211 ( idet + gain.ndet * iham ) ) ); 00212 c2 = *( gain.c2 + ibnd + gain.nbnd * ( igain + gain.ngain * 00213 ( idet + gain.ndet * iham ) ) ); 00214 /* 00215 * apply the gain, rvs to de cal 00216 */ 00217 if( c2 == 0 ) /* just linear fit */ 00218 { 00219 dn = ( lt * rvs_val - c0 ) / c1; 00220 } 00221 else /* a quadratic fit */ 00222 { 00223 dn = ( -c1 + sqrt( c1 * c1 - 4. * c2 * ( c0 - lt * rvs_val ) ) ) 00224 / ( 2. * c2 ); 00225 } 00226 big_dn = dn + dark; 00227 /* 00228 * if it saturates at low gain too, note it 00229 */ 00230 if( big_dn > dn_max ) 00231 { 00232 dn_sat = 1; 00233 big_dn = dn_max; 00234 } 00235 } 00236 else /* below, dual gain transition, set the raw dn */ 00237 { 00238 big_dn = dn + dark; 00239 } 00240 dn = big_dn - dark; 00241 } 00242 else /* single gain treatment */ 00243 { 00244 big_dn = dn + dark; 00245 if( big_dn > dn_max ) 00246 { 00247 dn_sat = 1; 00248 big_dn = dn_max; 00249 } 00250 dn = big_dn - dark; 00251 } 00252 /* 00253 * if DN is wanted use big_dn, else, remove the dark from big_dn 00254 */ 00255 if( ctl->count_dark_opt == 0 ) 00256 dn = big_dn - dark; 00257 else 00258 dn = big_dn; 00259 } 00260 else /* bad lt val, leave dn 0 */ 00261 { 00262 dn = 0.; 00263 dn_sat = -1; 00264 } 00265 /* 00266 * place the count... in the structure 00267 */ 00268 *( in_rec->dn[ibnd] + ipix + npix * ilin ) = dn; 00269 *( in_rec->dn_sat[ibnd] + ipix + npix * ilin ) = dn_sat; 00270 *( in_rec->gain_bit[ibnd] + ipix + npix * ilin ) = igain; 00271 } 00272 } 00273 } 00274 /* 00275 * report the rvs min, max 00276 */ 00277 if( rvs_note == 0 ) 00278 { 00279 rvs_note = 1; 00280 printf( "\n\n%s, %d:info\n", __FILE__, __LINE__ ); 00281 printf( "\nDe-calibration RVS min and max / band:\n" ); 00282 printf( "(file: %s\n", ctl->count_decal_rvs_file ); 00283 printf( "\nM band Min RVS Max RVS\n" ); 00284 for( ibnd = 0; ibnd < N_VNIR_BND; ibnd++ ) 00285 printf( " %3d %9.6f %9.6f\n", ibnd + 1, 00286 min_rvs[ibnd], max_rvs[ibnd] ); 00287 printf( "----------------------------------------\n\n" ); 00288 } 00289 /* 00290 * end loops and done 00291 */ 00292 return 0; 00293 }
1.7.6.1