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

ECC-1467: Debugging

This commit is contained in:
Shahram Najm 2022-12-29 17:30:01 +00:00
parent 78af61b989
commit 4dd790f2ff
5 changed files with 285 additions and 107 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
examples/C/grib_get_keys.c
View File

@ -1,19 +1,3 @@
/*
* (C) Copyright 2005- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
/*
* C Implementation: grib_get_keys
*
* Description: how to get values using keys from GRIB messages
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
@ -22,26 +6,14 @@
int main(int argc, char** argv)
{
int err = 0;
double* values = NULL;
float* values = NULL;
size_t values_len = 0;
size_t i = 0, len = 0;
double latitudeOfFirstGridPointInDegrees;
double longitudeOfFirstGridPointInDegrees;
double latitudeOfLastGridPointInDegrees;
double longitudeOfLastGridPointInDegrees;
double jDirectionIncrementInDegrees;
double iDirectionIncrementInDegrees;
long numberOfPointsAlongAParallel;
long numberOfPointsAlongAMeridian;
double average = 0;
char* packingType = NULL;
FILE* in = NULL;
const char* filename = "../../data/regular_latlon_surface.grib1";
const char* filename = "../../data/sample.grib2";
codes_handle* h = NULL;
in = fopen(filename, "rb");
@ -58,89 +30,25 @@ int main(int argc, char** argv)
}
fclose(in);
/* store the filename in the key "file" for this handle */
len = strlen(filename);
CODES_CHECK(codes_set_string(h, "file", filename, &len), 0);
/* get as a long*/
CODES_CHECK(codes_get_long(h, "Ni", &numberOfPointsAlongAParallel), 0);
printf("numberOfPointsAlongAParallel=%ld\n", numberOfPointsAlongAParallel);
/* get as a long*/
CODES_CHECK(codes_get_long(h, "Nj", &numberOfPointsAlongAMeridian), 0);
printf("numberOfPointsAlongAMeridian=%ld\n", numberOfPointsAlongAMeridian);
/* get as a double*/
CODES_CHECK(codes_get_double(h, "latitudeOfFirstGridPointInDegrees", &latitudeOfFirstGridPointInDegrees), 0);
printf("latitudeOfFirstGridPointInDegrees=%g\n", latitudeOfFirstGridPointInDegrees);
/* get as a double*/
CODES_CHECK(codes_get_double(h, "longitudeOfFirstGridPointInDegrees", &longitudeOfFirstGridPointInDegrees), 0);
printf("longitudeOfFirstGridPointInDegrees=%g\n", longitudeOfFirstGridPointInDegrees);
/* get as a double*/
CODES_CHECK(codes_get_double(h, "latitudeOfLastGridPointInDegrees", &latitudeOfLastGridPointInDegrees), 0);
printf("latitudeOfLastGridPointInDegrees=%g\n", latitudeOfLastGridPointInDegrees);
/* get as a double*/
CODES_CHECK(codes_get_double(h, "longitudeOfLastGridPointInDegrees", &longitudeOfLastGridPointInDegrees), 0);
printf("longitudeOfLastGridPointInDegrees=%g\n", longitudeOfLastGridPointInDegrees);
/* get as a double*/
CODES_CHECK(codes_get_double(h, "jDirectionIncrementInDegrees", &jDirectionIncrementInDegrees), 0);
printf("jDirectionIncrementInDegrees=%g\n", jDirectionIncrementInDegrees);
/* get as a double*/
CODES_CHECK(codes_get_double(h, "iDirectionIncrementInDegrees", &iDirectionIncrementInDegrees), 0);
printf("iDirectionIncrementInDegrees=%g\n", iDirectionIncrementInDegrees);
/* get as string */
CODES_CHECK(codes_get_length(h, "packingType", &len), 0);
packingType = (char*)malloc(len * sizeof(char));
codes_get_string(h, "packingType", packingType, &len);
printf("packingType=%s\n", packingType);
free(packingType);
/* get the size of the values array*/
CODES_CHECK(codes_get_size(h, "values", &values_len), 0);
values = (double*)malloc(values_len * sizeof(double));
values = (float*)malloc(values_len * sizeof(float));
/* get data values*/
CODES_CHECK(codes_get_double_array(h, "values", values, &values_len), 0);
CODES_CHECK(codes_get_float_array(h, "values", values, &values_len), 0);
average = 0;
for (i = 0; i < values_len; i++)
for (i = 0; i < values_len; i++){
//printf("%f\n",values[i]);
average += values[i];
}
average /= (double)values_len;
free(values);
printf("There are %d values, average is %g\n", (int)values_len, average);
{
int eq = 0;
/* now retrieve the value of the key "file" */
char file[256] = {0,};
CODES_CHECK(codes_get_length(h, "file", &len), 0);
assert(len == 1 + strlen(filename));
codes_get_string(h, "file", file, &len);
eq = strcmp(file, filename);
if (eq != 0) assert(!"file and filename not equal");
}
{
/* example of getting bytes */
const char* name = "reservedNeedNotBePresent";
unsigned char* byte_val = NULL;
size_t keySize = 0;
CODES_CHECK(codes_get_size(h, name, &keySize), 0);
byte_val = (unsigned char*)malloc(keySize * sizeof(unsigned char));
GRIB_CHECK(codes_get_bytes(h, name, byte_val, &keySize), name);
free(byte_val);
}
codes_handle_delete(h);
return 0;

1
src/eccodes_prototypes.h
View File

@ -1540,6 +1540,7 @@ int grib_encode_size_tb(unsigned char* p, size_t val, long* bitp, long nb);
/* grib_bits_any_endian_simple.c */
int grib_decode_long_array(const unsigned char* p, long* bitp, long bitsPerValue, size_t n_vals, long* val);
int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerValue, double reference_value, double s, double d, size_t n_vals, double* val);
int grib_decode_float_array(const unsigned char* p, long* bitp, long bitsPerValue, double reference_value, double s, double d, size_t n_vals, float* val);
int grib_decode_double_array_complex(const unsigned char* p, long* bitp, long nbits, double reference_value, double s, double* d, size_t size, double* val);
int grib_encode_long_array(size_t n_vals, const long* val, long bits_per_value, unsigned char* p, long* off);
int grib_encode_double_array(size_t n_vals, const double* val, long bits_per_value, double reference_value, double d, double divisor, unsigned char* p, long* off);

178
src/grib_accessor_class_data_simple_packing.c
View File

@ -23,7 +23,7 @@
CLASS = accessor
SUPER = grib_accessor_class_values
IMPLEMENTS = init
IMPLEMENTS = unpack_double
IMPLEMENTS = unpack_double;unpack_float
IMPLEMENTS = unpack_double_element;unpack_double_element_set
IMPLEMENTS = unpack_double_subarray
IMPLEMENTS = pack_double
@ -54,6 +54,7 @@ or edit "accessor.class" and rerun ./make_class.pl
static int pack_double(grib_accessor*, const double* val, size_t* len);
static int unpack_double(grib_accessor*, double* val, size_t* len);
static int unpack_float(grib_accessor*, float* val, size_t* len);
static int value_count(grib_accessor*, long*);
static void init(grib_accessor*, const long, grib_arguments*);
static void init_class(grib_accessor_class*);
@ -110,7 +111,7 @@ static grib_accessor_class _grib_accessor_class_data_simple_packing = {
&pack_double, /* grib_pack procedures double */
0, /* grib_pack procedures float */
&unpack_double, /* grib_unpack procedures double */
0, /* grib_unpack procedures float */
&unpack_float, /* grib_unpack procedures float */
0, /* grib_pack procedures string */
0, /* grib_unpack procedures string */
0, /* grib_pack array procedures string */
@ -150,7 +151,6 @@ static void init_class(grib_accessor_class* c)
c->pack_long = (*(c->super))->pack_long;
c->unpack_long = (*(c->super))->unpack_long;
c->pack_float = (*(c->super))->pack_float;
c->unpack_float = (*(c->super))->unpack_float;
c->pack_string = (*(c->super))->pack_string;
c->unpack_string = (*(c->super))->unpack_string;
c->pack_string_array = (*(c->super))->pack_string_array;
@ -317,6 +317,158 @@ static int unpack_double_element_set(grib_accessor* a, const size_t* index_array
return GRIB_SUCCESS;
}
// unpack an array of double or single precision real numbers. As doubles if dval!= NULL or floats if fval!=NULL
static int _unpack_real(grib_accessor* a, double* dval, float* fval, size_t* len, unsigned char* buf, long pos, size_t n_vals)
{
grib_accessor_data_simple_packing* self = (grib_accessor_data_simple_packing*)a;
grib_handle* gh = grib_handle_of_accessor(a);
size_t i = 0;
int err = 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;
Assert( ! (fval && dval) );
printf("DBG:: SIMPLE PACK unpack real.....\n");
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) {
if (dval) {
for (i = 0; i < n_vals; i++)
dval[i] = reference_value;
} else if (fval) {
for (i = 0; i < n_vals; i++)
fval[i] = reference_value;
}
*len = n_vals;
return GRIB_SUCCESS;
}
s = grib_power(binary_scale_factor, 2);
d = grib_power(-decimal_scale_factor, 10);
grib_context_log(a->context, GRIB_LOG_DEBUG,
"grib_accessor_data_simple_packing: unpack_double : creating %s, %d values",
a->name, n_vals);
offsetBeforeData = grib_byte_offset(a);
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,
"Data section size mismatch: offset before data=%ld, offset after data=%ld (num values=%ld, bits per value=%ld)",
offsetBeforeData, offsetAfterData, n_vals, bits_per_value);
return GRIB_DECODING_ERROR;
}
}
#if 0
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;
}
#endif
}
grib_context_log(a->context, GRIB_LOG_DEBUG,
"unpack_double: calling outline function : bpv %d, rv : %g, sf : %d, dsf : %d ",
bits_per_value, reference_value, binary_scale_factor, decimal_scale_factor);
if(dval)
grib_decode_double_array(buf, &pos, bits_per_value, reference_value, s, d, n_vals, dval);
if(fval)
grib_decode_float_array(buf, &pos, bits_per_value, reference_value, s, d, n_vals, fval);
*len = (long)n_vals;
if(dval) {
if (units_factor != 1.0) {
if (units_bias != 0.0)
for (i = 0; i < n_vals; i++)
dval[i] = dval[i] * units_factor + units_bias;
else
for (i = 0; i < n_vals; i++)
dval[i] *= units_factor;
}
else if (units_bias != 0.0)
for (i = 0; i < n_vals; i++)
dval[i] += units_bias;
}
if(fval) {
if (units_factor != 1.0) {
if (units_bias != 0.0)
for (i = 0; i < n_vals; i++)
fval[i] = fval[i] * units_factor + units_bias;
else
for (i = 0; i < n_vals; i++)
fval[i] *= units_factor;
}
else if (units_bias != 0.0)
for (i = 0; i < n_vals; i++)
fval[i] += units_bias;
}
return err;
}
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* self = (grib_accessor_data_simple_packing*)a;
@ -464,7 +616,7 @@ static int unpack_double_subarray(grib_accessor* a, double* val, size_t start, s
return _unpack_double(a, val, plen, buf, pos, nvals);
}
static int unpack_double(grib_accessor* a, double* val, size_t* len)
static int unpack_double(grib_accessor* a, double* dval, size_t* len)
{
unsigned char* buf = (unsigned char*)grib_handle_of_accessor(a)->buffer->data;
size_t nvals = 0;
@ -476,8 +628,24 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
if (err)
return err;
nvals = count;
printf("SIMPLE unpak double.....\n");
return _unpack_real(a, dval, NULL, len, buf, pos, nvals);
}
return _unpack_double(a, val, len, buf, pos, nvals);
static int unpack_float(grib_accessor* a, float* fval, size_t* len)
{
unsigned char* buf = (unsigned char*)grib_handle_of_accessor(a)->buffer->data;
size_t nvals = 0;
long pos = 0;
int err = 0;
long count = 0;
err = grib_value_count(a, &count);
if (err)
return err;
nvals = count;
printf("SIMPLE unpack float.....\n");
return _unpack_real(a, NULL, fval, len, buf, pos, nvals);
}
#if GRIB_IBMPOWER67_OPT

8
src/grib_accessor_class_gen.c
View File

@ -343,8 +343,10 @@ static int unpack_double(grib_accessor* a, double* v, size_t* len)
static int unpack_float(grib_accessor* a, float* v, size_t* len)
{
int type = GRIB_TYPE_UNDEFINED;
printf("====%s\n",a->name);
/* int type = GRIB_TYPE_UNDEFINED;
printf("grib_accessor_class_gen.c unpack_float a->name =%s\n",a->name);
printf("grib_accessor_class_gen.c unpack_float a->cclass->name=%s\n",a->cclass->name);
printf("grib_accessor_class_gen.c unpack_float alen=%zu\n",*len);
printf("DEBUG unpack_float:: a->cclass->unpack_double=%p\n", (void*)a->cclass->unpack_double);
printf("DEBUG unpack_float:: &unpack_double=%p\n", (void*)&unpack_double);
if (a->cclass->unpack_double && a->cclass->unpack_double != &unpack_double) {
@ -359,7 +361,7 @@ static int unpack_float(grib_accessor* a, float* v, size_t* len)
grib_context_log(a->context, GRIB_LOG_ERROR, "Cannot unpack %s as float", a->name);
if (grib_get_native_type(grib_handle_of_accessor(a), a->name, &type) == GRIB_SUCCESS) {
grib_context_log(a->context, GRIB_LOG_ERROR, "Hint: Try unpacking as %s", grib_get_type_name(type));
}
} */
return GRIB_NOT_IMPLEMENTED;
}

99
src/grib_bits_any_endian_simple.c
View File

@ -164,6 +164,105 @@ int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerVal
return 0;
}
/*
* TODO: First lame attempt at decoding a float array. This and the grib_decode_double_array function
* should be merged and refactored! Most probably using C++ templates
*/
/**
* decode an array of n_vals values from an octet-bitstream to float-representation
*
* @param p input bitstream, for technical reasons put into octets
* @param bitp current position in the bitstream
* @param bitsPerValue number of bits needed to build a number (e.g. 8=byte, 16=short, 32=int, but also other sizes allowed)
* @param n_vals number of values to decode
* @param val output, values encoded as 32/64bit numbers
*/
int grib_decode_float_array(const unsigned char* p, long* bitp, long bitsPerValue,
double reference_value, double s, double d,
size_t n_vals, float* val)
{
long i = 0;
unsigned long lvalue = 0;
double x;
printf("grib_decode_float_array\n");
#if 0
/* slow reference code */
int j=0;
for(i=0;i < n_vals;i++) {
lvalue=0;
for(j=0; j< bitsPerValue;j++){
lvalue <<= 1;
if(grib_get_bit( p, *bitp)) lvalue += 1;
*bitp += 1;
}
x=((lvalue*s)+reference_value)*d;
val[i] = (double)x;
}
#endif
if (bitsPerValue % 8 == 0) {
/* See ECC-386 */
int bc;
int l = bitsPerValue / 8;
size_t o = 0;
for (i = 0; i < n_vals; i++) {
lvalue = 0;
lvalue <<= 8;
lvalue |= p[o++];
for (bc = 1; bc < l; bc++) {
lvalue <<= 8;
lvalue |= p[o++];
}
x = ((lvalue * s) + reference_value) * d;
val[i] = (float)x;
/* *bitp += bitsPerValue * n_vals; */
}
}
else {
unsigned long mask = BIT_MASK1(bitsPerValue);
/* pi: position of bitp in p[]. >>3 == /8 */
long pi = *bitp / 8;
/* some bits might of the current byte at pi might be used */
/* by the previous number usefulBitsInByte gives remaining unused bits */
/* number of useful bits in current byte */
int usefulBitsInByte = 8 - (*bitp & 7);
for (i = 0; i < n_vals; i++) {
/* value read as long */
long bitsToRead = 0;
lvalue = 0;
bitsToRead = bitsPerValue;
/* read one byte after the other to lvalue until >= bitsPerValue are read */
while (bitsToRead > 0) {
lvalue <<= 8;
lvalue += p[pi];
pi++;
bitsToRead -= usefulBitsInByte;
usefulBitsInByte = 8;
}
*bitp += bitsPerValue;
/* bitsToRead is now <= 0, remove the last bits */
lvalue >>= -1 * bitsToRead;
/* set leading bits to 0 - removing bits used for previous number */
lvalue &= mask;
usefulBitsInByte = -1 * bitsToRead; /* prepare for next round */
if (usefulBitsInByte > 0) {
pi--; /* reread the current byte */
}
else {
usefulBitsInByte = 8; /* start with next full byte */
}
/* scaling and move value to output */
x = ((lvalue * s) + reference_value) * d;
val[i] = (float)x;
}
}
return 0;
}
int grib_decode_double_array_complex(const unsigned char* p, long* bitp, long nbits, double reference_value, double s, double* d, size_t size, double* val)
{
return GRIB_NOT_IMPLEMENTED;