1 : #include <stdio.h>
2 : #include <stdlib.h>
3 : #include "grib2.h"
4 :
5 : g2int getdim(unsigned char *,g2int *,g2int *,g2int *);
6 : g2int getpoly(unsigned char *,g2int *,g2int *,g2int *);
7 : void simpack(g2float *, g2int, g2int *, unsigned char *, g2int *);
8 : void cmplxpack(g2float *, g2int, g2int, g2int *, unsigned char *, g2int *);
9 : void specpack(g2float *,g2int,g2int,g2int,g2int,g2int *,unsigned char *,
10 : g2int *);
11 : void jpcpack(g2float *,g2int,g2int,g2int *,unsigned char *,g2int *);
12 : #ifdef USE_PNG
13 : void pngpack(g2float *,g2int,g2int,g2int *,unsigned char *,g2int *);
14 : #endif /* USE_PNG */
15 :
16 :
17 :
18 0 : g2int g2_addfield(unsigned char *cgrib,g2int ipdsnum,g2int *ipdstmpl,
19 : g2float *coordlist,g2int numcoord,g2int idrsnum,g2int *idrstmpl,
20 : g2float *fld,g2int ngrdpts,g2int ibmap,g2int *bmap)
21 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
22 : // . . . .
23 : // SUBPROGRAM: g2_addfield
24 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-11-05
25 : //
26 : // ABSTRACT: This routine packs up Sections 4 through 7 for a given field
27 : // and adds them to a GRIB2 message. They are Product Definition Section,
28 : // Data Representation Section, Bit-Map Section and Data Section,
29 : // respectively.
30 : // This routine is used with routines "g2_create", "g2_addlocal",
31 : // "g2_addgrid", and "g2_gribend" to create a complete GRIB2 message.
32 : // g2_create must be called first to initialize a new GRIB2 message.
33 : // Also, routine g2_addgrid must be called after g2_create and
34 : // before this routine to add the appropriate grid description to
35 : // the GRIB2 message. Also, a call to g2_gribend is required to complete
36 : // GRIB2 message after all fields have been added.
37 : //
38 : // PROGRAM HISTORY LOG:
39 : // 2002-11-05 Gilbert
40 : // 2002-12-23 Gilbert - Added complex spherical harmonic packing
41 : // 2003-08-27 Gilbert - Added support for new templates using
42 : // PNG and JPEG2000 algorithms/templates.
43 : // 2004-11-29 Gilbert - JPEG2000 now allowed to use WMO Template no. 5.40
44 : // PNG now allowed to use WMO Template no. 5.41
45 : // - Added check to determine if packing algorithm failed.
46 : // 2005-05-10 Gilbert - Imposed minimum size on cpack, used to hold encoded
47 : // bit string.
48 : //
49 : // USAGE: int g2_addfield(unsigned char *cgrib,g2int ipdsnum,g2int *ipdstmpl,
50 : // g2float *coordlist,g2int numcoord,g2int idrsnum,g2int *idrstmpl,
51 : // g2float *fld,g2int ngrdpts,g2int ibmap,g2int *bmap)
52 : // INPUT ARGUMENT LIST:
53 : // cgrib - Char array that contains the GRIB2 message to which sections
54 : // 4 through 7 should be added.
55 : // ipdsnum - Product Definition Template Number ( see Code Table 4.0)
56 : // ipdstmpl - Contains the data values for the specified Product Definition
57 : // Template ( N=ipdsnum ). Each element of this integer
58 : // array contains an entry (in the order specified) of Product
59 : // Defintion Template 4.N
60 : // coordlist- Array containg floating point values intended to document
61 : // the vertical discretisation associated to model data
62 : // on hybrid coordinate vertical levels.
63 : // numcoord - number of values in array coordlist.
64 : // idrsnum - Data Representation Template Number ( see Code Table 5.0 )
65 : // idrstmpl - Contains the data values for the specified Data Representation
66 : // Template ( N=idrsnum ). Each element of this integer
67 : // array contains an entry (in the order specified) of Data
68 : // Representation Template 5.N
69 : // Note that some values in this template (eg. reference
70 : // values, number of bits, etc...) may be changed by the
71 : // data packing algorithms.
72 : // Use this to specify scaling factors and order of
73 : // spatial differencing, if desired.
74 : // fld[] - Array of data points to pack.
75 : // ngrdpts - Number of data points in grid.
76 : // i.e. size of fld and bmap.
77 : // ibmap - Bitmap indicator ( see Code Table 6.0 )
78 : // 0 = bitmap applies and is included in Section 6.
79 : // 1-253 = Predefined bitmap applies
80 : // 254 = Previously defined bitmap applies to this field
81 : // 255 = Bit map does not apply to this product.
82 : // bmap[] - Integer array containing bitmap to be added. ( if ibmap=0 )
83 : //
84 : // OUTPUT ARGUMENT LIST:
85 : // cgrib - Character array to contain the updated GRIB2 message.
86 : // Must be allocated large enough to store the entire
87 : // GRIB2 message.
88 : //
89 : // RETURN VALUES:
90 : // ierr - Return code.
91 : // > 0 = Current size of updated GRIB2 message
92 : // -1 = GRIB message was not initialized. Need to call
93 : // routine g2_create first.
94 : // -2 = GRIB message already complete. Cannot add new section.
95 : // -3 = Sum of Section byte counts doesn't add to total byte count
96 : // -4 = Previous Section was not 3 or 7.
97 : // -5 = Could not find requested Product Definition Template.
98 : // -6 = Section 3 (GDS) not previously defined in message
99 : // -7 = Tried to use unsupported Data Representationi Template
100 : // -8 = Specified use of a previously defined bitmap, but one
101 : // does not exist in the GRIB message.
102 : // -9 = GDT of one of 5.50 through 5.53 required to pack field
103 : // using DRT 5.51.
104 : // -10 = Error packing data field.
105 : //
106 : // REMARKS: Note that the Sections 4 through 7 can only follow
107 : // Section 3 or Section 7 in a GRIB2 message.
108 : //
109 : // ATTRIBUTES:
110 : // LANGUAGE: C
111 : // MACHINE:
112 : //
113 : //$$$
114 : {
115 : g2int ierr;
116 : static unsigned char G=0x47; // 'G'
117 : static unsigned char R=0x52; // 'R'
118 : static unsigned char I=0x49; // 'I'
119 : static unsigned char B=0x42; // 'B'
120 : static unsigned char s7=0x37; // '7'
121 :
122 : unsigned char *cpack;
123 : static g2int zero=0,one=1,four=4,five=5,six=6,seven=7;
124 0 : const g2int minsize=50000;
125 : g2int iofst,ibeg,lencurr,len,nsize;
126 : g2int ilen,isecnum,i,nbits,temp,left;
127 : g2int ibmprev,j,lcpack,ioctet,newlen,ndpts;
128 : g2int lensec4,lensec5,lensec6,lensec7;
129 0 : g2int issec3,isprevbmap,lpos3=0,JJ,KK,MM;
130 : g2int *coordieee;
131 : g2int width,height,iscan,itemp;
132 : g2float *pfld;
133 : xxtemplate *mappds,*mapdrs;
134 0 : unsigned int allones=4294967295u;
135 :
136 0 : ierr=0;
137 : //
138 : // Check to see if beginning of GRIB message exists
139 : //
140 0 : if ( cgrib[0]!=G || cgrib[1]!=R || cgrib[2]!=I || cgrib[3]!=B ) {
141 0 : printf("g2_addfield: GRIB not found in given message.\n");
142 0 : printf("g2_addfield: Call to routine g2_create required to initialize GRIB messge.\n");
143 0 : ierr=-1;
144 0 : return(ierr);
145 : }
146 : //
147 : // Get current length of GRIB message
148 : //
149 0 : gbit(cgrib,&lencurr,96,32);
150 : //
151 : // Check to see if GRIB message is already complete
152 : //
153 0 : if ( cgrib[lencurr-4]==s7 && cgrib[lencurr-3]==s7 &&
154 0 : cgrib[lencurr-2]==s7 && cgrib[lencurr-1]==s7 ) {
155 0 : printf("g2_addfield: GRIB message already complete. Cannot add new section.\n");
156 0 : ierr=-2;
157 0 : return(ierr);
158 : }
159 : //
160 : // Loop through all current sections of the GRIB message to
161 : // find the last section number.
162 : //
163 0 : issec3=0;
164 0 : isprevbmap=0;
165 0 : len=16; // length of Section 0
166 : for (;;) {
167 : // Get number and length of next section
168 0 : iofst=len*8;
169 0 : gbit(cgrib,&ilen,iofst,32);
170 0 : iofst=iofst+32;
171 0 : gbit(cgrib,&isecnum,iofst,8);
172 0 : iofst=iofst+8;
173 : // Check if previous Section 3 exists
174 0 : if (isecnum == 3) {
175 0 : issec3=1;
176 0 : lpos3=len;
177 : }
178 : // Check if a previous defined bitmap exists
179 0 : if (isecnum == 6) {
180 0 : gbit(cgrib,&ibmprev,iofst,8);
181 0 : iofst=iofst+8;
182 0 : if ((ibmprev >= 0) && (ibmprev <= 253)) isprevbmap=1;
183 : }
184 0 : len=len+ilen;
185 : // Exit loop if last section reached
186 0 : if ( len == lencurr ) break;
187 : // If byte count for each section doesn't match current
188 : // total length, then there is a problem.
189 0 : if ( len > lencurr ) {
190 0 : printf("g2_addfield: Section byte counts don''t add to total.\n");
191 0 : printf("g2_addfield: Sum of section byte counts = %d\n",len);
192 0 : printf("g2_addfield: Total byte count in Section 0 = %d\n",lencurr);
193 0 : ierr=-3;
194 0 : return(ierr);
195 : }
196 0 : }
197 : //
198 : // Sections 4 through 7 can only be added after section 3 or 7.
199 : //
200 0 : if ( (isecnum != 3) && (isecnum != 7) ) {
201 0 : printf("g2_addfield: Sections 4-7 can only be added after Section 3 or 7.\n");
202 0 : printf("g2_addfield: Section ',isecnum,' was the last found in given GRIB message.\n");
203 0 : ierr=-4;
204 0 : return(ierr);
205 : //
206 : // Sections 4 through 7 can only be added if section 3 was previously defined.
207 : //
208 : }
209 0 : else if ( ! issec3) {
210 0 : printf("g2_addfield: Sections 4-7 can only be added if Section 3 was previously included.\n");
211 0 : printf("g2_addfield: Section 3 was not found in given GRIB message.\n");
212 0 : printf("g2_addfield: Call to routine addgrid required to specify Grid definition.\n");
213 0 : ierr=-6;
214 0 : return(ierr);
215 : }
216 : //
217 : // Add Section 4 - Product Definition Section
218 : //
219 0 : ibeg=lencurr*8; // Calculate offset for beginning of section 4
220 0 : iofst=ibeg+32; // leave space for length of section
221 0 : sbit(cgrib,&four,iofst,8); // Store section number ( 4 )
222 0 : iofst=iofst+8;
223 0 : sbit(cgrib,&numcoord,iofst,16); // Store num of coordinate values
224 0 : iofst=iofst+16;
225 0 : sbit(cgrib,&ipdsnum,iofst,16); // Store Prod Def Template num.
226 0 : iofst=iofst+16;
227 : //
228 : // Get Product Definition Template
229 : //
230 0 : mappds=getpdstemplate(ipdsnum);
231 0 : if (mappds == 0) { // undefined template
232 0 : ierr=-5;
233 0 : return(ierr);
234 : }
235 : //
236 : // Extend the Product Definition Template, if necessary.
237 : // The number of values in a specific template may vary
238 : // depending on data specified in the "static" part of the
239 : // template.
240 : //
241 0 : if ( mappds->needext ) {
242 0 : free(mappds);
243 0 : mappds=extpdstemplate(ipdsnum,ipdstmpl);
244 : }
245 : //
246 : // Pack up each input value in array ipdstmpl into the
247 : // the appropriate number of octets, which are specified in
248 : // corresponding entries in array mappds.
249 : //
250 0 : for (i=0;i<mappds->maplen;i++) {
251 0 : nbits=abs(mappds->map[i])*8;
252 0 : if ( (mappds->map[i] >= 0) || (ipdstmpl[i] >= 0) )
253 0 : sbit(cgrib,ipdstmpl+i,iofst,nbits);
254 : else {
255 0 : sbit(cgrib,&one,iofst,1);
256 0 : temp=abs(ipdstmpl[i]);
257 0 : sbit(cgrib,&temp,iofst+1,nbits-1);
258 : }
259 0 : iofst=iofst+nbits;
260 : }
261 : // Pack template extension, if appropriate
262 0 : j=mappds->maplen;
263 0 : if ( mappds->needext && (mappds->extlen > 0) ) {
264 0 : for (i=0;i<mappds->extlen;i++) {
265 0 : nbits=abs(mappds->ext[i])*8;
266 0 : if ( (mappds->ext[i] >= 0) || (ipdstmpl[j] >= 0) )
267 0 : sbit(cgrib,ipdstmpl+j,iofst,nbits);
268 : else {
269 0 : sbit(cgrib,&one,iofst,1);
270 0 : temp=abs(ipdstmpl[j]);
271 0 : sbit(cgrib,&temp,iofst+1,nbits-1);
272 : }
273 0 : iofst=iofst+nbits;
274 0 : j++;
275 : }
276 : }
277 0 : free(mappds);
278 : //
279 : // Add Optional list of vertical coordinate values
280 : // after the Product Definition Template, if necessary.
281 : //
282 0 : if ( numcoord != 0 ) {
283 0 : coordieee=(g2int *)calloc(numcoord,sizeof(g2int));
284 0 : mkieee(coordlist,coordieee,numcoord);
285 0 : sbits(cgrib,coordieee,iofst,32,0,numcoord);
286 0 : iofst=iofst+(32*numcoord);
287 0 : free(coordieee);
288 : }
289 : //
290 : // Calculate length of section 4 and store it in octets
291 : // 1-4 of section 4.
292 : //
293 0 : lensec4=(iofst-ibeg)/8;
294 0 : sbit(cgrib,&lensec4,ibeg,32);
295 : //
296 : // Pack Data using appropriate algorithm
297 : //
298 : //
299 : // Get Data Representation Template
300 : //
301 0 : mapdrs=getdrstemplate(idrsnum);
302 0 : if (mapdrs == 0) {
303 0 : ierr=-5;
304 0 : return(ierr);
305 : }
306 : //
307 : // contract data field, removing data at invalid grid points,
308 : // if bit-map is provided with field.
309 : //
310 0 : if ( ibmap == 0 || ibmap==254 ) {
311 0 : pfld=(g2float *)malloc(ngrdpts*sizeof(g2float));
312 0 : ndpts=0;
313 0 : for (j=0;j<ngrdpts;j++) {
314 0 : if ( bmap[j]==1 ) pfld[ndpts++]=fld[j];
315 : }
316 : }
317 : else {
318 0 : ndpts=ngrdpts;
319 0 : pfld=fld;
320 : }
321 0 : nsize=ndpts*4;
322 0 : if ( nsize < minsize ) nsize=minsize;
323 0 : cpack=malloc(nsize);
324 0 : if (idrsnum == 0) // Simple Packing
325 0 : simpack(pfld,ndpts,idrstmpl,cpack,&lcpack);
326 0 : else if (idrsnum==2 || idrsnum==3) // Complex Packing
327 0 : cmplxpack(pfld,ndpts,idrsnum,idrstmpl,cpack,&lcpack);
328 0 : else if (idrsnum == 50) { // Sperical Harmonic Simple Packing
329 0 : simpack(pfld+1,ndpts-1,idrstmpl,cpack,&lcpack);
330 0 : mkieee(pfld+0,idrstmpl+4,1); // ensure RE(0,0) value is IEEE format
331 : }
332 0 : else if (idrsnum == 51) { // Sperical Harmonic Complex Packing
333 0 : getpoly(cgrib+lpos3,&JJ,&KK,&MM);
334 0 : if ( JJ!=0 && KK!=0 && MM!=0 )
335 0 : specpack(pfld,ndpts,JJ,KK,MM,idrstmpl,cpack,&lcpack);
336 : else {
337 0 : printf("g2_addfield: Cannot pack DRT 5.51.\n");
338 0 : return (-9);
339 : }
340 : }
341 0 : else if (idrsnum == 40 || idrsnum == 40000) { /* JPEG2000 encoding */
342 0 : if (ibmap == 255) {
343 0 : getdim(cgrib+lpos3,&width,&height,&iscan);
344 0 : if ( width==0 || height==0 ) {
345 0 : width=ndpts;
346 0 : height=1;
347 : }
348 0 : else if ( width==allones || height==allones ) {
349 0 : width=ndpts;
350 0 : height=1;
351 : }
352 0 : else if ( (iscan&32) == 32) { /* Scanning mode: bit 3 */
353 0 : itemp=width;
354 0 : width=height;
355 0 : height=itemp;
356 : }
357 : }
358 : else {
359 0 : width=ndpts;
360 0 : height=1;
361 : }
362 0 : lcpack=nsize;
363 0 : jpcpack(pfld,width,height,idrstmpl,cpack,&lcpack);
364 : }
365 : #ifdef USE_PNG
366 : else if (idrsnum == 41 || idrsnum == 40010) { /* PNG encoding */
367 : if (ibmap == 255) {
368 : getdim(cgrib+lpos3,&width,&height,&iscan);
369 : if ( width==0 || height==0 ) {
370 : width=ndpts;
371 : height=1;
372 : }
373 : else if ( width==allones || height==allones ) {
374 : width=ndpts;
375 : height=1;
376 : }
377 : else if ( (iscan&32) == 32) { /* Scanning mode: bit 3 */
378 : itemp=width;
379 : width=height;
380 : height=itemp;
381 : }
382 : }
383 : else {
384 : width=ndpts;
385 : height=1;
386 : }
387 : pngpack(pfld,width,height,idrstmpl,cpack,&lcpack);
388 : }
389 : #endif /* USE_PNG */
390 : else {
391 0 : printf("g2_addfield: Data Representation Template 5.%d not yet implemented.\n",idrsnum);
392 0 : ierr=-7;
393 0 : return(ierr);
394 : }
395 0 : if ( ibmap == 0 || ibmap==254 ) { // free temp space
396 0 : if (fld != pfld) free(pfld);
397 : }
398 0 : if ( lcpack < 0 ) {
399 0 : if( cpack != 0 ) free(cpack);
400 0 : ierr=-10;
401 0 : return(ierr);
402 : }
403 :
404 : //
405 : // Add Section 5 - Data Representation Section
406 : //
407 0 : ibeg=iofst; // Calculate offset for beginning of section 5
408 0 : iofst=ibeg+32; // leave space for length of section
409 0 : sbit(cgrib,&five,iofst,8); // Store section number ( 5 )
410 0 : iofst=iofst+8;
411 0 : sbit(cgrib,&ndpts,iofst,32); // Store num of actual data points
412 0 : iofst=iofst+32;
413 0 : sbit(cgrib,&idrsnum,iofst,16); // Store Data Repr. Template num.
414 0 : iofst=iofst+16;
415 : //
416 : // Pack up each input value in array idrstmpl into the
417 : // the appropriate number of octets, which are specified in
418 : // corresponding entries in array mapdrs.
419 : //
420 0 : for (i=0;i<mapdrs->maplen;i++) {
421 0 : nbits=abs(mapdrs->map[i])*8;
422 0 : if ( (mapdrs->map[i] >= 0) || (idrstmpl[i] >= 0) )
423 0 : sbit(cgrib,idrstmpl+i,iofst,nbits);
424 : else {
425 0 : sbit(cgrib,&one,iofst,1);
426 0 : temp=abs(idrstmpl[i]);
427 0 : sbit(cgrib,&temp,iofst+1,nbits-1);
428 : }
429 0 : iofst=iofst+nbits;
430 : }
431 0 : free(mapdrs);
432 : //
433 : // Calculate length of section 5 and store it in octets
434 : // 1-4 of section 5.
435 : //
436 0 : lensec5=(iofst-ibeg)/8;
437 0 : sbit(cgrib,&lensec5,ibeg,32);
438 :
439 : //
440 : // Add Section 6 - Bit-Map Section
441 : //
442 0 : ibeg=iofst; // Calculate offset for beginning of section 6
443 0 : iofst=ibeg+32; // leave space for length of section
444 0 : sbit(cgrib,&six,iofst,8); // Store section number ( 6 )
445 0 : iofst=iofst+8;
446 0 : sbit(cgrib,&ibmap,iofst,8); // Store Bit Map indicator
447 0 : iofst=iofst+8;
448 : //
449 : // Store bitmap, if supplied
450 : //
451 0 : if (ibmap == 0) {
452 0 : sbits(cgrib,bmap,iofst,1,0,ngrdpts); // Store BitMap
453 0 : iofst=iofst+ngrdpts;
454 : }
455 : //
456 : // If specifying a previously defined bit-map, make sure
457 : // one already exists in the current GRIB message.
458 : //
459 0 : if ((ibmap==254) && ( ! isprevbmap)) {
460 0 : printf("g2_addfield: Requested previously defined bitmap,");
461 0 : printf(" but one does not exist in the current GRIB message.\n");
462 0 : ierr=-8;
463 0 : return(ierr);
464 : }
465 : //
466 : // Calculate length of section 6 and store it in octets
467 : // 1-4 of section 6. Pad to end of octect, if necessary.
468 : //
469 0 : left=8-(iofst%8);
470 0 : if (left != 8) {
471 0 : sbit(cgrib,&zero,iofst,left); // Pad with zeros to fill Octet
472 0 : iofst=iofst+left;
473 : }
474 0 : lensec6=(iofst-ibeg)/8;
475 0 : sbit(cgrib,&lensec6,ibeg,32);
476 :
477 : //
478 : // Add Section 7 - Data Section
479 : //
480 0 : ibeg=iofst; // Calculate offset for beginning of section 7
481 0 : iofst=ibeg+32; // leave space for length of section
482 0 : sbit(cgrib,&seven,iofst,8); // Store section number ( 7 )
483 0 : iofst=iofst+8;
484 : // Store Packed Binary Data values, if non-constant field
485 0 : if (lcpack != 0) {
486 0 : ioctet=iofst/8;
487 : //cgrib(ioctet+1:ioctet+lcpack)=cpack(1:lcpack)
488 0 : for (j=0;j<lcpack;j++) cgrib[ioctet+j]=cpack[j];
489 0 : iofst=iofst+(8*lcpack);
490 : }
491 : //
492 : // Calculate length of section 7 and store it in octets
493 : // 1-4 of section 7.
494 : //
495 0 : lensec7=(iofst-ibeg)/8;
496 0 : sbit(cgrib,&lensec7,ibeg,32);
497 :
498 0 : if( cpack != 0 ) free(cpack);
499 : //
500 : // Update current byte total of message in Section 0
501 : //
502 0 : newlen=lencurr+lensec4+lensec5+lensec6+lensec7;
503 0 : sbit(cgrib,&newlen,96,32);
504 :
505 0 : return(newlen);
506 :
507 : }
|
