1 /*
2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data.
3
4 Author: Uwe Schulzweida
5
6 */
7
8 /*
9 This module contains the following operators:
10
11 Seltime seltimestep Select timesteps
12 Seltime seltime Select times
13 Seltime selhour Select hours
14 Seltime selday Select days
15 Seltime selmonth Select months
16 Seltime selyear Select years
17 Seltime selseason Select seasons
18 Seltime seldate Select dates
19 Seltime selsmon Select single month
20 */
21 #include <cassert>
22 #include <algorithm> // sort
23 #include <cdi.h>
24
25 #include "cdo_options.h"
26 #include "process_int.h"
27 #include "param_conversion.h"
28 #include "util_string.h"
29
30 static std::vector<int>
get_season_list(const std::vector<std::string> & seasonString)31 get_season_list(const std::vector<std::string> &seasonString)
32 {
33 std::vector<int> listArrayInt;
34 int imon[13] = { 0 }; // 1-12 !
35
36 const int nsel = seasonString.size();
37 if (isdigit(seasonString[0][0]))
38 {
39 for (int i = 0; i < nsel; i++)
40 {
41 const auto ival = parameter_to_int(seasonString[i]);
42 // clang-format off
43 if (ival == 1) { imon[12]++; imon[ 1]++; imon[ 2]++; }
44 else if (ival == 2) { imon[ 3]++; imon[ 4]++; imon[ 5]++; }
45 else if (ival == 3) { imon[ 6]++; imon[ 7]++; imon[ 8]++; }
46 else if (ival == 4) { imon[ 9]++; imon[10]++; imon[11]++; }
47 else cdo_abort("Season %d not available!", ival);
48 // clang-format on
49 }
50 }
51 else
52 {
53 for (int i = 0; i < nsel; i++) season_to_months(seasonString[i].c_str(), imon);
54 }
55
56 for (int i = 1; i < 13; ++i)
57 if (imon[i]) listArrayInt.push_back(i);
58
59 return listArrayInt;
60 }
61
62 // last update: 2020-02-17 (Oliver Heidmann)
63 /*
64 * Input:
65 * std::vector<std::string> => vector of length 1 or 2 containing string representing a range of dates or a single date to be
66 * selected Output: std::vector<doube> => vector of length 1 or 2 containing the to double converted date or date range
67 *
68 * This function turns a date range string into a list of double values which
69 * represent the start and the end of the range.
70 * when only one argument is given and it contains no time string (marked with
71 * T e.g 2006-01-01T23:00:30) we set the second value of the retrun vector to
72 * the end of the day from the first argument. if only one argument is given
73 * AND it contains a time string we return a vector of length one which
74 * contains only the converted value of the first argument. This function
75 * accepts the string "-". This represents the possible maximum start or end
76 * date depending on where the string is in the arguments. If the first
77 * argument is "-" we set the start date to -999999999 and to 999999999 if the
78 * string is in the second.
79 */
80
81 static std::vector<double>
get_date_list(const std::vector<std::string> & arguments)82 get_date_list(const std::vector<std::string> &arguments)
83 {
84 std::vector<double> dateList;
85
86 if (arguments.size() < 1) cdo_abort("Too few arguments!");
87 if (arguments.size() > 2) cdo_abort("Too many arguments!");
88
89 bool containsTime = false;
90 for (size_t i = 0; i < arguments.size(); i++)
91 {
92 double fval = 0.0;
93 // if "-" set start date to maximum start (i == 0) or end (i == 1)
94 if (arguments[i] == "-")
95 {
96 fval = (i == 0) ? -99999999999. : 99999999999.;
97 }
98 // get the double value representing the date and check for time string
99 else
100 {
101 fval = date_str_to_double(arguments[i].c_str(), 0);
102 if (strchr(arguments[i].c_str(), 'T')) containsTime = true;
103 }
104 dateList.push_back(fval);
105 }
106 // if two arguments and no time:
107 // set second argument to end of day
108 if (arguments.size() == 2 && !containsTime) dateList[1] += 0.999;
109
110 // if time and only one argument:
111 // set second date to first and set second to end of day
112 if (arguments.size() == 1 && !containsTime) dateList.push_back(dateList[0] + 0.999);
113
114 return dateList;
115 }
116
117 static std::string
get_argv_string(const std::vector<std::string> & argvVec)118 get_argv_string(const std::vector<std::string> &argvVec)
119 {
120 std::string s = argvVec[0];
121 for (size_t i = 1; i < argvVec.size(); i++) s += "," + argvVec[i];
122 return s;
123 }
124
125 static bool
argv_has_negativ_values(const std::vector<std::string> & argv)126 argv_has_negativ_values(const std::vector<std::string> &argv)
127 {
128 for (const auto &argument : argv)
129 {
130 int first, last, inc;
131 split_intstring(argument, first, last, inc);
132 if (first < 0 || last < 0) return true;
133 }
134
135 return false;
136 }
137
138 static std::vector<int>
cdo_argv_to_int_timestep(const std::vector<std::string> & argv,int ntimesteps)139 cdo_argv_to_int_timestep(const std::vector<std::string> &argv, int ntimesteps)
140 {
141 std::vector<int> v;
142
143 for (const auto &argument : argv)
144 {
145 int first, last, inc;
146 split_intstring(argument, first, last, inc);
147
148 if (inc == -1 && (first < 0 || last < 0)) inc = 1;
149 if (first < 0) first += ntimesteps + 1;
150 if (last < 0) last += ntimesteps + 1;
151
152 if (inc >= 0)
153 {
154 for (auto ival = first; ival <= last; ival += inc) v.push_back(ival);
155 }
156 else
157 {
158 for (auto ival = first; ival >= last; ival += inc) v.push_back(ival);
159 }
160 }
161
162 return v;
163 }
164
165 void *Select(void *process);
166
167 void *
Seltime(void * process)168 Seltime(void *process)
169 {
170 CdoStreamID streamID2 = CDO_STREAM_UNDEF;
171 int varID, levelID;
172 int nsel = 0;
173 int ncts = 0, nts;
174 int nts1 = 0, nts2 = 0;
175 int its1 = 0, its2 = 0;
176 double selfval = 0;
177 bool lconstout = false;
178 bool process_nts1 = false, process_nts2 = false;
179 std::vector<int64_t> vdate_list;
180 std::vector<int> vtime_list;
181 std::vector<int> intarr;
182 std::vector<double> fltarr;
183
184 cdo_initialize(process);
185
186 const auto dataIsUnchanged = data_is_unchanged();
187
188 enum
189 {
190 func_time,
191 func_date,
192 func_step,
193 func_datetime
194 };
195
196 // clang-format off
197 const auto SELTIMESTEP = cdo_operator_add("seltimestep", func_step, 0, "timesteps");
198 const auto SELDATE = cdo_operator_add("seldate", func_datetime, 0, "start date and end date (format YYYY-MM-DDThh:mm:ss)");
199 const auto SELTIME = cdo_operator_add("seltime", func_time, 0, "times (format hh:mm:ss)");
200 const auto SELHOUR = cdo_operator_add("selhour", func_time, 0, "hours");
201 const auto SELDAY = cdo_operator_add("selday", func_date, 0, "days");
202 const auto SELMONTH = cdo_operator_add("selmonth", func_date, 0, "months");
203 const auto SELYEAR = cdo_operator_add("selyear", func_date, 0, "years");
204 const auto SELSEASON = cdo_operator_add("selseason", func_date, 0, "seasons");
205 const auto SELSMON = cdo_operator_add("selsmon", func_date, 0, "month[,nts1[,nts2]]");
206 // clang-format on
207
208 const auto operatorID = cdo_operator_id();
209 const auto operfunc = cdo_operator_f1(operatorID);
210
211 operator_input_arg(cdo_operator_enter(operatorID));
212 const auto argv = get_argv_string(cdo_get_oper_argv());
213 std::string newCommand = "select,";
214
215 auto redirectEnabled = false; // only for seperating prototype from code (redirects execution to Select.cc)
216 auto redirectFound = true; // for opers that are not redirectable for now
217 auto ID = operatorID;
218 if (redirectEnabled)
219 {
220 // clang-format off
221 if (ID == SELTIMESTEP ) { newCommand += "timestep=" + argv; }
222 else if (ID == SELDATE) { newCommand += "date=" + argv; }
223 //else if(ID == SELTIME) { newCommand += "time=" + argv; } // unimplemented
224 else if (ID == SELHOUR) { newCommand += "hour=" + argv; }
225 else if (ID == SELDAY) { newCommand += "day=" + argv; }
226 else if (ID == SELMONTH) { newCommand += "month=" + argv; }
227 else if (ID == SELYEAR) { newCommand += "year=" + argv; }
228 else if (ID == SELSEASON) { newCommand += "season=" + argv; }
229 else if (ID == SELSMON) { newCommand += "month=" + argv; }
230 else { redirectFound = false; }
231 // clang-format on
232 }
233
234 if (redirectFound && redirectEnabled)
235 {
236 Debug(Yellow("Redirecting to %s"), newCommand);
237 ((Process *) process)->init_process("select", {newCommand});
238 return Select(process);
239 // If a redirect was found the entire process is ended through this return!
240 }
241
242 if (operatorID == SELSEASON)
243 {
244 intarr = get_season_list(cdo_get_oper_argv());
245 nsel = intarr.size();
246 }
247 else if (operatorID == SELDATE)
248 {
249 fltarr = get_date_list(cdo_get_oper_argv());
250 nsel = fltarr.size();
251 }
252 else if (operatorID == SELTIME)
253 {
254 nsel = cdo_operator_argc();
255 if (nsel < 1) cdo_abort("Too few arguments!");
256 intarr.reserve(nsel);
257 for (int i = 0; i < nsel; i++)
258 {
259 auto currentArgument = cdo_operator_argv(i).c_str();
260 if (strchr(currentArgument, ':'))
261 {
262 int hour = 0, minute = 0, second = 0;
263 sscanf(currentArgument, "%d:%d:%d", &hour, &minute, &second);
264 intarr.push_back(cdiEncodeTime(hour, minute, second));
265 }
266 else
267 {
268 intarr.push_back(parameter_to_int(currentArgument));
269 }
270 }
271 }
272 else if (operatorID != SELTIMESTEP)
273 {
274 intarr = cdo_argv_to_int(cdo_get_oper_argv());
275 nsel = intarr.size();
276 }
277
278 if (operatorID != SELTIMESTEP && nsel < 1) cdo_abort("No timestep selected!");
279
280 if (operatorID == SELSMON)
281 {
282 if (nsel > 1) nts1 = intarr[1];
283 nts2 = (nsel > 2) ? intarr[2] : nts1;
284
285 if (nsel > 3) cdo_abort("Too many parameters");
286
287 if (Options::cdoVerbose) cdo_print("mon=%d nts1=%d nts2=%d", intarr[0], nts1, nts2);
288
289 nsel = 1;
290 }
291
292 const auto streamID1 = cdo_open_read(0);
293
294 const auto vlistID1 = cdo_stream_inq_vlist(streamID1);
295 const auto vlistID2 = vlistDuplicate(vlistID1);
296
297 Field field;
298
299 VarList varList1;
300 varListInit(varList1, vlistID1);
301
302 // add support for negative timestep values
303 if (operatorID == SELTIMESTEP)
304 {
305 auto ntsteps = vlistNtsteps(vlistID1);
306
307 if (argv_has_negativ_values(cdo_get_oper_argv()))
308 {
309 if (ntsteps < 0)
310 {
311 if (cdo_assert_files_only())
312 {
313 int tsID = 0;
314 while (cdo_stream_inq_timestep(streamID1, tsID)) tsID++;
315 ntsteps = tsID;
316 if (Options::cdoVerbose) cdo_print("Found %d timesteps", ntsteps);
317 }
318
319 cdo_abort("Negative timesteps not supported in CDO pipes!");
320 }
321
322 intarr = cdo_argv_to_int_timestep(cdo_get_oper_argv(), ntsteps);
323 }
324 else
325 {
326 intarr = cdo_argv_to_int(cdo_get_oper_argv());
327 }
328
329 nsel = intarr.size();
330
331 if (nsel < 1) cdo_abort("No timestep selected!");
332 }
333
334 if (operatorID == SELTIMESTEP)
335 {
336 std::sort(intarr.begin(), intarr.end());
337 const auto ip = std::unique(intarr.begin(), intarr.end());
338 intarr.resize(std::distance(intarr.begin(), ip));
339 nsel = intarr.size();
340 }
341
342 if (Options::cdoVerbose)
343 {
344 for (int i = 0; i < nsel; i++)
345 {
346 if (operatorID == SELDATE)
347 cdo_print("fltarr entry: %d %14.4f", i + 1, fltarr[i]);
348 else
349 cdo_print("intarr entry: %d %d", i + 1, intarr[i]);
350 }
351 }
352
353 vlistDefNtsteps(vlistID2, (nsel == 1 && operfunc == func_step) ? 1 : -1);
354
355 const auto taxisID1 = vlistInqTaxis(vlistID1);
356 const auto taxisID2 = taxisDuplicate(taxisID1);
357 vlistDefTaxis(vlistID2, taxisID2);
358
359 std::vector<bool> selfound(nsel, false);
360
361 int iselmax = -1;
362 if (operatorID != SELDATE)
363 for (int i = 0; i < nsel; i++) iselmax = std::max(iselmax, intarr[i]);
364
365 const auto nvars = vlistNvars(vlistID1);
366 int nconst = 0;
367 for (varID = 0; varID < nvars; varID++)
368 if (varList1[varID].timetype == TIME_CONSTANT) nconst++;
369
370 FieldVector3D vars;
371
372 const auto lnts1 = (operatorID == SELSMON) && (nts1 > 0);
373
374 if (lnts1 || nconst)
375 {
376 if (lnts1)
377 {
378 vdate_list.resize(nts1);
379 vtime_list.resize(nts1);
380 }
381 else
382 {
383 nts1 = 1;
384 }
385
386 vars.resize(nts1);
387
388 for (int tsID = 0; tsID < nts1; tsID++)
389 {
390 fields_from_vlist(vlistID1, vars[tsID]);
391
392 for (varID = 0; varID < nvars; varID++)
393 {
394 if (lnts1 || varList1[varID].timetype == TIME_CONSTANT)
395 {
396 const auto nlevels = varList1[varID].nlevels;
397 for (levelID = 0; levelID < nlevels; levelID++) vars[tsID][varID][levelID].resize(varList1[varID].gridsize);
398 }
399 }
400 }
401 }
402
403 int tsmax = -1;
404 if (operatorID == SELTIMESTEP)
405 for (int i = 0; i < nsel; i++) tsmax = std::max(tsmax, intarr[i]);
406
407 int indexNext = 0;
408 int tsID = 0;
409 int tsID2 = 0;
410 while (true)
411 {
412 const auto nrecs = cdo_stream_inq_timestep(streamID1, tsID);
413 if (nrecs == 0) break;
414
415 const auto vdate = taxisInqVdate(taxisID1);
416 const auto vtime = taxisInqVtime(taxisID1);
417
418 auto copytimestep = false;
419 int selival = -1;
420
421 if (operfunc == func_step)
422 {
423 selival = tsID + 1;
424 if (selival > iselmax) break;
425 }
426 else if (operfunc == func_date)
427 {
428 int year, month, day;
429 cdiDecodeDate(vdate, &year, &month, &day);
430 selival = (operatorID == SELYEAR) ? year : (operatorID == SELDAY) ? day : month;
431 }
432 else if (operfunc == func_time)
433 {
434 int hour, minute, second;
435 cdiDecodeTime(vtime, &hour, &minute, &second);
436 selival = (operatorID == SELHOUR) ? hour : vtime;
437 }
438 else if (operfunc == func_datetime)
439 {
440 selfval = vdate + vtime / 1000000.;
441 }
442
443 if (operatorID == SELDATE)
444 {
445 if (selfval >= fltarr[0] && selfval <= fltarr[nsel - 1])
446 {
447 copytimestep = true;
448 selfound[0] = true;
449 selfound[nsel - 1] = true;
450 }
451 else if (selfval > fltarr[nsel - 1])
452 {
453 break;
454 }
455 }
456 else if (operatorID == SELTIMESTEP)
457 {
458 if (tsID >= tsmax) break;
459
460 int index;
461 for (index = indexNext; index < nsel; index++)
462 if (selival == intarr[index]) break;
463 if (index < nsel)
464 {
465 copytimestep = true;
466 selfound[index] = true;
467 indexNext = index + 1;
468 }
469 }
470 else
471 {
472 for (int i = 0; i < nsel; i++)
473 if (selival == intarr[i])
474 {
475 copytimestep = true;
476 selfound[i] = true;
477 break;
478 }
479 }
480
481 auto copy_nts2 = false;
482 if (operatorID == SELSMON && !copytimestep)
483 {
484 copy_nts2 = false;
485
486 if (process_nts1)
487 {
488 process_nts2 = true;
489 its2 = 0;
490 process_nts1 = false;
491 }
492
493 if (process_nts2)
494 {
495 if (its2++ < nts2)
496 copy_nts2 = true;
497 else
498 process_nts2 = false;
499 }
500 }
501
502 if (copytimestep || copy_nts2)
503 {
504 cdo_taxis_copy_timestep(taxisID2, taxisID1);
505 if (tsID2 == 0)
506 {
507 streamID2 = cdo_open_write(1);
508 cdo_def_vlist(streamID2, vlistID2);
509 }
510
511 if (lnts1 && ncts == 0)
512 {
513 nts = nts1;
514 if (its1 < nts1)
515 {
516 nts = its1;
517 cdo_warning("%d timesteps missing before month %d!", nts1 - its1, intarr[0]);
518 }
519
520 for (int it = 0; it < nts; it++)
521 {
522 taxisDefVdate(taxisID2, vdate_list[it]);
523 taxisDefVtime(taxisID2, vtime_list[it]);
524 cdo_def_timestep(streamID2, tsID2++);
525
526 for (varID = 0; varID < nvars; varID++)
527 {
528 if (varList1[varID].timetype == TIME_CONSTANT && tsID2 > 1) continue;
529 const auto nlevels = varList1[varID].nlevels;
530 for (levelID = 0; levelID < nlevels; levelID++)
531 {
532 cdo_def_record(streamID2, varID, levelID);
533 auto single = vars[it][varID][levelID].vec_d.data();
534 const auto nmiss = vars[it][varID][levelID].nmiss;
535 cdo_write_record(streamID2, single, nmiss);
536 }
537 }
538 }
539
540 its1 = 0;
541 }
542
543 ncts++;
544 if (!process_nts2)
545 {
546 its2 = 0;
547 process_nts1 = true;
548 }
549
550 cdo_def_timestep(streamID2, tsID2++);
551
552 if (tsID > 0 && lconstout)
553 {
554 lconstout = false;
555 nts = nts1 - 1;
556 for (varID = 0; varID < nvars; varID++)
557 {
558 if (varList1[varID].timetype == TIME_CONSTANT)
559 {
560 const auto nlevel = varList1[varID].nlevels;
561 for (levelID = 0; levelID < nlevel; levelID++)
562 {
563 cdo_def_record(streamID2, varID, levelID);
564 auto single = vars[nts][varID][levelID].vec_d.data();
565 const auto nmiss = vars[nts][varID][levelID].nmiss;
566 cdo_write_record(streamID2, single, nmiss);
567 }
568 }
569 }
570 }
571
572 for (int recID = 0; recID < nrecs; recID++)
573 {
574 cdo_inq_record(streamID1, &varID, &levelID);
575 cdo_def_record(streamID2, varID, levelID);
576 if (dataIsUnchanged)
577 {
578 cdo_copy_record(streamID2, streamID1);
579 }
580 else
581 {
582 field.init(varList1[varID]);
583 cdo_read_record(streamID1, field);
584 cdo_write_record(streamID2, field);
585 }
586 }
587 }
588 else
589 {
590 ncts = 0;
591
592 if (lnts1 || tsID == 0)
593 {
594 if (tsID == 0 && nconst && (!lnts1)) lconstout = true;
595
596 nts = nts1 - 1;
597 if (lnts1)
598 {
599 if (its1 <= nts)
600 nts = its1;
601 else
602 for (int it = 0; it < nts; it++)
603 {
604 vdate_list[it] = vdate_list[it + 1];
605 vtime_list[it] = vtime_list[it + 1];
606 for (varID = 0; varID < nvars; varID++)
607 {
608 if (varList1[varID].timetype == TIME_CONSTANT) continue;
609 const auto nlevels = varList1[varID].nlevels;
610 for (levelID = 0; levelID < nlevels; levelID++)
611 {
612 vars[it][varID][levelID].vec_d = vars[it + 1][varID][levelID].vec_d;
613 vars[it][varID][levelID].nmiss = vars[it + 1][varID][levelID].nmiss;
614 }
615 }
616 }
617
618 vdate_list[nts] = taxisInqVdate(taxisID1);
619 vtime_list[nts] = taxisInqVtime(taxisID1);
620
621 its1++;
622 }
623
624 for (int recID = 0; recID < nrecs; recID++)
625 {
626 cdo_inq_record(streamID1, &varID, &levelID);
627 if (lnts1 || varList1[varID].timetype == TIME_CONSTANT)
628 {
629 auto single = vars[nts][varID][levelID].vec_d.data();
630 cdo_read_record(streamID1, single, &vars[nts][varID][levelID].nmiss);
631 }
632 }
633 }
634 }
635
636 tsID++;
637 }
638
639 if (streamID2 != CDO_STREAM_UNDEF) cdo_stream_close(streamID2);
640 cdo_stream_close(streamID1);
641
642 if (operatorID == SELSMON)
643 if (its2 < nts2) cdo_warning("%d timesteps missing after the last month!", nts2 - its2);
644
645 for (int isel = 0; isel < nsel; isel++)
646 {
647 if (selfound[isel] == false)
648 {
649 if (operatorID == SELTIMESTEP)
650 {
651 int isel2;
652 auto lcont = false;
653 for (isel2 = isel + 1; isel2 < nsel; isel2++)
654 if (selfound[isel2]) break;
655 if (isel2 == nsel && (nsel - isel) > 1) lcont = true;
656
657 auto lcont2 = false;
658 if (lcont)
659 {
660 for (isel2 = isel + 1; isel2 < nsel; isel2++)
661 if (intarr[isel2 - 1] != intarr[isel2] - 1) break;
662 if (isel2 == nsel) lcont2 = true;
663 }
664
665 if (lcont2)
666 {
667 cdo_warning("Timesteps %d-%d not found!", intarr[isel], intarr[nsel - 1]);
668 break;
669 }
670 else
671 cdo_warning("Timestep %d not found!", intarr[isel]);
672 }
673 else if (operatorID == SELDATE)
674 {
675 if (isel == 0) cdo_warning("Date between %14.4f and %14.4f not found!", fltarr[0], fltarr[nsel - 1]);
676 }
677 else if (operatorID == SELTIME)
678 {
679 cdo_warning("Time %d not found!", intarr[isel]);
680 }
681 else if (operatorID == SELHOUR)
682 {
683 cdo_warning("Hour %d not found!", intarr[isel]);
684 }
685 else if (operatorID == SELDAY)
686 {
687 cdo_warning("Day %d not found!", intarr[isel]);
688 }
689 else if (operatorID == SELMONTH)
690 {
691 cdo_warning("Month %d not found!", intarr[isel]);
692 }
693 else if (operatorID == SELYEAR)
694 {
695 cdo_warning("Year %d not found!", intarr[isel]);
696 }
697 else if (operatorID == SELSEASON)
698 {
699 if (isel < 3) cdo_warning("Month %d not found!", intarr[isel]);
700 }
701 }
702 }
703
704 vlistDestroy(vlistID2);
705
706 if (tsID2 == 0) cdo_abort("No timesteps selected!");
707
708 cdo_finish();
709
710 return nullptr;
711 }
712