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 <string.h>
  37 #include "gribfuncs.h"
  38
  39 #if defined(_WIN32)
  40 #include <io.h>
  41 #else
  42 #include <unistd.h>
  43 #endif
  44
  45 int rg_fwrite_grib(PDS_INPUT *pds, grid_desc_sec *gds, float **data, FILE *fid)
  46 {
  47
  48   GRIB_HDR *gh=NULL;
  49   DATA_INPUT data_input;
  50   GEOM_IN geom_in;
  51   USER_INPUT user_input;
  52   char errmsg[240];
  53   float *data_one_d;
  54   int status;
  55   int i,j;
  56
  57   if ((gds->head.usData_type == LATLON_PRJ) ||
  58       (gds->head.usData_type == GAUSS_PRJ)) {
  59     if (gds->head.usData_type == GAUSS_PRJ) {
  60       strcpy(geom_in.prjn_name,"gaussian");
  61     } else {
  62       strcpy(geom_in.prjn_name,"spherical");
  63     }
  64     geom_in.nx = gds->llg.usNi;
  65     geom_in.ny = gds->llg.usNj;
  66     geom_in.parm_1 = (gds->llg.iDj)/1000.;
  67     geom_in.parm_2 = (gds->llg.iDi)/1000.;
  68     geom_in.first_lat = gds->llg.lLat1/1000.;
  69     geom_in.first_lon = gds->llg.lLon1/1000.;
  70     geom_in.last_lat = gds->llg.lLat2/1000.;
  71     geom_in.last_lon = gds->llg.lLon2/1000.;
  72     geom_in.scan = gds->llg.usScan_mode;
  73     geom_in.usRes_flag = gds->llg.usRes_flag;
  74   } else if (gds->head.usData_type == LAMB_PRJ) {
  75     strcpy(geom_in.prjn_name,"lambert");
  76     geom_in.nx = gds->lam.iNx;
  77     geom_in.ny = gds->lam.iNy;
  78     geom_in.x_int_dis = gds->lam.ulDx/1000.;
  79     geom_in.y_int_dis = gds->lam.ulDy/1000.;
  80     geom_in.parm_1 = (gds->lam.lLat_cut2)/1000.;
  81     geom_in.parm_2 = (gds->lam.lLat_cut1)/1000.;
  82     geom_in.parm_3 = (gds->lam.lLon_orient)/1000.;
  83     geom_in.first_lat = gds->lam.lLat1/1000.;
  84     geom_in.first_lon = gds->lam.lLon1/1000.;
  85     geom_in.scan = gds->lam.usScan_mode;
  86     geom_in.usRes_flag = gds->lam.usRes_flag;
  87   } else {
  88     fprintf(stderr,"rg_fwrite_grib: invalid input projection %d\n",
  89             gds->head.usData_type);
  90     return -1;
  91   }
  92
  93   data_input.usProc_id = pds->usProc_id;
  94   data_input.usGrid_id = pds->usGrid_id;
  95   data_input.usParm_id = pds->usParm_id;
  96   data_input.usParm_sub_id = pds->usCenter_sub;
  97   data_input.usLevel_id = pds->usLevel_id;
  98   data_input.nLvl_1 = pds->usHeight1;
  99   data_input.nLvl_2 = pds->usHeight2;
 100   data_input.nYear = pds->usYear + (pds->usCentury-1)*100;
 101   data_input.nMonth = pds->usMonth;
 102   data_input.nDay = pds->usDay;
 103   data_input.nHour = pds->usHour;
 104   data_input.nMinute = pds->usMinute;
 105   data_input.nSecond = 0;
 106   data_input.usFcst_id = pds->usFcst_unit_id;
 107   data_input.usFcst_per1 = pds->usP1;
 108   data_input.usFcst_per2 = pds->usP2;
 109   data_input.usTime_range_id = pds->usTime_range;
 110   data_input.usTime_range_avg = pds->usTime_range_avg;
 111   data_input.usTime_range_mis = pds->usTime_range_mis;
 112   /* We add an extra digit here, because the grib library seems to cut off
 113    *   one more than I would prefer.
 114    */
 115   data_input.nDec_sc_fctr = pds->sDec_sc_fctr+1;
 116
 117   user_input.chCase_id='0';
 118   user_input.usParm_tbl=pds->usParm_tbl;
 119   user_input.usSub_tbl=pds->usSub_tbl;
 120   /*
 121   user_input.usCenter_id=190;
 122   */
 123   user_input.usCenter_id = pds->usCenter_id;
 124   user_input.usCenter_sub=pds->usCenter_sub;
 125   user_input.usTrack_num=0;
 126   user_input.usGds_bms_id = 128;
 127   user_input.usBDS_flag=0;
 128   user_input.usBit_pack_num=0;
 129
 130
 131   data_one_d = (float *)calloc(geom_in.nx*geom_in.ny,sizeof(float));
 132   if (data_one_d == NULL) {
 133     fprintf(stderr,"rg_fwrite_grib: could not allocate space for data_one_d\n");
 134     return -1;
 135   }
 136
 137   for (i=0; i<geom_in.nx; i++) {
 138     for (j=0; j<geom_in.ny; j++) {
 139       data_one_d[i+j*geom_in.nx] = data[j][i];
 140     }
 141   }
 142
 143   status = init_gribhdr(&gh,errmsg);
 144   if (status != 0) {
 145     fprintf (stderr,"rg_fwrite_grib: %s",errmsg);
 146     return -1;
 147   }
 148
 149   status = grib_enc(data_input,user_input,geom_in,data_one_d,gh,errmsg);
 150   if (status != 0) {
 151     fprintf (stderr,"rg_fwrite_grib: %s",errmsg);
 152     return -1;
 153   }
 154
 155   status = gribhdr2file(gh,fid,errmsg);
 156   if (status != 0) {
 157     fprintf (stderr,"rg_fwrite_grib: %s",errmsg);
 158     return -1;
 159   }
 160
 161   free_gribhdr(&gh);
 162
 163   return 1;
 164
 165 }
 166
 167 int rg_write_grib(PDS_INPUT *pds, grid_desc_sec *gds, char filename[],
 168                float **data)
 169 {
 170   char tmpfile[240];
 171   char tmpstring[240];
 172   FILE *fid;
 173   int status;
 174
 175   sprintf(tmpfile,"/tmp/tmpgribfile_%d",getpid());
 176   fid = fopen(tmpfile,"wb");
 177   if (fid == NULL) {
 178     fprintf(stderr,"rg_write_grib: Could not open %s\n",tmpfile);
 179     return -1;
 180   }
 181
 182   status = rg_fwrite_grib(pds,gds,data,fid);
 183   if (status != 1)
 184     {
 185       fprintf(stderr,"rg_write_grib: rg_fwrite_grib failed\n");
 186       return -1;
 187     }
 188
 189   /* append tmpfile to filename */
 190   sprintf(tmpstring,"cat %s >> %s",tmpfile,filename);
 191   system(tmpstring);
 192   unlink(tmpfile);
 193
 194   fclose(fid);
 195
 196   return(1);
 197 }