field.cc 11.8 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
/*
  This file is part of CDO. CDO is a collection of Operators to
  manipulate and analyse Climate model Data.

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14
15
16
17
  See COPYING file for copying and redistribution conditions.

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; version 2 of the License.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
*/

18
19
#include <cdi.h>

20
#include <cstdlib>  // qsort
Oliver Heidmann's avatar
Oliver Heidmann committed
21
#include <cassert>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
22
#include <algorithm>
23

24
#include "percentiles.h"
25
#include "array.h"
Oliver Heidmann's avatar
Oliver Heidmann committed
26
27
28
#include "field.h"
#include "functs.h"
#include "cdo_output.h"
29

30
Field::Field() : fpeRaised(0), nwpv(1), memType(MemType::Native), grid(-1), gridsize(0), size(0), nsamp(0), nmiss(0), missval(0) { m_count = 0; }
31
Field3D::Field3D() : nlevels(0) {}
32
33

void
34
Field::init(const CdoVar &var)
35
{
36
37
38
39
40
  grid = var.gridID;
  gridsize = var.gridsize;
  missval = var.missval;
  memType = var.memType;
  size = var.gridsize * var.nwpv;
41
42
43
44
45
  if (memType == MemType::Float)
    vec_f.resize(size);
  else
    vec_d.resize(size);
}
46
47
48
49
50

void
Field::resize(size_t count)
{
  m_count = count;
51
  vec_d.resize(m_count);
52
53
54
55
56
57
58
  if (!size) size = m_count;
}

void
Field::resize(size_t count, double value)
{
  m_count = count;
59
  vec_d.resize(m_count, value);
60
61
62
63
64
65
66
  if (!size) size = m_count;
}

void
Field::resizef(size_t count)
{
  m_count = count;
67
  vec_f.resize(m_count);
68
69
70
71
72
73
74
  if (!size) size = m_count;
}

void
Field::resizef(size_t count, float value)
{
  m_count = count;
75
  vec_f.resize(m_count, value);
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
  if (!size) size = m_count;
}

bool
Field::empty() const
{
  return m_count == 0;
}

void
Field::check_gridsize() const
{
  if (size == 0) fprintf(stderr, "Internal problem, size of field not set!");
  if (size > m_count) fprintf(stderr, "Internal problem, size of field is greater than allocated size of field!");
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91

92
void
93
Field3D::init(const CdoVar &var)
94
{
95
96
97
98
99
100
  nlevels = var.nlevels;
  grid = var.gridID;
  gridsize = var.gridsize;
  missval = var.missval;
  memType = var.memType;
  size = var.nlevels * var.gridsize * var.nwpv;
101
102
103
104
105
106
  if (memType == MemType::Float)
    vec_f.resize(size);
  else
    vec_d.resize(size);
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
108
void
fieldFill(Field &field, double value)
109
110
111
{
  field.check_gridsize();

112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
  if (field.memType == MemType::Float)
    std::fill(field.vec_f.begin(), field.vec_f.begin() + field.size, value);
  else
    std::fill(field.vec_d.begin(), field.vec_d.begin() + field.size, value);
}

void fieldCopy(const Field &field_src, Field &field_tgt)
{
  if (field_src.memType == MemType::Float)
    field_tgt.vec_f = field_src.vec_f;
  else
    field_tgt.vec_d = field_src.vec_d;
}

void fieldCopy(const Field3D &field_src, int levelID, Field &field_tgt)
{
  const auto size = field_src.gridsize * field_src.nwpv;
  const auto offset = levelID * size;
  if (field_src.memType == MemType::Float)
    std::copy(field_src.vec_f.begin() + offset, field_src.vec_f.begin() + offset + size, field_tgt.vec_f.begin());
  else
    std::copy(field_src.vec_d.begin() + offset, field_src.vec_d.begin() + offset + size, field_tgt.vec_d.begin());
134
135
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
137
138
139
140
141
142
143
144
145
146
void fieldAdd(Field &field1, const Field3D &field2, int levelID)
{
  const auto size = field1.gridsize * field1.nwpv;
  const auto offset = levelID * size;
  if (field1.memType == MemType::Float)
    for (size_t i = 0; i < size; i++) field1.vec_f[i] += field2.vec_f[offset + i];
  else
    for (size_t i = 0; i < size; i++) field1.vec_d[i] += field2.vec_d[offset + i];

}

147
// functor that returns true if value is equal to the value of the constructor parameter provided
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
149
class valueDblIsEqual
{
150
  double _missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151

152
public:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
155
156
  valueDblIsEqual(double missval) : _missval(missval) {}
  bool
  operator()(double value) const
  {
157
158
159
160
    return DBL_IS_EQUAL(value, _missval);
  }
};

Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
// functor that returns true if value is equal to the value of the constructor parameter provided
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
163
class valueIsEqual
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
  double _missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165

Uwe Schulzweida's avatar
Uwe Schulzweida committed
166
public:
Uwe Schulzweida's avatar
Uwe Schulzweida committed
167
168
169
170
  valueIsEqual(double missval) : _missval(missval) {}
  bool
  operator()(double value) const
  {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
171
172
173
174
    return IS_EQUAL(value, _missval);
  }
};

Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
176
size_t
fieldNumMiss(const Field &field)
177
178
179
{
  field.check_gridsize();

180
  if (std::isnan(field.missval))
181
    {
182
      return std::count_if(field.vec_d.begin(), field.vec_d.begin() + field.size, valueDblIsEqual(field.missval));
183
184
185
    }
  else
    {
186
      return std::count_if(field.vec_d.begin(), field.vec_d.begin() + field.size, valueIsEqual(field.missval));
187
188
189
    }
}

190
191
192
193
size_t
fieldNumMV(Field &field)
{
  if (field.memType == MemType::Float)
194
    field.nmiss = varrayNumMV(field.size, field.vec_f, (float)field.missval);
195
196
  else
    field.nmiss = varrayNumMV(field.size, field.vec_d, field.missval);
197
198
199
200

  return field.nmiss;
}

201
MinMaxVal
202
203
204
205
206
207
fieldMinMax(Field &field)
{
  if (field.memType == MemType::Float)
    return varrayMinMax(field.vec_f);
  else
    return varrayMinMax(field.vec_d);
208
209
}

210
static double
211
212
vfldrange(const Field &field)
{
213
  return field.nmiss ? varrayRangeMV(field.size, field.vec_d, field.missval) : varrayRange(field.size, field.vec_d);
214
215
}

216
217
218
double
vfldmin(const Field &field)
{
219
  return field.nmiss ? varrayMinMV(field.size, field.vec_d, field.missval) : varrayMin(field.size, field.vec_d);
220
221
222
223
224
}

double
vfldmax(const Field &field)
{
225
  return field.nmiss ? varrayMaxMV(field.size, field.vec_d, field.missval) : varrayMax(field.size, field.vec_d);
226
227
228
229
230
}

double
vfldsum(const Field &field)
{
231
  return field.nmiss ? varraySumMV(field.size, field.vec_d, field.missval) : varraySum(field.size, field.vec_d);
232
233
234
235
236
}

double
vfldmean(const Field &field)
{
237
  return field.nmiss ? varrayMeanMV(field.size, field.vec_d, field.missval) : varrayMean(field.size, field.vec_d);
238
239
240
241
242
}

double
vfldmeanw(const Field &field)
{
243
244
  return field.nmiss ? varrayWeightedMeanMV(field.size, field.vec_d, field.weightv, field.missval)
                     : varrayWeightedMean(field.size, field.vec_d, field.weightv, field.missval);
245
246
247
248
249
}

double
vfldavg(const Field &field)
{
250
  return field.nmiss ? varrayAvgMV(field.size, field.vec_d, field.missval) : varrayMean(field.size, field.vec_d);
251
252
253
254
255
}

double
vfldavgw(const Field &field)
{
256
257
  return field.nmiss ? varrayWeightedAvgMV(field.size, field.vec_d, field.weightv, field.missval)
                     : varrayWeightedMean(field.size, field.vec_d, field.weightv, field.missval);
258
259
260
261
262
}

double
vfldvar(const Field &field)
{
263
  return varrayVar(field.size, field.vec_d, field.nmiss, field.missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
}
265

266
267
268
double
vfldvar1(const Field &field)
{
269
  return varrayVar1(field.size, field.vec_d, field.nmiss, field.missval);
270
271
272
273
274
}

double
vfldkurt(const Field &field)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
  return varrayKurt(field.size, field.vec_d, field.nmiss, field.missval);
276
277
278
279
280
}

double
vfldskew(const Field &field)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
281
  return varraySkew(field.size, field.vec_d, field.nmiss, field.missval);
282
283
284
285
286
}

double
vfldvarw(const Field &field)
{
287
  return varrayWeightedVar(field.size, field.vec_d, field.weightv, field.nmiss, field.missval);
288
289
290
291
292
}

double
vfldvar1w(const Field &field)
{
293
  return varrayWeightedVar1(field.size, field.vec_d, field.weightv, field.nmiss, field.missval);
294
295
}

296
297
double
varToStd(double rvar, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298
{
299
  if (DBL_IS_EQUAL(rvar, missval) || rvar < 0)
300
    {
301
      return missval;
302
303
304
    }
  else
    {
305
      return IS_NOT_EQUAL(rvar, 0) ? std::sqrt(rvar) : 0;
306
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
307
308
}

309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
double
vfldstd(const Field &field)
{
  return varToStd(vfldvar(field), field.missval);
}

double
vfldstd1(const Field &field)
{
  return varToStd(vfldvar1(field), field.missval);
}

double
vfldstdw(const Field &field)
{
  return varToStd(vfldvarw(field), field.missval);
}

double
vfldstd1w(const Field &field)
{
  return varToStd(vfldvar1w(field), field.missval);
}

333
void
334
vfldrms(const Field &field, const Field &field2, Field &field3)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
{
336
  size_t i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
  size_t rnmiss = 0;
338
  auto grid1 = field.grid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
  //  size_t nmiss1   = field.nmiss;
340
  const auto array1 = field.vec_d.data();
341
  auto grid2 = field2.grid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
  //  size_t nmiss2   = field2.nmiss;
343
  const auto array2 = field2.vec_d.data();
344
345
  const auto missval1 = field.missval;
  const auto missval2 = field2.missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
346
  const auto &w = field.weightv;
347
  double rsum = 0, rsumw = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
348

349
  const auto len = gridInqSize(grid1);
350
  if (len != gridInqSize(grid2)) cdoAbort("fields have different size!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
352

  /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
  if ( nmiss1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
354
  */
355
356
357
358
  {
    for (i = 0; i < len; i++)
      if (!DBL_IS_EQUAL(w[i], missval1))
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
          rsum = ADDMN(rsum, MULMN(w[i], MULMN(SUBMN(array2[i], array1[i]), SUBMN(array2[i], array1[i]))));
360
361
362
363
364
365
366
367
368
369
370
371
372
          rsumw = ADDMN(rsumw, w[i]);
        }
  }
  /*
else
  {
    for ( i = 0; i < len; i++ )
      {
        rsum  += w[i] * array1[i];
        rsumw += w[i];
      }
  }
  */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373

374
  const double ravg = SQRTMN(DIVMN(rsum, rsumw));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
375

376
  if (DBL_IS_EQUAL(ravg, missval1)) rnmiss++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
377

378
  field3.vec_d[0] = ravg;
379
  field3.nmiss = rnmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
381
}

382
383
384
double
vfldpctl(Field &field, const double pn)
{
385
  auto pctl = field.missval;
386

Uwe Schulzweida's avatar
Uwe Schulzweida committed
387
  if ((field.size - field.nmiss) > 0)
388
389
390
    {
      if (field.nmiss)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391
          Varray<double> v(field.size - field.nmiss);
392
393

          size_t j = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
394
          for (size_t i = 0; i < field.size; i++)
395
            if (!DBL_IS_EQUAL(field.vec_d[i], field.missval)) v[j++] = field.vec_d[i];
396
397
398
399
400

          pctl = percentile(v.data(), j, pn);
        }
      else
        {
401
          pctl = percentile(field.vec_d.data(), field.size, pn);
402
403
404
405
406
407
        }
    }

  return pctl;
}

408
409
410
411
412
413
414
415
static int
compareDouble(const void *a, const void *b)
{
  const double *x = (const double *) a;
  const double *y = (const double *) b;
  return ((*x > *y) - (*x < *y)) * 2 + (*x > *y) - (*x < *y);
}

416
double
417
vfldrank(Field &field)
418
419
420
{
  double res = 0;
  // Using first value as reference (observation)
421
  double *array = &field.vec_d[1];
422
423
  const auto val = array[-1];
  const auto len = field.size - 1;
424

425
  if (field.nmiss) return field.missval;
426

427
  std::qsort(array, len, sizeof(double), compareDouble);
428

429
430
  if (val > array[len - 1])
    res = (double) len;
431
  else
432
    for (size_t j = 0; j < len; j++)
433
434
435
436
437
      if (array[j] >= val)
        {
          res = (double) j;
          break;
        }
438
439
440
441

  return res;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
double fldbrs(Field &field)  // Not used!
443
{
444
445
  const auto nmiss = field.nmiss;
  const auto len = field.size;
446
  const auto array = field.vec_d.data();
447
  const auto missval = field.missval;
448
449

  double brs = 0;
450
  size_t i, count = 0;
451
452

  // Using first value as reference
453
  if (nmiss == 0)
454
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
      for (i = 1; i < len; i++) brs += (array[i] - array[0]) * (array[i] - array[0]);
456
      count = i - 1;
457
458
459
    }
  else
    {
460
461
462
463
464
465
466
467
      if (DBL_IS_EQUAL(array[0], missval)) return missval;

      for (i = 1; i < len; i++)
        if (!DBL_IS_EQUAL(array[i], missval))
          {
            brs += (array[i] - array[0]) * (array[i] - array[0]);
            count++;
          }
468
469
    }

470
  return brs / count;
471
472
}

473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
double
vfldfun(const Field &field, int function)
{
  double rval = 0;

  // clang-format off
  switch (function)
    {
    case func_range:  rval = vfldrange(field);  break;
    case func_min:    rval = vfldmin(field);    break;
    case func_max:    rval = vfldmax(field);    break;
    case func_sum:    rval = vfldsum(field);    break;
    case func_mean:   rval = vfldmean(field);   break;
    case func_avg:    rval = vfldavg(field);    break;
    case func_std:    rval = vfldstd(field);    break;
    case func_std1:   rval = vfldstd1(field);   break;
    case func_var:    rval = vfldvar(field);    break;
    case func_var1:   rval = vfldvar1(field);   break;
    case func_meanw:  rval = vfldmeanw(field);  break;
    case func_avgw:   rval = vfldavgw(field);   break;
    case func_stdw:   rval = vfldstdw(field);   break;
    case func_std1w:  rval = vfldstd1w(field);  break;
    case func_varw:   rval = vfldvarw(field);   break;
    case func_var1w:  rval = vfldvar1w(field);  break;
    case func_skew:   rval = vfldskew(field);   break;
    case func_kurt:   rval = vfldkurt(field);   break;
    default: cdoAbort("%s: function %d not implemented!", __func__, function);
    }
  // clang-format on

  return rval;
}