Timstat2.cc 8.36 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
18
19
20
  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.
*/

/*
   This module contains the following operators:

Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
        Timstat2        timcor      correlates two data files on the same grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
22
23
*/

Ralf Müller's avatar
Ralf Müller committed
24
#include <cdi.h>
Oliver Heidmann's avatar
Oliver Heidmann committed
25

26
#include "functs.h"
27
#include "process_int.h"
28
#include "cdo_vlist.h"
29

30
// correlation in time
31
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
33
correlationInit(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals,
                double *work0, double *work1, double *work2, double *work3, double *work4)
34
{
35
  for (size_t i = 0; i < gridsize; ++i)
36
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
      if ((!DBL_IS_EQUAL(array1[i], missval1)) && (!DBL_IS_EQUAL(array2[i], missval2)))
38
39
40
        {
          work0[i] += array1[i];
          work1[i] += array2[i];
41
42
43
          work2[i] += array1[i] * array1[i];
          work3[i] += array2[i] * array2[i];
          work4[i] += array1[i] * array2[i];
44
45
          nofvals[i]++;
        }
46
    }
47
48
}

49
static size_t
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
correlation(size_t gridsize, double missval1, double missval2, size_t *nofvals, double *work0, double *work1, double *work2,
            double *work3, double *work4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
{
53
  size_t nmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54

55
56
  for (size_t i = 0; i < gridsize; ++i)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
      double cor;
58
      const auto nvals = nofvals[i];
59
60
      if (nvals > 0)
        {
61
62
63
64
65
66
67
          const auto temp0 = MULMN(work0[i], work1[i]);
          const auto temp1 = SUBMN(work4[i], DIVMN(temp0, nvals));
          const auto temp2 = MULMN(work0[i], work0[i]);
          const auto temp3 = MULMN(work1[i], work1[i]);
          const auto temp4 = SUBMN(work2[i], DIVMN(temp2, nvals));
          const auto temp5 = SUBMN(work3[i], DIVMN(temp3, nvals));
          const auto temp6 = MULMN(temp4, temp5);
68
69

          cor = DIVMN(temp1, SQRTMN(temp6));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
          cor = std::min(std::max(cor, -1.0), 1.0);
71
72
73

          if (DBL_IS_EQUAL(cor, missval1)) nmiss++;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
      else
75
76
77
78
        {
          nmiss++;
          cor = missval1;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79
80

      work0[i] = cor;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
82
83
84
85
    }

  return nmiss;
}

86
// covariance in time
87
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
89
covarianceInit(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals,
               double *work0, double *work1, double *work2)
90
{
91
  for (size_t i = 0; i < gridsize; ++i)
92
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
      if ((!DBL_IS_EQUAL(array1[i], missval1)) && (!DBL_IS_EQUAL(array2[i], missval2)))
94
95
96
        {
          work0[i] += array1[i];
          work1[i] += array2[i];
97
          work2[i] += array1[i] * array2[i];
98
99
          nofvals[i]++;
        }
100
    }
101
102
}

103
static size_t
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
covariance(size_t gridsize, double missval1, double missval2, size_t *nofvals, double *work0, double *work1, double *work2)
105
{
106
  size_t nmiss = 0;
107

108
109
  for (size_t i = 0; i < gridsize; ++i)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
110
      double covar;
111
      const auto nvals = nofvals[i];
112
113
      if (nvals > 0)
        {
114
          double dnvals = nvals;
115
          const auto temp = DIVMN(MULMN(work0[i], work1[i]), dnvals * dnvals);
116
117
118
          covar = SUBMN(DIVMN(work2[i], dnvals), temp);
          if (DBL_IS_EQUAL(covar, missval1)) nmiss++;
        }
119
      else
120
121
122
123
        {
          nmiss++;
          covar = missval1;
        }
124
125
126
127
128
129
130

      work0[i] = covar;
    }

  return nmiss;
}

131
132
void *
Timstat2(void *process)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
{
134
135
  int64_t vdate = 0;
  int vtime = 0;
136
  int varID, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
  size_t nmiss = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138

139
  cdoInitialize(process);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
140

Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142
143
  cdoOperatorAdd("timcor",   func_cor,   5, nullptr);
  cdoOperatorAdd("timcovar", func_covar, 3, nullptr);
144
  cdoOperatorAdd("timrms",   func_rms,   1, nullptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146

147
148
149
  const auto operatorID = cdoOperatorID();
  const auto operfunc = cdoOperatorF1(operatorID);
  const auto nwork = cdoOperatorF2(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
150

151
152
  operatorCheckArgc(0);

153
154
  const auto streamID1 = cdoOpenRead(0);
  const auto streamID2 = cdoOpenRead(1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
155

156
157
158
  const auto vlistID1 = cdoStreamInqVlist(streamID1);
  const auto vlistID2 = cdoStreamInqVlist(streamID2);
  const auto vlistID3 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159

Uwe Schulzweida's avatar
Uwe Schulzweida committed
160
  vlistCompare(vlistID1, vlistID2, CMP_ALL);
161

Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
163
164
165
  VarList varList1, varList2;
  varListInit(varList1, vlistID1);
  varListInit(varList2, vlistID2);

166
167
  const auto nvars = vlistNvars(vlistID1);
  auto nrecs = vlistNrecs(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
  const int nrecs3 = nrecs;
169
  std::vector<int> recVarID(nrecs), recLevelID(nrecs);
170

171
172
173
  const auto taxisID1 = vlistInqTaxis(vlistID1);
  // const auto taxisID2 = vlistInqTaxis(vlistID2);
  const auto taxisID3 = taxisDuplicate(taxisID1);
174

Uwe Schulzweida's avatar
Uwe Schulzweida committed
175
  vlistDefTaxis(vlistID3, taxisID3);
176
  const auto streamID3 = cdoOpenWrite(2);
177
  cdoDefVlist(streamID3, vlistID3);
178

179
  const auto gridsizemax = vlistGridsizeMax(vlistID1);
180
  Varray<double> array1(gridsizemax), array2(gridsizemax);
181

182
  Varray4D<double> work(nvars);
183
  Varray3D<size_t> nofvals(nvars);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184

185
  for (varID = 0; varID < nvars; varID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
    {
187
188
      const auto gridsize = varList1[varID].gridsize;
      const auto nlevs = varList1[varID].nlevels;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189

190
      work[varID].resize(nlevs);
191
      nofvals[varID].resize(nlevs);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192

193
194
      for (levelID = 0; levelID < nlevs; levelID++)
        {
195
          nofvals[varID][levelID].resize(gridsize, 0);
196
          work[varID][levelID].resize(nwork);
197
          for (int i = 0; i < nwork; i++) work[varID][levelID][i].resize(gridsize, 0.0);
198
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
    }
200

201
  int tsID = 0;
202
  while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
204
205
206
    {
      vdate = taxisInqVdate(taxisID1);
      vtime = taxisInqVtime(taxisID1);

207
      auto nrecs2 = cdoStreamInqTimestep(streamID2, tsID);
208
      if (nrecs != nrecs2) cdoWarning("Input streams have different number of records!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209

210
211
      for (int recID = 0; recID < nrecs; recID++)
        {
212
213
          cdoInqRecord(streamID1, &varID, &levelID);
          cdoInqRecord(streamID2, &varID, &levelID);
214
215
216
217
218
219
220

          if (tsID == 0)
            {
              recVarID[recID] = varID;
              recLevelID[recID] = levelID;
            }

221
222
223
          const auto gridsize = varList1[varID].gridsize;
          const auto missval1 = varList1[varID].missval;
          const auto missval2 = varList2[varID].missval;
224

225
226
          cdoReadRecord(streamID1, &array1[0], &nmiss);
          cdoReadRecord(streamID2, &array2[0], &nmiss);
227

Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
229
230
          auto &rwork = work[varID][levelID];
          auto &rnofvals = nofvals[varID][levelID];

231
232
          if (operfunc == func_cor)
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
234
              correlationInit(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(),
                              rwork[0].data(), rwork[1].data(), rwork[2].data(), rwork[3].data(), rwork[4].data());
235
236
237
            }
          else if (operfunc == func_covar)
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
239
              covarianceInit(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(),
                             rwork[0].data(), rwork[1].data(), rwork[2].data());
240
241
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242

Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
244
245
246
247
248
      tsID++;
    }

  tsID = 0;
  taxisDefVdate(taxisID3, vdate);
  taxisDefVtime(taxisID3, vtime);
249
  cdoDefTimestep(streamID3, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250

251
  for (int recID = 0; recID < nrecs3; recID++)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
    {
253
      varID = recVarID[recID];
254
      levelID = recLevelID[recID];
255

256
257
258
      const auto gridsize = varList1[varID].gridsize;
      const auto missval1 = varList1[varID].missval;
      const auto missval2 = varList2[varID].missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259

Uwe Schulzweida's avatar
Uwe Schulzweida committed
260
261
262
      auto &rwork = work[varID][levelID];
      auto &rnofvals = nofvals[varID][levelID];

263
264
      if (operfunc == func_cor)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
265
266
          nmiss = correlation(gridsize, missval1, missval2, rnofvals.data(), rwork[0].data(),
                              rwork[1].data(), rwork[2].data(), rwork[3].data(), rwork[4].data());
267
268
269
        }
      else if (operfunc == func_covar)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
270
271
          nmiss = covariance(gridsize, missval1, missval2, rnofvals.data(), rwork[0].data(),
                             rwork[1].data(), rwork[2].data());
272
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
273

274
      cdoDefRecord(streamID3, varID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
      cdoWriteRecord(streamID3, rwork[0].data(), nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
276
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277

278
279
280
  cdoStreamClose(streamID3);
  cdoStreamClose(streamID2);
  cdoStreamClose(streamID1);
281
282
283

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
284
  return nullptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
285
}