diff --git a/src/accessor/grib_accessor_class_data_apply_bitmap.cc b/src/accessor/grib_accessor_class_data_apply_bitmap.cc index 86e54ff63..cf3d5f4a5 100644 --- a/src/accessor/grib_accessor_class_data_apply_bitmap.cc +++ b/src/accessor/grib_accessor_class_data_apply_bitmap.cc @@ -11,13 +11,14 @@ #include "grib_accessor_class_data_apply_bitmap.h" -grib_accessor_class_data_apply_bitmap_t _grib_accessor_class_data_apply_bitmap{"data_apply_bitmap"}; +grib_accessor_class_data_apply_bitmap_t _grib_accessor_class_data_apply_bitmap{ "data_apply_bitmap" }; grib_accessor_class* grib_accessor_class_data_apply_bitmap = &_grib_accessor_class_data_apply_bitmap; -void grib_accessor_class_data_apply_bitmap_t::init(grib_accessor* a, const long v, grib_arguments* args){ +void grib_accessor_class_data_apply_bitmap_t::init(grib_accessor* a, const long v, grib_arguments* args) +{ grib_accessor_class_gen_t::init(a, v, args); - int n = 0; + int n = 0; grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; self->coded_values = grib_arguments_get_name(grib_handle_of_accessor(a), args, n++); @@ -29,14 +30,16 @@ void grib_accessor_class_data_apply_bitmap_t::init(grib_accessor* a, const long a->length = 0; } -void grib_accessor_class_data_apply_bitmap_t::dump(grib_accessor* a, grib_dumper* dumper){ +void grib_accessor_class_data_apply_bitmap_t::dump(grib_accessor* a, grib_dumper* dumper) +{ grib_dump_values(dumper, a); } -int grib_accessor_class_data_apply_bitmap_t::value_count(grib_accessor* a, long* count){ +int grib_accessor_class_data_apply_bitmap_t::value_count(grib_accessor* a, long* count) +{ grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; - size_t len = 0; - int ret = 0; + size_t len = 0; + int ret = GRIB_SUCCESS; if (grib_find_accessor(grib_handle_of_accessor(a), self->bitmap)) ret = grib_get_size(grib_handle_of_accessor(a), self->bitmap, &len); @@ -48,17 +51,18 @@ int grib_accessor_class_data_apply_bitmap_t::value_count(grib_accessor* a, long* return ret; } -int grib_accessor_class_data_apply_bitmap_t::unpack_double_element(grib_accessor* a, size_t idx, double* val){ +int grib_accessor_class_data_apply_bitmap_t::unpack_double_element(grib_accessor* a, size_t idx, double* val) +{ grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); - int err = 0, i = 0; - size_t cidx = 0; + grib_handle* gh = grib_handle_of_accessor(a); + size_t i = 0, cidx = 0; double missing_value = 0; double* bvals = NULL; size_t n_vals = 0; long nn = 0; - err = a->value_count(&nn); n_vals = nn; + int err = a->value_count(&nn); + n_vals = nn; if (err) return err; @@ -93,19 +97,21 @@ int grib_accessor_class_data_apply_bitmap_t::unpack_double_element(grib_accessor return grib_get_double_element_internal(gh, self->coded_values, cidx, val); } -int grib_accessor_class_data_apply_bitmap_t::unpack_double_element_set(grib_accessor* a, const size_t* index_array, size_t len, double* val_array){ +int grib_accessor_class_data_apply_bitmap_t::unpack_double_element_set(grib_accessor* a, const size_t* index_array, size_t len, double* val_array) +{ grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); + grib_handle* gh = grib_handle_of_accessor(a); int err = 0, all_missing = 1; - size_t cidx = 0; /* index into the coded_values array */ - size_t* cidx_array = NULL; /* array of indexes into the coded_values */ - double* cval_array = NULL; /* array of values of the coded_values */ + size_t cidx = 0; /* index into the coded_values array */ + size_t* cidx_array = NULL; /* array of indexes into the coded_values */ + double* cval_array = NULL; /* array of values of the coded_values */ double missing_value = 0; double* bvals = NULL; size_t n_vals = 0, i = 0, j = 0, idx = 0, count_1s = 0, ci = 0; long nn = 0; - err = a->value_count(&nn); n_vals = nn; + err = a->value_count(&nn); + n_vals = nn; if (err) return err; if (!grib_find_accessor(gh, self->bitmap)) @@ -119,7 +125,8 @@ int grib_accessor_class_data_apply_bitmap_t::unpack_double_element_set(grib_acce for (i = 0; i < len; i++) { if (val_array[i] == 0) { val_array[i] = missing_value; - } else { + } + else { all_missing = 0; count_1s++; } @@ -144,7 +151,7 @@ int grib_accessor_class_data_apply_bitmap_t::unpack_double_element_set(grib_acce ci = 0; for (i = 0; i < len; i++) { if (val_array[i] == 1) { - idx = index_array[i]; + idx = index_array[i]; cidx = 0; for (j = 0; j < idx; j++) { cidx += bvals[j]; @@ -171,17 +178,18 @@ int grib_accessor_class_data_apply_bitmap_t::unpack_double_element_set(grib_acce return GRIB_SUCCESS; } -int grib_accessor_class_data_apply_bitmap_t::pack_double(grib_accessor* a, const double* val, size_t* len){ +int grib_accessor_class_data_apply_bitmap_t::pack_double(grib_accessor* a, const double* val, size_t* len) +{ grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; - int err = 0; - size_t bmaplen = *len; - long coded_n_vals = 0; - double* coded_vals = NULL; - long i = 0; - long j = 0; - double missing_value = 0; - grib_handle* hand = grib_handle_of_accessor(a); - grib_context* ctxt = a->context; + int err = 0; + size_t bmaplen = *len; + long coded_n_vals = 0; + double* coded_vals = NULL; + long i = 0; + long j = 0; + double missing_value = 0; + grib_handle* hand = grib_handle_of_accessor(a); + grib_context* ctxt = a->context; if (*len == 0) return GRIB_NO_VALUES; @@ -231,22 +239,106 @@ int grib_accessor_class_data_apply_bitmap_t::pack_double(grib_accessor* a, const } -int grib_accessor_class_data_apply_bitmap_t::unpack_double(grib_accessor* a, double* val, size_t* len){ +template +static int unpack(grib_accessor* a, T* val, size_t* len) +{ + static_assert(std::is_floating_point::value, "Requires floating point numbers"); + grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; + + size_t i = 0; + size_t j = 0; + size_t n_vals = 0; + long nn = 0; + size_t coded_n_vals = 0; + T* coded_vals = NULL; + double missing_value = 0; + + int err = a->value_count(&nn); + n_vals = nn; + if (err) + return err; + + if (!grib_find_accessor(grib_handle_of_accessor(a), self->bitmap)) + return grib_get_array(grib_handle_of_accessor(a), self->coded_values, val, len); + + if ((err = grib_get_size(grib_handle_of_accessor(a), self->coded_values, &coded_n_vals)) != GRIB_SUCCESS) + return err; + + if ((err = grib_get_double_internal(grib_handle_of_accessor(a), self->missing_value, &missing_value)) != GRIB_SUCCESS) + return err; + + if (*len < n_vals) { + *len = n_vals; + return GRIB_ARRAY_TOO_SMALL; + } + + if (coded_n_vals == 0) { + for (i = 0; i < n_vals; i++) + val[i] = missing_value; + + *len = n_vals; + return GRIB_SUCCESS; + } + + if ((err = grib_get_array_internal(grib_handle_of_accessor(a), self->bitmap, val, &n_vals)) != GRIB_SUCCESS) + return err; + + coded_vals = (T*)grib_context_malloc(a->context, coded_n_vals * sizeof(T)); + if (coded_vals == NULL) + return GRIB_OUT_OF_MEMORY; + + if ((err = grib_get_array(grib_handle_of_accessor(a), self->coded_values, coded_vals, &coded_n_vals)) != GRIB_SUCCESS) { + grib_context_free(a->context, coded_vals); + return err; + } + + grib_context_log(a->context, GRIB_LOG_DEBUG, + "grib_accessor_class_data_apply_bitmap: %s : creating %s, %d values", + __func__, a->name, n_vals); + + for (i = 0; i < n_vals; i++) { + if (val[i] == 0) { + val[i] = missing_value; + } + else { + val[i] = coded_vals[j++]; + if (j > coded_n_vals) { + grib_context_free(a->context, coded_vals); + grib_context_log(a->context, GRIB_LOG_ERROR, + "grib_accessor_class_data_apply_bitmap [%s]:" + " %s : number of coded values does not match bitmap %ld %ld", + a->name, __func__, coded_n_vals, n_vals); + + return GRIB_ARRAY_TOO_SMALL; + } + } + } + + *len = n_vals; + + grib_context_free(a->context, coded_vals); + return err; +} + +int grib_accessor_class_data_apply_bitmap_t::unpack_double(grib_accessor* a, double* val, size_t* len) +{ return unpack(a, val, len); } -int grib_accessor_class_data_apply_bitmap_t::unpack_float(grib_accessor* a, float* val, size_t* len){ +int grib_accessor_class_data_apply_bitmap_t::unpack_float(grib_accessor* a, float* val, size_t* len) +{ return unpack(a, val, len); } -int grib_accessor_class_data_apply_bitmap_t::get_native_type(grib_accessor* a){ - //grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; - //return grib_accessor_get_native_type(grib_find_accessor(grib_handle_of_accessor(a),self->coded_values)); - +int grib_accessor_class_data_apply_bitmap_t::get_native_type(grib_accessor* a) +{ + // grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; + // return grib_accessor_get_native_type(grib_find_accessor(grib_handle_of_accessor(a),self->coded_values)); return GRIB_TYPE_DOUBLE; } -int grib_accessor_class_data_apply_bitmap_t::compare(grib_accessor* a, grib_accessor* b){ +int grib_accessor_class_data_apply_bitmap_t::compare(grib_accessor* a, grib_accessor* b) +{ int retval = GRIB_SUCCESS; double* aval = 0; double* bval = 0; @@ -256,11 +348,13 @@ int grib_accessor_class_data_apply_bitmap_t::compare(grib_accessor* a, grib_acce int err = 0; long count = 0; - err = a->value_count(&count); if (err) + err = a->value_count(&count); + if (err) return err; alen = count; - err = b->value_count(&count); if (err) + err = b->value_count(&count); + if (err) return err; blen = count; @@ -270,9 +364,10 @@ int grib_accessor_class_data_apply_bitmap_t::compare(grib_accessor* a, grib_acce aval = (double*)grib_context_malloc(a->context, alen * sizeof(double)); bval = (double*)grib_context_malloc(b->context, blen * sizeof(double)); - a->unpack_double(aval, &alen); b->unpack_double(bval, &blen); + a->unpack_double(aval, &alen); + b->unpack_double(bval, &blen); retval = GRIB_SUCCESS; - for (size_t i=0; i int unpack(grib_accessor*, T*, size_t*); }; - - -template -int grib_accessor_class_data_apply_bitmap_t::unpack(grib_accessor* a, T* val, size_t* len){ - static_assert(std::is_floating_point::value, "Requires floating point numbers"); - grib_accessor_data_apply_bitmap_t* self = (grib_accessor_data_apply_bitmap_t*)a; - - size_t i = 0; - size_t j = 0; - size_t n_vals = 0; - long nn = 0; - int err = 0; - size_t coded_n_vals = 0; - T* coded_vals = NULL; - double missing_value = 0; - - err = a->value_count(&nn); n_vals = nn; - if (err) - return err; - - if (!grib_find_accessor(grib_handle_of_accessor(a), self->bitmap)) - return grib_get_array(grib_handle_of_accessor(a), self->coded_values, val, len); - - if ((err = grib_get_size(grib_handle_of_accessor(a), self->coded_values, &coded_n_vals)) != GRIB_SUCCESS) - return err; - - if ((err = grib_get_double_internal(grib_handle_of_accessor(a), self->missing_value, &missing_value)) != GRIB_SUCCESS) - return err; - - if (*len < n_vals) { - *len = n_vals; - return GRIB_ARRAY_TOO_SMALL; - } - - if (coded_n_vals == 0) { - for (i = 0; i < n_vals; i++) - val[i] = missing_value; - - *len = n_vals; - return GRIB_SUCCESS; - } - - if ((err = grib_get_array_internal(grib_handle_of_accessor(a), self->bitmap, val, &n_vals)) != GRIB_SUCCESS) - return err; - - coded_vals = (T*)grib_context_malloc(a->context, coded_n_vals * sizeof(T)); - if (coded_vals == NULL) - return GRIB_OUT_OF_MEMORY; - - if ((err = grib_get_array(grib_handle_of_accessor(a), self->coded_values, coded_vals, &coded_n_vals)) != GRIB_SUCCESS) { - grib_context_free(a->context, coded_vals); - return err; - } - - grib_context_log(a->context, GRIB_LOG_DEBUG, - "grib_accessor_class_data_apply_bitmap: %s : creating %s, %d values", - __func__, - a->name, n_vals); - - for (i = 0; i < n_vals; i++) { - if (val[i] == 0) { - val[i] = missing_value; - } - else { - val[i] = coded_vals[j++]; - if (j > coded_n_vals) { - grib_context_free(a->context, coded_vals); - grib_context_log(a->context, GRIB_LOG_ERROR, - "grib_accessor_class_data_apply_bitmap [%s]:" - " %s : number of coded values does not match bitmap %ld %ld", - a->name, __func__, coded_n_vals, n_vals); - - return GRIB_ARRAY_TOO_SMALL; - } - } - } - - *len = n_vals; - - grib_context_free(a->context, coded_vals); - return err; -} - -#endif /* eccodes_accessor_data_apply_bitmap_h */ diff --git a/src/accessor/grib_accessor_class_data_complex_packing.cc b/src/accessor/grib_accessor_class_data_complex_packing.cc index 77a3ab822..81e402303 100644 --- a/src/accessor/grib_accessor_class_data_complex_packing.cc +++ b/src/accessor/grib_accessor_class_data_complex_packing.cc @@ -11,13 +11,14 @@ #include "grib_accessor_class_data_complex_packing.h" -grib_accessor_class_data_complex_packing_t _grib_accessor_class_data_complex_packing{"data_complex_packing"}; +grib_accessor_class_data_complex_packing_t _grib_accessor_class_data_complex_packing{ "data_complex_packing" }; grib_accessor_class* grib_accessor_class_data_complex_packing = &_grib_accessor_class_data_complex_packing; -void grib_accessor_class_data_complex_packing_t::init(grib_accessor* a, const long v, grib_arguments* args){ +void grib_accessor_class_data_complex_packing_t::init(grib_accessor* a, const long v, grib_arguments* args) +{ grib_accessor_class_data_simple_packing_t::init(a, v, args); grib_accessor_data_complex_packing_t* self = (grib_accessor_data_complex_packing_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); + grib_handle* gh = grib_handle_of_accessor(a); self->GRIBEX_sh_bug_present = grib_arguments_get_name(gh, args, self->carg++); self->ieee_floats = grib_arguments_get_name(gh, args, self->carg++); @@ -33,14 +34,16 @@ void grib_accessor_class_data_complex_packing_t::init(grib_accessor* a, const lo a->flags |= GRIB_ACCESSOR_FLAG_DATA; } -int grib_accessor_class_data_complex_packing_t::value_count(grib_accessor* a, long* count){ +int grib_accessor_class_data_complex_packing_t::value_count(grib_accessor* a, long* count) +{ grib_accessor_data_complex_packing_t* self = (grib_accessor_data_complex_packing_t*)a; - int ret = 0; + + int ret = GRIB_SUCCESS; grib_handle* gh = grib_handle_of_accessor(a); - long pen_j = 0; - long pen_k = 0; - long pen_m = 0; - *count = 0; + long pen_j = 0; + long pen_k = 0; + long pen_m = 0; + *count = 0; if (a->length == 0) return 0; @@ -62,7 +65,8 @@ int grib_accessor_class_data_complex_packing_t::value_count(grib_accessor* a, lo return ret; } -double calculate_pfactor(const grib_context* ctx, const double* spectralField, long fieldTruncation, long subsetTruncation){ +double calculate_pfactor(const grib_context* ctx, const double* spectralField, long fieldTruncation, long subsetTruncation) +{ /*long n_vals = ((fieldTruncation+1)*(fieldTruncation+2));*/ long loop, index, m, n = 0; double pFactor, zeps = 1.0e-15; @@ -95,7 +99,6 @@ double calculate_pfactor(const grib_context* ctx, const double* spectralField, l /* * Form norms for the rows which contain part of the unscaled subset. */ - index = -2; for (m = 0; m < subsetTruncation; m++) { for (n = m; n <= fieldTruncation; n++) { @@ -110,7 +113,6 @@ double calculate_pfactor(const grib_context* ctx, const double* spectralField, l /* * Form norms for the rows which do not contain part of the unscaled subset. */ - for (m = subsetTruncation; m <= fieldTruncation; m++) { for (n = m; n <= fieldTruncation; n++) { index += 2; @@ -146,7 +148,6 @@ double calculate_pfactor(const grib_context* ctx, const double* spectralField, l /* * Perform a least square fit for the equation */ - for (loop = ismin; loop <= ismax; loop++) { x = log((double)(loop * (loop + 1))); y = log(norms[loop]); @@ -168,10 +169,11 @@ double calculate_pfactor(const grib_context* ctx, const double* spectralField, l return pFactor; } -int grib_accessor_class_data_complex_packing_t::pack_double(grib_accessor* a, const double* val, size_t* len){ +int grib_accessor_class_data_complex_packing_t::pack_double(grib_accessor* a, const double* val, size_t* len) +{ grib_accessor_data_complex_packing_t* self = (grib_accessor_data_complex_packing_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); - const char* cclass_name = a->cclass->name; + grib_handle* gh = grib_handle_of_accessor(a); + const char* cclass_name = a->cclass->name; size_t i = 0; int ret = GRIB_SUCCESS; @@ -519,7 +521,7 @@ int grib_accessor_class_data_complex_packing_t::pack_double(grib_accessor* a, co grib_get_double_internal(gh, self->reference_value, &ref); if (ref != reference_value) { grib_context_log(a->context, GRIB_LOG_ERROR, "%s %s: %s (ref=%.10e != reference_value=%.10e)", - cclass_name, __func__, self->reference_value, ref, reference_value); + cclass_name, __func__, self->reference_value, ref, reference_value); return GRIB_INTERNAL_ERROR; } } @@ -536,34 +538,253 @@ int grib_accessor_class_data_complex_packing_t::pack_double(grib_accessor* a, co return ret; } +template +static int unpack_real(grib_accessor* a, T* val, size_t* len) +{ + static_assert(std::is_floating_point::value, "Requires floating point numbers"); + grib_accessor_data_complex_packing_t* self = (grib_accessor_data_complex_packing_t*)a; + grib_handle* gh = grib_handle_of_accessor(a); + const char* cclass_name = a->cclass->name; + + size_t i = 0; + int ret = GRIB_SUCCESS; + long hcount = 0; + long lcount = 0; + long hpos = 0; + long lup = 0; + long mmax = 0; + long n_vals = 0; + T* scals = NULL; + T *pscals = NULL, *pval = NULL; + + T s = 0; + T d = 0; + T laplacianOperator = 0; + unsigned char* buf = NULL; + unsigned char* hres = NULL; + unsigned char* lres = NULL; + unsigned long packed_offset; + long lpos = 0; + long maxv = 0; + long GRIBEX_sh_bug_present = 0; + long ieee_floats = 0; -int grib_accessor_class_data_complex_packing_t::unpack_double(grib_accessor* a, double* val, size_t* len){ - return unpack(a, val, len); + long offsetdata = 0; + long bits_per_value = 0; + T reference_value = 0; + long binary_scale_factor = 0; + long decimal_scale_factor = 0; + + long sub_j = 0; + long sub_k = 0; + long sub_m = 0; + long pen_j = 0; + long pen_k = 0; + long pen_m = 0; + + T operat = 0; + int bytes; + int err = 0; + double tmp; + + decode_float_proc decode_float = NULL; + + err = a->value_count(&n_vals); + if (err) + return err; + + if (*len < n_vals) { + *len = n_vals; + return GRIB_ARRAY_TOO_SMALL; + } + + if ((ret = grib_get_long_internal(gh, self->offsetdata, &offsetdata)) != GRIB_SUCCESS) + return ret; + if ((ret = grib_get_long_internal(gh, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS) + return ret; + if ((ret = grib_get_double_internal(gh, self->reference_value, &tmp)) != GRIB_SUCCESS) + return ret; + reference_value = tmp; + if ((ret = grib_get_long_internal(gh, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS) + return ret; + + if ((ret = grib_get_long_internal(gh, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS) + return ret; + + if ((ret = grib_get_long_internal(gh, self->GRIBEX_sh_bug_present, &GRIBEX_sh_bug_present)) != GRIB_SUCCESS) + return ret; + + /* ECC-774: don't use grib_get_long_internal */ + if ((ret = grib_get_long(gh, self->ieee_floats, &ieee_floats)) != GRIB_SUCCESS) + return ret; + + if ((ret = grib_get_double_internal(gh, self->laplacianOperator, &tmp)) != GRIB_SUCCESS) + return ret; + laplacianOperator = tmp; + + if ((ret = grib_get_long_internal(gh, self->sub_j, &sub_j)) != GRIB_SUCCESS) + return ret; + if ((ret = grib_get_long_internal(gh, self->sub_k, &sub_k)) != GRIB_SUCCESS) + return ret; + if ((ret = grib_get_long_internal(gh, self->sub_m, &sub_m)) != GRIB_SUCCESS) + return ret; + if ((ret = grib_get_long_internal(gh, self->pen_j, &pen_j)) != GRIB_SUCCESS) + return ret; + if ((ret = grib_get_long_internal(gh, self->pen_k, &pen_k)) != GRIB_SUCCESS) + return ret; + if ((ret = grib_get_long_internal(gh, self->pen_m, &pen_m)) != GRIB_SUCCESS) + return ret; + + self->dirty = 0; + + switch (ieee_floats) { + case 0: + decode_float = grib_long_to_ibm; + bytes = 4; + break; + case 1: + decode_float = grib_long_to_ieee; + bytes = 4; + break; + case 2: + decode_float = grib_long_to_ieee64; + bytes = 8; + break; + default: + return GRIB_NOT_IMPLEMENTED; + } + + if (sub_j != sub_k || sub_j != sub_m || pen_j != pen_k || pen_j != pen_m) { + grib_context_log(a->context, GRIB_LOG_ERROR, "%s: Invalid pentagonal resolution parameters", cclass_name); + return GRIB_DECODING_ERROR; + } + + buf = (unsigned char*)gh->buffer->data; + + maxv = pen_j + 1; + + buf += a->byte_offset(); + hres = buf; + lres = buf; + + if (pen_j == sub_j) { + n_vals = (pen_j + 1) * (pen_j + 2); + d = codes_power(-decimal_scale_factor, 10); + + grib_ieee_decode_array(a->context, buf, n_vals, bytes, val); + if (d) { + for (i = 0; i < n_vals; i++) + val[i] *= d; + } + return 0; + } + + packed_offset = a->byte_offset() + bytes * (sub_k + 1) * (sub_k + 2); + lpos = 8 * (packed_offset - offsetdata); + + s = codes_power(binary_scale_factor, 2); + d = codes_power(-decimal_scale_factor, 10); + + scals = (T*)grib_context_malloc(a->context, maxv * sizeof(T)); + if (!scals) return GRIB_OUT_OF_MEMORY; + + scals[0] = 0; + for (i = 1; i < maxv; i++) { + operat = pow(i * (i + 1), laplacianOperator); + if (operat != 0) + scals[i] = (1.0 / operat); + else { + grib_context_log(a->context, GRIB_LOG_WARNING, + "%s: Problem with operator div by zero at index %d of %d", cclass_name, i, maxv); + scals[i] = 0; + } + } + + /* + printf("UNPACKING LAPLACE=%.20f\n",laplacianOperator); + printf("packed offset=%ld\n",packed_offset); + for(i=0;i 0) { + lup = mmax; + if (sub_k >= 0) { + for (hcount = 0; hcount < sub_k + 1; hcount++) { + val[i++] = decode_float(grib_decode_unsigned_long(hres, &hpos, 8 * bytes)); + val[i++] = decode_float(grib_decode_unsigned_long(hres, &hpos, 8 * bytes)); + + if (GRIBEX_sh_bug_present && hcount == sub_k) { + /* bug in ecmwf data, last row (K+1)is scaled but should not */ + val[i - 2] *= scals[lup]; + val[i - 1] *= scals[lup]; + } + lup++; + } + sub_k--; + } + + pscals = scals + lup; + pval = val + i; +#if FAST_BIG_ENDIAN + grib_decode_double_array_complex(lres, + &lpos, bits_per_value, + reference_value, s, pscals, (maxv - hcount) * 2, pval); + i += (maxv - hcount) * 2; +#else + (void)pscals; /* suppress gcc warning */ + (void)pval; /* suppress gcc warning */ + for (lcount = hcount; lcount < maxv; lcount++) { + val[i++] = d * (T)((grib_decode_unsigned_long(lres, &lpos, bits_per_value) * s) + reference_value) * scals[lup]; + val[i++] = d * (T)((grib_decode_unsigned_long(lres, &lpos, bits_per_value) * s) + reference_value) * scals[lup]; + /* These values should always be zero, but as they are packed, + it is necessary to force them back to zero */ + if (mmax == 0) + val[i - 1] = 0; + lup++; + } +#endif + + maxv--; + hcount = 0; + mmax++; + } + + Assert(*len >= i); + *len = i; + + grib_context_free(a->context, scals); + + return ret; +} + +int grib_accessor_class_data_complex_packing_t::unpack_double(grib_accessor* a, double* val, size_t* len) +{ + return unpack_real(a, val, len); } -int grib_accessor_class_data_complex_packing_t::unpack_float(grib_accessor* a, float* val, size_t* len){ +int grib_accessor_class_data_complex_packing_t::unpack_float(grib_accessor* a, float* val, size_t* len) +{ // TODO(maee): See ECC-1579 // Investigate why results are not bit-identical // return unpack(a, val, len); - int err = 0; - size_t i = 0; - size_t size = *len; - double* val8 = NULL; - val8 = (double*)grib_context_malloc(a->context, size*(sizeof(double))); + size_t size = *len; + double* val8 = (double*)grib_context_malloc(a->context, size * (sizeof(double))); if (!val8) return GRIB_OUT_OF_MEMORY; - err = unpack(a, val8, len); + int err = unpack_real(a, val8, len); if (err) { - grib_context_free(a->context,val8); + grib_context_free(a->context, val8); return err; } - for(i=0; icontext,val8); + grib_context_free(a->context, val8); return GRIB_SUCCESS; } diff --git a/src/accessor/grib_accessor_class_data_complex_packing.h b/src/accessor/grib_accessor_class_data_complex_packing.h index d588eef38..a8276311a 100644 --- a/src/accessor/grib_accessor_class_data_complex_packing.h +++ b/src/accessor/grib_accessor_class_data_complex_packing.h @@ -9,8 +9,7 @@ * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. */ -#ifndef eccodes_accessor_data_complex_packing_h -#define eccodes_accessor_data_complex_packing_h +#pragma once #include "grib_accessor_class_data_simple_packing.h" #include "grib_scaling.h" @@ -45,228 +44,4 @@ class grib_accessor_class_data_complex_packing_t : public grib_accessor_class_da int unpack_float(grib_accessor*, float* val, size_t* len) override; int value_count(grib_accessor*, long*) override; void init(grib_accessor*, const long, grib_arguments*) override; - -private: - template int unpack(grib_accessor*, T*, size_t*); }; - -template -int grib_accessor_class_data_complex_packing_t::unpack(grib_accessor* a, T* val, size_t* len){ - static_assert(std::is_floating_point::value, "Requires floating point numbers"); - grib_accessor_data_complex_packing_t* self = (grib_accessor_data_complex_packing_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); - const char* cclass_name = a->cclass->name; - - size_t i = 0; - int ret = GRIB_SUCCESS; - long hcount = 0; - long lcount = 0; - long hpos = 0; - long lup = 0; - long mmax = 0; - long n_vals = 0; - T* scals = NULL; - T* pscals = NULL, *pval = NULL; - - T s = 0; - T d = 0; - T laplacianOperator = 0; - unsigned char* buf = NULL; - unsigned char* hres = NULL; - unsigned char* lres = NULL; - unsigned long packed_offset; - long lpos = 0; - - long maxv = 0; - long GRIBEX_sh_bug_present = 0; - long ieee_floats = 0; - - long offsetdata = 0; - long bits_per_value = 0; - T reference_value = 0; - long binary_scale_factor = 0; - long decimal_scale_factor = 0; - - long sub_j = 0; - long sub_k = 0; - long sub_m = 0; - long pen_j = 0; - long pen_k = 0; - long pen_m = 0; - - T operat = 0; - int bytes; - int err = 0; - double tmp; - - decode_float_proc decode_float = NULL; - - err = a->value_count(&n_vals); if (err) - return err; - - if (*len < n_vals) { - *len = n_vals; - return GRIB_ARRAY_TOO_SMALL; - } - - if ((ret = grib_get_long_internal(gh, self->offsetdata, &offsetdata)) != GRIB_SUCCESS) - return ret; - if ((ret = grib_get_long_internal(gh, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS) - return ret; - if ((ret = grib_get_double_internal(gh, self->reference_value, &tmp)) != GRIB_SUCCESS) - return ret; - reference_value = tmp; - if ((ret = grib_get_long_internal(gh, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS) - return ret; - - if ((ret = grib_get_long_internal(gh, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS) - return ret; - - if ((ret = grib_get_long_internal(gh, self->GRIBEX_sh_bug_present, &GRIBEX_sh_bug_present)) != GRIB_SUCCESS) - return ret; - - /* ECC-774: don't use grib_get_long_internal */ - if ((ret = grib_get_long(gh, self->ieee_floats, &ieee_floats)) != GRIB_SUCCESS) - return ret; - - if ((ret = grib_get_double_internal(gh, self->laplacianOperator, &tmp)) != GRIB_SUCCESS) - return ret; - laplacianOperator = tmp; - - if ((ret = grib_get_long_internal(gh, self->sub_j, &sub_j)) != GRIB_SUCCESS) - return ret; - if ((ret = grib_get_long_internal(gh, self->sub_k, &sub_k)) != GRIB_SUCCESS) - return ret; - if ((ret = grib_get_long_internal(gh, self->sub_m, &sub_m)) != GRIB_SUCCESS) - return ret; - if ((ret = grib_get_long_internal(gh, self->pen_j, &pen_j)) != GRIB_SUCCESS) - return ret; - if ((ret = grib_get_long_internal(gh, self->pen_k, &pen_k)) != GRIB_SUCCESS) - return ret; - if ((ret = grib_get_long_internal(gh, self->pen_m, &pen_m)) != GRIB_SUCCESS) - return ret; - - self->dirty = 0; - - switch (ieee_floats) { - case 0: - decode_float = grib_long_to_ibm; - bytes = 4; - break; - case 1: - decode_float = grib_long_to_ieee; - bytes = 4; - break; - case 2: - decode_float = grib_long_to_ieee64; - bytes = 8; - break; - default: - return GRIB_NOT_IMPLEMENTED; - } - - if (sub_j != sub_k || sub_j != sub_m || pen_j != pen_k || pen_j != pen_m) { - grib_context_log(a->context, GRIB_LOG_ERROR, "%s: Invalid pentagonal resolution parameters", cclass_name); - return GRIB_DECODING_ERROR; - } - - buf = (unsigned char*)gh->buffer->data; - - maxv = pen_j + 1; - - buf += a->byte_offset(); hres = buf; - lres = buf; - - if (pen_j == sub_j) { - n_vals = (pen_j + 1) * (pen_j + 2); - d = codes_power(-decimal_scale_factor, 10); - - grib_ieee_decode_array(a->context, buf, n_vals, bytes, val); - if (d) { - for (i = 0; i < n_vals; i++) - val[i] *= d; - } - return 0; - } - - packed_offset = a->byte_offset() + bytes * (sub_k + 1) * (sub_k + 2); - lpos = 8 * (packed_offset - offsetdata); - - s = codes_power(binary_scale_factor, 2); - d = codes_power(-decimal_scale_factor, 10); - - scals = (T*)grib_context_malloc(a->context, maxv * sizeof(T)); - if (!scals) return GRIB_OUT_OF_MEMORY; - - scals[0] = 0; - for (i = 1; i < maxv; i++) { - operat = pow(i * (i + 1), laplacianOperator); - if (operat != 0) - scals[i] = (1.0 / operat); - else { - grib_context_log(a->context, GRIB_LOG_WARNING, - "%s: Problem with operator div by zero at index %d of %d", cclass_name, i, maxv); - scals[i] = 0; - } - } - - /* - printf("UNPACKING LAPLACE=%.20f\n",laplacianOperator); - printf("packed offset=%ld\n",packed_offset); - for(i=0;i 0) { - lup = mmax; - if (sub_k >= 0) { - for (hcount = 0; hcount < sub_k + 1; hcount++) { - val[i++] = decode_float(grib_decode_unsigned_long(hres, &hpos, 8 * bytes)); - val[i++] = decode_float(grib_decode_unsigned_long(hres, &hpos, 8 * bytes)); - - if (GRIBEX_sh_bug_present && hcount == sub_k) { - /* bug in ecmwf data, last row (K+1)is scaled but should not */ - val[i - 2] *= scals[lup]; - val[i - 1] *= scals[lup]; - } - lup++; - } - sub_k--; - } - - pscals = scals + lup; - pval = val + i; -#if FAST_BIG_ENDIAN - grib_decode_double_array_complex(lres, - &lpos, bits_per_value, - reference_value, s, pscals, (maxv - hcount) * 2, pval); - i += (maxv - hcount) * 2; -#else - (void)pscals; /* suppress gcc warning */ - (void)pval; /* suppress gcc warning */ - for (lcount = hcount; lcount < maxv; lcount++) { - val[i++] = d * (T)((grib_decode_unsigned_long(lres, &lpos, bits_per_value) * s) + reference_value) * scals[lup]; - val[i++] = d * (T)((grib_decode_unsigned_long(lres, &lpos, bits_per_value) * s) + reference_value) * scals[lup]; - /* These values should always be zero, but as they are packed, - it is necessary to force them back to zero */ - if (mmax == 0) - val[i - 1] = 0; - lup++; - } -#endif - - maxv--; - hcount = 0; - mmax++; - } - - Assert(*len >= i); - *len = i; - - grib_context_free(a->context, scals); - - return ret; -} - -#endif /* eccodes_accessor_data_complex_packing_h */ diff --git a/src/accessor/grib_accessor_class_data_simple_packing.cc b/src/accessor/grib_accessor_class_data_simple_packing.cc index d1032808e..591e89526 100644 --- a/src/accessor/grib_accessor_class_data_simple_packing.cc +++ b/src/accessor/grib_accessor_class_data_simple_packing.cc @@ -11,26 +11,28 @@ #include "grib_accessor_class_data_simple_packing.h" #include "grib_optimize_decimal_factor.h" +#include "grib_bits_any_endian_simple.h" #include #include -grib_accessor_class_data_simple_packing_t _grib_accessor_class_data_simple_packing{"data_simple_packing"}; +grib_accessor_class_data_simple_packing_t _grib_accessor_class_data_simple_packing{ "data_simple_packing" }; grib_accessor_class* grib_accessor_class_data_simple_packing = &_grib_accessor_class_data_simple_packing; -void grib_accessor_class_data_simple_packing_t::init(grib_accessor* a, const long v, grib_arguments* args){ +void grib_accessor_class_data_simple_packing_t::init(grib_accessor* a, const long v, grib_arguments* args) +{ grib_accessor_class_values_t::init(a, v, args); grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); - self->units_factor = grib_arguments_get_name(gh, args, self->carg++); - self->units_bias = grib_arguments_get_name(gh, args, self->carg++); - self->changing_precision = grib_arguments_get_name(gh, args, self->carg++); - self->number_of_values = grib_arguments_get_name(gh, args, self->carg++); - self->bits_per_value = grib_arguments_get_name(gh, args, self->carg++); - self->reference_value = grib_arguments_get_name(gh, args, self->carg++); - self->binary_scale_factor = grib_arguments_get_name(gh, args, self->carg++); - self->decimal_scale_factor = grib_arguments_get_name(gh, args, self->carg++); - self->optimize_scaling_factor = grib_arguments_get_name(gh, args, self->carg++); + grib_handle* gh = grib_handle_of_accessor(a); + self->units_factor = grib_arguments_get_name(gh, args, self->carg++); + self->units_bias = grib_arguments_get_name(gh, args, self->carg++); + self->changing_precision = grib_arguments_get_name(gh, args, self->carg++); + self->number_of_values = grib_arguments_get_name(gh, args, self->carg++); + self->bits_per_value = grib_arguments_get_name(gh, args, self->carg++); + self->reference_value = grib_arguments_get_name(gh, args, self->carg++); + self->binary_scale_factor = grib_arguments_get_name(gh, args, self->carg++); + self->decimal_scale_factor = grib_arguments_get_name(gh, args, self->carg++); + self->optimize_scaling_factor = grib_arguments_get_name(gh, args, self->carg++); a->flags |= GRIB_ACCESSOR_FLAG_DATA; self->dirty = 1; } @@ -44,7 +46,8 @@ static const unsigned long nbits[32] = { 0x40000000, 0x80000000 }; -static int number_of_bits(unsigned long x, long* result){ +static int number_of_bits(unsigned long x, long* result) +{ const int count = sizeof(nbits) / sizeof(nbits[0]); const unsigned long* n = nbits; *result = 0; @@ -58,19 +61,21 @@ static int number_of_bits(unsigned long x, long* result){ return GRIB_SUCCESS; } -int grib_accessor_class_data_simple_packing_t::value_count(grib_accessor* a, long* number_of_values){ +int grib_accessor_class_data_simple_packing_t::value_count(grib_accessor* a, long* number_of_values) +{ grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; - *number_of_values = 0; + *number_of_values = 0; return grib_get_long_internal(grib_handle_of_accessor(a), self->number_of_values, number_of_values); } -int grib_accessor_class_data_simple_packing_t::unpack_double_element(grib_accessor* a, size_t idx, double* val){ +int grib_accessor_class_data_simple_packing_t::unpack_double_element(grib_accessor* a, size_t idx, double* val) +{ grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; - const char* cclass_name = a->cclass->name; + const char* cclass_name = a->cclass->name; - long n_vals; - int err = 0; + long n_vals = 0; + int err = 0; grib_handle* gh = grib_handle_of_accessor(a); double reference_value; @@ -78,12 +83,12 @@ int grib_accessor_class_data_simple_packing_t::unpack_double_element(grib_access long bits_per_value; long decimal_scale_factor; unsigned char* buf = (unsigned char*)gh->buffer->data; - double s = 0; - double d = 0; - long pos = 0; + double s = 0; + double d = 0; + long pos = 0; - n_vals = 0; - err = a->value_count(&n_vals); if (err) + err = a->value_count(&n_vals); + if (err) return err; if ((err = grib_get_long_internal(gh, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS) @@ -146,29 +151,176 @@ int grib_accessor_class_data_simple_packing_t::unpack_double_element(grib_access return err; } -int grib_accessor_class_data_simple_packing_t::unpack_double_element_set(grib_accessor* a, const size_t* index_array, size_t len, double* val_array){ - int err = 0; +int grib_accessor_class_data_simple_packing_t::unpack_double_element_set(grib_accessor* a, const size_t* index_array, size_t len, double* val_array) +{ + int err = 0; size_t i = 0; - for (i=0; i +static int unpack(grib_accessor* a, T* val, size_t* len) +{ + static_assert(std::is_floating_point::value, "Requires floating point numbers"); + + grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; + const char* cclass_name = a->cclass->name; + grib_handle* gh = grib_handle_of_accessor(a); + unsigned char* buf = (unsigned char*)grib_handle_of_accessor(a)->buffer->data; + + size_t i = 0; + int err = 0; + size_t n_vals = 0; + long pos = 0; + long count = 0; + + double reference_value; + long binary_scale_factor; + long bits_per_value; + long decimal_scale_factor; + long offsetBeforeData; + double s = 0; + double d = 0; + double units_factor = 1.0; + double units_bias = 0.0; + + err = a->value_count(&count); + if (err) + return err; + n_vals = count; + + if (*len < n_vals) { + *len = (long)n_vals; + return GRIB_ARRAY_TOO_SMALL; + } + + if ((err = grib_get_long_internal(gh, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS) + return err; + + /* + * check we don't decode bpv > max(ulong) as it is + * not currently supported by the algorithm + */ + if (bits_per_value > (sizeof(long) * 8)) { + return GRIB_INVALID_BPV; + } + + if (self->units_factor && + (grib_get_double_internal(gh, self->units_factor, &units_factor) == GRIB_SUCCESS)) { + grib_set_double_internal(gh, self->units_factor, 1.0); + } + + if (self->units_bias && + (grib_get_double_internal(gh, self->units_bias, &units_bias) == GRIB_SUCCESS)) { + grib_set_double_internal(gh, self->units_bias, 0.0); + } + + if (n_vals == 0) { + *len = 0; + return GRIB_SUCCESS; + } + + self->dirty = 0; + + if ((err = grib_get_double_internal(gh, self->reference_value, &reference_value)) != GRIB_SUCCESS) + return err; + + if ((err = grib_get_long_internal(gh, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS) + return err; + + if ((err = grib_get_long_internal(gh, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS) + return err; -int grib_accessor_class_data_simple_packing_t::unpack_double(grib_accessor* a, double* val, size_t* len){ + /* Special case */ + + if (bits_per_value == 0) { + for (i = 0; i < n_vals; i++) + val[i] = reference_value; + *len = n_vals; + return GRIB_SUCCESS; + } + + s = codes_power(binary_scale_factor, 2); + d = codes_power(-decimal_scale_factor, 10); + + grib_context_log(a->context, GRIB_LOG_DEBUG, + "%s %s: Creating %s, %zu values", cclass_name, __func__, a->name, n_vals); + + offsetBeforeData = a->byte_offset(); + buf += offsetBeforeData; + + /*Assert(((bits_per_value*n_vals)/8) < (1<<29));*/ /* See GRIB-787 */ + + /* ECC-941 */ + if (!a->context->ieee_packing) { + /* Must turn off this check when the environment variable ECCODES_GRIB_IEEE_PACKING is on */ + long offsetAfterData = 0; + err = grib_get_long(gh, "offsetAfterData", &offsetAfterData); + if (!err && offsetAfterData > offsetBeforeData) { + const long valuesSize = (bits_per_value * n_vals) / 8; /*in bytes*/ + if (offsetBeforeData + valuesSize > offsetAfterData) { + grib_context_log(a->context, GRIB_LOG_ERROR, + "%s: Data section size mismatch: " + "offset before data=%ld, offset after data=%ld (num values=%zu, bits per value=%ld)", + cclass_name, offsetBeforeData, offsetAfterData, n_vals, bits_per_value); + return GRIB_DECODING_ERROR; + } + } + // if (offsetBeforeData == offsetAfterData) { + // /* Crazy case: Constant field with bitsPerValue > 0 */ + // for (i = 0; i < n_vals; i++) + // val[i] = reference_value; + // *len = n_vals; + // return GRIB_SUCCESS; + // } + } + + grib_context_log(a->context, GRIB_LOG_DEBUG, + "%s %s: calling outline function: bpv: %ld, rv: %g, bsf: %ld, dsf: %ld", + cclass_name, __func__, bits_per_value, reference_value, binary_scale_factor, decimal_scale_factor); + grib_decode_array(buf, &pos, bits_per_value, reference_value, s, d, n_vals, val); + + *len = (long)n_vals; + + if (units_factor != 1.0) { + if (units_bias != 0.0) { + for (i = 0; i < n_vals; i++) { + val[i] = val[i] * units_factor + units_bias; + } + } + else { + for (i = 0; i < n_vals; i++) { + val[i] *= units_factor; + } + } + } + else if (units_bias != 0.0) { + for (i = 0; i < n_vals; i++) { + val[i] += units_bias; + } + } + return err; +} + +int grib_accessor_class_data_simple_packing_t::unpack_double(grib_accessor* a, double* val, size_t* len) +{ return unpack(a, val, len); } -int grib_accessor_class_data_simple_packing_t::unpack_float(grib_accessor* a, float* val, size_t* len){ +int grib_accessor_class_data_simple_packing_t::unpack_float(grib_accessor* a, float* val, size_t* len) +{ return unpack(a, val, len); } -int _unpack_double(grib_accessor* a, double* val, size_t* len, unsigned char* buf, long pos, size_t n_vals){ +static int _unpack_double(grib_accessor* a, double* val, size_t* len, unsigned char* buf, long pos, size_t n_vals) +{ grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); - const char* cclass_name = a->cclass->name; + grib_handle* gh = grib_handle_of_accessor(a); + const char* cclass_name = a->cclass->name; size_t i = 0; int err = 0; @@ -240,7 +392,8 @@ int _unpack_double(grib_accessor* a, double* val, size_t* len, unsigned char* bu grib_context_log(a->context, GRIB_LOG_DEBUG, "%s %s: Creating %s, %zu values", cclass_name, __func__, a->name, n_vals); - offsetBeforeData = a->byte_offset(); buf += offsetBeforeData; + offsetBeforeData = a->byte_offset(); + buf += offsetBeforeData; /*Assert(((bits_per_value*n_vals)/8) < (1<<29));*/ /* See GRIB-787 */ @@ -290,14 +443,16 @@ int _unpack_double(grib_accessor* a, double* val, size_t* len, unsigned char* bu return err; } -int grib_accessor_class_data_simple_packing_t::unpack_double_subarray(grib_accessor* a, double* val, size_t start, size_t len){ +int grib_accessor_class_data_simple_packing_t::unpack_double_subarray(grib_accessor* a, double* val, size_t start, size_t len) +{ grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; - unsigned char* buf = (unsigned char*)grib_handle_of_accessor(a)->buffer->data; - size_t nvals = len; - size_t* plen = &len; - long bits_per_value = 0; - long pos; - int err; + + unsigned char* buf = (unsigned char*)grib_handle_of_accessor(a)->buffer->data; + size_t nvals = len; + size_t* plen = &len; + long bits_per_value = 0; + long pos = 0; + int err = GRIB_SUCCESS; if ((err = grib_get_long_internal(grib_handle_of_accessor(a), self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS) @@ -308,14 +463,15 @@ int grib_accessor_class_data_simple_packing_t::unpack_double_subarray(grib_acces return _unpack_double(a, val, plen, buf, pos, nvals); } -int grib_accessor_class_data_simple_packing_t::pack_double(grib_accessor* a, const double* val, size_t* len){ +int grib_accessor_class_data_simple_packing_t::pack_double(grib_accessor* a, const double* val, size_t* len) +{ grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; - grib_handle* gh = grib_handle_of_accessor(a); - const char* cclass_name = a->cclass->name; + grib_handle* gh = grib_handle_of_accessor(a); + const char* cclass_name = a->cclass->name; - size_t i = 0; - size_t n_vals = *len; - int err = 0; + size_t i = 0; + size_t n_vals = *len; + int err = 0; double reference_value = 0; long binary_scale_factor = 0; long bits_per_value = 0; @@ -360,8 +516,10 @@ int grib_accessor_class_data_simple_packing_t::pack_double(grib_accessor* a, con max = val[0]; min = max; for (i = 1; i < n_vals; i++) { - if (val[i] > max) max = val[i]; - else if (val[i] < min) min = val[i]; + if (val[i] > max) + max = val[i]; + else if (val[i] < min) + min = val[i]; } if ((err = grib_check_data_values_minmax(gh, min, max)) != GRIB_SUCCESS) { @@ -460,7 +618,7 @@ int grib_accessor_class_data_simple_packing_t::pack_double(grib_accessor* a, con /* divisor=1; */ } else { - int last = 127; /* 'last' should be a parameter coming from a definitions file */ + int last = 127; /* 'last' should be a parameter coming from a definitions file */ if (c->gribex_mode_on && self->edition == 1) last = 99; /* bits_per_value is given and decimal_scale_factor and binary_scale_factor are calcualated */ diff --git a/src/accessor/grib_accessor_class_data_simple_packing.h b/src/accessor/grib_accessor_class_data_simple_packing.h index 92f355b3f..52761202d 100644 --- a/src/accessor/grib_accessor_class_data_simple_packing.h +++ b/src/accessor/grib_accessor_class_data_simple_packing.h @@ -9,11 +9,9 @@ * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction. */ -#ifndef eccodes_accessor_data_simple_packing_h -#define eccodes_accessor_data_simple_packing_h +#pragma once #include "grib_accessor_class_values.h" -#include "grib_bits_any_endian_simple.h" #include "grib_scaling.h" class grib_accessor_data_simple_packing_t : public grib_accessor_values_t @@ -45,151 +43,4 @@ class grib_accessor_class_data_simple_packing_t : public grib_accessor_class_val int unpack_double_element(grib_accessor*, size_t i, double* val) override; int unpack_double_element_set(grib_accessor*, const size_t* index_array, size_t len, double* val_array) override; int unpack_double_subarray(grib_accessor*, double* val, size_t start, size_t len) override; - -private: - template int unpack(grib_accessor*, T*, size_t*); }; - - -template -int grib_accessor_class_data_simple_packing_t::unpack(grib_accessor* a, T* val, size_t* len){ - static_assert(std::is_floating_point::value, "Requires floating point numbers"); - - grib_accessor_data_simple_packing_t* self = (grib_accessor_data_simple_packing_t*)a; - const char* cclass_name = a->cclass->name; - grib_handle* gh = grib_handle_of_accessor(a); - unsigned char* buf = (unsigned char*)grib_handle_of_accessor(a)->buffer->data; - - size_t i = 0; - int err = 0; - size_t n_vals = 0; - long pos = 0; - long count = 0; - - double reference_value; - long binary_scale_factor; - long bits_per_value; - long decimal_scale_factor; - long offsetBeforeData; - double s = 0; - double d = 0; - double units_factor = 1.0; - double units_bias = 0.0; - - err = a->value_count(&count); if (err) - return err; - n_vals = count; - - if (*len < n_vals) { - *len = (long)n_vals; - return GRIB_ARRAY_TOO_SMALL; - } - - if ((err = grib_get_long_internal(gh, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS) - return err; - - /* - * check we don't decode bpv > max(ulong) as it is - * not currently supported by the algorithm - */ - if (bits_per_value > (sizeof(long) * 8)) { - return GRIB_INVALID_BPV; - } - - if (self->units_factor && - (grib_get_double_internal(gh, self->units_factor, &units_factor) == GRIB_SUCCESS)) { - grib_set_double_internal(gh, self->units_factor, 1.0); - } - - if (self->units_bias && - (grib_get_double_internal(gh, self->units_bias, &units_bias) == GRIB_SUCCESS)) { - grib_set_double_internal(gh, self->units_bias, 0.0); - } - - if (n_vals == 0) { - *len = 0; - return GRIB_SUCCESS; - } - - self->dirty = 0; - - if ((err = grib_get_double_internal(gh, self->reference_value, &reference_value)) != GRIB_SUCCESS) - return err; - - if ((err = grib_get_long_internal(gh, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS) - return err; - - if ((err = grib_get_long_internal(gh, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS) - return err; - - /* Special case */ - - if (bits_per_value == 0) { - for (i = 0; i < n_vals; i++) - val[i] = reference_value; - *len = n_vals; - return GRIB_SUCCESS; - } - - s = codes_power(binary_scale_factor, 2); - d = codes_power(-decimal_scale_factor, 10); - - grib_context_log(a->context, GRIB_LOG_DEBUG, - "%s %s: Creating %s, %zu values", cclass_name, __func__, a->name, n_vals); - - offsetBeforeData = a->byte_offset(); buf += offsetBeforeData; - - /*Assert(((bits_per_value*n_vals)/8) < (1<<29));*/ /* See GRIB-787 */ - - /* ECC-941 */ - if (!a->context->ieee_packing) { - /* Must turn off this check when the environment variable ECCODES_GRIB_IEEE_PACKING is on */ - long offsetAfterData = 0; - err = grib_get_long(gh, "offsetAfterData", &offsetAfterData); - if (!err && offsetAfterData > offsetBeforeData) { - const long valuesSize = (bits_per_value * n_vals) / 8; /*in bytes*/ - if (offsetBeforeData + valuesSize > offsetAfterData) { - grib_context_log(a->context, GRIB_LOG_ERROR, - "%s: Data section size mismatch: " - "offset before data=%ld, offset after data=%ld (num values=%zu, bits per value=%ld)", - cclass_name, offsetBeforeData, offsetAfterData, n_vals, bits_per_value); - return GRIB_DECODING_ERROR; - } - } - // if (offsetBeforeData == offsetAfterData) { - // /* Crazy case: Constant field with bitsPerValue > 0 */ - // for (i = 0; i < n_vals; i++) - // val[i] = reference_value; - // *len = n_vals; - // return GRIB_SUCCESS; - // } - } - - grib_context_log(a->context, GRIB_LOG_DEBUG, - "%s %s: calling outline function: bpv: %ld, rv: %g, bsf: %ld, dsf: %ld", - cclass_name, __func__, bits_per_value, reference_value, binary_scale_factor, decimal_scale_factor); - grib_decode_array(buf, &pos, bits_per_value, reference_value, s, d, n_vals, val); - - *len = (long)n_vals; - - if (units_factor != 1.0) { - if (units_bias != 0.0) { - for (i = 0; i < n_vals; i++) { - val[i] = val[i] * units_factor + units_bias; - } - } else { - for (i = 0; i < n_vals; i++) { - val[i] *= units_factor; - } - } - } - else if (units_bias != 0.0) { - for (i = 0; i < n_vals; i++) { - val[i] += units_bias; - } - } - return err; -} - - -#endif /* eccodes_accessor_data_simple_packing_h */