core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Move function from header to local static. Use pragma once

This commit is contained in:
shahramn 2024-05-08 10:16:12 +01:00
parent 57b2698cc3
commit cf72dca971
6 changed files with 597 additions and 586 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

178
src/accessor/grib_accessor_class_data_apply_bitmap.cc
View File

@ -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 <typename T>
static int unpack(grib_accessor* a, T* val, size_t* len)
{
static_assert(std::is_floating_point<T>::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<T>(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<T>(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<T>(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<double>(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<float>(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<alen && retval == GRIB_SUCCESS; ++i) {
for (size_t i = 0; i < alen && retval == GRIB_SUCCESS; ++i) {
if (aval[i] != bval[i]) retval = GRIB_DOUBLE_VALUE_MISMATCH;
}
@ -281,4 +376,3 @@ int grib_accessor_class_data_apply_bitmap_t::compare(grib_accessor* a, grib_acce
return retval;
}

90
src/accessor/grib_accessor_class_data_apply_bitmap.h
View File

@ -9,8 +9,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#ifndef eccodes_accessor_data_apply_bitmap_h
#define eccodes_accessor_data_apply_bitmap_h
#pragma once
#include "grib_api_internal.h"
#include "grib_accessor_class_gen.h"
@ -42,91 +41,4 @@ public:
int compare(grib_accessor*, grib_accessor*) override;
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;
private:
template <typename T> int unpack(grib_accessor*, T*, size_t*);
};
template <typename T>
int grib_accessor_class_data_apply_bitmap_t::unpack(grib_accessor* a, T* val, size_t* len){
static_assert(std::is_floating_point<T>::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<T>(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<T>(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<T>(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 */

279
src/accessor/grib_accessor_class_data_complex_packing.cc
View File

@ -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 <typename T>
static int unpack_real(grib_accessor* a, T* val, size_t* len)
{
static_assert(std::is_floating_point<T>::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;
int grib_accessor_class_data_complex_packing_t::unpack_double(grib_accessor* a, double* val, size_t* len){
return unpack<double>(a, val, len);
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<T>(-decimal_scale_factor, 10);
grib_ieee_decode_array<T>(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<T>(binary_scale_factor, 2);
d = codes_power<T>(-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<maxv;i++)
printf("scals[%d]=%g\n",i,scals[i]);*/
i = 0;
while (maxv > 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_float(grib_accessor* a, float* val, size_t* len){
int grib_accessor_class_data_complex_packing_t::unpack_double(grib_accessor* a, double* val, size_t* len)
{
return unpack_real<double>(a, val, 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<float>(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<double>(a, val8, len);
int err = unpack_real<double>(a, val8, len);
if (err) {
grib_context_free(a->context,val8);
grib_context_free(a->context, val8);
return err;
}
for(i=0; i<size; i++)
for (size_t i = 0; i < size; i++)
val[i] = val8[i];
grib_context_free(a->context,val8);
grib_context_free(a->context, val8);
return GRIB_SUCCESS;
}

227
src/accessor/grib_accessor_class_data_complex_packing.h
View File

@ -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 @@ public:
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 <typename T> int unpack(grib_accessor*, T*, size_t*);
};
template <typename T>
int grib_accessor_class_data_complex_packing_t::unpack(grib_accessor* a, T* val, size_t* len){
static_assert(std::is_floating_point<T>::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<T>(-decimal_scale_factor, 10);
grib_ieee_decode_array<T>(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<T>(binary_scale_factor, 2);
d = codes_power<T>(-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<maxv;i++)
printf("scals[%d]=%g\n",i,scals[i]);*/
i = 0;
while (maxv > 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 */

258
src/accessor/grib_accessor_class_data_simple_packing.cc
View File

@ -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 <float.h>
#include <type_traits>
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<len; ++i) {
if ((err = unpack_double_element(a, index_array[i], val_array+i)) != GRIB_SUCCESS)
for (i = 0; i < len; ++i) {
if ((err = unpack_double_element(a, index_array[i], val_array + i)) != GRIB_SUCCESS)
return err;
}
return GRIB_SUCCESS;
}
template <typename T>
static int unpack(grib_accessor* a, T* val, size_t* len)
{
static_assert(std::is_floating_point<T>::value, "Requires floating point numbers");
int grib_accessor_class_data_simple_packing_t::unpack_double(grib_accessor* a, double* val, size_t* len){
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<T>(binary_scale_factor, 2);
d = codes_power<T>(-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<T>(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<double>(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<float>(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 */

151
src/accessor/grib_accessor_class_data_simple_packing.h
View File

@ -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 @@ public:
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 <typename T> int unpack(grib_accessor*, T*, size_t*);
};
template <typename T>
int grib_accessor_class_data_simple_packing_t::unpack(grib_accessor* a, T* val, size_t* len){
static_assert(std::is_floating_point<T>::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<T>(binary_scale_factor, 2);
d = codes_power<T>(-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<T>(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 */