Line data Source code
1 : #include <float.h>
2 : #include <memory.h>
3 : #include <stdio.h>
4 : #include <stdlib.h>
5 : #include <string.h>
6 : #include "grib2.h"
7 :
8 : #ifndef DoubleToFloatClamp_defined
9 : #define DoubleToFloatClamp_defined
10 26198 : static float DoubleToFloatClamp(double val) {
11 26198 : if (val >= FLT_MAX) return FLT_MAX;
12 26198 : if (val <= -FLT_MAX) return -FLT_MAX;
13 26198 : return (float)val;
14 : }
15 : #endif
16 :
17 163 : g2int g2_unpack7(unsigned char *cgrib,g2int cgrib_length,g2int *iofst,g2int igdsnum,g2int *igdstmpl,
18 : g2int idrsnum,g2int *idrstmpl,g2int ndpts,g2float **fld)
19 : //$$$ SUBPROGRAM DOCUMENTATION BLOCK
20 : // . . . .
21 : // SUBPROGRAM: g2_unpack7
22 : // PRGMMR: Gilbert ORG: W/NP11 DATE: 2002-10-31
23 : //
24 : // ABSTRACT: This subroutine unpacks Section 7 (Data Section)
25 : // as defined in GRIB Edition 2.
26 : //
27 : // PROGRAM HISTORY LOG:
28 : // 2002-10-31 Gilbert
29 : // 2002-12-20 Gilbert - Added GDT info to arguments
30 : // and added 5.51 processing.
31 : // 2003-08-29 Gilbert - Added support for new templates using
32 : // PNG and JPEG2000 algorithms/templates.
33 : // 2004-11-29 Gilbert - JPEG2000 now allowed to use WMO Template no. 5.40
34 : // PNG now allowed to use WMO Template no. 5.41
35 : // 2004-12-16 Taylor - Added check on comunpack return code.
36 : // 2008-12-23 Wesley - Initialize Number of data points unpacked
37 : //
38 : // USAGE: int g2_unpack7(unsigned char *cgrib,g2int *iofst,g2int igdsnum,
39 : // g2int *igdstmpl, g2int idrsnum,
40 : // g2int *idrstmpl, g2int ndpts,g2float **fld)
41 : // INPUT ARGUMENTS:
42 : // cgrib - char array containing Section 7 of the GRIB2 message
43 : // iofst - Bit offset of the beginning of Section 7 in cgrib.
44 : // igdsnum - Grid Definition Template Number ( see Code Table 3.0)
45 : // ( Only used for DRS Template 5.51 )
46 : // igdstmpl - Pointer to an integer array containing the data values for
47 : // the specified Grid Definition
48 : // Template ( N=igdsnum ). Each element of this integer
49 : // array contains an entry (in the order specified) of Grid
50 : // Definition Template 3.N
51 : // ( Only used for DRS Template 5.51 )
52 : // idrsnum - Data Representation Template Number ( see Code Table 5.0)
53 : // idrstmpl - Pointer to an integer array containing the data values for
54 : // the specified Data Representation
55 : // Template ( N=idrsnum ). Each element of this integer
56 : // array contains an entry (in the order specified) of Data
57 : // Representation Template 5.N
58 : // ndpts - Number of data points unpacked and returned.
59 : //
60 : // OUTPUT ARGUMENTS:
61 : // iofst - Bit offset at the end of Section 7, returned.
62 : // fld - Pointer to a float array containing the unpacked data field.
63 : //
64 : // RETURN VALUES:
65 : // ierr - Error return code.
66 : // 0 = no error
67 : // 2 = Not section 7
68 : // 4 = Unrecognized Data Representation Template
69 : // 5 = need one of GDT 3.50 through 3.53 to decode DRT 5.51
70 : // 6 = memory allocation error
71 : // 7 = corrupt section 7.
72 : //
73 : // REMARKS: None
74 : //
75 : // ATTRIBUTES:
76 : // LANGUAGE: C
77 : // MACHINE:
78 : //
79 : //$$$//
80 : {
81 : g2int ierr,isecnum;
82 : g2int ipos,lensec;
83 : g2float *lfld;
84 :
85 163 : ierr=0;
86 163 : *fld=0; //NULL
87 :
88 163 : gbit(cgrib,&lensec,*iofst,32); // Get Length of Section
89 163 : *iofst=*iofst+32;
90 163 : gbit(cgrib,&isecnum,*iofst,8); // Get Section Number
91 163 : *iofst=*iofst+8;
92 :
93 163 : if ( isecnum != 7 ) {
94 0 : ierr=2;
95 : //fprintf(stderr,"g2_unpack7: Not Section 7 data.\n");
96 0 : return(ierr);
97 : }
98 :
99 163 : ipos=(*iofst/8);
100 163 : if( ipos >= cgrib_length ) {
101 0 : return 7;
102 : }
103 163 : if (idrsnum == 40 || idrsnum == 40000) /* added by GDAL */
104 : {
105 18 : *fld= lfld = 0;
106 : }
107 : else
108 : {
109 145 : lfld=(g2float *)calloc(ndpts,sizeof(g2float));
110 145 : if (lfld == 0) {
111 0 : ierr=6;
112 0 : return(ierr);
113 : }
114 : else {
115 145 : *fld=lfld;
116 : }
117 : }
118 :
119 163 : if (idrsnum == 0)
120 36 : simunpack(cgrib+ipos,cgrib_length-ipos,idrstmpl,ndpts,lfld);
121 127 : else if (idrsnum == 2 || idrsnum == 3) {
122 34 : if (comunpack(cgrib+ipos,cgrib_length-ipos,lensec,idrsnum,idrstmpl,ndpts,lfld) != 0) {
123 0 : return 7;
124 : }
125 : }
126 93 : else if( idrsnum == 4 ) {
127 : // Grid point data - IEEE Floating Point Data
128 : static const int one = 1;
129 58 : int is_lsb = *((char*)&one) == 1;
130 58 : if (idrstmpl[0] == 1) {
131 : // IEEE754 single precision
132 8 : if( cgrib_length-ipos < 4 * ndpts )
133 0 : return 7;
134 8 : memcpy(lfld, cgrib+ipos, 4 * ndpts );
135 8 : if( is_lsb )
136 : {
137 : int i;
138 8 : unsigned char* ch_fld = (unsigned char*) lfld;
139 1616 : for(i=0;i<ndpts;i++)
140 : {
141 1608 : unsigned char temp = ch_fld[i*4];
142 1608 : ch_fld[i*4] = ch_fld[i*4+3];
143 1608 : ch_fld[i*4+3] = temp;
144 1608 : temp = ch_fld[i*4+1];
145 1608 : ch_fld[i*4+1] = ch_fld[i*4+2];
146 1608 : ch_fld[i*4+2] = temp;
147 : }
148 : }
149 : }
150 50 : else if( idrstmpl[0] == 2) {
151 : // IEEE754 double precision
152 : // FIXME? due to the interface: we downgrade it to float
153 : int i;
154 50 : unsigned char* src = cgrib+ipos;
155 50 : if( cgrib_length-ipos < 8 * ndpts )
156 0 : return 7;
157 50 : if( is_lsb )
158 : {
159 26248 : for(i=0;i<ndpts;i++)
160 : {
161 : unsigned char temp[8];
162 : double d;
163 : {
164 : int j;
165 235782 : for(j = 0; j < 8; j++ )
166 209584 : temp[j] = src[i * 8 + 7 - j];
167 : }
168 26198 : memcpy(&d, temp, 8);
169 26198 : lfld[i] = DoubleToFloatClamp(d);
170 : }
171 : }
172 : else
173 : {
174 0 : for(i=0;i<ndpts;i++)
175 : {
176 : double d;
177 0 : memcpy(&d, src + i * 8, 8);
178 0 : lfld[i] = DoubleToFloatClamp(d);
179 : }
180 : }
181 : }
182 : else
183 : {
184 0 : fprintf(stderr,"g2_unpack7: Invalid precision=%d for Data Section 5.4.\n", idrstmpl[0]);
185 : }
186 : }
187 35 : else if (idrsnum == 50) { // Spectral Simple
188 0 : if( ndpts > 0 )
189 : {
190 0 : simunpack(cgrib+ipos,cgrib_length-ipos,idrstmpl,ndpts-1,lfld+1);
191 0 : rdieee(idrstmpl+4,lfld+0,1);
192 : }
193 : }
194 35 : else if (idrsnum == 51) // Spectral complex
195 0 : if ( igdsnum>=50 && igdsnum <=53 )
196 0 : specunpack(cgrib+ipos,idrstmpl,ndpts,igdstmpl[0],igdstmpl[2],igdstmpl[2],lfld);
197 : else {
198 0 : fprintf(stderr,"g2_unpack7: Cannot use GDT 3.%d to unpack Data Section 5.51.\n",(int)igdsnum);
199 0 : ierr=5;
200 0 : if ( lfld != 0 ) free(lfld);
201 0 : *fld=0; //NULL
202 0 : return(ierr);
203 : }
204 35 : else if (idrsnum == 40 || idrsnum == 40000) {
205 18 : if( jpcunpack(cgrib+ipos,lensec-5,idrstmpl,ndpts,fld) != 0 )
206 : {
207 0 : ierr=7;
208 0 : if ( *fld != 0 ) free(*fld);
209 0 : *fld=0; //NULL
210 0 : return(ierr);
211 : }
212 : }
213 : #ifdef USE_PNG
214 17 : else if (idrsnum == 41 || idrsnum == 40010) {
215 16 : pngunpack(cgrib+ipos,lensec-5,idrstmpl,ndpts,lfld);
216 : }
217 : #endif /* USE_PNG */
218 1 : else if (idrsnum == 42) {
219 : #ifdef USE_AEC
220 1 : aecunpack(cgrib+ipos,lensec-5,idrstmpl,ndpts,lfld);
221 : #else
222 : fprintf(stderr,"g2_unpack7: Data Representation Template 5.42 decoding requires building against libaec.\n");
223 : ierr=4;
224 : if ( lfld != 0 ) free(lfld);
225 : *fld=0; //NULL
226 : return(ierr);
227 : #endif
228 : }
229 : else {
230 0 : fprintf(stderr,"g2_unpack7: Data Representation Template 5.%d not yet implemented.\n",(int)idrsnum);
231 0 : ierr=4;
232 0 : if ( lfld != 0 ) free(lfld);
233 0 : *fld=0; //NULL
234 0 : return(ierr);
235 : }
236 :
237 163 : *iofst=*iofst+(8*lensec);
238 :
239 163 : return(ierr); // End of Section 7 processing
240 :
241 : }
|
