Single-precision: Support for second order packing (GRIB1)

This commit is contained in:
Shahram Najm 2023-05-10 18:02:09 +01:00 committed by shahramn
parent 1a4eaeaa61
commit 1e68497461
4 changed files with 112 additions and 49 deletions

View File

@ -45,7 +45,8 @@
MEMBERS=const char* orderOfSPD
MEMBERS=const char* numberOfPoints
MEMBERS=const char* dataFlag
MEMBERS=double* values
MEMBERS=double* dvalues
MEMBERS=float* fvalues
MEMBERS=size_t size
@ -117,7 +118,10 @@ typedef struct grib_accessor_data_g1second_order_general_extended_packing
const char* orderOfSPD;
const char* numberOfPoints;
const char* dataFlag;
double* values;
double* dvalues;
float* fvalues;
int double_dirty;
int float_dirty;
size_t size;
} grib_accessor_data_g1second_order_general_extended_packing;
@ -250,7 +254,9 @@ static void init(grib_accessor* a, const long v, grib_arguments* args)
self->dataFlag = grib_arguments_get_name(handle, args, self->carg++);
self->edition = 1;
self->dirty = 1;
self->values = NULL;
self->dvalues = NULL;
self->fvalues = NULL;
self->double_dirty = self->float_dirty = 1;
self->size = 0;
a->flags |= GRIB_ACCESSOR_FLAG_DATA;
}
@ -344,12 +350,7 @@ static int unpack_double_element_set(grib_accessor* a, const size_t* index_array
return GRIB_SUCCESS;
}
static int unpack_float(grib_accessor*, float* val, size_t* len)
{
return GRIB_NOT_IMPLEMENTED;
}
static int unpack_double(grib_accessor* a, double* values, size_t* len)
static int unpack(grib_accessor* a, double* dvalues, float* fvalues, size_t* len)
{
grib_accessor_data_g1second_order_general_extended_packing* self = (grib_accessor_data_g1second_order_general_extended_packing*)a;
int ret = 0;
@ -372,19 +373,33 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
long bias = 0;
long y = 0, z = 0, w = 0;
size_t k, ngroups;
Assert(!(dvalues && fvalues));
if (!self->dirty) {
if (dvalues) {
if (!self->double_dirty) {
if (*len < self->size) {
return GRIB_ARRAY_TOO_SMALL;
}
for (k = 0; k < self->size; k++)
values[k] = self->values[k];
dvalues[k] = self->dvalues[k];
*len = self->size;
return GRIB_SUCCESS;
}
self->double_dirty = 0;
}
self->dirty = 0;
if (fvalues) {
if (!self->float_dirty) {
if (*len < self->size) {
return GRIB_ARRAY_TOO_SMALL;
}
for (k = 0; k < self->size; k++)
fvalues[k] = self->fvalues[k];
*len = self->size;
return GRIB_SUCCESS;
}
self->float_dirty = 0;
}
buf += grib_byte_offset(a);
ret = value_count(a, &numberOfValues);
@ -506,21 +521,42 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
break;
}
if (self->values) {
if (dvalues) { // double-precision
if (self->dvalues) {
if (numberOfValues != self->size) {
grib_context_free(a->context, self->values);
self->values = (double*)grib_context_malloc_clear(a->context, sizeof(double) * numberOfValues);
grib_context_free(a->context, self->dvalues);
self->dvalues = (double*)grib_context_malloc_clear(a->context, sizeof(double) * numberOfValues);
}
}
else {
self->values = (double*)grib_context_malloc_clear(a->context, sizeof(double) * numberOfValues);
self->dvalues = (double*)grib_context_malloc_clear(a->context, sizeof(double) * numberOfValues);
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
for (i = 0; i < numberOfValues; i++) {
values[i] = (double)(((X[i] * s) + reference_value) * d);
self->values[i] = values[i];
dvalues[i] = (double)(((X[i] * s) + reference_value) * d);
self->dvalues[i] = dvalues[i];
}
}
else {
// single-precision
if (self->fvalues) {
if (numberOfValues != self->size) {
grib_context_free(a->context, self->fvalues);
self->fvalues = (float*)grib_context_malloc_clear(a->context, sizeof(float) * numberOfValues);
}
}
else {
self->fvalues = (float*)grib_context_malloc_clear(a->context, sizeof(float) * numberOfValues);
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
for (i = 0; i < numberOfValues; i++) {
fvalues[i] = (float)(((X[i] * s) + reference_value) * d);
self->fvalues[i] = fvalues[i];
}
}
*len = numberOfValues;
@ -536,6 +572,16 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
return ret;
}
static int unpack_float(grib_accessor* a, float* values, size_t* len)
{
return unpack(a, NULL, values, len);
}
static int unpack_double(grib_accessor* a, double* values, size_t* len)
{
return unpack(a, values, NULL, len);
}
static void grib_split_long_groups(grib_handle* hand, grib_context* c, long* numberOfGroups, long* lengthOfSecondOrderValues,
long* groupLengths, long* widthOfLengths,
long* groupWidths, long widthOfWidths,
@ -1318,7 +1364,7 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
grib_handle* handle = grib_handle_of_accessor(a);
long optimize_scaling_factor = 0;
self->dirty = 1;
self->double_dirty = 1;
numberOfValues = *len;
@ -1933,8 +1979,12 @@ static int pack_double(grib_accessor* a, const double* val, size_t* len)
static void destroy(grib_context* context, grib_accessor* a)
{
grib_accessor_data_g1second_order_general_extended_packing* self = (grib_accessor_data_g1second_order_general_extended_packing*)a;
if (self->values != NULL) {
grib_context_free(context, self->values);
self->values = NULL;
if (self->dvalues != NULL) {
grib_context_free(context, self->dvalues);
self->dvalues = NULL;
}
if (self->fvalues != NULL) {
grib_context_free(context, self->fvalues);
self->fvalues = NULL;
}
}

View File

@ -195,13 +195,10 @@ static int value_count(grib_accessor* a, long* numberOfSecondOrderPackedValues)
return err;
}
static int unpack_float(grib_accessor*, float* val, size_t* len)
{
return GRIB_NOT_IMPLEMENTED;
}
static int unpack_double(grib_accessor* a, double* values, size_t* len)
template <typename T>
static int unpack(grib_accessor* a, T* values, size_t* len)
{
static_assert(std::is_floating_point<T>::value, "Requires floating point numbers");
grib_accessor_data_g1second_order_general_packing* self = (grib_accessor_data_g1second_order_general_packing*)a;
int ret = 0;
long numberOfGroups, numberOfSecondOrderPackedValues;
@ -291,7 +288,7 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
for (i = 0; i < numberOfSecondOrderPackedValues; i++) {
values[i] = (double)(((X[i] * s) + reference_value) * d);
values[i] = (T)(((X[i] * s) + reference_value) * d);
}
*len = numberOfSecondOrderPackedValues;
@ -303,6 +300,16 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
return ret;
}
static int unpack_float(grib_accessor* a, float* values, size_t* len)
{
return unpack<float>(a, values, len);
}
static int unpack_double(grib_accessor* a, double* values, size_t* len)
{
return unpack<double>(a, values, len);
}
static int pack_double(grib_accessor* a, const double* cval, size_t* len)
{
/* return GRIB_NOT_IMPLEMENTED; */

View File

@ -259,12 +259,8 @@ static int value_count(grib_accessor* a, long* count)
return ret;
}
static int unpack_float(grib_accessor*, float* val, size_t* len)
{
return GRIB_NOT_IMPLEMENTED;
}
static int unpack_double(grib_accessor* a, double* values, size_t* len)
template <typename T>
static int unpack(grib_accessor* a, T* values, size_t* len)
{
grib_accessor_data_g1second_order_row_by_row_packing* self = (grib_accessor_data_g1second_order_row_by_row_packing*)a;
grib_handle* gh = grib_handle_of_accessor(a);
@ -418,7 +414,7 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
for (i = 0; i < n; i++) {
values[i] = (double)(((X[i] * s) + reference_value) * d);
values[i] = (T)(((X[i] * s) + reference_value) * d);
}
grib_context_free(a->context, firstOrderValues);
grib_context_free(a->context, X);
@ -431,6 +427,16 @@ static int unpack_double(grib_accessor* a, double* values, size_t* len)
return ret;
}
static int unpack_float(grib_accessor* a, float* values, size_t* len)
{
return unpack<float>(a, values, len);
}
static int unpack_double(grib_accessor* a, double* values, size_t* len)
{
return unpack<double>(a, values, len);
}
static int pack_double(grib_accessor* a, const double* cval, size_t* len)
{
int err = 0;

View File

@ -30,7 +30,7 @@ gfiles="$gfiles spherical_pressure_level.grib2" # spectral_complex, edition=2
gfiles="$gfiles gfs.complex.mvmu.grib2" # grid_complex, edition=2, g22order_packing
# Second order
gfiles="$gfiles lfpw.grib1" # grid_second_order, edition=1
gfiles="$gfiles lfpw.grib1 gen.grib gen_ext.grib second_ord_rbr.grib1" # grid_second_order, edition=1
#gfiles="$gfiles " # grid_second_order, edition=2
# CCSDS