grb_write.c 8.53 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_LIBGRIB

#include "dmemory.h"
#include "cdi.h"
#include "cdi_int.h"
#include "stream_cgribex.h"
#include "stream_grb.h"
#include "stream_gribapi.h"
#include "file.h"
#include "cgribex.h"  /* gribZip gribGetZip gribGinfo */
#include "namespace.h"


static
size_t grbEncode(int filetype, int memtype, int varID, int levelID, int vlistID, int gridID, int zaxisID,
		 int date, int time, int tsteptype, int numavg,
21
		 size_t datasize, const void *data, int nmiss, void **gribbuffer,
22
23
24
25
26
		 int comptype, void *gribContainer)
{
  size_t nbytes = 0;

#ifdef HAVE_LIBCGRIBEX
27
  if ( filetype == CDI_FILETYPE_GRB )
28
29
    {
      size_t gribbuffersize = datasize*4+3000;
30
      *gribbuffer = Malloc(gribbuffersize);
31
32
33

      nbytes = cgribexEncode(memtype, varID, levelID, vlistID, gridID, zaxisID,
			     date, time, tsteptype, numavg,
34
			     (long) datasize, data, nmiss, *gribbuffer, gribbuffersize);
35
36
37
38
39
    }
  else
#endif
#ifdef HAVE_LIBGRIB_API
    {
40
41
42
43
44
45
46
47
      const void *datap = data;
      if ( memtype == MEMTYPE_FLOAT )
        {
          const float *dataf = (const float*) data;
          double *datad = (double*) Malloc(datasize*sizeof(double));
          for ( size_t i = 0; i < datasize; ++i ) datad[i] = (double) dataf[i];
          datap = (const void*) datad;
        }
48
49
50
51

      size_t gribbuffersize;
      nbytes = gribapiEncode(varID, levelID, vlistID, gridID, zaxisID,
			     date, time, tsteptype, numavg,
52
			     (long) datasize, datap, nmiss, gribbuffer, &gribbuffersize,
53
			     comptype, gribContainer);
54
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
      if ( memtype == MEMTYPE_FLOAT ) Free((void*)datap);
56
57
    }
#else
58
59
60
61
62
    {
      Error("GRIB_API support not compiled in!");
      (void)gribContainer;
      (void)comptype;
    }
63
64
65
66
67
68
#endif

  return nbytes;
}

static
69
size_t grbSzip(int filetype, void *gribbuffer, size_t gribbuffersize)
70
71
{
  size_t buffersize = gribbuffersize + 1000; /* compressed record can be greater than source record */
72
  void *buffer = Malloc(buffersize);
73
74
75

  /*  memcpy(buffer, gribbuffer, gribbuffersize); */

76
  size_t nbytes = 0;
77
  if ( filetype == CDI_FILETYPE_GRB )
78
    {
79
      nbytes = (size_t)gribZip((unsigned char *)gribbuffer, (long) gribbuffersize, (unsigned char *)buffer, (long) buffersize);
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    }
  else
    {
      static int lszip_warn = 1;
      if ( lszip_warn ) Warning("Szip compression of GRIB2 records not implemented!");
      lszip_warn = 0;
      nbytes = gribbuffersize;
    }

  Free(buffer);

  return nbytes;
}


void grbCopyRecord(stream_t * streamptr2, stream_t * streamptr1)
{
  int filetype = streamptr1->filetype;
  int fileID1 = streamptr1->fileID;
  int fileID2 = streamptr2->fileID;
  int tsID    = streamptr1->curTsID;
  int vrecID  = streamptr1->tsteps[tsID].curRecID;
  int recID   = streamptr1->tsteps[tsID].recIDs[vrecID];
  off_t recpos  = streamptr1->tsteps[tsID].records[recID].position;
  size_t recsize = streamptr1->tsteps[tsID].records[recID].size;

  fileSetPos(fileID1, recpos, SEEK_SET);

  /* round up recsize to next multiple of 8 */
  size_t gribbuffersize = ((recsize + 7U) & ~7U);

  unsigned char *gribbuffer = (unsigned char *) Malloc(gribbuffersize);

  if (fileRead(fileID1, gribbuffer, recsize) != recsize)
    Error("Could not read GRIB record for copying!");

  size_t nbytes = recsize;

118
119
120
121
122
123
124
125
126
127
128
129
130
131
  if ( filetype == CDI_FILETYPE_GRB )
    {
      extern int GribChangeParameterIdentification;
      if ( GribChangeParameterIdentification == 1 )
        {
          extern int GribChangeParID_code;
          extern int GribChangeParID_ltype;
          extern int GribChangeParID_lev;
          // Even if you are stream-copy records you might need to change a bit of grib-header !
#if defined HAVE_LIBCGRIBEX
          void *gh = cgribex_handle_new_from_meassage((void*) gribbuffer, recsize);
          cgribexChangeParameterIdentification(gh, GribChangeParID_code, GribChangeParID_ltype, GribChangeParID_lev);
          cgribex_handle_delete(gh);
#elif defined HAVE_LIBGRIB_API
132
          void *gh = (void*)grib_handle_new_from_message(NULL, (void*) gribbuffer, recsize);
133
134
135
136
137
138
139
          gribapiChangeParameterIdentification(gh, GribChangeParID_code, GribChangeParID_ltype, GribChangeParID_lev);
          grib_handle_delete(gh);
#endif
          GribChangeParameterIdentification = -1; // after grib attributes have been changed turn it off again
        }
    }

140
  if ( filetype == CDI_FILETYPE_GRB )
141
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
143
      size_t unzipsize;
      int izip = gribGetZip(recsize, gribbuffer, &unzipsize);
144

145
      if ( izip == 0 && streamptr2->comptype == CDI_COMPRESS_SZIP )
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
          nbytes = grbSzip(filetype, gribbuffer, nbytes);
    }

  while ( nbytes & 7 ) gribbuffer[nbytes++] = 0;

  size_t nwrite = fileWrite(fileID2, gribbuffer, nbytes);
  if ( nwrite != nbytes )
    {
      perror(__func__);
      Error("Could not write record for copying!");
    }

  Free(gribbuffer);
}


void grb_write_var_slice(stream_t *streamptr, int varID, int levelID, int memtype, const void *data, int nmiss)
{
164
  void *gribbuffer = NULL;
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
  void *gc = NULL;

  int filetype  = streamptr->filetype;
  int fileID    = streamptr->fileID;
  int vlistID   = streamptr->vlistID;
  int gridID    = vlistInqVarGrid(vlistID, varID);
  int zaxisID   = vlistInqVarZaxis(vlistID, varID);
  int tsteptype = vlistInqVarTsteptype(vlistID, varID);
  int comptype  = streamptr->comptype;
  int tsID      = streamptr->curTsID;
  int date      = streamptr->tsteps[tsID].taxis.vdate;
  int time      = streamptr->tsteps[tsID].taxis.vtime;
  int numavg    = 0;
  if ( vlistInqVarTimave(vlistID, varID) )
    numavg = streamptr->tsteps[tsID].taxis.numavg;

  if ( CDI_Debug )
    Message("gridID = %d zaxisID = %d", gridID, zaxisID);

  size_t datasize = (size_t)gridInqSize(gridID);
185

186
#ifdef HAVE_LIBCGRIBEX
187
  if ( filetype == CDI_FILETYPE_GRB )
188
189
190
191
192
193
194
195
196
197
198
199
200
201
    {
    }
  else
#endif
    {
#ifdef GRIBCONTAINER2D
      gribContainer_t **gribContainers =  (gribContainer_t **) streamptr->gribContainers;
      gc = (void *) &gribContainers[varID][levelID];
#else
      gribContainer_t *gribContainers =  (gribContainer_t *) streamptr->gribContainers;
      gc = (void *) &gribContainers[varID];
#endif
    }

202
  if ( comptype != CDI_COMPRESS_JPEG && comptype != CDI_COMPRESS_SZIP ) comptype = CDI_COMPRESS_NONE;
203

204
  if ( filetype == CDI_FILETYPE_GRB && comptype == CDI_COMPRESS_JPEG )
205
206
207
208
209
210
211
    {
      static int ljpeg_warn = 1;
      if ( ljpeg_warn ) Warning("JPEG compression of GRIB1 records not available!");
      ljpeg_warn = 0;
    }

  size_t nbytes = grbEncode(filetype, memtype, varID, levelID, vlistID, gridID, zaxisID, date, time, tsteptype, numavg,
212
                            datasize, data, nmiss, &gribbuffer, comptype, gc);
213

214
  if ( filetype == CDI_FILETYPE_GRB && streamptr->comptype == CDI_COMPRESS_SZIP )
215
216
    nbytes = grbSzip(filetype, gribbuffer, nbytes);

217
  size_t (*myFileWrite)(int fileID, const void *restrict buffer, size_t len, int tsID)
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
    = (size_t (*)(int, const void *restrict, size_t, int))
    namespaceSwitchGet(NSSWITCH_FILE_WRITE).func;
  size_t nwrite = myFileWrite(fileID, gribbuffer, nbytes, tsID);

  if ( nwrite != nbytes )
    {
      perror(__func__);
      Error("Failed to write GRIB slice!");
    }

  if ( gribbuffer ) Free(gribbuffer);
}


void grb_write_var(stream_t *streamptr, int varID, int memtype, const void *data, int nmiss)
{
  int vlistID  = streamptr->vlistID,
    gridID   = vlistInqVarGrid(vlistID, varID),
    gridsize = gridInqSize(gridID),
    zaxisID  = vlistInqVarZaxis(vlistID, varID),
    nlevs    = zaxisInqSize(zaxisID);
  double missval = vlistInqVarMissval(vlistID, varID);

  size_t chunkLen = (size_t)gridsize;
  if ( memtype == MEMTYPE_FLOAT )
    for ( int levelID = 0; levelID < nlevs; levelID++ )
      {
        const float *restrict fdata = ((const float *)data)+levelID*gridsize;
246
247
248
249
250
251
        
        int nmiss_slice = 0;
        if ( nmiss )
          for ( size_t i = 0; i < chunkLen; ++i )
            nmiss_slice += DBL_IS_EQUAL(fdata[i], missval);

252
253
254
255
256
257
        grb_write_var_slice(streamptr, varID, levelID, memtype, fdata, nmiss_slice);
      }
  else
    for ( int levelID = 0; levelID < nlevs; levelID++ )
      {
        const double *restrict ddata = ((const double *)data)+levelID*gridsize;
258
259
260
261
262
263
        
        int nmiss_slice = 0;
        if ( nmiss )
          for ( size_t i = 0; i < chunkLen; ++i )
            nmiss_slice += DBL_IS_EQUAL(ddata[i], missval);

264
265
266
267
268
        grb_write_var_slice(streamptr, varID, levelID, memtype, ddata, nmiss_slice);
      }
}


269
void grb_write_record(stream_t *streamptr, int memtype, const void *data, int nmiss)
270
271
272
273
274
275
276
{
  int varID   = streamptr->record->varID;
  int levelID = streamptr->record->levelID;

  grb_write_var_slice(streamptr, varID, levelID, memtype, data, nmiss);
}

277

278
#endif