1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Modali Kameswarrao
5
6 */
7
8 #ifdef HAVE_CONFIG_H
9 #include "config.h" /* HAVE_LIBMAGICS */
10 #endif
11
12 #include <climits>
13 #include <cdi.h>
14 #include <cctype>
15
16 #include "process_int.h"
17 #include <mpim_grid.h>
18
19 #ifdef HAVE_LIBMAGICS
20
21 #include "magics_api.h"
22 #include "magics_template_parser.h"
23 #include "results_template_parser.h"
24 #include "string_utilities.h"
25 #include "util_string.h"
26 #include "cdo_vlist.h"
27 #include "printinfo.h"
28
29 #define DBG 0
30
31 const char *line_colours[] = {
32 "red",
33 "green",
34 "blue",
35 "yellow",
36 "cyan",
37 "magenta",
38 "avocado",
39 "beige",
40 "brick",
41 "brown",
42 "burgundy",
43 "charcoal",
44 "chestnut",
45 "coral",
46 "cream",
47 "evergreen",
48 "gold",
49 "khaki",
50 "kellygreen",
51 "lavender",
52 "mustard",
53 "navy",
54 "ochre",
55 "olive",
56 "peach",
57 "pink",
58 "rose",
59 "rust",
60 "sky",
61 "tan",
62 "tangerine",
63 "turquoise",
64 "violet",
65 "reddishpurple",
66 "purplered",
67 "purplishred",
68 "orangishred",
69 "redorange",
70 "reddishorange",
71 "orange",
72 "yellowishorange",
73 "orangeyellow",
74 "orangishyellow",
75 "greenishyellow",
76 "yellowgreen",
77 "yellowishgreen",
78 "bluishgreen",
79 "bluegreen",
80 "greenishblue",
81 "purplishblue",
82 "bluepurple",
83 "bluishpurple",
84 "purple",
85 };
86
87 const char *graph_params[] = { "ymin", "ymax", "sigma", "stat", "obsv", "device", "linewidth" };
88
89 int graph_param_count = sizeof(graph_params) / sizeof(char *);
90 int num_colours = sizeof(line_colours) / sizeof(char *);
91
92 extern int checkdevice(char *device_in);
93
94 extern char *DEVICE;
95 extern char *DEVICE_TABLE;
96 extern int DEVICE_COUNT;
97
98 static int
compareDate(int64_t date1,int64_t date2)99 compareDate(int64_t date1, int64_t date2)
100 {
101 int c1[3], c2[3];
102 cdiDecodeDate(date1, c1, c1 + 1, c1 + 2);
103 cdiDecodeDate(date2, c2, c2 + 1, c2 + 2);
104
105 for (int i = 0; i < 3; i++)
106 {
107 const auto flag = c1[i] - c2[i];
108 if (flag > 0)
109 return 1;
110 else if (flag < 0)
111 return -1;
112 else
113 continue;
114 }
115
116 return 0;
117 }
118
119 static int
compareTime(int time1,int time2)120 compareTime(int time1, int time2)
121 {
122 int c1[3], c2[3];
123 cdiDecodeTime(time1, c1, c1 + 1, c1 + 2);
124 cdiDecodeTime(time2, c2, c2 + 1, c2 + 2);
125
126 for (int i = 0; i < 3; i++)
127 {
128 const auto flag = c1[i] - c2[i];
129 if (flag > 0)
130 return 1;
131 else if (flag < 0)
132 return -1;
133 else
134 continue;
135 }
136
137 return 0;
138 }
139
140 static void
maggraph(const char * plotfile,const char * varname,const char * varunits,long nfiles,std::vector<long> nts,std::vector<std::vector<int64_t>> vdates,std::vector<std::vector<int>> vtimes,std::vector<std::vector<double>> datatab,int nparam,std::vector<std::string> & params)141 maggraph(const char *plotfile, const char *varname, const char *varunits, long nfiles, std::vector<long> nts,
142 std::vector<std::vector<int64_t>> vdates, std::vector<std::vector<int>> vtimes, std::vector<std::vector<double>> datatab,
143 int nparam, std::vector<std::string> ¶ms)
144 {
145 char min_date_time_str[1024], max_date_time_str[1024];
146 int min_index = 0, max_index = 0;
147 char legend_text_data[256];
148 int num_sigma = 2;
149 bool stat = false, obsv = false;
150 int file_begin = 0;
151 int count;
152 int ret;
153 long tsID, fileID, i, ntime_steps = 0;
154 double min_val = 1.0e+200, max_val = -1.0e+200;
155 double y_min_val = 1.0e+200, y_max_val = -1.0e+200;
156 int linewidth_val = 8;
157
158 if (DBG)
159 {
160 fprintf(stderr, "Num params %d\n", nparam);
161
162 for (i = 0; i < nparam; i++) fprintf(stderr, "Param %s\n", params[i].c_str());
163 }
164
165 std::string temp_str;
166 for (i = 0; i < nparam; ++i)
167 {
168 const auto splitStrings = cstr_split_with_seperator(params[i].c_str(), "=");
169 const auto &key = splitStrings[0];
170 const auto &value = splitStrings[1];
171
172 if (key == "obsv")
173 {
174 temp_str = string_to_lower(value);
175 if (temp_str == "true")
176 {
177 obsv = true;
178 file_begin = 1;
179 if (DBG) fprintf(stderr, "OBSV true\n");
180 }
181 }
182
183 if (key == "stat")
184 {
185 temp_str = string_to_lower(value);
186 if (temp_str == "true")
187 {
188 stat = true;
189 if (DBG) fprintf(stderr, "STAT true\n");
190 }
191 }
192
193 if (key == "ymin")
194 {
195 y_min_val = std::stod(value);
196 if (DBG) fprintf(stderr, "Y min Val %g\n", y_min_val);
197 }
198
199 if (key == "ymax")
200 {
201 y_max_val = std::stod(value);
202 if (DBG) fprintf(stderr, "Y max Val %g\n", y_max_val);
203 }
204
205 if (key == "linewidth")
206 {
207 linewidth_val = std::stoi(value);
208 if (DBG) fprintf(stderr, "linewidth Val %d\n", linewidth_val);
209 }
210
211 if (key == "sigma")
212 {
213 num_sigma = std::stod(value);
214 if (DBG) fprintf(stderr, "SIGMA %d\n", num_sigma);
215 }
216
217 if (key == "device")
218 {
219 temp_str = string_to_upper(value);
220 DEVICE = strdup(temp_str.c_str());
221 if (DBG) fprintf(stderr, "DEVICE %s\n", DEVICE);
222
223 mag_setc("output_format", DEVICE);
224 }
225 }
226
227 if (DBG)
228 {
229 ntime_steps = nts[0];
230 fprintf(stderr, " nfiles=%ld ntime_steps=%ld\n", nfiles, ntime_steps);
231 fprintf(stderr, "STAT %d\n", (int) stat);
232 }
233
234 if (stat)
235 {
236 ntime_steps = nts[0];
237
238 for (fileID = 1; fileID < nfiles; fileID++)
239 {
240 if (nts[fileID] != ntime_steps)
241 {
242 cdo_warning(" Unequal number of time steps! Statistics disabled.");
243 stat = false;
244 break;
245 }
246
247 // First date & time of the present file
248 if (compareDate(vdates[0][0], vdates[fileID][0]))
249 {
250 cdo_warning(" Incosistent start date! Statistics disabled.");
251 stat = false;
252 break;
253 }
254
255 // First time of the present file
256 if (compareTime(vtimes[0][0], vtimes[fileID][0]))
257 {
258 cdo_warning(" Incosistent start time! Statistics disabled.");
259 stat = false;
260 break;
261 }
262
263 // Last date of the present file
264 if (compareDate(vdates[fileID][nts[fileID] - 1], vdates[0][nts[0] - 1]))
265 {
266 cdo_warning(" Incosistent end date! Statistics disabled.");
267 stat = false;
268 break;
269 }
270
271 // Last time of the present file
272 if (compareTime(vtimes[fileID][nts[fileID] - 1], vtimes[0][nts[0] - 1]))
273 {
274 cdo_warning(" Incosistent end time! Statistics disabled.");
275 stat = false;
276 break;
277 }
278 }
279 }
280
281 if (DBG) fprintf(stderr, "STAT %d\n", (int) stat);
282
283 char ***date_time_str = (char ***) malloc(nfiles * sizeof(char **));
284
285 std::vector<double> date_time;
286 std::vector<double> mean_val, std_dev_val;
287 std::vector<double> spread_min, spread_max;
288
289 if (stat)
290 {
291 // if all files are of same number of steps, only one date_time_str array is being used
292 date_time_str[0] = (char **) malloc(ntime_steps * sizeof(char *));
293
294 date_time.resize(ntime_steps);
295 mean_val.resize(ntime_steps);
296 std_dev_val.resize(ntime_steps);
297 spread_min.resize(ntime_steps);
298 spread_max.resize(ntime_steps);
299
300 for (tsID = 0; tsID < ntime_steps; ++tsID)
301 {
302 date_time[tsID] = tsID + 1;
303 date_time_str[0][tsID] = (char *) malloc(256);
304 sprintf(date_time_str[0][tsID], "%s %s", date_to_string(vdates[0][tsID]).c_str(),
305 time_to_string(vtimes[0][tsID]).c_str());
306 mean_val[tsID] = 0.;
307 std_dev_val[tsID] = 0.;
308
309 if (DBG)
310 {
311 fprintf(stderr, "tsID=%ld: %s\n", tsID, date_time_str[0][tsID]);
312 fprintf(stderr, "%6ld %6d", (long) vdates[0][tsID], vtimes[0][tsID]);
313 }
314
315 for (fileID = 0; fileID < nfiles; ++fileID)
316 {
317 if (DBG) fprintf(stderr, "fileID=%ld\n", fileID);
318
319 if (datatab[fileID][tsID] < min_val) min_val = datatab[fileID][tsID];
320 if (datatab[fileID][tsID] > max_val) max_val = datatab[fileID][tsID];
321
322 mean_val[tsID] += datatab[fileID][tsID];
323 std_dev_val[tsID] = 0.;
324 spread_min[tsID] = 0.;
325 spread_max[tsID] = 0.;
326
327 if (DBG)
328 {
329 fprintf(stderr, " %6g", datatab[fileID][tsID]);
330 fprintf(stderr, "\n");
331 }
332 }
333 }
334
335 for (tsID = 0; tsID < ntime_steps; ++tsID)
336 {
337 mean_val[tsID] /= (double) nfiles;
338 spread_min[tsID] = mean_val[tsID];
339 spread_max[tsID] = mean_val[tsID];
340
341 for (fileID = 0; fileID < nfiles; ++fileID)
342 {
343 std_dev_val[tsID] += (datatab[fileID][tsID] - mean_val[tsID]) * (datatab[fileID][tsID] - mean_val[tsID]);
344 }
345 std_dev_val[tsID] /= (double) nfiles;
346 std_dev_val[tsID] = std::pow(std_dev_val[tsID], 0.5);
347
348 if (DBG) fprintf(stderr, " Mean : %g Std Dev: %g\n", mean_val[tsID], std_dev_val[tsID]);
349
350 spread_min[tsID] = mean_val[tsID] - num_sigma * std_dev_val[tsID];
351 spread_max[tsID] = mean_val[tsID] + num_sigma * std_dev_val[tsID];
352
353 if (DBG) fprintf(stderr, " Min : %g Max: %g\n", spread_min[tsID], spread_max[tsID]);
354 }
355
356 for (tsID = 0; tsID < ntime_steps; ++tsID)
357 {
358 if (spread_min[tsID] < min_val) min_val = spread_min[tsID];
359 if (spread_max[tsID] > max_val) max_val = spread_max[tsID];
360 }
361
362 if (DBG)
363 {
364 fprintf(stderr, " %6g %6g\n", min_val, max_val);
365 fprintf(stderr, " %s %s\n", date_time_str[0][0], date_time_str[0][ntime_steps - 1]);
366 fprintf(stderr, "\n");
367 }
368
369 strcpy(min_date_time_str, date_time_str[0][0]);
370 strcpy(max_date_time_str, date_time_str[0][ntime_steps - 1]);
371 }
372 else
373 {
374 /* Find the min_date_time_str from the min's of nfiles
375 Find the max_date_time_str from the max's of nfiles
376 Construct the date_time_str array
377 */
378
379 if (DBG) fprintf(stderr, "STAT %d\n", (int) stat);
380
381 for (fileID = 0; fileID < nfiles; fileID++)
382 {
383 if (DBG) fprintf(stderr, "FILE %ld\n", fileID);
384 date_time.resize(nts[fileID]);
385 date_time_str[fileID] = (char **) malloc(nts[fileID] * sizeof(char *));
386
387 for (tsID = 0; tsID < nts[fileID]; ++tsID)
388 {
389 date_time[tsID] = tsID + 1;
390
391 date_time_str[fileID][tsID] = (char *) malloc(256);
392 sprintf(date_time_str[fileID][tsID], "%s %s", date_to_string(vdates[fileID][tsID]).c_str(),
393 time_to_string(vtimes[fileID][tsID]).c_str());
394 if (DBG && (tsID == 0 || tsID == nts[fileID] - 1))
395 fprintf(stderr, "%s %s %s\n", date_to_string(vdates[fileID][tsID]).c_str(),
396 time_to_string(vtimes[fileID][tsID]).c_str(), date_time_str[fileID][tsID]);
397
398 if (datatab[fileID][tsID] < min_val) min_val = datatab[fileID][tsID];
399 if (datatab[fileID][tsID] > max_val) max_val = datatab[fileID][tsID];
400 }
401
402 if (fileID == 0)
403 {
404 if (DBG) fprintf(stderr, "\n %s %s\n", date_time_str[fileID][0], date_time_str[fileID][nts[0] - 1]);
405 min_index = 0;
406 max_index = 0;
407 }
408 else
409 {
410 ret = compareDate(vdates[min_index][0], vdates[fileID][0]);
411 if (ret == 1)
412 min_index = fileID;
413 else if (!ret)
414 {
415 ret = compareTime(vtimes[min_index][0], vtimes[fileID][0]);
416 if (ret == -999)
417 cdo_abort("Error in input Date Time");
418 else if (ret == 1)
419 min_index = fileID;
420 }
421 if (DBG) fprintf(stderr, "Min File ID %d\n", min_index);
422
423 if (DBG) fprintf(stderr, "compareDateOrTime %s\n", date_time_str[fileID][nts[fileID] - 1]);
424
425 ret = compareDate(vdates[max_index][nts[max_index] - 1], vdates[fileID][nts[fileID] - 1]);
426 if (ret == -1)
427 max_index = fileID;
428 else if (!ret)
429 {
430 ret = compareTime(vtimes[max_index][nts[max_index] - 1], vtimes[fileID][nts[fileID] - 1]);
431 if (ret == -999)
432 cdo_abort("Error in input Date Time");
433 else if (ret == -1)
434 max_index = fileID;
435 }
436
437 if (DBG) fprintf(stderr, "Max File ID %d\n", max_index);
438 }
439 }
440
441 strcpy(min_date_time_str, date_time_str[min_index][0]);
442 strcpy(max_date_time_str, date_time_str[max_index][nts[max_index] - 1]);
443 if (DBG) fprintf(stderr, "%s %s\n", min_date_time_str, max_date_time_str);
444 }
445
446 if (DBG) fprintf(stderr, "%s %s\n", min_date_time_str, max_date_time_str);
447
448 auto splitStrings = cstr_split_with_seperator(max_date_time_str, "-");
449
450 int num_years = std::stoi(splitStrings[0]);
451 int num_months = std::stoi(splitStrings[1]);
452 int num_days = std::stoi(splitStrings[2]);
453
454 splitStrings = cstr_split_with_seperator(min_date_time_str, "-");
455 num_years -= std::stoi(splitStrings[0]);
456
457 if (num_years <= 1)
458 {
459 if (num_years == 1)
460 num_months += (12 - std::stoi(splitStrings[1]));
461 else
462 num_months -= (std::stoi(splitStrings[1]));
463
464 if (!num_months)
465 num_days -= std::stoi(splitStrings[2]);
466 else if (num_months == 1)
467 num_days += (31 - std::stoi(splitStrings[2]));
468 }
469
470 if (DBG) fprintf(stderr, " num_years=%d num_months=%d num_days=%d\n", num_years, num_months, num_days);
471
472 /*
473 1. Loop over the Files
474 2. Loop over the number of time steps
475 3. Set the attributes for the magics data and plot
476 */
477
478 // magics_template_parser( magics_node );
479
480 mag_setc("output_name", plotfile);
481 mag_setc("subpage_map_projection", "cartesian");
482 mag_setr("subpage_y_length", 14.);
483 mag_setr("subpage_y_position", 1.5);
484
485 // Horizontal Axis attributes
486 mag_setc("axis_orientation", "horizontal");
487 mag_setc("axis_grid", "on");
488 mag_setc("axis_grid_colour", "grey");
489 mag_seti("axis_grid_thickness", 1);
490 mag_setc("axis_grid_line_style", "dot");
491 mag_setc("axis_type", "date");
492
493 const char *dateType = (num_years > 1) ? "years" : (num_months > 1) ? "months" : (num_months == 1 || num_days) ? "days" : "hours";
494 mag_setc("axis_date_type", dateType);
495
496 mag_setc("axis_date_min_value", min_date_time_str);
497 mag_setc("axis_date_max_value", max_date_time_str);
498 mag_setc("axis_title_text", "Time");
499 mag_setc("axis_title_orientation", "horizontal");
500
501 mag_seti("axis_tick_label_frequency", 2);
502 mag_setr("axis_years_label_height", 0.4);
503
504 mag_axis();
505
506 // Vertical Axis attributes
507 mag_setc("axis_orientation", "vertical");
508 mag_setc("axis_grid", "on");
509 mag_setc("axis_type", "regular");
510 mag_setc("axis_grid_colour", "grey");
511 mag_seti("axis_grid_thickness", 1);
512 mag_setc("axis_grid_line_style", "dot");
513
514 // To redefine the y- axis scale based on user input in .xml file
515
516 // min & max values from the input data files
517 mag_setr("axis_min_value", min_val);
518 mag_setr("axis_max_value", max_val);
519
520 // min & max values specified by the user in the command line args
521 if (y_min_val < 1.0e+200) mag_setr("axis_min_value", y_min_val);
522 if (y_max_val > -1.0e+200) mag_setr("axis_max_value", y_max_val);
523
524 mag_setc("axis_title_text", varname);
525
526 mag_setc("axis_title_orientation", "vertical");
527
528 mag_seti("axis_tick_label_frequency", 2);
529 mag_setr("axis_tick_label_height", 0.5);
530
531 mag_axis();
532
533 // Legend
534 mag_setc("legend", "on");
535 mag_setc("legend_text_colour", "black");
536
537 mag_setc("graph_symbol", "off");
538 mag_seti("graph_line_thickness", linewidth_val);
539
540 if (DBG) fprintf(stderr, "FILE BEGIN %d\n", file_begin);
541
542 for (i = file_begin; i < nfiles; ++i)
543 {
544 count = i;
545 if (obsv) count = i - 1;
546 if (DBG) fprintf(stderr, "Current File %ld\n", i);
547 // sprintf(legend_text_data, "ens_%d", count + 1);
548 sprintf(legend_text_data, "data_%d", count + 1);
549 mag_setc("graph_line_colour", line_colours[count % num_colours]);
550 mag_setc("legend_user_text", legend_text_data);
551 if (stat)
552 mag_set1c("graph_curve_date_x_values", (const char **) date_time_str[0], ntime_steps);
553 else
554 mag_set1c("graph_curve_date_x_values", (const char **) date_time_str[i], nts[i]);
555
556 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * begin
557 mag_setr("graph_x_suppress_below", (double) LLONG_MIN);
558 mag_setr("graph_x_suppress_above", (double) LLONG_MAX);
559 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * end
560
561 mag_set1r("graph_curve_y_values", datatab[i].data(), nts[i]);
562 mag_graph();
563 }
564
565 if (obsv)
566 {
567 mag_setc("graph_line_colour", "black");
568 sprintf(legend_text_data, "%s", "Obsv");
569 mag_setc("legend_user_text", legend_text_data);
570 mag_set1c("graph_curve_date_x_values", (const char **) date_time_str[0], nts[0]);
571
572 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * begin
573 mag_setr("graph_x_suppress_below", (double) LLONG_MIN);
574 mag_setr("graph_x_suppress_above", (double) LLONG_MAX);
575 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * end
576
577 mag_set1r("graph_curve_y_values", datatab[0].data(), nts[0]);
578 mag_setc("graph_line_style", "dot");
579 mag_seti("graph_line_thickness", linewidth_val + 2);
580 mag_graph();
581 }
582
583 if (DBG) fprintf(stderr, "NTIME STEPS %ld\n", ntime_steps);
584
585 if (stat)
586 {
587 if (DBG) fprintf(stderr, "NTIME STEPS %ld\n", ntime_steps);
588
589 mag_seti("graph_line_thickness", linewidth_val);
590 mag_setc("graph_line_colour", "grey");
591 mag_setc("graph_line_style", "dash");
592 mag_set1c("graph_curve_date_x_values", (const char **) date_time_str[0], ntime_steps);
593
594 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * begin
595 mag_setr("graph_x_suppress_below", (double) LLONG_MIN);
596 mag_setr("graph_x_suppress_above", (double) LLONG_MAX);
597 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * end
598
599 mag_set1r("graph_curve_y_values", mean_val.data(), ntime_steps);
600 sprintf(legend_text_data, "Mean");
601 mag_setc("legend_user_text", legend_text_data);
602 mag_graph();
603
604 mag_reset("graph_type");
605 mag_setc("graph_type", "area");
606 mag_seti("graph_line_thickness", 1);
607 mag_setc("graph_shade_style", "dot");
608 mag_setr("graph_shade_dot_size", 1.);
609 mag_set1c("graph_curve2_date_x_values", (const char **) date_time_str[0], ntime_steps);
610 mag_set1r("graph_curve2_y_values", spread_max.data(), ntime_steps);
611 mag_set1c("graph_curve_date_x_values", (const char **) date_time_str[0], ntime_steps);
612 mag_set1r("graph_curve_y_values", spread_min.data(), ntime_steps);
613 mag_setc("graph_shade_colour", "grey");
614 sprintf(legend_text_data, "%dSigma", num_sigma);
615 mag_setc("legend_user_text", legend_text_data);
616
617 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * begin
618 mag_setr("graph_x_suppress_below", (double) LLONG_MIN);
619 mag_setr("graph_x_suppress_above", (double) LLONG_MAX);
620 // TEMPORARY FIX, UNITL NEW MAGICS LIBRARY RELEASE * end
621
622 mag_graph();
623 }
624
625 char *lines[1];
626 lines[0] = (char *) malloc(1024);
627 // To be obtained from Meta Data
628 // sprintf( lines[0],"%s","ExpID : " );
629 // sprintf( lines[0],"%sxxxx Variable : %s[%s]",lines[0], varname, varunits );
630 // sprintf( lines[0],"Variable : %s[%s]",varname, varunits );
631 // sprintf( lines[0],"%s Date : %s --%s",lines[0], min_date_time_str, max_date_time_str );
632 sprintf(lines[0], "Variable : %s[%s] Date : %s --%s", varname, varunits, min_date_time_str, max_date_time_str);
633 mag_set1c("text_lines", (const char **) lines, 1);
634
635 mag_setc("text_html", "true");
636 mag_setc("text_colour", "black");
637 mag_setr("text_font_size", 0.6);
638 mag_setc("text_mode", "positional");
639 mag_setr("text_box_x_position", 1.5);
640 mag_setr("text_box_y_position", 16.5);
641 mag_setr("text_box_x_length", 20.);
642 mag_setr("text_box_y_length", 2.5);
643 mag_setc("text_border", "off");
644 mag_setc("text_justification", "left");
645 mag_text();
646
647 free(date_time_str);
648
649 if (DBG) fprintf(stderr, "lines=%s\n", lines[0]);
650
651 free(lines[0]);
652 }
653
654 static void
init_MAGICS()655 init_MAGICS()
656 {
657 setenv("MAGPLUS_QUIET", "1", 1); /* To suppress magics messages */
658 mag_open();
659
660 // Some standard parameters affectng the magics environment, moved from the xml file ** begin **
661 mag_setc("page_id_line", "off");
662 // Some standard parameters affectng the magics environment, moved from the xml file ** end **
663 }
664
665 static void
quit_MAGICS()666 quit_MAGICS()
667 {
668 mag_close();
669 if (DBG) fprintf(stdout, "Exiting From MAGICS\n");
670 }
671
672 static void
VerifyGraphParameters(int num_param,std::vector<std::string> & param_names)673 VerifyGraphParameters(int num_param, std::vector<std::string> ¶m_names)
674 {
675 int i, j;
676 bool found = false, syntax = true, halt_flag = false;
677 char *temp_str;
678
679 for (i = 0; i < num_param; ++i)
680 {
681 found = false;
682 syntax = true;
683 const auto splitStrings = cstr_split_with_seperator(param_names[i].c_str(), "=");
684
685 if (splitStrings.size() > 1)
686 {
687 const auto &key = splitStrings[0];
688 const auto &value = splitStrings[1];
689
690 for (j = 0; j < graph_param_count; ++j)
691 {
692 if (key == graph_params[j])
693 {
694 found = true;
695 if (key == "obsv" || key == "stat")
696 {
697 if (cstr_is_numeric(value.c_str()))
698 syntax = false;
699 else
700 {
701 temp_str = strdup(value.c_str());
702 cstr_to_lower_case(temp_str);
703 if (strcmp(temp_str, "true") && strcmp(temp_str, "false")) syntax = false;
704 }
705 }
706
707 if (key == "ymin" || key == "ymax" || key == "sigma" || key == "linewidth")
708 {
709 if (!cstr_is_numeric(value.c_str())) syntax = false;
710 }
711
712 if (key == "device")
713 {
714 if (cstr_is_numeric(value.c_str()))
715 syntax = false;
716 else
717 {
718 if (DBG) fprintf(stderr, "Parameter value '%s'\n", value.c_str());
719 char *deviceCstr = strdup(value.c_str());
720 if (checkdevice(deviceCstr)) syntax = false;
721
722 // Graph not supported in google earth format
723 if (value == "GIF_ANIMATION" || value == "gif_animation")
724 {
725 syntax = false;
726 fprintf(stderr, "Animation not supported for Graph!\n");
727 if (DBG) fprintf(stderr, "Parameter value '%s'\n", value.c_str());
728 }
729 if (value == "KML" || value == "kml")
730 {
731 syntax = false;
732 fprintf(stderr, " 'kml' format not supported for Graph!\n");
733 if (DBG) fprintf(stderr, "Parameter value '%s'\n", value.c_str());
734 }
735 }
736 }
737
738 /*
739 if(key == "xml")
740 {
741 if( ( fp = fopen(value.c_str(),"r") ) == nullptr )
742 {
743 fprintf( stderr,"Input XML File not found in specified path '%s'\n", value.c_str() );
744 halt_flag = true;
745 }
746 else
747 {
748 // HARDCODED THE FILE NAME .. TO BE SENT AS COMMAND LINE ARGUMENT FOR THE MAGICS OPERATOR
749 fclose(fp);
750 init_XML_template_parser( value.c_str() ); updatemagics_and_results_nodes( );
751 }
752 }
753 */
754 }
755 }
756 }
757 else
758 {
759 syntax = false;
760 }
761
762 if (!found)
763 {
764 halt_flag = true;
765 fprintf(stderr, "Unknown parameter '%s'!\n", param_names[i].c_str());
766 }
767 if (found && !syntax)
768 {
769 halt_flag = true;
770 fprintf(stderr, "Invalid parameter specification '%s'!\n", param_names[i].c_str());
771 }
772 }
773
774 if (halt_flag) exit(0);
775 }
776 #endif
777
778 void *
Maggraph(void * process)779 Maggraph(void *process)
780 {
781 cdo_initialize(process);
782
783 #ifdef HAVE_LIBMAGICS
784 char varname[CDI_MAX_NAME], units[CDI_MAX_NAME];
785 int gridID;
786 int vlistID0 = -1;
787 size_t nts_alloc;
788
789 int nparam = cdo_operator_argc();
790 auto pnames = cdo_get_oper_argv();
791
792 if (nparam) VerifyGraphParameters(nparam, pnames);
793
794 int nfiles = cdo_stream_cnt() - 1;
795 const char *ofilename = cdo_get_stream_name(nfiles);
796
797 if (DBG)
798 {
799 fprintf(stderr, " Num of files %d\n", nfiles);
800 fprintf(stderr, " files %s\n", ofilename);
801 }
802
803 std::vector<std::vector<double>> datatab(nfiles);
804 std::vector<std::vector<int64_t>> vdates(nfiles);
805 std::vector<std::vector<int>> vtimes(nfiles);
806 std::vector<long> nts(nfiles, 0);
807
808 for (int fileID = 0; fileID < nfiles; fileID++)
809 {
810 if (DBG) fprintf(stderr, " file %d is %s\n", fileID, cdo_get_stream_name(fileID));
811
812 const auto streamID = cdo_open_read(fileID);
813 const auto vlistID = cdo_stream_inq_vlist(streamID);
814 const auto taxisID = vlistInqTaxis(vlistID);
815
816 vlistInqVarUnits(vlistID, 0, units);
817 if (DBG) fprintf(stderr, " units=%s\n", units);
818 if (fileID == 0)
819 {
820 vlistInqVarName(vlistID, 0, varname);
821
822 gridID = vlistInqVarGrid(vlistID, 0);
823 const auto zaxisID = vlistInqVarZaxis(vlistID, 0);
824 const auto nvars = vlistNvars(vlistID);
825
826 if (nvars > 1) cdo_abort("Input stream has more than on variable!");
827 if (gridInqSize(gridID) != 1) cdo_abort("Variable has more than one grid point!");
828 if (zaxisInqSize(zaxisID) != 1) cdo_abort("Variable has more than one level!");
829
830 vlistID0 = vlistDuplicate(vlistID);
831 }
832 else
833 {
834 vlist_compare(vlistID0, vlistID, CMP_ALL);
835 }
836
837 int tsID = 0;
838 nts_alloc = 0;
839 while (true)
840 {
841 const auto nrecs = cdo_stream_inq_timestep(streamID, tsID);
842 if (nrecs == 0) break;
843
844 if (nrecs != 1) cdo_abort("Input stream has more than one point in time!");
845
846 if ((size_t) tsID >= nts_alloc)
847 {
848 constexpr size_t NALLOC_INC = 1024;
849 nts_alloc += NALLOC_INC;
850 datatab[fileID].resize(nts_alloc);
851 vdates[fileID].resize(nts_alloc);
852 vtimes[fileID].resize(nts_alloc);
853 }
854
855 nts[fileID]++;
856
857 int varID, levelID;
858 cdo_inq_record(streamID, &varID, &levelID);
859 size_t nmiss;
860 double val;
861 cdo_read_record(streamID, &val, &nmiss);
862 datatab[fileID][tsID] = val;
863 vdates[fileID][tsID] = taxisInqVdate(taxisID);
864 vtimes[fileID][tsID] = taxisInqVtime(taxisID);
865
866 tsID++;
867 }
868
869 cdo_stream_close(streamID);
870 }
871
872 /* HARDCODED THE FILE NAME .. TO BE SENT AS COMMAND LINE ARGUMENT FOR THE MAGICS OPERATOR */
873 // init_XML_template_parser( Filename );
874 // updatemagics_and_results_nodes( );
875
876 init_MAGICS();
877
878 cdo_print(" Creating PLOT for %s", varname);
879 if (DBG)
880 {
881 fprintf(stderr, "Num params %d\n", nparam);
882
883 for (int i = 0; i < nparam; i++) fprintf(stderr, "Param %s\n", pnames[i].c_str());
884 }
885
886 maggraph(ofilename, varname, units, nfiles, nts, vdates, vtimes, datatab, nparam, pnames);
887
888 // quit_XML_template_parser( );
889
890 quit_MAGICS();
891
892 if (vlistID0 != -1) vlistDestroy(vlistID0);
893
894 #else
895
896 cdo_abort("MAGICS support not compiled in!");
897
898 #endif
899
900 cdo_finish();
901
902 return nullptr;
903 }
904