1 /*
2 * (C) Copyright 2005- ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence Version 2.0
5 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
8 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
9 */
10
11 /*
12 *
13 * Description: Check the geometry of a global GRIB field
14 * which has a Gaussian Grid (reduced or regular)
15 *
16 */
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <assert.h>
20 #include <stdarg.h>
21 #include <sys/stat.h>
22
23 #include "grib_api_internal.h"
24
25 #define STR_EQUAL(s1, s2) (strcmp((s1), (s2)) == 0)
26
27 int exit_on_error = 1; /* By default exit if any check fails */
28 int verbose = 0; /* By default quiet unless errors */
29 int error_count = 0;
30
DBL_EQUAL(double d1,double d2,double tolerance)31 int DBL_EQUAL(double d1, double d2, double tolerance)
32 {
33 return fabs(d1 - d2) <= tolerance;
34 }
35
usage(const char * prog)36 static void usage(const char* prog)
37 {
38 printf("Usage: %s [-f] [-v] grib_file grib_file ...\n\n", prog);
39 printf("Check geometry of GRIB fields with a Gaussian Grid.\n");
40 printf("(The grid is assumed to be GLOBAL)\n\n");
41 printf("Options:\n");
42 printf(" -f Do not exit on first error\n");
43 printf(" -v Verbose\n");
44 printf("\n");
45 exit(1);
46 }
47
48 /* Print an error message and optionally die */
error(const char * filename,int msg_num,const char * fmt,...)49 void error(const char* filename, int msg_num, const char* fmt, ...)
50 {
51 char buf[1024] = {0,};
52 va_list list;
53 va_start(list, fmt);
54 if (verbose)
55 sprintf(buf, " Error: %s", fmt); /* indent a bit */
56 else
57 sprintf(buf, "Error: %s #%d: %s", filename, msg_num, fmt);
58 vfprintf(stderr, buf, list);
59 va_end(list);
60
61 ++error_count;
62 if (exit_on_error) {
63 exit(1);
64 }
65 }
66
get_precision(long edition)67 double get_precision(long edition)
68 {
69 if (edition == 1)
70 return 1.0 / 1000.0; /* milli degrees */
71 if (edition == 2)
72 return 1.0 / 1000000.0; /* micro degrees */
73 assert(!"Invalid edition");
74 return 0.0;
75 }
76
process_file(const char * filename)77 int process_file(const char* filename)
78 {
79 int err = 0, msg_num = 0;
80 grib_handle* h = NULL;
81 FILE* in = NULL;
82
83 if (!path_is_regular_file(filename)) {
84 if (verbose)
85 printf(" WARNING: '%s' not a regular file! Ignoring\n", filename);
86 return GRIB_IO_PROBLEM;
87 }
88
89 in = fopen(filename, "r");
90 if (!in) {
91 fprintf(stderr, "ERROR: unable to open input file '%s'\n", filename);
92 exit(1);
93 }
94
95 if (verbose)
96 printf("Checking file %s\n", filename);
97
98 while ((h = grib_handle_new_from_file(0, in, &err)) != NULL) {
99 int is_reduced = 0, is_regular = 0, grid_ok = 0;
100 long edition = 0, N = 0, Nj = 0, numberOfDataPoints;
101 size_t len = 0, sizeOfValuesArray = 0;
102 double* lats = NULL;
103 long* pl = NULL;
104 char gridType[128] = {0,};
105 double angular_tolerance, lat1, lon1, lat2, lon2, expected_lon2;
106 double iDirectionIncrementInDegrees;
107
108 if (err != GRIB_SUCCESS)
109 GRIB_CHECK(err, 0);
110 ++msg_num;
111 GRIB_CHECK(grib_get_long(h, "edition", &edition), 0);
112 if (verbose)
113 printf(" Processing GRIB message #%d (edition=%ld)\n", msg_num, edition);
114
115 len = 32;
116 GRIB_CHECK(grib_get_string(h, "gridType", gridType, &len), 0);
117 is_regular = STR_EQUAL(gridType, "regular_gg");
118 is_reduced = STR_EQUAL(gridType, "reduced_gg");
119 grid_ok = is_regular || is_reduced;
120 if (!grid_ok) {
121 /*error("ERROR: gridType should be Reduced or Regular Gaussian Grid!\n");*/
122 if (verbose)
123 printf(" WARNING: gridType=%s. It should be Reduced or Regular Gaussian Grid! Ignoring\n", gridType);
124 grib_handle_delete(h);
125 continue;
126 }
127
128 GRIB_CHECK(grib_get_long(h, "N", &N), 0);
129 GRIB_CHECK(grib_get_long(h, "Nj", &Nj), 0);
130 GRIB_CHECK(grib_get_long(h, "numberOfDataPoints", &numberOfDataPoints), 0);
131 GRIB_CHECK(grib_get_double(h, "latitudeOfFirstGridPointInDegrees", &lat1), 0);
132 GRIB_CHECK(grib_get_double(h, "longitudeOfFirstGridPointInDegrees", &lon1), 0);
133 GRIB_CHECK(grib_get_double(h, "latitudeOfLastGridPointInDegrees", &lat2), 0);
134 GRIB_CHECK(grib_get_double(h, "longitudeOfLastGridPointInDegrees", &lon2), 0);
135 GRIB_CHECK(grib_get_double(h, "iDirectionIncrementInDegrees", &iDirectionIncrementInDegrees), 0);
136
137 angular_tolerance = get_precision(edition);
138
139 if (N <= 0) {
140 error(filename, msg_num, "N should be > 0\n", N);
141 }
142 if (Nj != 2 * N) {
143 error(filename, msg_num, "Nj is %ld but should be 2*N (%ld)\n", Nj, 2 * N);
144 }
145
146 if (lon1 != 0) {
147 error(filename, msg_num, "latitudeOfFirstGridPointInDegrees=%f but should be 0\n", lon1);
148 }
149 expected_lon2 = 360.0 - 360.0 / (4 * N);
150
151 /* Check first and last latitudes */
152 if (lat1 != -lat2) {
153 error(filename, msg_num, "First latitude must be = last latitude but opposite in sign: lat1=%f, lat2=%f\n",
154 lat1, lat2);
155 }
156 /* Note: grib_get_gaussian_latitudes() assumes the 'lats' array has 2N elements! */
157 /* So do not allocate Nj */
158 lats = (double*)malloc(sizeof(double) * 2 * N);
159 GRIB_CHECK(grib_get_gaussian_latitudes(N, lats), 0);
160
161 if (!DBL_EQUAL(lats[0], lat1, angular_tolerance)) {
162 error(filename, msg_num, "latitudeOfFirstGridPointInDegrees=%f but should be %f\n", lat1, lats[0]);
163 }
164 if (!DBL_EQUAL(lats[Nj - 1], lat2, angular_tolerance)) {
165 error(filename, msg_num, "latitudeOfLastGridPointInDegrees=%f but should be %f\n", lat2, lats[Nj - 1]);
166 }
167
168 if (is_reduced) {
169 int pl_sum = 0, max_pl = 0, is_missing_Ni = 0, is_missing_Di = 0;
170 size_t i = 0, pl_len = 0;
171 long is_octahedral = 0;
172 is_missing_Ni = grib_is_missing(h, "Ni", &err);
173 assert(err == GRIB_SUCCESS);
174 is_missing_Di = grib_is_missing(h, "iDirectionIncrement", &err);
175 assert(err == GRIB_SUCCESS);
176 if (!is_missing_Ni) {
177 error(filename, msg_num, "For a reduced gaussian grid Ni should be missing\n");
178 }
179 if (!is_missing_Di) {
180 error(filename, msg_num, "For a reduced gaussian grid iDirectionIncrement should be missing\n");
181 }
182
183 GRIB_CHECK(grib_get_size(h, "pl", &pl_len), 0);
184 assert(pl_len > 0);
185 if (pl_len != 2 * N) {
186 error(filename, msg_num, "Length of pl array is %ld but should be 2*N (%ld)\n", pl_len, 2 * N);
187 }
188 pl = (long*)malloc(pl_len * sizeof(long));
189 assert(pl);
190 GRIB_CHECK(grib_get_long_array(h, "pl", pl, &pl_len), 0);
191 max_pl = pl[0];
192
193 /* Check pl is symmetric */
194 for (i = 0; i < pl_len / 2; ++i) {
195 const long pl_start = pl[i];
196 const long pl_end = pl[pl_len - 1 - i];
197 if (pl_start != pl_end) {
198 error(filename, msg_num, "pl array is not symmetric: pl[%ld]=%ld, pl[%ld]=%ld\n",
199 i, pl_start, pl_len - 1 - i, pl_end);
200 }
201 }
202
203 /* Check sum of pl array and total number of points */
204 for (i = 0; i < pl_len; ++i) {
205 pl_sum += pl[i];
206 if (pl[i] > max_pl)
207 max_pl = pl[i];
208 }
209 if (pl_sum != numberOfDataPoints) {
210 error(filename, msg_num, "Sum of pl array %ld does not match numberOfDataPoints %ld\n", pl_sum, numberOfDataPoints);
211 }
212 GRIB_CHECK(grib_get_long(h, "isOctahedral", &is_octahedral), 0);
213 if (is_octahedral) {
214 if (verbose)
215 printf(" This is an Octahedral Gaussian grid\n");
216 expected_lon2 = 360.0 - 360.0 / max_pl;
217 }
218 free(pl);
219 }
220
221 if (fabs(lon2 - expected_lon2) > angular_tolerance) {
222 error(filename, msg_num, "longitudeOfLastGridPointInDegrees=%f but should be %f\n", lon2, expected_lon2);
223 }
224
225 GRIB_CHECK(grib_get_size(h, "values", &sizeOfValuesArray), 0);
226 if (sizeOfValuesArray != numberOfDataPoints) {
227 error(filename, msg_num, "Number of data points %d different from size of values array %d\n",
228 numberOfDataPoints, sizeOfValuesArray);
229 }
230
231 free(lats);
232 grib_handle_delete(h);
233 }
234 fclose(in);
235 if (verbose)
236 printf("\n");
237 return GRIB_SUCCESS;
238 }
239
main(int argc,char ** argv)240 int main(int argc, char** argv)
241 {
242 int i = 0;
243
244 if (argc < 2) {
245 usage(argv[0]);
246 return 1;
247 }
248
249 for (i = 1; i < argc; ++i) {
250 const char* arg = argv[i];
251 if (STR_EQUAL(arg, "-f")) {
252 if (argc < 3) {
253 usage(argv[0]);
254 return 1;
255 }
256 exit_on_error = 0;
257 }
258 else if (STR_EQUAL(arg, "-v")) {
259 if (argc < 3) {
260 usage(argv[0]);
261 return 1;
262 }
263 verbose = 1;
264 }
265 else {
266 /* We have a file (not an option) */
267 process_file(arg);
268 }
269 }
270
271 if (verbose)
272 printf("###############\n");
273 if (error_count == 0) {
274 if (verbose)
275 printf("ALL OK\n");
276 }
277 else {
278 printf("Error count: %d\n", error_count);
279 return 1;
280 }
281
282 return 0;
283 }
284