|
ocssw
1.0
|
00001 #include "viirs_sim_sdr.h" 00002 #include <math.h> 00003 #include <stdlib.h> 00004 00005 int scan_cvt( in_rec_struc *in_rec, out_rec_struc *out_rec ) 00006 /*----------------------------------------------------------------------------- 00007 scan_cvt 00008 00009 purpose: make proper output record from the input record 00010 based on the scan formats for both 00011 00012 Returns type: int - 0 if good 00013 00014 Parameters (in calling order): 00015 Type Name I/O Description 00016 ---- ---- --- ----------- 00017 in_rec_struc * in_rec I/O controls for input record reading 00018 out_rec_struc * out_rec I/O controls for output file writing 00019 00020 Modification history: 00021 00022 W. Robinson, SAIC 29 Apr 2010 Original development 00023 W. Robinson, SAIC 18 Nov 2010 create the QF1_VIIRSMBANDSDR value in 00024 out_rec array qual1_m from saturation and gain info in in_rec 00025 00026 ----------------------------------------------------------------------------*/ 00027 { 00028 static int cvt_mode = -1; /* conversion mode for this run: -1 not started, 00029 0 - no conversion, 1 - make aggregated, 00030 2 - make unaggregated */ 00031 static int npixin, npixout; 00032 static float *cvt_lat, *cvt_lon, *cvt_senz, *cvt_sena, *cvt_solz, *cvt_sola, 00033 *cvt_bnd_lt[MAX_BND]; 00034 int ibnd, ilin, olin, opix, nxfr, nxfrb, ipix_st, ag_rng, num, iag, ipix; 00035 float theta1, theta2, theta, phi1, phi2, phi, sum; 00036 static unsigned char *cvt_bnd_q[MAX_BND]; 00037 unsigned char uc_val, ua_gain_vals[] = { 4, 0 }; 00038 unsigned char ua_sat_vals[] = { 0, 8 }; 00039 int v_unsat, v_sat, ag_sat_vals[] = { 8, 0, 8, 4 }; 00040 00041 /* 00042 * if in_rec->lat is null, this is signal to for end, so free local storage 00043 */ 00044 if( in_rec->lat == NULL ) 00045 { 00046 free( cvt_lat ); 00047 free( cvt_lon ); 00048 free( cvt_senz ); 00049 free( cvt_sena ); 00050 free( cvt_solz ); 00051 free( cvt_sola ); 00052 for( ibnd = 0; ibnd < out_rec->nbnd; ibnd++ ) 00053 { 00054 free( cvt_bnd_lt[ibnd] ); 00055 free( cvt_bnd_q[ibnd] ); 00056 } 00057 return 0; 00058 } 00059 /* 00060 * initialization phase 00061 */ 00062 if( cvt_mode == -1 ) 00063 { 00064 /* 00065 * determine the conversion mode based on the input scan format and 00066 * requested output scan format change 00067 */ 00068 if( ( in_rec->scn_fmt == 1 ) && ( out_rec->scn_fmt == 0 ) ) 00069 cvt_mode = 1; 00070 else if( ( in_rec->scn_fmt == 2 ) && ( out_rec->scn_fmt == 0 ) ) 00071 cvt_mode = 1; 00072 else if( ( in_rec->scn_fmt ==2 ) && ( out_rec->scn_fmt == 1 ) ) 00073 cvt_mode = 2; 00074 else 00075 { 00076 cvt_mode = 0; 00077 } 00078 if( cvt_mode != 0 ) 00079 { 00080 /* 00081 * A conversion is needed for all scans. allocate space for 00082 * the converted data and re-assign it to the out_rec locations 00083 * (from default done in init_sdr with in_rec) 00084 */ 00085 npixin = in_rec->npix; 00086 npixout = out_rec->npix; 00087 if( ( cvt_lat = (float *) 00088 malloc( out_rec->npix * out_rec->ndet_scan * sizeof (float) ) ) 00089 == NULL ) 00090 { 00091 printf( "%s, %d: Unable to allocate storage for lat conversion\n", 00092 __FILE__, __LINE__ ); 00093 return 1; 00094 } 00095 out_rec->lat = cvt_lat; 00096 if( ( cvt_lon = (float *) 00097 malloc( out_rec->npix * out_rec->ndet_scan * sizeof (float) ) ) 00098 == NULL ) 00099 { 00100 printf( "%s, %d: Unable to allocate storage for lon conversion\n", 00101 __FILE__, __LINE__ ); 00102 return 1; 00103 } 00104 out_rec->lon = cvt_lon; 00105 if( ( cvt_senz = (float *) 00106 malloc( out_rec->npix * out_rec->ndet_scan * sizeof (float) ) ) 00107 == NULL ) 00108 { 00109 printf( "%s, %d: Unable to allocate storage for senz conversion\n", 00110 __FILE__, __LINE__ ); 00111 return 1; 00112 } 00113 out_rec->senz = cvt_senz; 00114 if( ( cvt_sena = (float *) 00115 malloc( out_rec->npix * out_rec->ndet_scan * sizeof (float) ) ) 00116 == NULL ) 00117 { 00118 printf( "%s, %d: Unable to allocate storage for sena conversion\n", 00119 __FILE__, __LINE__ ); 00120 return 1; 00121 } 00122 out_rec->sena = cvt_sena; 00123 if( ( cvt_solz = (float *) 00124 malloc( out_rec->npix * out_rec->ndet_scan * sizeof (float) ) ) 00125 == NULL ) 00126 { 00127 printf( "%s, %d: Unable to allocate storage for solz conversion\n", 00128 __FILE__, __LINE__ ); 00129 return 1; 00130 } 00131 out_rec->solz = cvt_solz; 00132 if( ( cvt_sola = (float *) 00133 malloc( out_rec->npix * out_rec->ndet_scan * sizeof (float) ) ) 00134 == NULL ) 00135 { 00136 printf( "%s, %d: Unable to allocate storage for sola conversion\n", 00137 __FILE__, __LINE__ ); 00138 return 1; 00139 } 00140 out_rec->sola = cvt_sola; 00141 for( ibnd = 0; ibnd < out_rec->nbnd; ibnd++ ) 00142 { 00143 if( ( cvt_bnd_lt[ibnd] = (float *) 00144 malloc( out_rec->npix * out_rec->ndet_scan * sizeof (float) ) ) 00145 == NULL ) 00146 { 00147 printf( 00148 "%s, %d: Unable to allocate storage for bnd_lt[%d] conversion\n", 00149 __FILE__, __LINE__, ibnd ); 00150 return 1; 00151 } 00152 out_rec->bnd_lt[ibnd] = cvt_bnd_lt[ibnd]; 00153 /* 00154 */ 00155 if( ( cvt_bnd_q[ibnd] = (unsigned char *) 00156 malloc( out_rec->npix * out_rec->ndet_scan * 00157 sizeof (unsigned char) ) ) == NULL ) 00158 { 00159 printf( 00160 "%s, %d: Unable to allocate storage for bnd_q[%d] conversion\n", 00161 __FILE__, __LINE__, ibnd ); 00162 return 1; 00163 } 00164 out_rec->bnd_q[ibnd] = cvt_bnd_q[ibnd]; 00165 } 00166 } 00167 } 00168 /* 00169 * work on different conversion modes 00170 */ 00171 if( cvt_mode == 0 ) 00172 { 00173 /* 00174 * if the data is aggregated, set up qual1_m from the dn_sat * 4 00175 * (most likely, dn_sat always 0 in this case, but if not, we'll 00176 * assume the 3 state format for the standard saturation field) 00177 */ 00178 if( out_rec->scn_fmt == 0 ) 00179 { 00180 for( ibnd = 0; ibnd < in_rec->nbnd; ibnd++ ) 00181 { 00182 for( ilin = 0; ilin < in_rec->ndet_scan; ilin++ ) 00183 { 00184 for( ipix = 0; ipix < in_rec->npix; ipix++ ) 00185 { 00186 *( out_rec->qual1_m[ibnd] + ipix + in_rec->npix * ilin ) = 00187 *( in_rec->dn_sat[ibnd] + ipix + in_rec->npix * ilin ) * 4; 00188 } 00189 } 00190 } 00191 } 00192 /* 00193 * if the data is un-aggregated, fold the gain and sat into qual1_m 00194 * bit 2 = gain 0 hi, 1 lo, bit 3 = sat 0 not, 1 saturated 00195 */ 00196 else 00197 { 00198 for( ibnd = 0; ibnd < in_rec->nbnd; ibnd++ ) 00199 { 00200 for( ilin = 0; ilin < in_rec->ndet_scan; ilin++ ) 00201 { 00202 for( ipix = 0; ipix < in_rec->npix; ipix++ ) 00203 { 00204 uc_val = ua_gain_vals[ *( in_rec->gain_bit[ibnd] + ipix + 00205 in_rec->npix * ilin ) ]; 00206 uc_val = uc_val | ua_sat_vals[ (int) *( in_rec->dn_sat[ibnd] + ipix + 00207 in_rec->npix * ilin ) ]; 00208 *( out_rec->qual1_m[ibnd] + ipix + in_rec->npix * ilin ) = uc_val; 00209 } 00210 } 00211 } 00212 } 00213 return 0; 00214 } 00215 else if( cvt_mode == 2 ) 00216 { 00217 /* 00218 * mode 2 will go from margin to unaggregated, so move to proper storage 00219 */ 00220 nxfr = npixout * sizeof(float); 00221 nxfrb = npixout * sizeof(unsigned char); 00222 for( olin = 0, ilin = in_rec->margin[0]; olin < NDET; olin++, ilin++ ) 00223 { 00224 memcpy( (void *)( cvt_lat + npixout * olin ), 00225 (void *)( in_rec->lat + npixin * ilin + in_rec->margin[1] ), nxfr ); 00226 memcpy( (void *)( cvt_lon + npixout * olin ), 00227 (void *)( in_rec->lon + npixin * ilin + in_rec->margin[1] ), nxfr ); 00228 memcpy( (void *)( cvt_sena + npixout * olin ), 00229 (void *)( in_rec->sena + npixin * ilin + in_rec->margin[1] ), nxfr ); 00230 memcpy( (void *)( cvt_senz + npixout * olin ), 00231 (void *)( in_rec->senz + npixin * ilin + in_rec->margin[1] ), nxfr ); 00232 memcpy( (void *)( cvt_sola + npixout * olin ), 00233 (void *)( in_rec->sola + npixin * ilin + in_rec->margin[1] ), nxfr ); 00234 memcpy( (void *)( cvt_solz + npixout * olin ), 00235 (void *)( in_rec->solz + npixin * ilin + in_rec->margin[1] ), nxfr ); 00236 for( ibnd = 0; ibnd < out_rec->nbnd; ibnd++ ) 00237 { 00238 memcpy( (void *)( cvt_bnd_lt[ibnd] + npixout * olin ), 00239 (void *)( in_rec->bnd_lt[ibnd] + npixin * ilin + in_rec->margin[1] ), 00240 nxfr ); 00241 memcpy( (void *)( cvt_bnd_q[ibnd] + npixout * olin ), 00242 (void *)( in_rec->bnd_q[ibnd] + npixin * ilin + in_rec->margin[1] ), 00243 nxfrb ); 00244 /* 00245 * handle the conversion and transfer of saturation and gain 00246 */ 00247 for( ipix = in_rec->margin[1], opix = 0; 00248 ipix < in_rec->npix - in_rec->margin[1]; ipix++, opix++ ) 00249 { 00250 uc_val = ua_gain_vals[ *( in_rec->gain_bit[ibnd] + ipix + 00251 in_rec->npix * ilin ) ]; 00252 uc_val = uc_val | ua_sat_vals[ (int) *( in_rec->dn_sat[ibnd] + ipix + 00253 in_rec->npix * ilin ) ]; 00254 *( out_rec->qual1_m[ibnd] + opix + npixout * olin ) = uc_val; 00255 } 00256 } 00257 } 00258 } 00259 else if( cvt_mode == 1 ) 00260 { 00261 /* 00262 * mode 1 will involve aggregating and possibly trimming off the margins 00263 * Method for aggregating depends on the agg zone: 00264 * ZONE (scan location) 00265 * PARM 1:1 agg (edge) 2:1 agg 3:1 agg (center) 00266 * ------- 00267 * lat, lon, Just 'Angle Use center 00268 * view angles Transfer Average' angle 00269 * 00270 * Lt data values Average Average 00271 * 2 values 3 values 00272 */ 00273 for( olin = 0, ilin = in_rec->margin[0]; olin < out_rec->ndet_scan; 00274 olin++, ilin++ ) 00275 { 00276 for( opix = 0; opix < npixout; opix++ ) 00277 { 00278 if( ( opix >= 0 ) && ( opix <= 639 ) ) 00279 { 00280 ipix_st = opix + in_rec->margin[1]; 00281 ag_rng = 1; 00282 } 00283 else if( ( opix >= 640 ) && ( opix <= 1007 ) ) 00284 { 00285 ipix_st = 640 + in_rec->margin[1] + ( opix - 640 ) * 2; 00286 ag_rng = 2; 00287 } 00288 else if( ( opix >= 1008 ) && ( opix <= 2191 ) ) 00289 { 00290 ipix_st = 1376 + in_rec->margin[1] + ( opix - 1008 ) * 3; 00291 ag_rng = 3; 00292 } 00293 else if( ( opix >= 2192 ) && ( opix <= 2559 ) ) 00294 { 00295 ipix_st = 4928 + in_rec->margin[1] + ( opix - 2192 ) * 2; 00296 ag_rng = 2; 00297 } 00298 else if( ( opix >= 2560 ) && ( opix <= 3199 ) ) 00299 { 00300 ipix_st = 5664 + in_rec->margin[1] + opix - 2560; 00301 ag_rng = 1; 00302 } 00303 /* 00304 * Handle each agg range for the parameters 00305 */ 00306 switch( ag_rng ) 00307 { 00308 case 1: 00309 /* 00310 * at 1:1 zones, just transfer (to proper locations) 00311 */ 00312 *( cvt_lat + npixout * olin + opix ) = 00313 *( in_rec->lat + ipix_st + npixin * ilin ); 00314 *( cvt_lon + npixout * olin + opix ) = 00315 *( in_rec->lon + ipix_st + npixin * ilin ); 00316 *( cvt_senz + npixout * olin + opix ) = 00317 *( in_rec->senz + ipix_st + npixin * ilin ); 00318 *( cvt_sena + npixout * olin + opix ) = 00319 *( in_rec->sena + ipix_st + npixin * ilin ); 00320 *( cvt_solz + npixout * olin + opix ) = 00321 *( in_rec->solz + ipix_st + npixin * ilin ); 00322 *( cvt_sola + npixout * olin + opix ) = 00323 *( in_rec->sola + ipix_st + npixin * ilin ); 00324 for( ibnd = 0; ibnd < out_rec->nbnd; ibnd++ ) 00325 { 00326 *( cvt_bnd_lt[ibnd] + npixout * olin + opix ) = 00327 *( in_rec->bnd_lt[ibnd] + ipix_st + npixin * ilin ); 00328 *( cvt_bnd_q[ibnd] + npixout * olin + opix ) = 00329 *( in_rec->bnd_q[ibnd] + ipix_st + npixin * ilin ); 00330 /* 00331 * saturation aggregation - the unagg saturation is either 0 or 1 00332 * it xlates to 0 or 2 (all sat) and at bit location it 00333 * translates to 0 or 8 (x4) 00334 */ 00335 *( out_rec->qual1_m[ibnd] + opix + npixout * olin ) = 00336 *( in_rec->dn_sat[ibnd] + ipix_st + npixin * ilin ); 00337 } 00338 break; 00339 case 2: 00340 /* 00341 * at 2:1 zones, average the radiances and treat lat, lon and 00342 * view angles as vectors to find the average of 00343 * 00344 * Lt 00345 */ 00346 for( ibnd = 0; ibnd < out_rec->nbnd; ibnd++ ) 00347 { 00348 sum = 0.; 00349 num = 0; 00350 v_unsat = 0; 00351 v_sat = 0; 00352 for( iag = 0; iag < ag_rng; iag++ ) 00353 { 00354 if( 00355 ( *( in_rec->bnd_lt[ibnd] + ilin * npixin + ipix_st + iag ) > 0 ) 00356 && ( *( in_rec->bnd_q[ibnd] + ilin * npixin + ipix_st + iag ) == 0 ) ) 00357 { 00358 sum += 00359 *( in_rec->bnd_lt[ibnd] + ilin * npixin + ipix_st + iag ); 00360 num++; 00361 /* 00362 * see if saturated or unsaturated or both - the v_unsat and 00363 * v_sat work together with ag_sat_vals to set proper value 00364 */ 00365 if( *( in_rec->dn_sat[ibnd] + ilin + npixin + ipix_st + iag ) 00366 == 0 ) v_unsat = 1; 00367 if( *( in_rec->dn_sat[ibnd] + ilin + npixin + ipix_st + iag ) 00368 == 1 ) v_sat = 2; 00369 } 00370 } 00371 if( num > 0 ) 00372 { 00373 *( cvt_bnd_lt[ibnd] + npixout * olin + opix ) = sum / num; 00374 *( cvt_bnd_q[ibnd] + npixout * olin + opix ) = 0; 00375 } 00376 else 00377 { 00378 *( cvt_bnd_lt[ibnd] + npixout * olin + opix ) = -32767.; 00379 *( cvt_bnd_q[ibnd] + npixout * olin + opix ) = 2; 00380 } 00381 /* for saturation setting */ 00382 *( out_rec->qual1_m[ibnd] + npixout * olin + opix ) = 00383 ag_sat_vals[ v_sat + v_unsat ]; 00384 } 00385 /* 00386 * lat, lon, and view angles 00387 */ 00388 theta1 = *( in_rec->lat + ilin * npixin + ipix_st ); 00389 theta2 = *( in_rec->lat + ilin * npixin + ipix_st + 1 ); 00390 phi1 = *( in_rec->lon + ilin * npixin + ipix_st ); 00391 phi2 = *( in_rec->lon + ilin * npixin + ipix_st + 1 ); 00392 ang_avg( theta1, phi1, theta2, phi2, &theta, &phi ); 00393 *( cvt_lat + npixout * olin + opix ) = theta; 00394 *( cvt_lon + npixout * olin + opix ) = phi; 00395 00396 theta1 = *( in_rec->senz + ilin * npixin + ipix_st ); 00397 theta2 = *( in_rec->senz + ilin * npixin + ipix_st + 1 ); 00398 phi1 = *( in_rec->sena + ilin * npixin + ipix_st ); 00399 phi2 = *( in_rec->sena + ilin * npixin + ipix_st + 1 ); 00400 ang_avg( theta1, phi1, theta2, phi2, &theta, &phi ); 00401 *( cvt_senz + npixout * olin + opix ) = theta; 00402 *( cvt_sena + npixout * olin + opix ) = phi; 00403 00404 theta1 = *( in_rec->solz + ilin * npixin + ipix_st ); 00405 theta2 = *( in_rec->solz + ilin * npixin + ipix_st + 1 ); 00406 phi1 = *( in_rec->sola + ilin * npixin + ipix_st ); 00407 phi2 = *( in_rec->sola + ilin * npixin + ipix_st + 1 ); 00408 ang_avg( theta1, phi1, theta2, phi2, &theta, &phi ); 00409 *( cvt_solz + npixout * olin + opix ) = theta; 00410 *( cvt_sola + npixout * olin + opix ) = phi; 00411 break; 00412 case 3: 00413 /* 00414 * at 3:1 zones, average 3 lt and take center angle 00415 */ 00416 for( ibnd = 0; ibnd < out_rec->nbnd; ibnd++ ) 00417 { 00418 sum = 0.; 00419 num = 0; 00420 v_unsat = 0; 00421 v_sat = 0; 00422 for( iag = 0; iag < ag_rng; iag++ ) 00423 { 00424 if( 00425 ( *( in_rec->bnd_lt[ibnd] + ilin * npixin + ipix_st + iag ) > 0 ) 00426 && ( *( in_rec->bnd_q[ibnd] + ilin * npixin + ipix_st + iag ) == 0 ) ) 00427 { 00428 sum += 00429 *( in_rec->bnd_lt[ibnd] + ilin * npixin + ipix_st + iag ); 00430 num++; 00431 /* 00432 * see if saturated or unsaturated or both - the v_unsat and 00433 * v_sat work together with ag_sat_vals to set proper value 00434 */ 00435 if( *( in_rec->dn_sat[ibnd] + ilin + npixin + ipix_st + iag ) 00436 == 0 ) v_unsat = 1; 00437 if( *( in_rec->dn_sat[ibnd] + ilin + npixin + ipix_st + iag ) 00438 == 1 ) v_sat = 2; 00439 } 00440 } 00441 if( num > 0 ) 00442 { 00443 *( cvt_bnd_lt[ibnd] + npixout * olin + opix ) = sum / num; 00444 *( cvt_bnd_q[ibnd] + npixout * olin + opix ) = 0; 00445 } 00446 else 00447 { 00448 *( cvt_bnd_lt[ibnd] + npixout * olin + opix ) = -32767.; 00449 *( cvt_bnd_q[ibnd] + npixout * olin + opix ) = 2; 00450 } 00451 /* for saturation setting */ 00452 *( out_rec->qual1_m[ibnd] + npixout * olin + opix ) = 00453 ag_sat_vals[ v_sat + v_unsat ]; 00454 } 00455 00456 *( cvt_lat + npixout * olin + opix ) = 00457 *( in_rec->lat + ilin * npixin + ipix_st + 1 ); 00458 *( cvt_lon + npixout * olin + opix ) = 00459 *( in_rec->lon + ilin * npixin + ipix_st + 1 ); 00460 *( cvt_senz + npixout * olin + opix ) = 00461 *( in_rec->senz + ilin * npixin + ipix_st + 1 ); 00462 *( cvt_sena + npixout * olin + opix ) = 00463 *( in_rec->sena + ilin * npixin + ipix_st + 1 ); 00464 *( cvt_solz + npixout * olin + opix ) = 00465 *( in_rec->solz + ilin * npixin + ipix_st + 1 ); 00466 *( cvt_sola + npixout * olin + opix ) = 00467 *( in_rec->sola + ilin * npixin + ipix_st + 1 ); 00468 break; 00469 } 00470 } 00471 } 00472 } 00473 return 0; 00474 } 00475 00476 int ang_avg( float theta1, float phi1, float theta2, float phi2, 00477 float *theta, float *phi ) 00478 /*----------------------------------------------------------------------------- 00479 ang_avg 00480 00481 purpose: get resultant vector that is the average of 2 vectors 00482 00483 Returns type: int - 0 if good 00484 00485 Parameters (in calling order): 00486 Type Name I/O Description 00487 ---- ---- --- ----------- 00488 float theta1 I first zenith or lat angle in 00489 degrees 00490 float phi1 I first azimuth or lon angle in 00491 degrees 00492 float theta2 I second zenith or lat angle in 00493 degrees 00494 float phi2 I second azimuth or lon angle in 00495 degrees 00496 float * theta O resultant zenith or lat angle in 00497 degrees 00498 float * phi O result azimuth or lon angle in 00499 degrees 00500 00501 Modification history: 00502 00503 W. Robinson, SAIC 4 May 2010 Original development 00504 00505 ----------------------------------------------------------------------------*/ 00506 { 00507 float v1[2], v2[2], v[2]; 00508 float rad2deg = 180. / M_PI; 00509 /* 00510 * convert to just the x, y parts of a vector - this assumes all are 00511 * unit length 00512 *** don't know why but sinf, cosf... do not work, so use less compatible 00513 with floats sin, cos (gcc?) 00514 */ 00515 v1[0] = sin( theta1 / rad2deg ) * cos( phi1 / rad2deg ); 00516 v1[1] = sin( theta1 / rad2deg ) * sin( phi1 / rad2deg ); 00517 00518 v2[0] = sin( theta2 / rad2deg ) * cos( phi2 / rad2deg ); 00519 v2[1] = sin( theta2 / rad2deg ) * sin( phi2 / rad2deg ); 00520 00521 /* average the vectors and return to unit length */ 00522 v[0] = ( v1[0] + v2[0] ) / 2.; 00523 v[1] = ( v1[1] + v2[1] ) / 2.; 00524 00525 /* go back to angles */ 00526 *theta = asin( sqrt( v[0] * v[0] + v[1] * v[1] ) ) * rad2deg; 00527 *phi = atan2( v[1], v[0] ) * rad2deg; 00528 00529 return 0; 00530 }
1.7.6.1