1 : #include <stdlib.h>
2 : #include <math.h>
3 : #include "grib2.h"
4 :
5 :
6 0 : void simpack(g2float *fld,g2int ndpts,g2int *idrstmpl,unsigned char *cpack,g2int *lcpack)
7 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
8 : // . . . .
9 : // SUBPROGRAM: simpack
10 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-11-06
11 : //
12 : // ABSTRACT: This subroutine packs up a data field using the simple
13 : // packing algorithm as defined in the GRIB2 documention. It
14 : // also fills in GRIB2 Data Representation Template 5.0 with the
15 : // appropriate values.
16 : //
17 : // PROGRAM HISTORY LOG:
18 : // 2002-11-06 Gilbert
19 : //
20 : // USAGE: CALL simpack(fld,ndpts,idrstmpl,cpack,lcpack)
21 : // INPUT ARGUMENT LIST:
22 : // fld[] - Contains the data values to pack
23 : // ndpts - The number of data values in array fld[]
24 : // idrstmpl - Contains the array of values for Data Representation
25 : // Template 5.0
26 : // [0] = Reference value - ignored on input
27 : // [1] = Binary Scale Factor
28 : // [2] = Decimal Scale Factor
29 : // [3] = Number of bits used to pack data, if value is
30 : // > 0 and <= 31.
31 : // If this input value is 0 or outside above range
32 : // then the num of bits is calculated based on given
33 : // data and scale factors.
34 : // [4] = Original field type - currently ignored on input
35 : // Data values assumed to be reals.
36 : //
37 : // OUTPUT ARGUMENT LIST:
38 : // idrstmpl - Contains the array of values for Data Representation
39 : // Template 5.0
40 : // [0] = Reference value - set by simpack routine.
41 : // [1] = Binary Scale Factor - unchanged from input
42 : // [2] = Decimal Scale Factor - unchanged from input
43 : // [3] = Number of bits used to pack data, unchanged from
44 : // input if value is between 0 and 31.
45 : // If this input value is 0 or outside above range
46 : // then the num of bits is calculated based on given
47 : // data and scale factors.
48 : // [4] = Original field type - currently set = 0 on output.
49 : // Data values assumed to be reals.
50 : // cpack - The packed data field
51 : // lcpack - length of packed field starting at cpack.
52 : //
53 : // REMARKS: None
54 : //
55 : // ATTRIBUTES:
56 : // LANGUAGE: C
57 : // MACHINE:
58 : //
59 : //$$$
60 : {
61 :
62 : static g2int zero=0;
63 : g2int *ifld;
64 : g2int j,nbits,imin,imax,maxdif,nbittot,left;
65 : g2float bscale,dscale,rmax,rmin,temp;
66 : double maxnum;
67 : static g2float alog2=0.69314718; // ln(2.0)
68 :
69 0 : bscale=int_power(2.0,-idrstmpl[1]);
70 0 : dscale=int_power(10.0,idrstmpl[2]);
71 0 : if (idrstmpl[3] <= 0 || idrstmpl[3] > 31)
72 0 : nbits=0;
73 : else
74 0 : nbits=idrstmpl[3];
75 : //
76 : // Find max and min values in the data
77 : //
78 0 : rmax=fld[0];
79 0 : rmin=fld[0];
80 0 : for (j=1;j<ndpts;j++) {
81 0 : if (fld[j] > rmax) rmax=fld[j];
82 0 : if (fld[j] < rmin) rmin=fld[j];
83 : }
84 :
85 0 : ifld=calloc(ndpts,sizeof(g2int));
86 : //
87 : // If max and min values are not equal, pack up field.
88 : // If they are equal, we have a constant field, and the reference
89 : // value (rmin) is the value for each point in the field and
90 : // set nbits to 0.
91 : //
92 0 : if (rmin != rmax) {
93 : //
94 : // Determine which algorithm to use based on user-supplied
95 : // binary scale factor and number of bits.
96 : //
97 0 : if (nbits==0 && idrstmpl[1]==0) {
98 : //
99 : // No binary scaling and calculate minumum number of
100 : // bits in which the data will fit.
101 : //
102 0 : imin=(g2int)RINT(rmin*dscale);
103 0 : imax=(g2int)RINT(rmax*dscale);
104 0 : maxdif=imax-imin;
105 0 : temp=log((double)(maxdif+1))/alog2;
106 0 : nbits=(g2int)ceil(temp);
107 0 : rmin=(g2float)imin;
108 : // scale data
109 0 : for(j=0;j<ndpts;j++)
110 0 : ifld[j]=(g2int)RINT(fld[j]*dscale)-imin;
111 : }
112 0 : else if (nbits!=0 && idrstmpl[1]==0) {
113 : //
114 : // Use minimum number of bits specified by user and
115 : // adjust binary scaling factor to accomodate data.
116 : //
117 0 : rmin=rmin*dscale;
118 0 : rmax=rmax*dscale;
119 0 : maxnum=int_power(2.0,nbits)-1;
120 0 : temp=log(maxnum/(rmax-rmin))/alog2;
121 0 : idrstmpl[1]=(g2int)ceil(-1.0*temp);
122 0 : bscale=int_power(2.0,-idrstmpl[1]);
123 : // scale data
124 0 : for (j=0;j<ndpts;j++)
125 0 : ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
126 : }
127 0 : else if (nbits==0 && idrstmpl[1]!=0) {
128 : //
129 : // Use binary scaling factor and calculate minumum number of
130 : // bits in which the data will fit.
131 : //
132 0 : rmin=rmin*dscale;
133 0 : rmax=rmax*dscale;
134 0 : maxdif=(g2int)RINT((rmax-rmin)*bscale);
135 0 : temp=log((double)(maxdif+1))/alog2;
136 0 : nbits=(g2int)ceil(temp);
137 : // scale data
138 0 : for (j=0;j<ndpts;j++)
139 0 : ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
140 : }
141 0 : else if (nbits!=0 && idrstmpl[1]!=0) {
142 : //
143 : // Use binary scaling factor and use minumum number of
144 : // bits specified by user. Dangerous - may loose
145 : // information if binary scale factor and nbits not set
146 : // properly by user.
147 : //
148 0 : rmin=rmin*dscale;
149 : // scale data
150 0 : for (j=0;j<ndpts;j++)
151 0 : ifld[j]=(g2int)RINT(((fld[j]*dscale)-rmin)*bscale);
152 : }
153 : //
154 : // Pack data, Pad last octet with Zeros, if necessary,
155 : // and calculate the length of the packed data in bytes
156 : //
157 0 : sbits(cpack,ifld+0,0,nbits,0,ndpts);
158 0 : nbittot=nbits*ndpts;
159 0 : left=8-(nbittot%8);
160 0 : if (left != 8) {
161 0 : sbit(cpack,&zero,nbittot,left); // Pad with zeros to fill Octet
162 0 : nbittot=nbittot+left;
163 : }
164 0 : *lcpack=nbittot/8;
165 : }
166 : else {
167 0 : nbits=0;
168 0 : *lcpack=0;
169 : }
170 :
171 : //
172 : // Fill in ref value and number of bits in Template 5.0
173 : //
174 : //printf("SAGmkieee %f\n",rmin);
175 0 : mkieee(&rmin,idrstmpl+0,1); // ensure reference value is IEEE format
176 : //printf("SAGmkieee %ld\n",idrstmpl[0]);
177 0 : idrstmpl[3]=nbits;
178 0 : idrstmpl[4]=0; // original data were reals
179 :
180 0 : free(ifld);
181 0 : }
|
