1 /*
2 * Copyright 2005-2018 ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence Version 2.0
5 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
8 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
9 */
10
11 #include "grib_api_internal.h"
12 #include <math.h>
13 /*
14 This is used by make_class.pl
15
16 START_CLASS_DEF
17 CLASS = accessor
18 SUPER = grib_accessor_class_data_simple_packing
19 IMPLEMENTS = unpack_double
20 IMPLEMENTS = value_count
21 IMPLEMENTS = init
22 MEMBERS= const char* GRIBEX_sh_bug_present
23 MEMBERS= const char* ieee_floats
24 MEMBERS= const char* laplacianOperatorIsSet
25 MEMBERS= const char* laplacianOperator
26 MEMBERS= const char* sub_j
27 MEMBERS= const char* sub_k
28 MEMBERS= const char* sub_m
29 MEMBERS= const char* pen_j
30 MEMBERS= const char* pen_k
31 MEMBERS= const char* pen_m
32 END_CLASS_DEF
33
34 */
35
36
37 /* START_CLASS_IMP */
38
39 /*
40
41 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
42 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
43 or edit "accessor.class" and rerun ./make_class.pl
44
45 */
46
47 static int unpack_double(grib_accessor*, double* val,size_t *len);
48 static int value_count(grib_accessor*,long*);
49 static void init(grib_accessor*,const long, grib_arguments* );
50 static void init_class(grib_accessor_class*);
51
52 typedef struct grib_accessor_data_sh_unpacked {
53 grib_accessor att;
54 /* Members defined in gen */
55 /* Members defined in values */
56 int carg;
57 const char* seclen;
58 const char* offsetdata;
59 const char* offsetsection;
60 int dirty;
61 /* Members defined in data_simple_packing */
62 int edition;
63 const char* units_factor;
64 const char* units_bias;
65 const char* changing_precision;
66 const char* number_of_values;
67 const char* bits_per_value;
68 const char* reference_value;
69 const char* binary_scale_factor;
70 const char* decimal_scale_factor;
71 /* Members defined in data_sh_unpacked */
72 const char* GRIBEX_sh_bug_present;
73 const char* ieee_floats;
74 const char* laplacianOperatorIsSet;
75 const char* laplacianOperator;
76 const char* sub_j;
77 const char* sub_k;
78 const char* sub_m;
79 const char* pen_j;
80 const char* pen_k;
81 const char* pen_m;
82 } grib_accessor_data_sh_unpacked;
83
84 extern grib_accessor_class* grib_accessor_class_data_simple_packing;
85
86 static grib_accessor_class _grib_accessor_class_data_sh_unpacked = {
87 &grib_accessor_class_data_simple_packing, /* super */
88 "data_sh_unpacked", /* name */
89 sizeof(grib_accessor_data_sh_unpacked), /* size */
90 0, /* inited */
91 &init_class, /* init_class */
92 &init, /* init */
93 0, /* post_init */
94 0, /* free mem */
95 0, /* describes himself */
96 0, /* get length of section */
97 0, /* get length of string */
98 &value_count, /* get number of values */
99 0, /* get number of bytes */
100 0, /* get offset to bytes */
101 0, /* get native type */
102 0, /* get sub_section */
103 0, /* grib_pack procedures long */
104 0, /* grib_pack procedures long */
105 0, /* grib_pack procedures long */
106 0, /* grib_unpack procedures long */
107 0, /* grib_pack procedures double */
108 &unpack_double, /* grib_unpack procedures double */
109 0, /* grib_pack procedures string */
110 0, /* grib_unpack procedures string */
111 0, /* grib_pack procedures bytes */
112 0, /* grib_unpack procedures bytes */
113 0, /* pack_expression */
114 0, /* notify_change */
115 0, /* update_size */
116 0, /* preferred_size */
117 0, /* resize */
118 0, /* nearest_smaller_value */
119 0, /* next accessor */
120 0, /* compare vs. another accessor */
121 0, /* unpack only ith value */
122 0, /* unpack a subarray */
123 0, /* clear */
124 };
125
126
127 grib_accessor_class* grib_accessor_class_data_sh_unpacked = &_grib_accessor_class_data_sh_unpacked;
128
129
init_class(grib_accessor_class * c)130 static void init_class(grib_accessor_class* c)
131 {
132 c->dump = (*(c->super))->dump;
133 c->next_offset = (*(c->super))->next_offset;
134 c->string_length = (*(c->super))->string_length;
135 c->byte_count = (*(c->super))->byte_count;
136 c->byte_offset = (*(c->super))->byte_offset;
137 c->get_native_type = (*(c->super))->get_native_type;
138 c->sub_section = (*(c->super))->sub_section;
139 c->pack_missing = (*(c->super))->pack_missing;
140 c->is_missing = (*(c->super))->is_missing;
141 c->pack_long = (*(c->super))->pack_long;
142 c->unpack_long = (*(c->super))->unpack_long;
143 c->pack_double = (*(c->super))->pack_double;
144 c->pack_string = (*(c->super))->pack_string;
145 c->unpack_string = (*(c->super))->unpack_string;
146 c->pack_bytes = (*(c->super))->pack_bytes;
147 c->unpack_bytes = (*(c->super))->unpack_bytes;
148 c->pack_expression = (*(c->super))->pack_expression;
149 c->notify_change = (*(c->super))->notify_change;
150 c->update_size = (*(c->super))->update_size;
151 c->preferred_size = (*(c->super))->preferred_size;
152 c->resize = (*(c->super))->resize;
153 c->nearest_smaller_value = (*(c->super))->nearest_smaller_value;
154 c->next = (*(c->super))->next;
155 c->compare = (*(c->super))->compare;
156 c->unpack_double_element = (*(c->super))->unpack_double_element;
157 c->unpack_double_subarray = (*(c->super))->unpack_double_subarray;
158 c->clear = (*(c->super))->clear;
159 }
160
161 /* END_CLASS_IMP */
162
163 typedef unsigned long (*encode_float_proc)(double);
164 typedef double (*decode_float_proc)(unsigned long);
165
init(grib_accessor * a,const long v,grib_arguments * args)166 static void init(grib_accessor* a,const long v, grib_arguments* args)
167 {
168 grib_accessor_data_sh_unpacked *self =(grib_accessor_data_sh_unpacked*)a;
169
170 self->GRIBEX_sh_bug_present = grib_arguments_get_name(a->parent->h,args,self->carg++);
171 self->ieee_floats = grib_arguments_get_name(a->parent->h,args,self->carg++);
172 self->laplacianOperatorIsSet = grib_arguments_get_name(a->parent->h,args,self->carg++);
173 self->laplacianOperator = grib_arguments_get_name(a->parent->h,args,self->carg++);
174 self->sub_j = grib_arguments_get_name(a->parent->h,args,self->carg++);
175 self->sub_k = grib_arguments_get_name(a->parent->h,args,self->carg++);
176 self->sub_m = grib_arguments_get_name(a->parent->h,args,self->carg++);
177 self->pen_j = grib_arguments_get_name(a->parent->h,args,self->carg++);
178 self->pen_k = grib_arguments_get_name(a->parent->h,args,self->carg++);
179 self->pen_m = grib_arguments_get_name(a->parent->h,args,self->carg++);
180
181 a->flags |= GRIB_ACCESSOR_FLAG_DATA;
182 a->length=0;
183 }
184
185
value_count(grib_accessor * a,long * count)186 static int value_count(grib_accessor* a,long* count)
187 {
188 grib_accessor_data_sh_unpacked *self =(grib_accessor_data_sh_unpacked*)a;
189 int ret = 0;
190
191 long sub_j= 0;
192 long sub_k= 0;
193 long sub_m= 0;
194
195 if((ret = grib_get_long_internal(a->parent->h,self->sub_j,&sub_j)) != GRIB_SUCCESS)
196 return ret;
197 if((ret = grib_get_long_internal(a->parent->h,self->sub_k,&sub_k)) != GRIB_SUCCESS)
198 return ret;
199 if((ret = grib_get_long_internal(a->parent->h,self->sub_m,&sub_m)) != GRIB_SUCCESS)
200 return ret;
201
202 if (sub_j != sub_k || sub_j!=sub_m ) {
203 grib_context_log(a->parent->h->context,GRIB_LOG_ERROR,"sub_j=%ld, sub_k=%ld, sub_m=%ld\n",sub_j,sub_k,sub_m);
204 Assert ((sub_j == sub_k) && (sub_j == sub_m));
205 }
206 *count=(sub_j+1)*(sub_j+2);
207 return ret;
208 }
209
210
unpack_double(grib_accessor * a,double * val,size_t * len)211 static int unpack_double(grib_accessor* a, double* val, size_t *len)
212 {
213 grib_accessor_data_sh_unpacked* self = (grib_accessor_data_sh_unpacked*)a;
214
215 size_t i = 0;
216 int ret = GRIB_SUCCESS;
217 long hcount = 0;
218 long lcount = 0;
219 long hpos = 0;
220 long lup = 0;
221 long mmax = 0;
222 long n_vals = 0;
223 double *scals = NULL;
224 /* double *pscals=NULL; */
225 double dummy=0;
226
227 double s = 0;
228 double d = 0;
229 double laplacianOperator = 0;
230 unsigned char* buf = NULL;
231 unsigned char* hres = NULL;
232 unsigned char* lres = NULL;
233 unsigned long packed_offset;
234 long lpos = 0;
235
236 long maxv = 0;
237 long GRIBEX_sh_bug_present =0;
238 long ieee_floats = 0;
239
240 long offsetdata = 0;
241 long bits_per_value = 0;
242 double reference_value = 0;
243 long binary_scale_factor = 0;
244 long decimal_scale_factor = 0;
245
246
247 long sub_j= 0;
248 long sub_k= 0;
249 long sub_m= 0;
250 long pen_j= 0;
251 long pen_k= 0;
252 long pen_m= 0;
253
254 double operat= 0;
255 int err=0;
256
257 decode_float_proc decode_float = NULL;
258
259 n_vals = 0;
260 err=grib_value_count(a,&n_vals);
261 if (err) return err;
262
263 if(*len < n_vals){
264 *len = n_vals;
265 return GRIB_ARRAY_TOO_SMALL;
266 }
267
268 if((ret = grib_get_long_internal(a->parent->h,self->offsetdata,&offsetdata))
269 != GRIB_SUCCESS) return ret;
270
271 if((ret = grib_get_long_internal(a->parent->h,self->GRIBEX_sh_bug_present,&GRIBEX_sh_bug_present))
272 != GRIB_SUCCESS)
273 return ret;
274
275 if((ret = grib_get_long_internal(a->parent->h,self->ieee_floats,&ieee_floats)) != GRIB_SUCCESS)
276 return ret;
277
278 if((ret = grib_get_long_internal(a->parent->h,self->sub_j,&sub_j)) != GRIB_SUCCESS)
279 return ret;
280 if((ret = grib_get_long_internal(a->parent->h,self->sub_k,&sub_k)) != GRIB_SUCCESS)
281 return ret;
282 if((ret = grib_get_long_internal(a->parent->h,self->sub_m,&sub_m)) != GRIB_SUCCESS)
283 return ret;
284 if((ret = grib_get_long_internal(a->parent->h,self->pen_j,&pen_j)) != GRIB_SUCCESS)
285 return ret;
286 if((ret = grib_get_long_internal(a->parent->h,self->pen_k,&pen_k)) != GRIB_SUCCESS)
287 return ret;
288 if((ret = grib_get_long_internal(a->parent->h,self->pen_m,&pen_m)) != GRIB_SUCCESS)
289 return ret;
290
291 self->dirty=0;
292
293 switch (ieee_floats) {
294 case 0:
295 decode_float=grib_long_to_ibm;
296 break;
297 case 1:
298 decode_float=grib_long_to_ieee;
299 break;
300 case 2:
301 decode_float=grib_long_to_ieee64;
302 break;
303 default:
304 return GRIB_NOT_IMPLEMENTED;
305 }
306
307 Assert (sub_j == sub_k);
308 Assert (sub_j == sub_m);
309 Assert (pen_j == pen_k);
310 Assert (pen_j == pen_m);
311
312 buf = (unsigned char*)a->parent->h->buffer->data;
313
314 maxv = pen_j+1;
315
316 buf += offsetdata;
317 hres = buf;
318 lres = buf;
319
320 packed_offset = offsetdata + 4*(sub_k+1)*(sub_k+2);
321
322 lpos = 8*(packed_offset-offsetdata);
323
324 s = grib_power(binary_scale_factor,2);
325 d = grib_power(-decimal_scale_factor,10) ;
326
327 scals = (double*)grib_context_malloc(a->parent->h->context,maxv*sizeof(double));
328 Assert(scals);
329 if((ret = grib_get_double_internal(a->parent->h,self->laplacianOperator,&laplacianOperator))
330 != GRIB_SUCCESS)
331 return ret;
332
333 scals[0] = 0;
334 for(i=1;i<maxv;i++){
335 operat = pow(i*(i+1),laplacianOperator);
336 if(operat != 0)
337 scals[i] = (1.0/operat);
338 else{
339 scals[i] = 0;
340 }
341 }
342
343
344 i=0;
345
346 while(maxv>0)
347 {
348 lup=mmax;
349 if(sub_k>=0)
350 {
351 for(hcount=0;hcount<sub_k+1;hcount++)
352 {
353 val[i++] = decode_float(grib_decode_unsigned_long(hres,&hpos,32))*d;
354 val[i++] = decode_float(grib_decode_unsigned_long(hres,&hpos,32))*d;
355
356 if (GRIBEX_sh_bug_present && hcount==sub_k){
357 /* bug in ecmwf data, last row (K+1)is scaled but should not */
358 val[i-2] *= scals[lup];
359 val[i-1] *= scals[lup];
360 }
361 lup++;
362 }
363 sub_k--;
364 }
365
366 /* pscals=scals+lup; */
367 for(lcount=hcount; lcount < maxv ; lcount++)
368 {
369 dummy = (double) ((grib_decode_unsigned_long(lres, &lpos,
370 bits_per_value)*s)+reference_value);
371 dummy = (double) ((grib_decode_unsigned_long(lres, &lpos,
372 bits_per_value)*s)+reference_value);
373 lup++;
374 }
375
376 maxv--;
377 hcount=0;
378 mmax++;
379 }
380
381 Assert(*len >= i);
382 *len = n_vals;
383
384 if(d != 1) {
385 for(i=0;i<*len;i++)
386 val[i++] *= d;
387 }
388
389 (void)dummy; /* suppress gcc warning */
390 grib_context_free(a->parent->h->context,scals);
391
392 return ret;
393
394 }
395
396