NCEPLIBS-g2c  1.8.0
g2_addfield.c
Go to the documentation of this file.
1 
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include "grib2_int.h"
10 
98 g2int
99 g2_addfield(unsigned char *cgrib, g2int ipdsnum, g2int *ipdstmpl,
100  float *coordlist, g2int numcoord, g2int idrsnum, g2int *idrstmpl,
101  float *fld, g2int ngrdpts, g2int ibmap, g2int *bmap)
102 {
103  unsigned char *cpack;
104  static g2int zero = 0, one = 1, four = 4, five = 5, six = 6, seven = 7;
105  const g2int minsize = 50000;
106  g2int iofst, ibeg, lencurr, len, nsize;
107  g2int ilen, isecnum, i, nbits, temp, left;
108  g2int ibmprev, j, lcpack, ioctet, newlen, ndpts;
109  g2int lensec4, lensec5, lensec6, lensec7;
110  g2int issec3 = 0, isprevbmap = 0, lpos3 = 0, JJ, KK, MM;
111  g2int *coordieee;
112  float *pfld;
113  gtemplate *mappds, *mapdrs;
114 #if defined USE_PNG || defined USE_JPEG2000 || defined USE_OPENJPEG || defined USE_AEC
115  unsigned int allones = 4294967295u;
116  g2int width, height, iscan, itemp;
117 #endif
118  int ret;
119 
120  /* Check for GRIB header and terminator. Translate the error codes
121  * to the legacy G2 error codes. */
122  if ((ret = g2c_check_msg(cgrib, &lencurr, 1)))
123  {
124  if (ret == G2C_ENOTGRIB)
125  return G2_ADD_MSG_INIT;
126  if (ret == G2C_EMSGCOMPLETE)
127  return G2_ADD_MSG_COMPLETE;
128  }
129 
130  /* Loop through all current sections of the GRIB message to find
131  * the last section number. */
132  len = 16; /* length of Section 0 */
133  for (;;)
134  {
135  /* Get number and length of next section. */
136  iofst = len * 8;
137  gbit(cgrib, &ilen, iofst, 32);
138  iofst = iofst + 32;
139  gbit(cgrib, &isecnum, iofst, 8);
140  iofst = iofst + 8;
141 
142  /* Check if previous Section 3 exists. */
143  if (isecnum == 3)
144  {
145  issec3 = 1;
146  lpos3 = len;
147  }
148  /* Check if a previous defined bitmap exists. */
149  if (isecnum == 6)
150  {
151  gbit(cgrib, &ibmprev, iofst, 8);
152  iofst = iofst + 8;
153  if (ibmprev >= 0 && ibmprev <= 253)
154  isprevbmap = 1;
155  }
156  len = len + ilen;
157 
158  /* Exit loop if last section reached. */
159  if (len == lencurr)
160  break;
161 
162  /* If byte count for each section doesn't match current */
163  /* total length, then there is a problem. */
164  if (len > lencurr)
165  {
166  printf("g2_addfield: Section byte counts don''t add to total.\n");
167  printf("g2_addfield: Sum of section byte counts = %ld\n", len);
168  printf("g2_addfield: Total byte count in Section 0 = %ld\n", lencurr);
169  return G2_BAD_SEC_COUNTS;
170  }
171  }
172 
173  /* Sections 4 through 7 can only be added after section 3 or 7. */
174  if (isecnum != 3 && isecnum != 7)
175  {
176  printf("g2_addfield: Sections 4-7 can only be added after Section 3 or 7.\n");
177  printf("g2_addfield: Section %ld was the last found in given GRIB message.\n",
178  isecnum);
179  return G2_BAD_SEC;
180  }
181  else if (!issec3)
182  {
183  /* Sections 4 through 7 can only be added if section 3 was previously defined. */
184  printf("g2_addfield: Sections 4-7 can only be added if Section 3 was previously included.\n");
185  printf("g2_addfield: Section 3 was not found in given GRIB message.\n");
186  printf("g2_addfield: Call to routine addgrid required to specify Grid definition.\n");
187  return G2_ADDFIELD_BAD_GDS;
188  }
189 
190  /* Add Section 4 - Product Definition Section. */
191  ibeg = lencurr * 8; /* Calculate offset for beginning of section 4 */
192  iofst = ibeg + 32; /* leave space for length of section */
193  sbit(cgrib, &four, iofst, 8); /* Store section number (4) */
194  iofst = iofst + 8;
195  sbit(cgrib, &numcoord, iofst, 16); /* Store num of coordinate values */
196  iofst = iofst + 16;
197  sbit(cgrib, &ipdsnum, iofst, 16); /* Store Prod Def Template num. */
198  iofst = iofst + 16;
199 
200  /* Get Product Definition Template. */
201  if (!(mappds = getpdstemplate(ipdsnum)))
202  return G2_ADDFIELD_BAD_PDT;
203 
204  /* Extend the Product Definition Template, if necessary. The */
205  /* number of values in a specific template may vary depending on */
206  /* data specified in the "static" part of the template. */
207  if (mappds->needext)
208  {
209  free(mappds);
210  mappds = extpdstemplate(ipdsnum, ipdstmpl);
211  }
212 
213  /* Pack up each input value in array ipdstmpl into the the */
214  /* appropriate number of octets, which are specified in */
215  /* corresponding entries in array mappds. */
216  for (i = 0; i < mappds->maplen; i++)
217  {
218  nbits = abs(mappds->map[i]) * 8;
219  if ((mappds->map[i] >= 0) || (ipdstmpl[i] >= 0))
220  sbit(cgrib, ipdstmpl + i, iofst, nbits);
221  else
222  {
223  sbit(cgrib, &one, iofst, 1);
224  temp = abs(ipdstmpl[i]);
225  sbit(cgrib, &temp, iofst + 1, nbits-1);
226  }
227  iofst = iofst + nbits;
228  }
229 
230  /* Pack template extension, if appropriate. */
231  j = mappds->maplen;
232  if (mappds->needext && (mappds->extlen > 0))
233  {
234  for (i = 0; i < mappds->extlen; i++)
235  {
236  nbits = abs(mappds->ext[i]) * 8;
237  if (mappds->ext[i] >= 0 || ipdstmpl[j] >= 0)
238  sbit(cgrib, ipdstmpl + j, iofst, nbits);
239  else
240  {
241  sbit(cgrib, &one, iofst, 1);
242  temp = abs(ipdstmpl[j]);
243  sbit(cgrib, &temp, iofst + 1, nbits - 1);
244  }
245  iofst = iofst + nbits;
246  j++;
247  }
248  }
249  if (mappds->ext)
250  free(mappds->ext);
251  free(mappds);
252 
253  /* Add Optional list of vertical coordinate values after the */
254  /* Product Definition Template, if necessary. */
255  if (numcoord != 0)
256  {
257  coordieee = calloc(numcoord, sizeof(g2int));
258  mkieee(coordlist, coordieee, numcoord);
259  sbits(cgrib, coordieee, iofst, 32, 0, numcoord);
260  iofst = iofst + (32 * numcoord);
261  free(coordieee);
262  }
263 
264  /* Calculate length of section 4 and store it in octets 1-4 of */
265  /* section 4. */
266  lensec4 = (iofst - ibeg) / 8;
267  sbit(cgrib, &lensec4, ibeg, 32);
268 
269  /* Pack Data using appropriate algorithm Get Data Representation */
270  /* Template */
271  if (!(mapdrs = getdrstemplate(idrsnum)))
272  return G2_ADDFIELD_BAD_PDT;
273 
274  /* Contract data field, removing data at invalid grid points, if */
275  /* bit-map is provided with field. */
276  if (ibmap == 0 || ibmap == 254)
277  {
278  pfld = malloc(ngrdpts * sizeof(float));
279  ndpts = 0;
280  for (j = 0; j < ngrdpts; j++)
281  {
282  if (bmap[j] == 1)
283  pfld[ndpts++] = fld[j];
284  }
285  }
286  else
287  {
288  ndpts = ngrdpts;
289  pfld = fld;
290  }
291 
292  /* Allocate storage for the packed data. */
293  nsize = ndpts * 4;
294  if (nsize < minsize)
295  nsize = minsize;
296  cpack = malloc(nsize);
297 
298  /* Call packing function based on idrsnum. */
299  if (idrsnum == 0) /* Simple Packing */
300  simpack(pfld, ndpts, idrstmpl, cpack, &lcpack);
301  else if (idrsnum == 2 || idrsnum == 3) /* Complex Packing */
302  cmplxpack(pfld, ndpts, idrsnum, idrstmpl, cpack, &lcpack);
303  else if (idrsnum == 50) /* Sperical Harmonic Simple Packing */
304  {
305  simpack(pfld + 1, ndpts - 1, idrstmpl, cpack, &lcpack);
306  mkieee(pfld, idrstmpl + 4, 1); /* ensure RE(0, 0) value is IEEE format */
307  }
308  else if (idrsnum == 51) /* Sperical Harmonic Complex Packing */
309  {
310  getpoly(cgrib + lpos3, &JJ, &KK, &MM);
311  if (JJ != 0 && KK != 0 && MM != 0)
312  specpack(pfld, ndpts, JJ, KK, MM, idrstmpl, cpack, &lcpack);
313  else
314  {
315  printf("g2_addfield: Cannot pack DRT 5.51.\n");
316  return G2_ADDFIELD_BAD_GDT;
317  }
318  }
319 #if defined USE_JPEG2000 || defined USE_OPENJPEG
320  else if (idrsnum == 40 || idrsnum == 40000)
321  { /* JPEG2000 encoding */
322  if (ibmap == 255)
323  {
324  getdim(cgrib + lpos3, &width, &height, &iscan);
325  if (width == 0 || height == 0)
326  {
327  width = ndpts;
328  height = 1;
329  }
330  else if (width == allones || height == allones)
331  {
332  width = ndpts;
333  height = 1;
334  }
335  else if ((iscan & 32) == 32)
336  { /* Scanning mode: bit 3 */
337  itemp = width;
338  width = height;
339  height = itemp;
340  }
341  }
342  else
343  {
344  width = ndpts;
345  height = 1;
346  }
347  lcpack = nsize;
348  jpcpack(pfld, width, height, idrstmpl, cpack, &lcpack);
349  }
350 #endif /* USE_JPEG2000 */
351 #ifdef USE_PNG
352  else if (idrsnum == 41 || idrsnum == 40010)
353  { /* PNG encoding */
354  if (ibmap == 255)
355  {
356  getdim(cgrib + lpos3, &width, &height, &iscan);
357  if (width == 0 || height == 0)
358  {
359  width = ndpts;
360  height = 1;
361  }
362  else if (width == allones || height == allones)
363  {
364  width = ndpts;
365  height = 1;
366  }
367  else if ((iscan & 32) == 32)
368  { /* Scanning mode: bit 3 */
369  itemp = width;
370  width = height;
371  height = itemp;
372  }
373  }
374  else
375  {
376  width = ndpts;
377  height = 1;
378  }
379  pngpack(pfld, width, height, idrstmpl, cpack, &lcpack);
380  }
381 #endif /* USE_PNG */
382 #ifdef USE_AEC
383  else if (idrsnum == 42)
384  { /* AEC Encoding */
385  if (ibmap == 255)
386  {
387  getdim(cgrib + lpos3, &width, &height, &iscan);
388  if (width == 0 || height == 0)
389  {
390  width = ndpts;
391  height = 1;
392  }
393  else if (width == allones || height == allones)
394  {
395  width = ndpts;
396  height = 1;
397  }
398  else if ((iscan & 32) == 32)
399  { /* Scanning mode: bit 3 */
400  itemp = width;
401  width = height;
402  height = itemp;
403  }
404  }
405  else
406  {
407  width = ndpts;
408  height = 1;
409  }
410  //printf("IN G2C, in g2_addfield, using AEC compression...\n");
411  lcpack = nsize;
412  aecpack(pfld, width, height, idrstmpl, cpack, &lcpack);
413  }
414 #endif /* USE_AEC */
415  else
416  {
417  printf("g2_addfield: Data Representation Template 5.%ld not yet implemented.\n", idrsnum);
418  return G2_ADDFIELD_BAD_DRT;
419  }
420 
421  /* Free temp space. */
422  if (ibmap == 0 || ibmap == 254)
423  {
424  if (fld != pfld)
425  free(pfld);
426  }
427 
428  /* The packing functions return an lcpack of -1 if there was an
429  * error packing. */
430  if (lcpack < 0)
431  {
432  if (cpack)
433  free(cpack);
434  return G2_ADDFIELD_ERR;
435  }
436 
437  /* Add Section 5 - Data Representation Section */
438  ibeg = iofst; /* Calculate offset for beginning of section 5 */
439  iofst = ibeg + 32; /* leave space for length of section */
440  sbit(cgrib, &five, iofst, 8); /* Store section number (5) */
441  iofst = iofst + 8;
442  sbit(cgrib, &ndpts, iofst, 32); /* Store num of actual data points */
443  iofst = iofst + 32;
444  sbit(cgrib, &idrsnum, iofst, 16); /* Store Data Repr. Template num. */
445  iofst = iofst + 16;
446 
447  /* Pack up each input value in array idrstmpl into the */
448  /* the appropriate number of octets, which are specified in */
449  /* corresponding entries in array mapdrs. */
450  for (i = 0; i < mapdrs->maplen; i++)
451  {
452  nbits = abs(mapdrs->map[i]) * 8;
453  if (mapdrs->map[i] >= 0 || idrstmpl[i] >= 0)
454  sbit(cgrib, idrstmpl + i, iofst, nbits);
455  else
456  {
457  sbit(cgrib, &one, iofst, 1);
458  temp = abs(idrstmpl[i]);
459  sbit(cgrib, &temp, iofst + 1, nbits - 1);
460  }
461  iofst = iofst + nbits;
462  }
463  free(mapdrs);
464 
465  /* Calculate length of section 5 and store it in octets */
466  /* 1-4 of section 5. */
467  lensec5 = (iofst - ibeg) / 8;
468  sbit(cgrib, &lensec5, ibeg, 32);
469 
470  /* Add Section 6 - Bit-Map Section */
471  ibeg = iofst; /* Calculate offset for beginning of section 6 */
472  iofst = ibeg + 32; /* leave space for length of section */
473  sbit(cgrib, &six, iofst, 8); /* Store section number (6) */
474  iofst = iofst + 8;
475  sbit(cgrib, &ibmap, iofst, 8); /* Store Bit Map indicator */
476  iofst = iofst + 8;
477 
478  /* Store bitmap, if supplied */
479  if (ibmap == 0)
480  {
481  sbits(cgrib, bmap, iofst, 1, 0, ngrdpts); /* Store BitMap */
482  iofst = iofst + ngrdpts;
483  }
484 
485  /* If specifying a previously defined bit-map, make sure */
486  /* one already exists in the current GRIB message. */
487  if (ibmap == 254 && !isprevbmap)
488  {
489  printf("g2_addfield: Requested previously defined bitmap,");
490  printf(" but one does not exist in the current GRIB message.\n");
491  return G2_ADDFIELD_BAD_BITMAP;
492  }
493 
494  /* Calculate length of section 6 and store it in octets */
495  /* 1-4 of section 6. Pad to end of octect, if necessary. */
496  left = 8 - (iofst % 8);
497  if (left != 8)
498  {
499  sbit(cgrib, &zero, iofst, left); /* Pad with zeros to fill Octet */
500  iofst = iofst + left;
501  }
502  lensec6 = (iofst - ibeg) / 8;
503  sbit(cgrib, &lensec6, ibeg, 32);
504 
505  /* Add Section 7 - Data Section */
506  ibeg = iofst; /* Calculate offset for beginning of section 7 */
507  iofst = ibeg + 32; /* leave space for length of section */
508  sbit(cgrib, &seven, iofst, 8); /* Store section number (7) */
509  iofst = iofst + 8;
510 
511  /* Store Packed Binary Data values, if non-constant field. */
512  if (lcpack != 0)
513  {
514  ioctet = iofst / 8;
515  /*cgrib(ioctet + 1:ioctet + lcpack)=cpack(1:lcpack) */
516  for (j = 0; j < lcpack; j++)
517  cgrib[ioctet + j] = cpack[j];
518  iofst = iofst + (8 * lcpack);
519  }
520 
521  /* Calculate length of section 7 and store it in octets */
522  /* 1-4 of section 7. */
523  lensec7 = (iofst - ibeg) / 8;
524  sbit(cgrib, &lensec7, ibeg, 32);
525 
526  if (cpack)
527  free(cpack);
528 
529  /* Update current byte total of message in Section 0 */
530  newlen = lencurr + lensec4 + lensec5 + lensec6 + lensec7;
531  sbit(cgrib, &newlen, 96, 32);
532 
533  return newlen;
534 }
void aecpack(float *fld, g2int width, g2int height, g2int *idrstmpl, unsigned char *cpack, g2int *lcpack)
This function packs up a float array into a AEC code stream.
Definition: aecpack.c:297
void cmplxpack(float *fld, g2int ndpts, g2int idrsnum, g2int *idrstmpl, unsigned char *cpack, g2int *lcpack)
This subroutine packs up a data field using a complex packing algorithm as defined in the GRIB2 docum...
Definition: cmplxpack.c:36
gtemplate * getdrstemplate(g2int number)
This subroutine returns DRS template information for a specified Data Representation Template.
Definition: drstemplates.c:167
g2int g2_addfield(unsigned char *cgrib, g2int ipdsnum, g2int *ipdstmpl, float *coordlist, g2int numcoord, g2int idrsnum, g2int *idrstmpl, float *fld, g2int ngrdpts, g2int ibmap, g2int *bmap)
This routine packs up Sections 4 through 7 for a given field and adds them to a GRIB2 message.
Definition: g2_addfield.c:99
void gbit(unsigned char *in, g2int *iout, g2int iskip, g2int nbits)
Get bits - unpack bits: Extract arbitrary size values from a packed bit string, right justifying each...
Definition: gbits.c:20
void sbits(unsigned char *out, g2int *in, g2int iskip, g2int nbits, g2int nskip, g2int n)
Store bits - put arbitrary size values into a packed bit string, taking the low order bits from each ...
Definition: gbits.c:180
void sbit(unsigned char *out, g2int *in, g2int iskip, g2int nbits)
Store bits - put arbitrary size values into a packed bit string, taking the low order bits from each ...
Definition: gbits.c:38
g2int getdim(unsigned char *csec3, g2int *width, g2int *height, g2int *iscan)
This subroutine returns the dimensions and scanning mode of a grid definition packed in GRIB2 Grid De...
Definition: getdim.c:29
g2int getpoly(unsigned char *csec3, g2int *jj, g2int *kk, g2int *mm)
This subroutine returns the J, K, and M pentagonal resolution parameters specified in a GRIB Grid Def...
Definition: getpoly.c:40
#define G2_BAD_SEC_COUNTS
Sum of Section byte counts doesn't add to total byte count.
Definition: grib2.h:475
#define G2_ADDFIELD_BAD_GDS
In g2_addfield() section 3 (GDS) not previously defined in message.
Definition: grib2.h:477
#define G2_BAD_SEC
Previous Section was unexpected.
Definition: grib2.h:463
#define G2_ADDFIELD_BAD_BITMAP
In g2_addfield() no bitmap in the GRIB message.
Definition: grib2.h:479
#define G2_ADDFIELD_BAD_DRT
In g2_addfield() unsupported Data Representationi Template.
Definition: grib2.h:478
#define G2_ADD_MSG_COMPLETE
GRIB message already complete.
Definition: grib2.h:474
#define G2_ADDFIELD_BAD_PDT
In g2_addfield() could not find requested Product Definition Template.
Definition: grib2.h:476
#define G2_ADD_MSG_INIT
GRIB message was not initialized - call g2_create() first.
Definition: grib2.h:473
#define G2C_ENOTGRIB
GRIB header not found.
Definition: grib2.h:493
#define G2_ADDFIELD_BAD_GDT
In g2_addfield() GDT of one of 5.50 through 5.53 required when using DRT 5.51.
Definition: grib2.h:480
#define G2C_EMSGCOMPLETE
GRIB message already complete.
Definition: grib2.h:494
#define G2_ADDFIELD_ERR
In g2_addfield() error packing data field.
Definition: grib2.h:481
int64_t g2int
Long integer type.
Definition: grib2.h:33
Header file with internal function prototypes NCEPLIBS-g2c library.
void pngpack(float *fld, g2int width, g2int height, g2int *idrstmpl, unsigned char *cpack, g2int *lcpack)
This subroutine packs up a float data field into PNG image format.
Definition: pngpack.c:242
void mkieee(float *a, g2int *rieee, g2int num)
This subroutine stores a list of real values in 32-bit IEEE floating point format.
Definition: mkieee.c:22
g2int * ext
Number of octets of each entry in the extension part of the template.
Definition: grib2_int.h:300
int g2c_check_msg(unsigned char *cgrib, g2int *lencurr, int verbose)
Check for 'GRIB' at the beginning of a GRIB message, and check to see if the message is already termi...
Definition: util.c:26
gtemplate * extpdstemplate(g2int number, g2int *list)
This subroutine generates the remaining octet map for a given Product Definition Template,...
Definition: pdstemplates.c:325
g2int extlen
Number of entries in the template extension.
Definition: grib2_int.h:296
g2int * map
Number of octets of each entry in the static part of the template.
Definition: grib2_int.h:290
g2int needext
Indicates whether or not the template needs to be extended.
Definition: grib2_int.h:293
void specpack(float *fld, g2int ndpts, g2int JJ, g2int KK, g2int MM, g2int *idrstmpl, unsigned char *cpack, g2int *lcpack)
This subroutine packs a spectral data field using the complex packing algorithm for spherical harmoni...
Definition: specpack.c:33
void simpack(float *fld, g2int ndpts, g2int *idrstmpl, unsigned char *cpack, g2int *lcpack)
This subroutine packs up a data field using the simple packing algorithm as defined in the GRIB2 docu...
Definition: simpack.c:34
void jpcpack(float *fld, g2int width, g2int height, g2int *idrstmpl, unsigned char *cpack, g2int *lcpack)
This function packs up a float array into a JPEG2000 code stream.
Definition: jpcpack.c:270
gtemplate * getpdstemplate(g2int number)
This subroutine returns PDS template information for a specified Product Definition Template.
Definition: pdstemplates.c:277
g2int maplen
Number of entries in the static part of the template.
Definition: grib2_int.h:286
Struct for GRIB template, returned by getgridtemplate().
Definition: grib2_int.h:276