|
ocssw
1.0
|
00001 /*----------------------------------------------------------------------------- 00002 Program: anc_seaice 00003 00004 Description: convert binary file of sea ice from wgrib2 into a 00005 stabndard file of sea ice 00006 00007 Arguments: 00008 Type Name I/O Description 00009 ---- ---- --- ----------- 00010 int argc I count of command line args - 3 00011 char *[] argv I command line arguments: 00012 see check_usage for all options 00013 00014 Modification history: 00015 00016 W. Robinson, SAIC 9 Dec 2009 initial routine 00017 W. Robinson, SAIC 14 Jan 2011 Add smoothed icemask 00018 00019 ----------------------------------------------------------------------------*/ 00020 #include "ancil.h" 00021 #include <time.h> 00022 #include <limits.h> 00023 #include <unistd.h> 00024 #include "ancnrt_proto.h" 00025 00026 /* 00027 * Met data specific settings 00028 */ 00029 00030 #define VGROUPNAME "Geophysical Data" 00031 00032 #define VSIZE 0.0833 00033 #define HSIZE 0.0833 00034 #define MAX_NORTH 90.0 00035 #define MAX_SOUTH -90.0 00036 #define MAX_WEST -180.0 00037 #define MAX_EAST 180.0 00038 #define GRIBLATSZ 2160 /* GRIB file size in lat */ 00039 #define GRIBLONSZ 4320 /* GRIB file size in lon */ 00040 #define OUTLATSZ 2160 /* output # lines */ 00041 #define OUTLONSZ 4320 /* output # pixels */ 00042 #define GRIB_MODE_1 1 /* the GRIB modes, GRIB 1 old, GRIB 2 post jan 08 */ 00043 #define GRIB_MODE_2 2 00044 00045 /* #define VSIZE 2.5 */ 00046 /* #define HSIZE 2.5 */ 00047 /* #define MAX_EAST 177.5 */ 00048 /* #define WPHLATSZ 73 */ /* latitude of GRIB files */ 00049 /* #define WPHLONSZ 144 */ /* long of MODIFIED GRIB files */ 00050 00051 /* fortran prototype */ 00052 void mk_smooth_ice_map_(char* char_ice, float* frac_ice_smoothed); 00053 00054 00055 int main(int argc, char *argv[]) 00056 { 00057 int i, j, k, l, anctyp, int_char_min, int_char_max; 00058 int rank, SDSref = 0; 00059 int result = 0; 00060 int array_size = 0; 00061 int numannarr = 0; 00062 int one = 1, rgmode = 1; 00063 int year = 0, month = 0, day = 0, hour = 0, msec = 0, 00064 julval = 0, pyear=0; 00065 int32 shape[2], shape_smooth[2]; 00066 00067 float32 *datarr1, dat; 00068 00069 static float32 latscl[GRIBLATSZ], lonscl[GRIBLONSZ], dlat, dlon; 00070 char message[MAXNAMELNG], inz_wind[MAXNAMELNG], inm_wind[MAXNAMELNG]; 00071 char inpress[MAXNAMELNG], inp_water[MAXNAMELNG], inrel_hum[MAXNAMELNG]; 00072 char in_oz[MAXNAMELNG], longname[MAXNAMELNG], *sdsname; 00073 char *units, *dataattr, *longnamelabel, *longnamevalue; 00074 char *unitslabel, *unitsvalue, *datafmt, outfile[MAXNAMELNG]; 00075 char outfilename[MAXNAMELNG], outdir[MAXNAMELNG]; 00076 char annotfile[MAXNAMELNG], source_name[100]; 00077 int n_opt_arg; 00078 float latstep, lonstep; /* lat/lon incr. */ 00079 int latsz, lonsz; /* lat/lon size */ 00080 float nmostlat, smostlat, wmostlon, emostlon; /* lat/lon corners */ 00081 time_t t; /* processing time */ 00082 struct tm *localtm; /* processing time */ 00083 struct annotation *xannot; 00084 int grib_mode, npix = 4320, nlin = 2160; 00085 char grib2_t_str[300]; 00086 00087 /* 00088 * HDF datafile variables 00089 */ 00090 00091 int32 sdfid, fid, gridid, dimid, sdsid, sdsref; 00092 int32 datatype; 00093 00094 /* 00095 * data type array pointers 00096 */ 00097 00098 float32 *float_SDSz_wind, *float_SDSm_wind, *float_SDSpress; 00099 float32 *float_SDSp_water, *float_SDSrel_hum; 00100 int16 *float_SDS_oz, *oz_newsiz; 00101 int8 *int8_SDSdataQC, *seaice8; 00102 00103 /* 00104 * external functions used 00105 */ 00106 00107 int readgrib(), readgrib1(), julian_(), wrtattr(), startHDF(); 00108 int addAttr(), setSDSref(); 00109 int32 wrtsds(); 00110 void deattachHDFGrid(), pexit (); 00111 int8 check_usage(); 00112 struct annotation *fillenv_annot(); 00113 00114 /* 00115 * ------- check command line arguments and set args ------------------ 00116 */ 00117 00118 if( ( check_usage( argc, argv, &anctyp, &n_opt_arg, source_name, 00119 &grib_mode, grib2_t_str ) ) != 0 ) 00120 exit (-1); 00121 strcpy(annotfile, argv[ n_opt_arg ]); 00122 strcpy(in_oz, argv[ 1 + n_opt_arg ]); 00123 if (argc - n_opt_arg < 3) strcpy(outdir, "./"); 00124 else strcpy(outdir, argv[ 2 + n_opt_arg ]); 00125 00126 if( ( datarr1 = malloc( GRIBLATSZ * GRIBLONSZ * sizeof( float ) ) ) 00127 == NULL ) 00128 { 00129 printf( "%s - %d: Unable to allocate space for datarr1 array\n", __FILE__, 00130 __LINE__ ); 00131 return 1; 00132 } 00133 /* 00134 * read the data 00135 */ 00136 if( grib_mode == GRIB_MODE_2 ) 00137 result = readgrib2( in_oz, grib2_t_str, npix, nlin, rgmode, datarr1, 00138 &year, &month, &day, &hour ); 00139 else 00140 result = readgrib1(in_oz, npix, nlin, datarr1, &year, &month, &day, &hour ); 00141 00142 if (result == 1) { 00143 strcpy(message, "readgrib file: "); 00144 strcat(message, in_oz ); 00145 pexit (message); 00146 } 00147 00148 /* 00149 * Create outfile name 00150 * Note that for grib 2 files, whole year is available already 00151 */ 00152 if ( grib_mode == GRIB_MODE_1 ) 00153 { 00154 if (year < 90) year = year + 2000; /* 2 digit 19xx or 20xx yrs to 4 */ 00155 else year = year + 1900; 00156 } 00157 00158 julval = julian_(&day, &month, &year); 00159 sprintf( outfile, "%s/N%04d%03d00_SEAICE_%s_24h.hdf", 00160 outdir, year, julval, source_name ); 00161 00162 sprintf(outfilename, "N%04d%03d00_SEAICE_%s_24h.hdf", 00163 year, julval, source_name ); 00164 /* For SDPS' benefit! */ 00165 printf("%s+%04d%03d000000000+%04d%03d235959999\n", outfile, year, julval, 00166 year, julval); 00167 00168 00169 /* 00170 * Write HDF annotations, lat and lon SDSs 00171 */ 00172 00173 /**** check for existing HDF ***/ 00174 00175 if (!access(outfile, F_OK)) pexit ("....output file exists"); 00176 00177 if ((numannarr = count_annot(annotfile)) == 0) 00178 pexit ("no annotations found"); 00179 00180 if ((xannot = (struct annotation *) 00181 malloc (sizeof(struct annotation) * numannarr)) 00182 == NULL) pexit ("malloc Annotation"); 00183 00184 xannot = fillenv_annot(annotfile); 00185 00186 /* 00187 * Determine processing time 00188 */ 00189 00190 (void) time(&t); 00191 localtm = localtime(&t); 00192 pyear = localtm->tm_year; 00193 if (pyear < 90) pyear = pyear + 2000; /* 2 digit 19xx or 20xx yrs to 4 */ 00194 else pyear = pyear + 1900; 00195 00196 if (hour > 0) msec = hour * 60 * 60 * 1000; 00197 else msec = 0; 00198 00199 /* 00200 * ------- assign metadata values to local descriptive variables ---- 00201 * ------- insert dates and other values to metadata array --------- 00202 * ------- follow order of fillenv data file ----------------------- 00203 */ 00204 00205 for(i = 0; i < numannarr; i++) { 00206 00207 /* Insert values to descr array element */ 00208 00209 if (!strcmp(xannot[i].label, "Product Name")) 00210 sprintf(xannot[i].descr, "%s", outfilename); 00211 00212 else if (!strcmp(xannot[i].label, "Processing Time")) 00213 sprintf(xannot[i].descr, "%04d%03d%02d%02d%02d%03d", 00214 pyear, ( localtm->tm_yday + 1 ), localtm->tm_hour, localtm->tm_min, 00215 localtm->tm_sec, 0); 00216 00217 else if (!strcmp(xannot[i].label, "Input Files")) 00218 sprintf(xannot[i].descr, "%s", in_oz ); 00219 00220 else if (!strcmp(xannot[i].label, "Processing Control")) 00221 sprintf(xannot[i].descr, "%s %s %s %s", 00222 argv[0], annotfile, in_oz, outdir); 00223 00224 else if (!strcmp(xannot[i].label, "Start Time")) 00225 sprintf(xannot[i].descr, "%04d%03d000000000", year, julval ); 00226 00227 else if (!strcmp(xannot[i].label, "End Time")) 00228 sprintf(xannot[i].descr, "%04d%03d230000000", year, julval ); 00229 00230 else if (!strcmp(xannot[i].label, "Start Year")) 00231 sprintf(xannot[i].descr, "%04d", year); 00232 00233 else if (!strcmp(xannot[i].label, "Start Day")) 00234 sprintf(xannot[i].descr, "%03d", julval); 00235 00236 else if (!strcmp(xannot[i].label, "Start Millisec")) 00237 sprintf(xannot[i].descr, "00000000" ); 00238 00239 else if (!strcmp(xannot[i].label, "End Year")) 00240 sprintf(xannot[i].descr, "%04d", year); 00241 00242 else if (!strcmp(xannot[i].label, "End Day")) 00243 sprintf(xannot[i].descr, "%03d", julval); 00244 00245 else if (!strcmp(xannot[i].label, "End Millisec")) 00246 sprintf(xannot[i].descr, "82800000" ); 00247 00248 /* Extract values from descr array element for program use */ 00249 00250 else if (!strcmp(xannot[i].label, "Northernmost Latitude")) 00251 sscanf(xannot[i].descr, "%f", &nmostlat); 00252 00253 else if (!strcmp(xannot[i].label, "Southernmost Latitude")) 00254 sscanf(xannot[i].descr, "%f", &smostlat); 00255 00256 else if (!strcmp(xannot[i].label, "Westernmost Longitude")) 00257 sscanf(xannot[i].descr, "%f", &wmostlon); 00258 00259 else if (!strcmp(xannot[i].label, "Easternmost Longitude")) 00260 sscanf(xannot[i].descr, "%f", &emostlon); 00261 00262 else if (!strcmp(xannot[i].label, "Latitude Step")) 00263 sscanf(xannot[i].descr, "%f", &latstep); 00264 } 00265 00266 /* 00267 * Create HDF file 00268 */ 00269 00270 if ((result = startHDF(outfile, &sdfid, &fid, DFACC_CREATE)) != 0) 00271 pexit ("Fatal error starting HDF file"); 00272 00273 /* 00274 * Write attribute array to HDF file 00275 */ 00276 00277 if ((result = wrtattr(sdfid, xannot, numannarr)) != 0) pexit ("wrtattr"); 00278 free(xannot); 00279 00280 /* 00281 * -------- Allocate space for 2D data arrays ----------------------- 00282 */ 00283 00284 rank = 2; 00285 00286 shape[0] = OUTLATSZ; /* lat */ 00287 shape[1] = OUTLONSZ; /* lon */ 00288 array_size = shape[0] * shape[1]; 00289 00290 if ( ( int8_SDSdataQC = 00291 ( int8 * ) calloc( array_size, sizeof(int8) ) ) == NULL ) 00292 pexit ("calloc int8_SDSdataQC"); 00293 00294 if ((float_SDS_oz = 00295 (int16 *) malloc (sizeof(int16) * GRIBLATSZ * GRIBLONSZ )) == NULL) 00296 pexit ("malloc float_SDS_oz"); 00297 00298 if ((oz_newsiz = 00299 (int16 *) malloc (sizeof(int16) * array_size )) == NULL) 00300 pexit ("malloc oz_newsiz"); 00301 if( ( seaice8 = (int8 *) malloc( sizeof( int8 ) * array_size ) ) == NULL )\ 00302 pexit ("malloc seaice8" ); 00303 00304 /* 00305 * Process ozone the same way here 00306 */ 00307 l = 0; 00308 for (j = 0; j < GRIBLATSZ; j++) 00309 { 00310 for (k = ( GRIBLONSZ / 2 ); k < GRIBLONSZ; k++) 00311 { 00312 dat = *( datarr1 + j * GRIBLONSZ + k ); 00313 if( dat < 0 || dat > SHRT_MAX ) 00314 float_SDS_oz[l] = 0; 00315 else 00316 float_SDS_oz[l] = (int16)( 100. * dat ); 00317 l++; /* index counter */ 00318 } /* for k */ 00319 00320 for (k = 0; k < (GRIBLONSZ/2); k++) 00321 { /* lons */ 00322 dat = *( datarr1 + j * GRIBLONSZ + k ); 00323 if( dat< 0 || dat > SHRT_MAX ) 00324 float_SDS_oz[l] = 0; 00325 else 00326 float_SDS_oz[l] = (int16)( 100. * dat ); 00327 l++; /* index counter */ 00328 } /* for k */ 00329 } /* for j */ 00330 00331 /* 00332 * interpolate the array to the traditional TOVS ozone size 00333 */ 00334 resize_oz( (short *) float_SDS_oz, GRIBLONSZ, GRIBLATSZ, 00335 OUTLONSZ, OUTLATSZ, (short *) oz_newsiz ); 00336 int_char_min = 0; 00337 int_char_max = 100; 00338 for( j = 0; j < array_size; j++ ) 00339 { 00340 if( *( oz_newsiz + j ) < int_char_min ) 00341 { 00342 *( int8_SDSdataQC + j ) = 1; 00343 *( seaice8 + j ) = -1; 00344 } 00345 else if( *( oz_newsiz + j ) > int_char_max ) 00346 { 00347 *( int8_SDSdataQC + j ) = 1; 00348 *( seaice8 + j ) = 101; 00349 } 00350 else 00351 *( seaice8 + j ) = *( oz_newsiz + j ); 00352 } 00353 /* 00354 * ----------------- Create Vgroup ---------------------------- 00355 */ 00356 00357 gridid = setupGrid(fid, VGROUPNAME); 00358 00359 /* 00360 * ----------------- Write ozone SDS ---------------------------- 00361 */ 00362 00363 sdsname = "seaice"; 00364 strcpy(longname, "Sea ice concentration"); 00365 units = "percent coverage"; 00366 datafmt = "int8"; 00367 datatype = DFNT_INT8; 00368 00369 if ((SDSinFile (sdsname, longname, units, datafmt, 00370 datatype, sdfid, rank, shape, 00371 seaice8, gridid)) != 0) 00372 pexit ("SDSinFile ozone"); 00373 00374 00375 // Smoothed 1/2 degree icemask for Aquarius processing 00376 float frac_ice_smooth[720*360]; 00377 mk_smooth_ice_map_( seaice8, frac_ice_smooth); 00378 00379 sdsname = "seaice_smoothed"; 00380 strcpy(longname, "Sea ice concentration (smoothed)"); 00381 units = "percent coverage"; 00382 datafmt = "float32"; 00383 datatype = DFNT_FLOAT32; 00384 shape_smooth[0] = 360; /* lat */ 00385 shape_smooth[1] = 720; /* lon */ 00386 SDSinFile (sdsname, longname, units, datafmt, 00387 datatype, sdfid, rank, shape_smooth, 00388 frac_ice_smooth, gridid); 00389 00390 free(float_SDS_oz); 00391 00392 /* 00393 * ---- Write QC flag SDS, units attribute, and add to grid Vgroup ----- 00394 */ 00395 00396 sdsname = "seaice_QC"; 00397 strcpy(longname, "Sea ice Q/C flag"); 00398 units = ""; 00399 datafmt = "int8"; 00400 datatype = DFNT_INT8; 00401 00402 if ((SDSinFile (sdsname, longname, units, datafmt, 00403 datatype, sdfid, rank, shape, 00404 int8_SDSdataQC, gridid)) != 0) 00405 pexit ("SDSinFile ozone_QC"); 00406 00407 /* 00408 * free storage for QC array (used for all products) 00409 * also free the input array 00410 */ 00411 00412 free(int8_SDSdataQC); 00413 free( datarr1 ); 00414 free( seaice8 ); 00415 free( oz_newsiz ); 00416 /* 00417 * deattach HDF grid (used for all products) 00418 */ 00419 00420 deattachHDFgrid(gridid); 00421 00422 /* 00423 * close HDF structures 00424 */ 00425 00426 if ((result = closeHDFstructs(sdfid, fid)) != 0) pexit ("closeHDFstructs"); 00427 00428 return 0; 00429 } /* main */ 00430 00431 00432 int8 check_usage(int argc, char *argv[], int *anctyp, int *n_opt_arg, 00433 char *source_name, int *grib_mode, char *grib2_t_str ) 00434 /***************************************************************** 00435 File: check_usage 00436 00437 Purpose: check command line arguments 00438 00439 Description: check command line arguements and report proper usage 00440 on argument count error or no argumnts. 00441 00442 Returns: success (0) or failure(-1) 00443 00444 Parameters: (in calling order) 00445 Type Name I/O Description 00446 ---- ---- --- ----------- 00447 int argc I arg count 00448 char *[] argv I input arg list 00449 int * anctyp O type of anc data, either 00450 ANC_MET or ANC_OZONE 00451 int * n_opt_arg O # option args 00452 char * source_name O An alternate source (as in gbl 00453 attrib source) name 00454 int * grib_mode O Type of grib file the binary 00455 file came from: GRIB_MODE_1 for 00456 a GRIB 1 file, GRIB_MODE_2 for a 00457 GRIB 2 file (NCEP std as of 00458 Jan 2008) 00459 char * grib2_t_str O string with time of the GRIB 2 00460 data - used for setting file 00461 name and attributes 00462 00463 Note: This routine is custom for each program that uses it since 00464 the arguments are different of course. 00465 00466 Author: Brian D. Schieber, GSC, 10/93 00467 00468 Mods: 10/4/93 BDS, new HDF inputs handle all 3 NCEP parms 00469 W. Robinson, GSC, 4 Feb 97 changes for 5th param:rel_hum 00470 W. Robinson, GSC, 18 Jun 98 add flexability to either accept 00471 7 or 8 inputs for met data or 3 or 4 for ozone data 00472 W. Robinson, SAIC, 12 May 2005 have optional argument -s 00473 for an optional source name 00474 W. Robinson, SAIC, 24 Jan 2008 add a grib mode -2 for grib 2 00475 format files 00476 J. Gales, Futuretech, 20 Jan 2010 Add version number 00477 *****************************************************************/ 00478 { 00479 extern char *optarg; 00480 extern int optind, optopt; 00481 int c, opt_source_found, errflg; 00482 00483 /* 00484 * get the options first, one for the source name 00485 */ 00486 *grib_mode = GRIB_MODE_1; 00487 opt_source_found = 0; 00488 errflg = 0; 00489 *n_opt_arg = 0; 00490 while ((c = getopt(argc, argv, ":s:2:")) != -1) 00491 { 00492 switch(c) 00493 { 00494 case 's': 00495 strcpy( source_name, optarg ); 00496 opt_source_found = 1; 00497 break; 00498 case '2': 00499 strcpy( grib2_t_str, optarg ); 00500 *grib_mode = GRIB_MODE_2; 00501 break; 00502 case ':': /* -s without operand */ 00503 fprintf(stderr, 00504 "Option -%c requires an operand\n", optopt); 00505 errflg++; 00506 break; 00507 case '?': 00508 fprintf(stderr, "Unrecognized option: -%c\n", optopt); 00509 errflg++; 00510 } 00511 } 00512 if( !errflg ) 00513 { 00514 *n_opt_arg = optind; 00515 /* 00516 * for regular arguments 00517 */ 00518 if (argc - optind == 3 || argc - optind == 4 ) 00519 { 00520 if( !opt_source_found ) strcpy( source_name, "NCEP" ); 00521 } 00522 else 00523 { 00524 errflg++; 00525 } 00526 } 00527 if( errflg ) 00528 { 00529 printf("\nVersion 1.0 (%s %s)",__DATE__,__TIME__); 00530 printf("\n\nUsage:\n"); 00531 printf("\t%s [-s<sname>] <metadata> <file> [outdir]\n", argv[0]); 00532 printf("\nWhere:\n"); 00533 printf("\tmetadata: DAAC-style HDF metadata (ASCII file)\n"); 00534 printf("\tfile: GRIB converted binary file to process\n"); 00535 printf("\t(this can be either from a GRIB or a GRIB2 file\n" ); 00536 printf("\toutdir: Optional output directory ('.' default)\n"); 00537 printf("\t Options:\n" ); 00538 printf("\t-s <sname> Optional source name to output file with:\n" ); 00539 printf("\t file name form: NYYYYDDD00_SEAICE_<sname>_24h.hdf\n" ); 00540 printf("\t default is NCEP\n" ); 00541 printf( 00542 "\t-2 <GRIB str> Specify the -2 with the <GRIB str> to process\n" ); 00543 printf( 00544 "\t binary files derived from a GRIB 2 format file. The\n" ); 00545 printf( 00546 "\t binary files are derived using wgrib2 with the -bin \n" ); 00547 printf( 00548 "\t option. The <GRIB str> is a required date / time\n" ); 00549 printf( 00550 "\t identification string gotten using wgrib2 -t\n" ); 00551 printf("\n"); 00552 exit(-1); 00553 } 00554 return (0); 00555 }
1.7.6.1