/* Program : make_grib_log (was printer.c) Programmer : Todd J. Kienitz, SAIC Date : January 10, 1996 Purpose : To produce the information file output of the GRIB message. Revisions : 04/17/96 Steve Lowe, SAIC: modified data print-out 04/22/96 Alice Nakajima (ATN), SAIC: added BMS summary 12/12/96 ATN: implement combined Decoder/Encdoer structs replaced (table2 tab2[], table3 tab3[], tables mgotab); 06/14/97 ATN: print upto encoded pricision. 02/22/98 ATN: replace projection id with constants, add printing for prjns: Rotated Lat/Lon, Stretched Lat/Lon, Stretched Rotated Lat/Lon, Rotated Gaussian, Stretched Gauss, Stretched Rotated Gaussian , Oblique Lambert, and for Albers equal-area. 09/10/98 ATN: extension flag printing. */ #include #include #include #include "grib_lookup.h" /* combined encoder/decoder structs */ #include "dprints.h" /* for debug printing */ #include "gribfuncs.h" /* prototypes */ /* ********************************************************************** * A. FUNCTION: make_grib_log * Produces debug file GRIB.log from the GRIB message in the Grib Header * * INTERFACE: * int make_grib_log (input_fn, lookup_fn, msg_length, offset, * pds, gds, bds, bms, grib_data, errmsg) * * ARGUMENTS (I=input, O=output, I&O=input and output): * (I) char *input_fn; name of input GRIB file * (I) char *lookup_fn; name of Lookup file, nil if not used * (I) unsigned long msg_length; total length of GRIB message * (I) long offset; starting location of GRIB message in bytes * (I) PDS_INPUT pds; product definition section structure * (I) grid_desc_sec gds; grid description section structure * (I) BDS_HEAD_INPUT bds; binary data section header structure * (I) BMS_INPUT bms; bit map definition section structure * (I) float *grib_data; array of decoded data * * ACCESSES GLOBAL VARS: * int UseTables; * set to one if lookup table used * CTRS_DEFN db_ctr_tbl[NCTRS]; * predefined array holding Originating Center info * PARM_DEFN db_parm_tbl [MAX_PARM_TBLS * NPARM]; * predefined arr of Parameter info * LVL_DEFN db_lvl_tbl [NLVL]; * predefined arr of Level info struct * MODEL_DEFN db_mdl_tbl [NMODEL]; * predefined arr of Model info struct * GEOM_DEFN db_geom_tbl [NGEOM]; * predefined arr of Geometry info struct * * RETURN CODE: * 0> no errors; file GRIB.log has been created; * 1> error, errmsg filled; ********************************************************************** */ extern int UseTables; /* set means use lookup tbl defns */ extern PARM_DEFN db_parm_tbl[]; /* parameter conversion info */ extern MODEL_DEFN db_mdl_tbl[]; /* model conversion info */ extern LVL_DEFN db_lvl_tbl[]; /* level conversion info */ extern GEOM_DEFN db_geom_tbl[]; /* Geom conversion info */ extern CTR_DEFN db_ctr_tbl[]; /* Ctr conversion info */ #if PROTOTYPE_NEEDED int make_grib_log ( char *input_fn, char *lookup_fn, unsigned long msg_length, long offset, PDS_INPUT pds, grid_desc_sec gds, BDS_HEAD_INPUT bds, BMS_INPUT bms, float *grib_data, char *errmsg) #else int make_grib_log (input_fn, lookup_fn, msg_length, offset, pds, gds, bds, bms, grib_data,errmsg) char *input_fn; char *lookup_fn; unsigned long msg_length; long offset; PDS_INPUT pds; grid_desc_sec gds; BDS_HEAD_INPUT bds; BMS_INPUT bms; float *grib_data; char *errmsg; #endif { char *func="make_grib_log"; int i, indx, k, fd, numpts=100; float dsf, res, min, max; FILE *fp; /* * * A.0 DEBUG printing */ DPRINT1 ("Entering %s\n", func); /* * * A.1 OPEN file "GRIB.log" in APPEND mode */ fp=fopen ("GRIB.log", "a+"); if (!fp) { DPRINT1("%s: failed to open 'GRIB.log' for appending, skip logfile\n", func); sprintf (errmsg, "%s: failed to open 'GRIB.log'\n", func); return (1); } /* * * A.2 WRITE Indicator Section information to file * !message length * !GRIB Edition number */ fseek(fp, 0L, 2); if (ftell(fp) == 0L) fprintf (fp, "%s: InFile= %s\n%s: Lookup=%d, fn='%s'\n\n" , func, input_fn, func,UseTables, lookup_fn); fprintf (fp, "**** VALID MESSAGE FOUND AT %ld BYTES ****\n" , offset); fprintf(fp, "\n********* SECTION 0 IDS *********\n" ); fprintf(fp, "Total Message length = %ld\nEdition Number = %d\n", msg_length, pds.usEd_num); /* * * A.3 WRITE Product Definition Section information to file * !Section length * !Parameter Table version * !Parameter Sub-Table version if defined and flagged by Extension flag * !Tracking id if defined and flagged by Extension flag */ fprintf(fp, "\n********* SECTION 1 PDS *********\n" \ "Section length = %d\nTable version = %d\n", pds.uslength, pds.usParm_tbl); if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG ) { if (pds.usSub_tbl != 0) fprintf(fp,"Local Table version = %d\n",pds.usSub_tbl); if(pds.usTrack_num != 0) fprintf(fp,"Tracking ID = %d\n",pds.usTrack_num); } /* * !Originating Center id * !IF (using tables) Name of Originating Center */ fprintf(fp,"Originating Center id = %d\n",pds.usCenter_id); if (UseTables) if ( db_ctr_tbl[pds.usCenter_id].ctr_dsc[0] ) fprintf(fp,"Originating Center = %s\n", db_ctr_tbl[pds.usCenter_id].ctr_dsc); else fprintf(fp,"Originating Center ID %d not defined in current table.\n", pds.usCenter_id); /* * !Sub-Table Entry for Originating Center if non-zero and if * !extension flag is set */ if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG && pds.usCenter_sub != 0) fprintf(fp,"Sub-Table Entry Originating Center = %d\n",pds.usCenter_sub); /* * !Extension flag */ fprintf(fp,"Extension flag = %d (extensions %s)\n", pds.usExt_flag, (pds.usExt_flag == (unsigned short)EXTENSION_FLAG ? "used" : "not used")); /* * !Model Identification * !IF (using tables) Model Description */ fprintf(fp,"Model id = %d\n",pds.usProc_id); if (UseTables) if ( db_mdl_tbl[pds.usProc_id].grib_dsc[0] ) fprintf(fp,"Model Description = %s\n", db_mdl_tbl[pds.usProc_id].grib_dsc); else fprintf(fp,"Model ID %d not defined in current table.\n", pds.usProc_id); /* * !Grid Identification * !IF (using tables) Grid Description */ fprintf(fp,"Grid id = %d\n",pds.usGrid_id); if (UseTables) if ( db_geom_tbl[pds.usGrid_id].grib_dsc[0] ) fprintf(fp,"Grid Description = %s\n", db_geom_tbl[pds.usGrid_id].grib_dsc); else fprintf(fp,"Grid ID %d not defined in current table.\n", pds.usGrid_id); /* * !Parameter Identification */ fprintf(fp,"Parameter id = %d\n",pds.usParm_id); /* * !IF (usExt_flag is set AND * ! (Parm id between 250 and 254) AND (Sub Parm ID defined))) * ! PRINT Parm_sub * !ENDIF */ if (pds.usExt_flag == (unsigned short)EXTENSION_FLAG && pds.usParm_id>=250 && pds.usParm_id<=254 && pds.usParm_sub!=0) fprintf(fp,"Parameter sub-id = %d\n",pds.usParm_sub); /* * !IF (using lookup table) THEN * ! CALCULATE index in Parm Conversion Array to use * ! let index= (usParm_Id - 249)*256 + usParm_sub; * ! * ! IF this index in Parm Conversion Array is defined THEN * ! PRINT its grib_dsc and grib_unit_dsc * ! ELSE * ! PRINT it's not defined mesage * ! ENDIF * !ENDIF */ if(UseTables) { indx = PARMTBL_INDX (pds.usParm_id, pds.usParm_sub); if ( db_parm_tbl[indx].grib_dsc[0] ) { fprintf(fp,"Parameter name = %s\n",db_parm_tbl[indx].grib_dsc); fprintf(fp,"Parameter units = %s\n",db_parm_tbl[indx].grib_unit_dsc); } else fprintf(fp,"Parameter ID %d not defined in current table.\n", pds.usParm_id); } /* * !Level Id * !IF (using tables) * ! Level description * ! SWITCH (number of octets to store Height1) * ! 2: Level = Height1 * ! 1: Bottom of Layer = Height1 * ! Top of Layer = Height2 * ! 0: (no Height value required) * ! default: (corrupt table entry or message) * ! ENDSWITCH * !ELSE (not using tables) * ! Level = Height1 (Level assumed) * !ENDIF */ fprintf(fp,"Level_type = %d\n",pds.usLevel_id); if(UseTables) { if ( db_lvl_tbl[pds.usLevel_id].grib_dsc[0] ) { fprintf(fp,"Level description = %s\n", db_lvl_tbl[pds.usLevel_id].grib_dsc); switch(db_lvl_tbl[pds.usLevel_id].num_octets){ case 2: fprintf(fp,"%s = %u\n", db_lvl_tbl[pds.usLevel_id].lvl_name_1, pds.usHeight1); break; case 1: fprintf(fp,"%s = %u\n%s = %u\n", db_lvl_tbl[pds.usLevel_id].lvl_name_1, pds.usHeight1, db_lvl_tbl[pds.usLevel_id].lvl_name_2, pds.usHeight2); break; case 0: break; default: fprintf(fp,"***Number of octets for table 3 undefined - possibly " "corrupt dataset.***\n"); } }else fprintf(fp,"Level ID %d not defined in current table.\n", pds.usLevel_id); } /* end UseTables 'if' statement */ else fprintf(fp,"Level = %u\n",pds.usHeight1); /* * !Reference Date/Time: * ! Century * ! Year * ! Month * ! Day * ! Hour * ! Minute * ! Second if defined */ fprintf(fp, "Reference Date/Time of Data Set:\n" \ " Century = %d\n Year = %d\n Month = %d\n Day = %d\n"\ " Hour = %d\n Minute = %d\n", pds.usCentury,pds.usYear,pds.usMonth,pds.usDay,pds.usHour,pds.usMinute); if(pds.usExt_flag == (unsigned short)EXTENSION_FLAG) fprintf(fp," Second = %d\n",pds.usSecond); /* * !Forecast Time Unit * ! Forecast Period 1 * ! Forecast Period 2 */ switch(pds.usFcst_unit_id){ case 0: fprintf(fp,"Forecast Time Unit = Minute\n"); break; case 1: fprintf(fp,"Forecast Time Unit = Hour\n"); break; case 2: fprintf(fp,"Forecast Time Unit = Day\n"); break; case 3: fprintf(fp,"Forecast Time Unit = Month\n"); break; case 4: fprintf(fp,"Forecast Time Unit = Year\n"); break; case 5: fprintf(fp,"Forecast Time Unit = Decade (10 years)\n"); break; case 6: fprintf(fp,"Forecast Time Unit = Normal (30 years)\n"); break; case 7: fprintf(fp,"Forecast Time Unit = Century (100 years)\n"); break; case 254: fprintf(fp,"Forecast Time Unit = Second\n"); break; default: fprintf(fp,"Forecast Time Unit = UNDEFINED!!\n"); } fprintf(fp," Forecast Period 1 = %d\n",pds.usP1); fprintf(fp," Forecast Period 2 = %d\n",pds.usP2); /* * !Time Range Indicator * !Number in Average * !Number Missing */ fprintf(fp,"Time Range = %d\n",pds.usTime_range); fprintf(fp,"Number in Average = %d\n",pds.usTime_range_avg); fprintf(fp,"Number Missing = %d\n",pds.usTime_range_mis); /* * !Decimal Scale Factor */ fprintf(fp,"Decimal Scale Factor = %d\n",pds.sDec_sc_fctr); /* * * A.4 IF (GDS included) THEN * A.4.1 WRITE Grid Definition Section information to file * !Section length * !Parm_nv * !Parm_pv_pl * !Data type */ if(pds.usGds_bms_id >> 7 & 1) { fprintf(fp,"\n********* SECTION 2 GDS *********\n"); fprintf(fp,"Section length = %d\n",gds.head.uslength); fprintf(fp,"Parm_nv = %d\n",gds.head.usNum_v); fprintf(fp,"Parm_pv_pl = %d\n",gds.head.usPl_Pv); fprintf(fp,"Data_type = %d\n",gds.head.usData_type); /* * A.4.2 SWITCH (Data Type, Table 6) * ! For each Data Type, write the following to file: * ! Number of points along rows/columns of grid * ! Reference Lat/Lon information * ! Resolution and Component Flags (Table 7) * ! Direction increments if given * ! Assumption of Earth shape * ! U&V component orientation * ! Scanning mode flags (Table 8) * Default: Projection not supported, exit; */ fprintf(fp,"Projection = %s\n", prjn_name[gds.head.usData_type]); switch(gds.head.usData_type) { /* * Case 0: Lat/Lon projection * Case 10: Rotated Lat/Lon projection * Case 20: Stretched Lat/Lon projection * Case 30: Stretched Rotated Lat/Lon projection */ case LATLON_PRJ: /* Lat/Lon Grid */ case ROT_LATLON_PRJ: /* Rotated Lat/Lon */ case STR_LATLON_PRJ: /* Stretched Lat/Lon */ case STR_ROT_LATLON_PRJ : /* Stretched and Rotated Lat/Lon */ fprintf(fp,"Number of points along a parallel = %d\n",gds.llg.usNi); fprintf(fp,"Number of points along a meridian = %d\n",gds.llg.usNj); fprintf(fp,"Latitude of first grid point = %.3f deg\n", ((float)gds.llg.lLat1)/1000.); fprintf(fp,"Longitude of first grid point = %.3f deg\n", ((float)gds.llg.lLon1)/1000.); fprintf(fp,"Latitude of last grid point = %.3f deg\n", ((float)gds.llg.lLat2)/1000.); fprintf(fp,"Longitude of last grid point = %.3f deg\n", ((float)gds.llg.lLon2)/1000.); fprintf(fp,"Resolution and Component Flags: \n"); if ((gds.llg.usRes_flag >> 7) & 1) { fprintf(fp," Longitudinal increment = %f deg\n", ((float)gds.llg.iDi)/1000.); fprintf(fp," Latitudinal increment = %f deg\n", ((float)gds.llg.iDj)/1000.); }else fprintf(fp," Direction increments not given.\n"); if ((gds.llg.usRes_flag >> 6) & 1) fprintf(fp," Earth assumed oblate spherical.\n"); else fprintf(fp," Earth assumed spherical.\n"); if ((gds.llg.usRes_flag >> 3) & 1) fprintf(fp," U&V components resolved relative to +I and " "+J\n"); else fprintf(fp," U&V components resolved relative to east " "and north.\n"); fprintf(fp,"Scanning Mode Flags: \n"); if ((gds.llg.usScan_mode >> 7) & 1) fprintf(fp," Points scan in -I direction.\n"); else fprintf(fp," Points scan in +I direction.\n"); if ((gds.llg.usScan_mode >> 6) & 1) fprintf(fp," Points scan in +J direction.\n"); else fprintf(fp," Points scan in -J direction.\n"); if ((gds.llg.usScan_mode >> 5) & 1) fprintf(fp," Adjacent points in J direction are " "consecutive.\n"); else fprintf(fp," Adjacent points in I direction are " "consecutive.\n"); /* added 02/98 This code pertains only to the Stretch/Rotate grids, so skip over if it's a LATLON_PRJN type; */ if (gds.head.usData_type != LATLON_PRJ) { fprintf(fp," Latitude of southern pole = %.3f deg\n", ((float)gds.llg.lLat_southpole)/1000.); fprintf(fp," Longitude of southern pole = %.3f deg\n", ((float)gds.llg.lLon_southpole)/1000.); /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation a single precision floating point value */ fprintf(fp," Angle of rotation = %.3f\n", grib_ibm_local ((unsigned long) gds.llg.lRotate)); fprintf(fp," Latitude of pole of stretching = %.3f deg\n", ((float)gds.llg.lPole_lat)/1000.); fprintf(fp," Longitude of pole of stretching = %.3f deg\n", ((float)gds.llg.lPole_lon)/1000.); /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation a single precision floating point value */ fprintf(fp," Stretching factor = %.3f\n", grib_ibm_local ((unsigned long)gds.llg.lStretch)); } break; /* * Case 1: Mercator Projection */ case MERC_PRJ: /* Mercator Projection Grid */ fprintf(fp,"Number of points along a parallel = %d\n",gds.merc.cols); fprintf(fp,"Number of points along a meridian = %d\n",gds.merc.rows); fprintf(fp,"Latitude of first grid point = %.3f deg\n", ((float)gds.merc.first_lat)/1000.); fprintf(fp,"Longitude of first grid point = %.3f deg\n", ((float)gds.merc.first_lon)/1000.); fprintf(fp,"Latitude of last grid point = %.3f deg\n", ((float)gds.merc.La2)/1000.); fprintf(fp,"Longitude of last grid point = %.3f deg\n", ((float)gds.merc.Lo2)/1000.); fprintf(fp,"Latitude of intersection with Earth = %.3f deg\n", ((float)gds.merc.latin)/1000.); fprintf(fp,"Resolution and Component Flags: \n"); if ((gds.merc.usRes_flag >> 7) & 1) { fprintf(fp," Longitudinal increment = %f deg\n", ((float)gds.merc.lon_inc)/1000.); fprintf(fp," Latitudinal increment = %f deg\n", ((float)gds.merc.lat_inc)/1000.); }else fprintf(fp," Direction increments not given.\n"); if ((gds.merc.usRes_flag >> 6) & 1) fprintf(fp," Earth assumed oblate spherical.\n"); else fprintf(fp," Earth assumed spherical.\n"); if ((gds.merc.usRes_flag >> 3) & 1) fprintf(fp," U&V components resolved relative to +I and " "+J\n"); else fprintf(fp," U&V components resolved relative to east " "and north.\n"); fprintf(fp,"Scanning Mode Flags: \n"); if ((gds.merc.usScan_mode >> 7) & 1) fprintf(fp," Points scan in -I direction.\n"); else fprintf(fp," Points scan in +I direction.\n"); if ((gds.merc.usScan_mode >> 6) & 1) fprintf(fp," Points scan in +J direction.\n"); else fprintf(fp," Points scan in -J direction.\n"); if ((gds.merc.usScan_mode >> 5) & 1) fprintf(fp," Adjacent points in J direction are " "consecutive.\n"); else fprintf(fp," Adjacent points in I direction are " "consecutive.\n"); break; /* * Case 3: Lambert Conformal Projection * Case 13: Oblique Lambert Conformal Projection * Case 8: Alberts equal-area secant/tangent conic/bipolar Prj */ case LAMB_PRJ: /* Lambert Conformal */ case OBLIQ_LAMB_PRJ: /* Oblique Lambert Conformal */ case ALBERS_PRJ: /* Albers equal-area */ fprintf(fp,"Number of points along X-axis = %d\n",gds.lam.iNx); fprintf(fp,"Number of points along Y-axis = %d\n",gds.lam.iNy); fprintf(fp,"Latitude of first grid point = %.3f deg\n", ((float)gds.lam.lLat1)/1000.); fprintf(fp,"Longitude of first grid point = %.3f deg\n", ((float)gds.lam.lLon1)/1000.); fprintf(fp,"Orientation of grid = %.3f deg\n", ((float)gds.lam.lLon_orient)/1000.); fprintf(fp,"First Latitude Cut = %.3f deg\n", ((float)gds.lam.lLat_cut1)/1000.); fprintf(fp,"Second Latitude Cut = %.3f deg\n", ((float)gds.lam.lLat_cut2)/1000.); fprintf(fp,"Resolution and Component Flags: \n"); if ((gds.lam.usRes_flag >> 7) & 1) { fprintf(fp," X-direction increment = %d meters\n", gds.lam.ulDx); fprintf(fp," Y-direction increment = %d meters\n", gds.lam.ulDy); }else fprintf(fp," Direction increments not given.\n"); if ((gds.lam.usRes_flag >> 6) & 1) fprintf(fp," Earth assumed oblate spherical.\n"); else fprintf(fp," Earth assumed spherical.\n"); if ((gds.lam.usRes_flag >> 3) & 1) fprintf(fp," U&V components resolved relative to +I and " "+J\n"); else fprintf(fp," U&V components resolved relative to east " "and north.\n"); fprintf(fp,"Scanning Mode Flags: \n"); if ((gds.lam.usScan_mode >> 7) & 1) fprintf(fp," Points scan in -I direction.\n"); else fprintf(fp," Points scan in +I direction.\n"); if ((gds.lam.usScan_mode >> 6) & 1) fprintf(fp," Points scan in +J direction.\n"); else fprintf(fp," Points scan in -J direction.\n"); if ((gds.lam.usScan_mode >> 5) & 1) fprintf(fp," Adjacent points in J direction are " "consecutive.\n"); else fprintf(fp," Adjacent points in I direction are " "consecutive.\n"); /* 02/98 This code pertains only to the Albers projection */ if (gds.head.usData_type == ALBERS_PRJ) { fprintf(fp," Latitude of the southern pole = %.3f\n", ((float)gds.lam.lLat_southpole)/1000.); fprintf(fp," Longitude of the southern pole = %.3f\n", ((float)gds.lam.lLon_southpole)/1000.); } break; /* * Case 4: Gaussian Lat/Lon Projection * Case 14: Rotated Gaussian Lat/Lon Projection * Case 24: Stretched Gaussian Lat/Lon Projection * Case 34: Stretched Rotated Gaussian Lat/Lon Projection */ case GAUSS_PRJ: /* Gaussian Latitude/Longitude Grid */ case ROT_GAUSS_PRJ: /* Rotated Gaussian */ case STR_GAUSS_PRJ : /* Stretched Gaussian */ case STR_ROT_GAUSS_PRJ : /* Stretched and Rotated Gaussian */ fprintf(fp,"Number of points along a parallel = %d\n",gds.llg.usNi); fprintf(fp,"Number of points along a meridian = %d\n",gds.llg.usNj); fprintf(fp,"Latitude of first grid point = %.3f deg\n", ((float)gds.llg.lLat1)/1000.); fprintf(fp,"Longitude of first grid point = %.3f deg\n", ((float)gds.llg.lLon1)/1000.); fprintf(fp,"Latitude of last grid point = %.3f deg\n", ((float)gds.llg.lLat2)/1000.); fprintf(fp,"Longitude of last grid point = %.3f deg\n", ((float)gds.llg.lLon2)/1000.); fprintf(fp,"Resolution and Component Flags: \n"); if ((gds.llg.usRes_flag >> 7) & 1) { fprintf(fp," i direction increment = %f deg\n", ((float)gds.llg.iDi)/1000.); fprintf(fp, " Number of parallels between pole and equator = %d\n", gds.llg.iDj); }else fprintf(fp," Direction increments not given.\n"); if ((gds.llg.usRes_flag >> 6) & 1) fprintf(fp," Earth assumed oblate spherical.\n"); else fprintf(fp," Earth assumed spherical.\n"); if ((gds.llg.usRes_flag >> 3) & 1) fprintf(fp," U&V components resolved relative to +I and " "+J\n"); else fprintf(fp," U&V components resolved relative to east " "and north.\n"); fprintf(fp,"Scanning Mode Flags: \n"); if ((gds.llg.usScan_mode >> 7) & 1) fprintf(fp," Points scan in -I direction.\n"); else fprintf(fp," Points scan in +I direction.\n"); if ((gds.llg.usScan_mode >> 6) & 1) fprintf(fp," Points scan in +J direction.\n"); else fprintf(fp," Points scan in -J direction.\n"); if ((gds.llg.usScan_mode >> 5) & 1) fprintf(fp," Adjacent points in J direction are " "consecutive.\n"); else fprintf(fp," Adjacent points in I direction are " "consecutive.\n"); /* added 02/98 This code pertains only to the Stretch/Rotate grids */ if (gds.head.usData_type != GAUSS_PRJ) { fprintf(fp," Latitude of southern pole = %.3f deg\n", ((float)gds.llg.lLat_southpole)/1000.); fprintf(fp," Longitude of southern pole = %.3f deg\n", ((float)gds.llg.lLon_southpole)/1000.); /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation a single precision floating point value */ fprintf(fp," Angle of rotation = %.3f\n", grib_ibm_local ((unsigned long) gds.llg.lRotate)); fprintf(fp," Latitude of pole of stretching = %.3f deg\n", (float)gds.llg.lPole_lat); fprintf(fp," Longitude of pole of stretching = %.3f deg\n", (float)gds.llg.lPole_lon); /* conv from 'saaaaaaa bbbbbbbb bbbbbbbb bbbbbbbb' representation a single precision floating point value */ fprintf(fp," Stretching factor = %.3f\n", grib_ibm_local ((unsigned long)gds.llg.lStretch)); } break; /* * Case 5: Polar Sterographic Projection */ case POLAR_PRJ: /* Polar Stereographic Projection Grid */ fprintf(fp,"Number of points along X-axis = %d\n",gds.pol.usNx); fprintf(fp,"Number of points along Y-axis = %d\n",gds.pol.usNy); fprintf(fp,"Latitude of first grid point = %.3f deg\n", ((float)gds.pol.lLat1)/1000.); fprintf(fp,"Longitude of first grid point = %.3f deg\n", ((float)gds.pol.lLon1)/1000.); fprintf(fp,"Orientation of grid = %.3f deg\n", ((float)gds.pol.lLon_orient)/1000.); fprintf(fp,"Projection Center: "); if ((gds.pol.usProj_flag >> 7) & 1) fprintf(fp,"South Pole\n"); else fprintf(fp,"North Pole\n"); fprintf(fp,"Resolution and Component Flags: \n"); if ((gds.pol.usRes_flag >> 7) & 1) { fprintf(fp," X-direction grid length = %d meters\n",gds.pol.ulDx); fprintf(fp," Y-direction grid length = %d meters\n",gds.pol.ulDy); }else fprintf(fp," Direction increments not given.\n"); if ((gds.pol.usRes_flag >> 6) & 1) fprintf(fp," Earth assumed oblate spherical.\n"); else fprintf(fp," Earth assumed spherical.\n"); if ((gds.pol.usRes_flag >> 3) & 1) fprintf(fp," U&V components resolved relative to +I and " "+J\n"); else fprintf(fp," U&V components resolved relative to east " "and north.\n"); fprintf(fp,"Scanning Mode Flags: \n"); if ((gds.pol.usScan_mode >> 7) & 1) fprintf(fp," Points scan in -I direction.\n"); else fprintf(fp," Points scan in +I direction.\n"); if ((gds.pol.usScan_mode >> 6) & 1) fprintf(fp," Points scan in +J direction.\n"); else fprintf(fp," Points scan in -J direction.\n"); if ((gds.pol.usScan_mode >> 5) & 1) fprintf(fp," Adjacent points in J direction are " "consecutive.\n"); else fprintf(fp," Adjacent points in I direction are " "consecutive.\n"); break; default: /* Bad projection: ignore & continue */ fprintf(stdout, "\n\n***%s WARNING ***:\nProjection %d is INVALID;\n\n", func, gds.head.usData_type); fprintf(fp,"================================================\n"\ "%s: projection %d is not currently implemented\n"\ "================================================\n", func, gds.head.usData_type); break; /* * A.4.2 ENDSWITCH (Data Type) */ } /* Switch */ } /* gds included */ /* * * A.4 ELSE * PRINT no Gds message * A.4 ENDIF */ else fprintf(fp,"\n******* NO SECTION 2 GDS *********\n" ); /* * * A.5 IF (Bitmap Section is present) * THEN * WRITE Bitmap Section information to file * ELSE * PRINT no bms mesg * ENDIF */ if(pds.usGds_bms_id >> 6 & 1) { fprintf(fp,"\n********* SECTION 3 BMS **********\n" ); fprintf(fp,"Section length = %ld\n", bms.uslength); if (bms.uslength <= 6) fprintf(fp,"Bitmap is predefined (Not in message).\n"); else fprintf(fp,"Bitmap is included with message.\n"); fprintf(fp,"Bitmap ID = %d \n", bms.usBMS_id); fprintf(fp,"Number of unused bits = %d\n", bms.usUnused_bits); fprintf(fp,"Number of datapoints set = %ld\n", bms.ulbits_set); }else{ fprintf(fp,"\n******* NO SECTION 3 BMS *********\n" ); } /* * * A.6 WRITE out Binary Data Section Information to file * !Section Length */ fprintf(fp,"\n********* SECTION 4 BDS *********\n" ); fprintf(fp,"Section length = %ld\n",bds.length); /* * !Table 11 Flags */ fprintf(fp,"Table 11 Flags:\n"); if ((bds.usBDS_flag >> 7) & 1) fprintf(fp," Spherical harmonic coefficients.\n"); else fprintf(fp," Grid-point data.\n"); if ((bds.usBDS_flag >> 6) & 1) fprintf(fp," Second-order packing.\n"); else fprintf(fp," Simple Packing.\n"); if ((bds.usBDS_flag >> 5) & 1) fprintf(fp," Integer values.\n"); else fprintf(fp," Floating point values.\n"); if ((bds.usBDS_flag >> 4) & 1) fprintf(fp," Octet 14 contains additional flag bits.\n"); else fprintf(fp," No additional flags at octet 14.\n"); /* * !Decimal Scale Factor (Repeated from PDS) */ fprintf(fp,"\nDecimal Scale Factor = %d\n",pds.sDec_sc_fctr); /* * !Binary Scale Factor * !Bit Width * !Number of Data Points */ fprintf(fp,"Binary scale factor = %d\n", bds.Bin_sc_fctr); fprintf(fp,"Bit width = %d\n", bds.usBit_pack_num); fprintf(fp,"Number of data points = %ld\n",bds.ulGrid_size); /* * A.6.1 WRITE Data Summary to file * !Compute Data Min/Max and Resolution */ dsf = (float) pow( (double) 10, (double) pds.sDec_sc_fctr); res = (float) pow((double)2,(double)bds.Bin_sc_fctr) / dsf; min = bds.fReference / dsf; max = (float) (pow((double)2, (double)bds.usBit_pack_num) - 1); max = min + max * res; fprintf(fp,"Data Minimum = %f\n", min ); fprintf(fp,"Data Maximum = %f\n", max ); fprintf(fp,"Resolution = %f\n",res ); /* * !Compute Format Specifier for printing Data */ fd = (int) -1 * (float) log10((double) res) + .5; if (fd <= 0) { fd = 0; fprintf(fp,"DATA will be displayed as integers (res > 0.1).\n"); } /* * !WRITE First 100 Data Points to file up to Encoded Precision */ if (bds.ulGrid_size > 1) { if (bds.ulGrid_size < 100) numpts = bds.ulGrid_size; fprintf(fp,"\nDATA ARRAY: (first %d)\n",numpts); if (fd > 0) { for (i=0; i