Timselstat.cc 9.97 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-2017 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
20
/*
   This module contains the following operators:

Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
22
23
24
25
26
27
28
29
30
      Timselstat    timselrange        Time selection range
      Timselstat    timselmin          Time selection minimum
      Timselstat    timselmax          Time selection maximum
      Timselstat    timselsum          Time selection sum
      Timselstat    timselmean         Time selection mean
      Timselstat    timselavg          Time selection average
      Timselstat    timselvar          Time selection variance
      Timselstat    timselvar1         Time selection variance [Normalize by (n-1)]
      Timselstat    timselstd          Time selection standard deviation
      Timselstat    timselstd1         Time selection standard deviation [Normalize by (n-1)]
31
32
33
*/


Ralf Mueller's avatar
Ralf Mueller committed
34
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
36
37
38
39
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"


40
41
42
43
44
45
typedef struct {
  short varID;
  short levelID;
} recinfo_t;


46
void *Timselstat(void *argument)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
{
48
  int timestat_date = TIMESTAT_MEAN;
49
  int nrecs = 0;
50
  int varID, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
52
53
54
55
56
  int tsID;
  int nsets;
  int nmiss;

  cdoInitialize(argument);

57
58
59
60
61
62
63
64
65
66
  cdoOperatorAdd("timselrange", func_range, 0, NULL);
  cdoOperatorAdd("timselmin",   func_min,   0, NULL);
  cdoOperatorAdd("timselmax",   func_max,   0, NULL);
  cdoOperatorAdd("timselsum",   func_sum,   0, NULL);
  cdoOperatorAdd("timselmean",  func_mean,  0, NULL);
  cdoOperatorAdd("timselavg",   func_avg,   0, NULL);
  cdoOperatorAdd("timselvar",   func_var,   0, NULL);
  cdoOperatorAdd("timselvar1",  func_var1,  0, NULL);
  cdoOperatorAdd("timselstd",   func_std,   0, NULL);
  cdoOperatorAdd("timselstd1",  func_std1,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67

68
69
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
71
72

  operatorInputArg("nsets <noffset <nskip>>");

73
  int nargc  = operatorArgc();
74
  int ndates = parameter2int(operatorArgv()[0]);
75
76
  int noffset = (nargc > 1) ? parameter2int(operatorArgv()[1]) : 0;
  int nskip   = (nargc > 2) ? parameter2int(operatorArgv()[2]) : 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
78
79

  if ( cdoVerbose ) cdoPrint("nsets = %d, noffset = %d, nskip = %d", ndates, noffset, nskip);

80
81
82
83
84
  bool lrange  = operfunc == func_range;
  bool lmean   = operfunc == func_mean || operfunc == func_avg;
  bool lstd    = operfunc == func_std || operfunc == func_std1;
  bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
  int  divisor = operfunc == func_std1 || operfunc == func_var1;
85

86
  int streamID1 = streamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
87

88
89
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90

91
92
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
93
94
  vlistDefTaxis(vlistID2, taxisID2);

95
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
97
98

  streamDefVlist(streamID2, vlistID2);

99
  int maxrecs = vlistNrecs(vlistID1);
100

101
  recinfo_t *recinfo = (recinfo_t *) Malloc(maxrecs*sizeof(recinfo_t));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102

103
104
105
  dtlist_type *dtlist = dtlist_new();
  dtlist_set_stat(dtlist, timestat_date);
  dtlist_set_calendar(dtlist, taxisInqCalendar(taxisID1));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106

107
  int gridsize = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108

109
  field_type field;
110
  field_init(&field);
111
  field.ptr = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112

113
114
115
  field_type **samp1 = field_malloc(vlistID1, FIELD_NONE);
  field_type **vars1 = field_malloc(vlistID1, FIELD_PTR);
  field_type **vars2 = NULL;
116
  if ( lvarstd || lrange ) vars2 = field_malloc(vlistID1, FIELD_PTR);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
118
119
120
121

  for ( tsID = 0; tsID < noffset; tsID++ )
    {
      nrecs = streamInqTimestep(streamID1, tsID);
      if ( nrecs == 0 ) break;
122

123
      for ( int recID = 0; recID < nrecs; recID++ )
124
125
126
127
128
	{
	  streamInqRecord(streamID1, &varID, &levelID);

	  if ( tsID == 0 )
	    {
129
130
              recinfo[recID].varID   = varID;
              recinfo[recID].levelID = levelID;
131
132
	    }
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
134
    }

135
  int otsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
137
  if ( tsID < noffset )
    {
138
      cdoWarning("noffset is larger than number of timesteps!");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
141
142
143
144
145
146
147
148
      goto LABEL_END;
    }

  while ( TRUE )
    {
      for ( nsets = 0; nsets < ndates; nsets++ )
	{
	  nrecs = streamInqTimestep(streamID1, tsID);
	  if ( nrecs == 0 ) break;

149
	  dtlist_taxisInqTimestep(dtlist, taxisID1, nsets);
150

151
	  for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
153
154
	    {
	      streamInqRecord(streamID1, &varID, &levelID);

155
156
	      if ( tsID == 0 )
		{
157
158
                  recinfo[recID].varID   = varID;
                  recinfo[recID].levelID = levelID;
159
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
160

161
162
163
              field_type *pvars1 = &vars1[varID][levelID];
              field_type *pvars2 = vars2 ? &vars2[varID][levelID] : NULL;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
165
166
167
	      gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));

	      if ( nsets == 0 )
		{
168
169
170
171
172
173
174
175
		  streamReadRecord(streamID1, pvars1->ptr, &nmiss);
		  pvars1->nmiss = (size_t)nmiss;
                  if ( lrange )
                    {
                      pvars2->nmiss = (size_t)nmiss;
		      for ( int i = 0; i < gridsize; i++ )
                        pvars2->ptr[i] = pvars1->ptr[i];
                    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
176
177
178
179

		  if ( nmiss > 0 || samp1[varID][levelID].ptr )
		    {
		      if ( samp1[varID][levelID].ptr == NULL )
180
			samp1[varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181

182
		      for ( int i = 0; i < gridsize; i++ )
183
                        samp1[varID][levelID].ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
185
186
187
		    }
		}
	      else
		{
188
189
		  streamReadRecord(streamID1, field.ptr, &nmiss);
                  field.nmiss   = (size_t)nmiss;
190
191
		  field.grid    = pvars1->grid;
		  field.missval = pvars1->missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192
193
194
195
196

		  if ( field.nmiss > 0 || samp1[varID][levelID].ptr )
		    {
		      if ( samp1[varID][levelID].ptr == NULL )
			{
197
			  samp1[varID][levelID].ptr = (double*) Malloc(gridsize*sizeof(double));
198
			  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
200
201
			    samp1[varID][levelID].ptr[i] = nsets;
			}

202
		      for ( int i = 0; i < gridsize; i++ )
203
			if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
204
205
206
			  samp1[varID][levelID].ptr[i]++;
		    }

207
		  if ( lvarstd )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
208
		    {
209
210
		      farsumq(pvars2, field);
		      farsum(pvars1, field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
		    }
212
213
214
215
216
                  else if ( lrange )
                    {
                      farmin(pvars2, field);
                      farmax(pvars1, field);
                    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
218
		  else
		    {
219
		      farfun(pvars1, field, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
220
221
222
223
		    }
		}
	    }

224
	  if ( nsets == 0 && lvarstd )
225
226
227
228
229
230
231
            for ( int recID = 0; recID < maxrecs; recID++ )
              {
                int varID   = recinfo[recID].varID;
                int levelID = recinfo[recID].levelID;
                field_type *pvars1 = &vars1[varID][levelID];
                field_type *pvars2 = &vars2[varID][levelID];

Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
		if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
233
234

                farmoq(pvars2, *pvars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
236
237
238
239
240
241
	      }

	  tsID++;
	}

      if ( nrecs == 0 && nsets == 0 ) break;

242
      if ( lmean )
243
244
245
246
247
248
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
            field_type *pvars1 = &vars1[varID][levelID];

Uwe Schulzweida's avatar
Uwe Schulzweida committed
249
	    if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
250
251
252
253
254

            if ( samp1[varID][levelID].ptr == NULL )
              farcdiv(pvars1, (double)nsets);
            else
              fardiv(pvars1, samp1[varID][levelID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
	  }
256
      else if ( lvarstd )
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
            field_type *pvars1 = &vars1[varID][levelID];
            field_type *pvars2 = &vars2[varID][levelID];

            if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;

            if ( samp1[varID][levelID].ptr == NULL )
              {
                if ( lstd ) farcstd(pvars1, *pvars2, nsets, divisor);
                else        farcvar(pvars1, *pvars2, nsets, divisor);
              }
            else
              {
                if ( lstd ) farstd(pvars1, *pvars2, samp1[varID][levelID], divisor);
                else        farvar(pvars1, *pvars2, samp1[varID][levelID], divisor);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
275
276
	      }
	  }
277
278
279
280
281
282
283
284
285
286
287
288
      else if ( lrange )
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
            field_type *pvars1 = &vars1[varID][levelID];
            field_type *pvars2 = &vars2[varID][levelID];

            if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;

            farsub(pvars1, *pvars2);
	  }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289

290
      dtlist_stat_taxisDefTimestep(dtlist, taxisID2, nsets);
291
      streamDefTimestep(streamID2, otsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
292

293
      for ( int recID = 0; recID < maxrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
	{
295
296
297
          int varID   = recinfo[recID].varID;
          int levelID = recinfo[recID].levelID;
          field_type *pvars1 = &vars1[varID][levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298

Uwe Schulzweida's avatar
Uwe Schulzweida committed
299
	  if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
300
301

	  streamDefRecord(streamID2, varID, levelID);
302
	  streamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
303
304
305
	}

      if ( nrecs == 0 ) break;
306
      otsID++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
307

308
      for ( int i = 0; i < nskip; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309
310
311
312
313
314
315
316
317
318
319
320
	{
	  nrecs = streamInqTimestep(streamID1, tsID);
	  if ( nrecs == 0 ) break;
	  tsID++;
	}

      if ( nrecs == 0 ) break;
    }

 LABEL_END:


321
322
323
  field_free(vars1, vlistID1);
  field_free(samp1, vlistID1);
  if ( lvarstd ) field_free(vars2, vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324

325
326
  dtlist_delete(dtlist);

327
  Free(recinfo);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328

329
  if ( field.ptr ) Free(field.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
330
331
332
333
334
335

  streamClose(streamID2);
  streamClose(streamID1);

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
}