external/io_grib1/grib1_util/write_grib.c

   1 /*
   2  *****************************************************************************
   3  * rg_write_grib - function which encapsulates parts of the MEL grib library
   4  *              writing routines.
   5  *
   6  * Todd Hutchinson
   7  * TASC
   8  * tahutchinson@tasc.com
   9  * (781)942-2000 x3108
  10  * 7/1/99
  11  *
  12  * Interface:
  13  *   Input:
  14  *     pds -      a pointer to a structure which stores information about the
  15  *                grib product definition section.
  16  *     gds -      a pointer to a structure which stores information about the
  17  *                grib grid definition section.
  18  *     filename - This is the output grib file which will be created if
  19  *                it does not already exist.  If this file already exists,
  20  *                grib records will be appended to it.
  21  *     data -     a 2-d array holding the values at grid points for the grid
  22  *                that is to be output.
  23  *   Return:
  24  *     1 for success
  25  *    -1 for failure
  26  *
  27  *   Caveats:   This function only supports a "latlon" grid
  28  *              currently, this is equivalent to a cylindrical equidistant
  29  *              projection.
  30  *
  31  ***************************************************************************
  32  */
  33
  34 #include <stdio.h>
  35 #include <stdlib.h>
  36 #include "gribfuncs.h"
  37
  38 int rg_write_grib(PDS_INPUT *pds, grid_desc_sec *gds, char filename[],
  39                float **data)
  40 {
  41   char tmpfile[240];
  42   char tmpstring[240];
  43   FILE *fid;
  44   int status;
  45
  46   sprintf(tmpfile,"/tmp/tmpgribfile_%d",getpid());
  47   fid = fopen(tmpfile,"wb");
  48   if (fid == NULL) {
  49     fprintf(stderr,"rg_write_grib: Could not open %s\n",tmpfile);
  50     return -1;
  51   }
  52
  53   status = rg_fwrite_grib(pds,gds,data,fid);
  54   if (status != 1)
  55     {
  56       fprintf(stderr,"rg_write_grib: rg_fwrite_grib failed\n");
  57       return -1;
  58     }
  59
  60   /* append tmpfile to filename */
  61   sprintf(tmpstring,"cat %s >> %s",tmpfile,filename);
  62   system(tmpstring);
  63   unlink(tmpfile);
  64
  65   close(fid);
  66
  67   return(1);
  68 }
  69
  70 int rg_fwrite_grib(PDS_INPUT *pds, grid_desc_sec *gds, float **data, FILE *fid)
  71 {
  72
  73   GRIB_HDR *gh=NULL;
  74   DATA_INPUT data_input;
  75   GEOM_IN geom_in;
  76   USER_INPUT user_input;
  77   char errmsg[240];
  78   float *data_one_d;
  79   int status;
  80   int i,j;
  81
  82   if ((gds->head.usData_type == LATLON_PRJ) ||
  83       (gds->head.usData_type == GAUSS_PRJ)) {
  84     if (gds->head.usData_type == GAUSS_PRJ) {
  85       strcpy(geom_in.prjn_name,"gaussian");
  86     } else {
  87       strcpy(geom_in.prjn_name,"spherical");
  88     }
  89     geom_in.nx = gds->llg.usNi;
  90     geom_in.ny = gds->llg.usNj;
  91     geom_in.parm_1 = (gds->llg.iDj)/1000.;
  92     geom_in.parm_2 = (gds->llg.iDi)/1000.;
  93     geom_in.first_lat = gds->llg.lLat1/1000.;
  94     geom_in.first_lon = gds->llg.lLon1/1000.;
  95     geom_in.last_lat = gds->llg.lLat2/1000.;
  96     geom_in.last_lon = gds->llg.lLon2/1000.;
  97     geom_in.scan = gds->llg.usScan_mode;
  98     geom_in.usRes_flag = gds->llg.usRes_flag;
  99   } else if (gds->head.usData_type == LAMB_PRJ) {
 100     strcpy(geom_in.prjn_name,"lambert");
 101     geom_in.nx = gds->lam.iNx;
 102     geom_in.ny = gds->lam.iNy;
 103     geom_in.x_int_dis = gds->lam.ulDx/1000.;
 104     geom_in.y_int_dis = gds->lam.ulDy/1000.;
 105     geom_in.parm_1 = (gds->lam.lLat_cut2)/1000.;
 106     geom_in.parm_2 = (gds->lam.lLat_cut1)/1000.;
 107     geom_in.parm_3 = (gds->lam.lLon_orient)/1000.;
 108     geom_in.first_lat = gds->lam.lLat1/1000.;
 109     geom_in.first_lon = gds->lam.lLon1/1000.;
 110     geom_in.scan = gds->lam.usScan_mode;
 111     geom_in.usRes_flag = gds->lam.usRes_flag;
 112   } else {
 113     fprintf(stderr,"rg_fwrite_grib: invalid input projection %d\n",
 114             gds->head.usData_type);
 115     return -1;
 116   }
 117
 118   data_input.usProc_id = pds->usProc_id;
 119   data_input.usGrid_id = pds->usGrid_id;
 120   data_input.usParm_id = pds->usParm_id;
 121   data_input.usParm_sub_id = pds->usCenter_sub;
 122   data_input.usLevel_id = pds->usLevel_id;
 123   data_input.nLvl_1 = pds->usHeight1;
 124   data_input.nLvl_2 = pds->usHeight2;
 125   data_input.nYear = pds->usYear + (pds->usCentury-1)*100;
 126   data_input.nMonth = pds->usMonth;
 127   data_input.nDay = pds->usDay;
 128   data_input.nHour = pds->usHour;
 129   data_input.nMinute = pds->usMinute;
 130   data_input.nSecond = 0;
 131   data_input.usFcst_id = pds->usFcst_unit_id;
 132   data_input.usFcst_per1 = pds->usP1;
 133   data_input.usFcst_per2 = pds->usP2;
 134   data_input.usTime_range_id = pds->usTime_range;
 135   data_input.usTime_range_avg = pds->usTime_range_avg;
 136   data_input.usTime_range_mis = pds->usTime_range_mis;
 137   /* We add an extra digit here, because the grib library seems to cut off
 138    *   one more than I would prefer.
 139    */
 140   data_input.nDec_sc_fctr = pds->sDec_sc_fctr+1;
 141
 142   user_input.chCase_id='0';
 143   user_input.usParm_tbl=pds->usParm_tbl;
 144   user_input.usSub_tbl=pds->usSub_tbl;
 145   /*
 146   user_input.usCenter_id=190;
 147   */
 148   user_input.usCenter_id = pds->usCenter_id;
 149   user_input.usCenter_sub=pds->usCenter_sub;
 150   user_input.usTrack_num=0;
 151   user_input.usGds_bms_id = 128;
 152   user_input.usBDS_flag=0;
 153   user_input.usBit_pack_num=0;
 154
 155
 156   data_one_d = (float *)calloc(geom_in.nx*geom_in.ny,sizeof(float));
 157   if (data_one_d == NULL) {
 158     fprintf(stderr,"rg_fwrite_grib: could not allocate space for data_one_d\n");
 159     return -1;
 160   }
 161
 162   for (i=0; i<geom_in.nx; i++) {
 163     for (j=0; j<geom_in.ny; j++) {
 164       data_one_d[i+j*geom_in.nx] = data[j][i];
 165     }
 166   }
 167
 168   status = init_gribhdr(&gh,errmsg);
 169   if (status != 0) {
 170     fprintf (stderr,"rg_fwrite_grib: %s",errmsg);
 171     return -1;
 172   }
 173
 174   status = grib_enc(data_input,user_input,geom_in,data_one_d,gh,errmsg);
 175   if (status != 0) {
 176     fprintf (stderr,"rg_fwrite_grib: %s",errmsg);
 177     return -1;
 178   }
 179
 180   status = gribhdr2file(gh,fid,errmsg);
 181   if (status != 0) {
 182     fprintf (stderr,"rg_fwrite_grib: %s",errmsg);
 183     return -1;
 184   }
 185
 186   free_gribhdr(&gh);
 187
 188   return 1;
 189
 190 }