1 #include <stdio.h>
2 #include <stdlib.h>
3 #include "grib2.h"
4
5
comunpack(unsigned char * cpack,g2int lensec,g2int idrsnum,g2int * idrstmpl,g2int ndpts,g2float * fld)6 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 g2int nbitsd=0,isign;
50 g2int j,iofst,ival1,ival2,minsd,itemp,l,k,n,non=0;
51 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 rdieee(idrstmpl+0,&ref,1);
60 // printf("SAGTref: %f\n",ref);
61 bscale = (g2float)int_power(2.0,idrstmpl[1]);
62 dscale = (g2float)int_power(10.0,-idrstmpl[2]);
63 nbitsgref = idrstmpl[3];
64 itype = idrstmpl[4];
65 ngroups = idrstmpl[9];
66 nbitsgwidth = idrstmpl[11];
67 nbitsglen = idrstmpl[15];
68 if (idrsnum == 3)
69 nbitsd=idrstmpl[17]*8;
70
71 // Constant field
72
73 if (ngroups == 0) {
74 for (j=0;j<ndpts;j++) fld[j]=ref;
75 return(0);
76 }
77
78 iofst=0;
79 ifld=(g2int *)calloc(ndpts,sizeof(g2int));
80 //printf("ALLOC ifld: %d %x\n",(int)ndpts,ifld);
81 gref=(g2int *)calloc(ngroups,sizeof(g2int));
82 //printf("ALLOC gref: %d %x\n",(int)ngroups,gref);
83 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 if ( idrstmpl[6] == 1 ) {
89 if (itype == 0)
90 rdieee(idrstmpl+7,&rmiss1,1);
91 else
92 rmiss1=(g2float)idrstmpl[7];
93 }
94 if ( idrstmpl[6] == 2 ) {
95 if (itype == 0) {
96 rdieee(idrstmpl+7,&rmiss1,1);
97 rdieee(idrstmpl+8,&rmiss2,1);
98 }
99 else {
100 rmiss1=(g2float)idrstmpl[7];
101 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 if (idrsnum == 3) {
110 if (nbitsd != 0) {
111 gbit(cpack,&isign,iofst,1);
112 iofst=iofst+1;
113 gbit(cpack,&ival1,iofst,nbitsd-1);
114 iofst=iofst+nbitsd-1;
115 if (isign == 1) ival1=-ival1;
116 if (idrstmpl[16] == 2) {
117 gbit(cpack,&isign,iofst,1);
118 iofst=iofst+1;
119 gbit(cpack,&ival2,iofst,nbitsd-1);
120 iofst=iofst+nbitsd-1;
121 if (isign == 1) ival2=-ival2;
122 }
123 gbit(cpack,&isign,iofst,1);
124 iofst=iofst+1;
125 gbit(cpack,&minsd,iofst,nbitsd-1);
126 iofst=iofst+nbitsd-1;
127 if (isign == 1) minsd=-minsd;
128 }
129 else {
130 ival1=0;
131 ival2=0;
132 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 if (nbitsgref != 0) {
141 gbits(cpack,gref+0,iofst,nbitsgref,0,ngroups);
142 itemp=nbitsgref*ngroups;
143 iofst=iofst+itemp;
144 if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
145 }
146 else {
147 for (j=0;j<ngroups;j++)
148 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 if (nbitsgwidth != 0) {
155 gbits(cpack,gwidth+0,iofst,nbitsgwidth,0,ngroups);
156 itemp=nbitsgwidth*ngroups;
157 iofst=iofst+itemp;
158 if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
159 }
160 else {
161 for (j=0;j<ngroups;j++)
162 gwidth[j]=0;
163 }
164
165 for (j=0;j<ngroups;j++)
166 gwidth[j]=gwidth[j]+idrstmpl[10];
167
168 //
169 // Extract Each Group's length (number of values in each group)
170 //
171 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 if (nbitsglen != 0) {
175 gbits(cpack,glen,iofst,nbitsglen,0,ngroups);
176 itemp=nbitsglen*ngroups;
177 iofst=iofst+itemp;
178 if (itemp%8 != 0) iofst=iofst+(8-(itemp%8));
179 }
180 else {
181 for (j=0;j<ngroups;j++)
182 glen[j]=0;
183 }
184 for (j=0;j<ngroups;j++)
185 glen[j]=(glen[j]*idrstmpl[13])+idrstmpl[12];
186 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 totBit = 0;
192 totLen = 0;
193 for (j=0;j<ngroups;j++) {
194 totBit += (gwidth[j]*glen[j]);
195 totLen += glen[j];
196 }
197 if (totLen != ndpts) {
198 return 1;
199 }
200 if (totBit / 8. > lensec) {
201 return 1;
202 }
203 //
204 // For each group, unpack data values
205 //
206 if ( idrstmpl[6] == 0 ) { // no missing values
207 n=0;
208 for (j=0;j<ngroups;j++) {
209 if (gwidth[j] != 0) {
210 gbits(cpack,ifld+n,iofst,gwidth[j],0,glen[j]);
211 for (k=0;k<glen[j];k++) {
212 ifld[n]=ifld[n]+gref[j];
213 n=n+1;
214 }
215 }
216 else {
217 for (l=n;l<n+glen[j];l++) ifld[l]=gref[j];
218 n=n+glen[j];
219 }
220 iofst=iofst+(gwidth[j]*glen[j]);
221 }
222 }
223 else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
224 // missing values included
225 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 n=0;
229 non=0;
230 for (j=0;j<ngroups;j++) {
231 //printf(" SAGNGP %d %d %d %d\n",j,gwidth[j],glen[j],gref[j]);
232 if (gwidth[j] != 0) {
233 msng1=(g2int)int_power(2.0,gwidth[j])-1;
234 msng2=msng1-1;
235 gbits(cpack,ifld+n,iofst,gwidth[j],0,glen[j]);
236 iofst=iofst+(gwidth[j]*glen[j]);
237 for (k=0;k<glen[j];k++) {
238 if (ifld[n] == msng1) {
239 ifldmiss[n]=1;
240 //ifld[n]=0;
241 }
242 else if (idrstmpl[6]==2 && ifld[n]==msng2) {
243 ifldmiss[n]=2;
244 //ifld[n]=0;
245 }
246 else {
247 ifldmiss[n]=0;
248 ifld[non++]=ifld[n]+gref[j];
249 }
250 n++;
251 }
252 }
253 else {
254 msng1=(g2int)int_power(2.0,nbitsgref)-1;
255 msng2=msng1-1;
256 if (gref[j] == msng1) {
257 for (l=n;l<n+glen[j];l++) ifldmiss[l]=1;
258 }
259 else if (idrstmpl[6]==2 && gref[j]==msng2) {
260 for (l=n;l<n+glen[j];l++) ifldmiss[l]=2;
261 }
262 else {
263 for (l=n;l<n+glen[j];l++) ifldmiss[l]=0;
264 for (l=non;l<non+glen[j];l++) ifld[l]=gref[j];
265 non += glen[j];
266 }
267 n=n+glen[j];
268 }
269 }
270 }
271
272 if ( gref != 0 ) free(gref);
273 if ( gwidth != 0 ) free(gwidth);
274 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 if (idrsnum == 3) { // spatial differencing
281 if (idrstmpl[16] == 1) { // first order
282 ifld[0]=ival1;
283 if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
284 else itemp=non;
285 for (n=1;n<itemp;n++) {
286 ifld[n]=ifld[n]+minsd;
287 ifld[n]=ifld[n]+ifld[n-1];
288 }
289 }
290 else if (idrstmpl[16] == 2) { // second order
291 ifld[0]=ival1;
292 ifld[1]=ival2;
293 if ( idrstmpl[6] == 0 ) itemp=ndpts; // no missing values
294 else itemp=non;
295 for (n=2;n<itemp;n++) {
296 ifld[n]=ifld[n]+minsd;
297 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 if ( idrstmpl[6] == 0 ) { // no missing values
306 for (n=0;n<ndpts;n++) {
307 fld[n]=(((g2float)ifld[n]*bscale)+ref)*dscale;
308 }
309 }
310 else if ( idrstmpl[6]==1 || idrstmpl[6]==2 ) {
311 // missing values included
312 non=0;
313 for (n=0;n<ndpts;n++) {
314 if ( ifldmiss[n] == 0 ) {
315 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 else if ( ifldmiss[n] == 1 )
319 fld[n]=rmiss1;
320 else if ( ifldmiss[n] == 2 )
321 fld[n]=rmiss2;
322 }
323 if ( ifldmiss != 0 ) free(ifldmiss);
324 }
325
326 if ( ifld != 0 ) free(ifld);
327
328 return(0);
329
330 }
331