1 : #include <stdio.h>
2 : #include <stdlib.h>
3 : #include "grib2.h"
4 :
5 :
6 2 : int comunpack(unsigned char *cpack,g2int lensec,g2int idrsnum,g2int *idrstmpl,g2int ndpts,g2float *fld)
7 : ////$$$ SUBPROGRAM DOCUMENTATION BLOCK
8 : // . . . .
9 : // SUBPROGRAM: comunpack
10 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-10-29
11 : //
12 : // ABSTRACT: This subroutine unpacks a data field that was packed using a
13 : // complex packing algorithm as defined in the GRIB2 documention,
14 : // using info from the GRIB2 Data Representation Template 5.2 or 5.3.
15 : // Supports GRIB2 complex packing templates with or without
16 : // spatial differences (i.e. DRTs 5.2 and 5.3).
17 : //
18 : // PROGRAM HISTORY LOG:
19 : // 2002-10-29 Gilbert
20 : // 2004-12-16 Gilbert - Added test ( provided by Arthur Taylor/MDL )
21 : // to verify that group widths and lengths are
22 : // consistent with section length.
23 : //
24 : // USAGE: int comunpack(unsigned char *cpack,g2int lensec,g2int idrsnum,
25 : // g2int *idrstmpl, g2int ndpts,g2float *fld)
26 : // INPUT ARGUMENT LIST:
27 : // cpack - pointer to the packed data field.
28 : // lensec - length of section 7 (used for error checking).
29 : // idrsnum - Data Representation Template number 5.N
30 : // Must equal 2 or 3.
31 : // idrstmpl - pointer to the array of values for Data Representation
32 : // Template 5.2 or 5.3
33 : // ndpts - The number of data values to unpack
34 : //
35 : // OUTPUT ARGUMENT LIST:
36 : // fld - Contains the unpacked data values. fld must be allocated
37 : // with at least ndpts*sizeof(g2float) bytes before
38 : // calling this routine.
39 : //
40 : // REMARKS: None
41 : //
42 : // ATTRIBUTES:
43 : // LANGUAGE: C
44 : // MACHINE:
45 : //
46 : //$$$//
47 : {
48 :
49 2 : g2int nbitsd=0,isign;
50 2 : g2int j,iofst,ival1,ival2,minsd,itemp,l,k,n,non=0;
51 2 : g2int *ifld,*ifldmiss=0;
52 : g2int *gref,*gwidth,*glen;
53 : g2int itype,ngroups,nbitsgref,nbitsgwidth,nbitsglen;
54 : g2int msng1,msng2;
55 : g2float ref,bscale,dscale,rmiss1,rmiss2;
56 : g2int totBit, totLen;
57 :
58 : //printf('IDRSTMPL: ',(idrstmpl(j),j=1,16)
59 2 : rdieee(idrstmpl+0,&ref,1);
60 : // printf("SAGTref: %f\n",ref);
61 2 : bscale = (g2float)int_power(2.0,idrstmpl[1]);
62 2 : dscale = (g2float)int_power(10.0,-idrstmpl[2]);
63 2 : nbitsgref = idrstmpl[3];
64 2 : itype = idrstmpl[4];
65 2 : ngroups = idrstmpl[9];
66 2 : nbitsgwidth = idrstmpl[11];
67 2 : nbitsglen = idrstmpl[15];
68 2 : if (idrsnum == 3)
69 2 : nbitsd=idrstmpl[17]*8;
70 :
71 : // Constant field
72 :
73 2 : if (ngroups == 0) {
74 0 : for (j=0;j<ndpts;j++) fld[j]=ref;
75 0 : return(0);
76 : }
77 :
78 2 : iofst=0;
79 2 : ifld=(g2int *)calloc(ndpts,sizeof(g2int));
80 : //printf("ALLOC ifld: %d %x\n",(int)ndpts,ifld);
81 2 : gref=(g2int *)calloc(ngroups,sizeof(g2int));
82 : //printf("ALLOC gref: %d %x\n",(int)ngroups,gref);
83 2 : gwidth=(g2int *)calloc(ngroups,sizeof(g2int));
84 : //printf("ALLOC gwidth: %d %x\n",(int)ngroups,gwidth);
85 : //
86 : // Get missing values, if supplied
87 : //
88 2 : if ( idrstmpl[6] == 1 ) {
89 2 : if (itype == 0)
90 2 : rdieee(idrstmpl+7,&rmiss1,1);
91 : else
92 0 : rmiss1=(g2float)idrstmpl[7];
93 : }
94 2 : if ( idrstmpl[6] == 2 ) {
95 0 : if (itype == 0) {
96 0 : rdieee(idrstmpl+7,&rmiss1,1);
97 0 : rdieee(idrstmpl+8,&rmiss2,1);
98 : }
99 : else {
100 0 : rmiss1=(g2float)idrstmpl[7];
101 0 : rmiss2=(g2float)idrstmpl[8];
102 : }
103 : }
104 :
105 : //printf("RMISSs: %f %f %f \n",rmiss1,rmiss2,ref);
106 : //
107 : // Extract Spatial differencing values, if using DRS Template 5.3
108 : //
109 2 : if (idrsnum == 3) {
110 2 : if (nbitsd != 0) {
111 2 : gbit(cpack,&isign,iofst,1);
112 2 : iofst=iofst+1;
113 2 : gbit(cpack,&ival1,iofst,nbitsd-1);
114 2 : iofst=iofst+nbitsd-1;
115 2 : if (isign == 1) ival1=-ival1;
116 2 : if (idrstmpl[16] == 2) {
117 2 : gbit(cpack,&isign,iofst,1);
118 2 : iofst=iofst+1;
119 2 : gbit(cpack,&ival2,iofst,nbitsd-1);
120 2 : iofst=iofst+nbitsd-1;
121 2 : if (isign == 1) ival2=-ival2;
122 : }
123 2 : gbit(cpack,&isign,iofst,1);
124 2 : iofst=iofst+1;
125 2 : gbit(cpack,&minsd,iofst,nbitsd-1);
126 2 : iofst=iofst+nbitsd-1;
127 2 : if (isign == 1) minsd=-minsd;
128 : }
129 : else {
130 0 : ival1=0;
131 0 : ival2=0;
132 0 : minsd=0;
133 : }
134 : //printf("SDu %ld %ld %ld %ld \n",ival1,ival2,minsd,nbitsd);
135 : }
136 : //
137 : // Extract Each Group's reference value
138 : //
139 : //printf("SAG1: %ld %ld %ld \n",nbitsgref,ngroups,iofst);
140 2 : if (nbitsgref != 0) {
141 2 : gbits(cpack,gref+0,iofst,nbitsgref,0,ngroups);
142 2 : itemp=nbitsgref*ngroups;
143 2 : iofst=iofst+itemp;
144 2 : if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
145 : }
146 : else {
147 0 : for (j=0;j<ngroups;j++)
148 0 : gref[j]=0;
149 : }
150 : //
151 : // Extract Each Group's bit width
152 : //
153 : //printf("SAG2: %ld %ld %ld %ld \n",nbitsgwidth,ngroups,iofst,idrstmpl[10]);
154 2 : if (nbitsgwidth != 0) {
155 2 : gbits(cpack,gwidth+0,iofst,nbitsgwidth,0,ngroups);
156 2 : itemp=nbitsgwidth*ngroups;
157 2 : iofst=iofst+itemp;
158 2 : if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
159 : }
160 : else {
161 0 : for (j=0;j<ngroups;j++)
162 0 : gwidth[j]=0;
163 : }
164 :
165 1017 : for (j=0;j<ngroups;j++)
166 1015 : gwidth[j]=gwidth[j]+idrstmpl[10];
167 :
168 : //
169 : // Extract Each Group's length (number of values in each group)
170 : //
171 2 : glen=(g2int *)calloc(ngroups,sizeof(g2int));
172 : //printf("ALLOC glen: %d %x\n",(int)ngroups,glen);
173 : //printf("SAG3: %ld %ld %ld %ld %ld \n",nbitsglen,ngroups,iofst,idrstmpl[13],idrstmpl[12]);
174 2 : if (nbitsglen != 0) {
175 2 : gbits(cpack,glen,iofst,nbitsglen,0,ngroups);
176 2 : itemp=nbitsglen*ngroups;
177 2 : iofst=iofst+itemp;
178 2 : if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
179 : }
180 : else {
181 0 : for (j=0;j<ngroups;j++)
182 0 : glen[j]=0;
183 : }
184 1017 : for (j=0;j<ngroups;j++)
185 1015 : glen[j]=(glen[j]*idrstmpl[13])+idrstmpl[12];
186 2 : glen[ngroups-1]=idrstmpl[14];
187 : //
188 : // Test to see if the group widths and lengths are consistent with number of
189 : // values, and length of section 7.
190 : //
191 2 : totBit = 0;
192 2 : totLen = 0;
193 1017 : for (j=0;j<ngroups;j++) {
194 1015 : totBit += (gwidth[j]*glen[j]);
195 1015 : totLen += glen[j];
196 : }
197 2 : if (totLen != ndpts) {
198 0 : return 1;
199 : }
200 2 : if (totBit / 8. > lensec) {
201 0 : return 1;
202 : }
203 : //
204 : // For each group, unpack data values
205 : //
206 2 : if ( idrstmpl[6] == 0 ) { // no missing values
207 0 : n=0;
208 0 : for (j=0;j<ngroups;j++) {
209 0 : if (gwidth[j] != 0) {
210 0 : gbits(cpack,ifld+n,iofst,gwidth[j],0,glen[j]);
211 0 : for (k=0;k<glen[j];k++) {
212 0 : ifld[n]=ifld[n]+gref[j];
213 0 : n=n+1;
214 : }
215 : }
216 : else {
217 0 : for (l=n;l<n+glen[j];l++) ifld[l]=gref[j];
218 0 : n=n+glen[j];
219 : }
220 0 : iofst=iofst+(gwidth[j]*glen[j]);
221 : }
222 : }
223 2 : else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
224 : // missing values included
225 2 : ifldmiss=(g2int *)malloc(ndpts*sizeof(g2int));
226 : //printf("ALLOC ifldmiss: %d %x\n",(int)ndpts,ifldmiss);
227 : //for (j=0;j<ndpts;j++) ifldmiss[j]=0;
228 2 : n=0;
229 2 : non=0;
230 1017 : for (j=0;j<ngroups;j++) {
231 : //printf(" SAGNGP %d %d %d %d\n",j,gwidth[j],glen[j],gref[j]);
232 1015 : if (gwidth[j] != 0) {
233 899 : msng1=(g2int)int_power(2.0,gwidth[j])-1;
234 899 : msng2=msng1-1;
235 899 : gbits(cpack,ifld+n,iofst,gwidth[j],0,glen[j]);
236 899 : iofst=iofst+(gwidth[j]*glen[j]);
237 39511 : for (k=0;k<glen[j];k++) {
238 38612 : if (ifld[n] == msng1) {
239 458 : ifldmiss[n]=1;
240 : //ifld[n]=0;
241 : }
242 38154 : else if (idrstmpl[6]==2 && ifld[n]==msng2) {
243 0 : ifldmiss[n]=2;
244 : //ifld[n]=0;
245 : }
246 : else {
247 38154 : ifldmiss[n]=0;
248 38154 : ifld[non++]=ifld[n]+gref[j];
249 : }
250 38612 : n++;
251 : }
252 : }
253 : else {
254 116 : msng1=(g2int)int_power(2.0,nbitsgref)-1;
255 116 : msng2=msng1-1;
256 116 : if (gref[j] == msng1) {
257 116 : for (l=n;l<n+glen[j];l++) ifldmiss[l]=1;
258 : }
259 0 : else if (idrstmpl[6]==2 && gref[j]==msng2) {
260 0 : for (l=n;l<n+glen[j];l++) ifldmiss[l]=2;
261 : }
262 : else {
263 0 : for (l=n;l<n+glen[j];l++) ifldmiss[l]=0;
264 0 : for (l=non;l<non+glen[j];l++) ifld[l]=gref[j];
265 0 : non += glen[j];
266 : }
267 116 : n=n+glen[j];
268 : }
269 : }
270 : }
271 :
272 2 : if ( gref != 0 ) free(gref);
273 2 : if ( gwidth != 0 ) free(gwidth);
274 2 : if ( glen != 0 ) free(glen);
275 : //
276 : // If using spatial differences, add overall min value, and
277 : // sum up recursively
278 : //
279 : //printf("SAGod: %ld %ld\n",idrsnum,idrstmpl[16]);
280 2 : if (idrsnum == 3) { // spatial differencing
281 2 : if (idrstmpl[16] == 1) { // first order
282 0 : ifld[0]=ival1;
283 0 : if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
284 0 : else itemp=non;
285 0 : for (n=1;n<itemp;n++) {
286 0 : ifld[n]=ifld[n]+minsd;
287 0 : ifld[n]=ifld[n]+ifld[n-1];
288 : }
289 : }
290 2 : else if (idrstmpl[16] == 2) { // second order
291 2 : ifld[0]=ival1;
292 2 : ifld[1]=ival2;
293 2 : if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
294 2 : else itemp=non;
295 38152 : for (n=2;n<itemp;n++) {
296 38150 : ifld[n]=ifld[n]+minsd;
297 38150 : ifld[n]=ifld[n]+(2*ifld[n-1])-ifld[n-2];
298 : }
299 : }
300 : }
301 : //
302 : // Scale data back to original form
303 : //
304 : //printf("SAGT: %f %f %f\n",ref,bscale,dscale);
305 2 : if ( idrstmpl[6] == 0 ) { // no missing values
306 0 : for (n=0;n<ndpts;n++) {
307 0 : fld[n]=(((g2float)ifld[n]*bscale)+ref)*dscale;
308 : }
309 : }
310 2 : else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
311 : // missing values included
312 2 : non=0;
313 45668 : for (n=0;n<ndpts;n++) {
314 45666 : if ( ifldmiss[n] == 0 ) {
315 38154 : fld[n]=(((g2float)ifld[non++]*bscale)+ref)*dscale;
316 : //printf(" SAG %d %f %d %f %f %f\n",n,fld[n],ifld[non-1],bscale,ref,dscale);
317 : }
318 7512 : else if ( ifldmiss[n] == 1 )
319 7512 : fld[n]=rmiss1;
320 0 : else if ( ifldmiss[n] == 2 )
321 0 : fld[n]=rmiss2;
322 : }
323 2 : if ( ifldmiss != 0 ) free(ifldmiss);
324 : }
325 :
326 2 : if ( ifld != 0 ) free(ifld);
327 :
328 2 : return(0);
329 :
330 : }
|
