1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
|
/*
* (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_precision
*
* Description: how to control decimal precision when packing fields.
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include "eccodes.h"
#define EPSILON 1e-6
#define EXPECTED_MIN 270.466796875
#define EXPECTED_MAX 311.096796875
int main(int argc, char** argv)
{
int err = 0;
size_t size = 0;
FILE* in = NULL;
const char* infile = "../../data/regular_latlon_surface.grib1";
FILE* out = NULL;
const char* outfile = "out.precision.grib1";
codes_handle* h = NULL;
const void* buffer = NULL;
double* values1 = NULL;
double* values2 = NULL;
double maxa = 0;
double maxv = 0, minv = 0;
double maxr = 0, r = 0;
long decimalPrecision;
long bitsPerValue1 = 0, bitsPerValue2 = 0;
int i = 0;
in = fopen(infile, "rb");
if (!in) {
fprintf(stderr, "Error: unable to open input file %s\n", infile);
return 1;
}
out = fopen(outfile, "wb");
if (!out) {
fprintf(stderr, "Error: unable to open output file %s\n", outfile);
fclose(in);
return 1;
}
/* create a new handle from a message in a file */
h = codes_handle_new_from_file(0, in, PRODUCT_GRIB, &err);
if (h == NULL) {
fprintf(stderr, "Error: unable to create handle from file %s\n", infile);
fclose(out);
return 1;
}
/* bitsPerValue before changing the packing parameters */
CODES_CHECK(codes_get_long(h, "bitsPerValue", &bitsPerValue1), 0);
assert(bitsPerValue1 == 16);
/* get the size of the values array*/
CODES_CHECK(codes_get_size(h, "values", &size), 0);
assert(size == 496);
values1 = (double*)malloc(size * sizeof(double));
/* get data values before changing the packing parameters*/
CODES_CHECK(codes_get_double_array(h, "values", values1, &size), 0);
/* changing decimal precision to 2 means that 2 decimal digits
are preserved when packing. */
decimalPrecision = 2;
CODES_CHECK(codes_set_long(h, "changeDecimalPrecision", decimalPrecision), 0);
/* bitsPerValue after changing the packing parameters */
CODES_CHECK(codes_get_long(h, "bitsPerValue", &bitsPerValue2), 0);
assert(bitsPerValue2 == 12);
values2 = (double*)malloc(size * sizeof(double));
/* get data values after changing the packing parameters*/
CODES_CHECK(codes_get_double_array(h, "values", values2, &size), 0);
/* computing error */
maxa = 0;
maxr = 0;
maxv = values2[0];
minv = maxv;
for (i = 0; i < size; i++) {
double a = fabs(values2[i] - values1[i]);
if (values2[i] > maxv) maxv = values2[i];
if (values2[i] < minv) minv = values2[i];
if (values2[i] != 0) r = fabs((values2[i] - values1[i]) / values2[i]);
if (a > maxa) maxa = a;
if (r > maxr) maxr = r;
}
printf("max absolute error = %g\n", maxa);
printf("max relative error = %g\n", maxr);
printf("min value = %g\n", minv);
printf("max value = %g\n", maxv);
assert(fabs(minv - EXPECTED_MIN) < EPSILON);
assert(fabs(maxv - EXPECTED_MAX) < EPSILON);
/* get the coded message in a buffer */
CODES_CHECK(codes_get_message(h, &buffer, &size), 0);
/* write the buffer in a file*/
if (fwrite(buffer, 1, size, out) != size) {
perror(outfile);
exit(1);
}
codes_handle_delete(h);
fclose(in);
fclose(out);
free(values1);
free(values2);
return 0;
}
|
