Yseasstat.cc 9.29 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.
*/

Uwe Schulzweida's avatar
Uwe Schulzweida committed
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
      Yseasstat  yseasrange      Multi-year seasonal range
      Yseasstat  yseasmin        Multi-year seasonal minimum
      Yseasstat  yseasmax        Multi-year seasonal maximum
      Yseasstat  yseassum        Multi-year seasonal sum
      Yseasstat  yseasmean       Multi-year seasonal mean
      Yseasstat  yseasavg        Multi-year seasonal average
      Yseasstat  yseasvar        Multi-year seasonal variance
      Yseasstat  yseasvar1       Multi-year seasonal variance [Normalize by (n-1)]
      Yseasstat  yseasstd        Multi-year seasonal standard deviation
      Yseasstat  yseasstd1       Multi-year seasonal standard deviation [Normalize by (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
32
*/

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


#define  NSEAS       4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
typedef struct {
43
44
  int vdate;
  int vtime;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
46
47
}
date_time_t;

48
49
50
51
52
53
typedef struct {
  short varID;
  short levelID;
} recinfo_t;


Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
55
void set_date(int vdate_new, int vtime_new, date_time_t *datetime)
{
56
57
58
59
60
61
62
63
64
  int year, month, day;
  cdiDecodeDate(vdate_new, &year, &month, &day);
  if ( month == 12 ) vdate_new = cdiEncodeDate(year-1, month, day);

  if ( vdate_new > datetime->vdate )
    {
      datetime->vdate = vdate_new;
      datetime->vtime = vtime_new;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65
66
67
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
69
70
void *Yseasstat(void *argument)
{
  int varID;
71
  int year, month, day;
72
  int nrecs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
  int levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
  int seas_nsets[NSEAS];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
  int nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
76
  date_time_t datetime[NSEAS];
77
  field_type **vars1[NSEAS], **vars2[NSEAS], **samp1[NSEAS];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
79
80

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
82
83
84
85
86
87
88
89
90
  cdoOperatorAdd("yseasrange", func_range, 0, NULL);
  cdoOperatorAdd("yseasmin",   func_min,   0, NULL);
  cdoOperatorAdd("yseasmax",   func_max,   0, NULL);
  cdoOperatorAdd("yseassum",   func_sum,   0, NULL);
  cdoOperatorAdd("yseasmean",  func_mean,  0, NULL);
  cdoOperatorAdd("yseasavg",   func_avg,   0, NULL);
  cdoOperatorAdd("yseasvar",   func_var,   0, NULL);
  cdoOperatorAdd("yseasvar1",  func_var1,  0, NULL);
  cdoOperatorAdd("yseasstd",   func_std,   0, NULL);
  cdoOperatorAdd("yseasstd1",  func_std1,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91

92
93
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
94

95
  for ( int seas = 0; seas < NSEAS; seas++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
98
99
      vars1[seas]  = NULL;
      vars2[seas]  = NULL;
      samp1[seas]  = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
      seas_nsets[seas]  = 0;
101
102
      datetime[seas].vdate = 0;
      datetime[seas].vtime = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
104
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
105
  bool lrange  = operfunc == func_range;
106
  bool lmean   = operfunc == func_mean || operfunc == func_avg;
107
108
  bool lstd    = operfunc == func_std || operfunc == func_std1;
  bool lvarstd = operfunc == func_std || operfunc == func_var || operfunc == func_std1 || operfunc == func_var1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
  int  divisor = operfunc == func_std1 || operfunc == func_var1;
110
111

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

113
114
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115

116
117
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
118
  if ( taxisHasBounds(taxisID2) ) taxisDeleteBounds(taxisID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119
120
  vlistDefTaxis(vlistID2, taxisID2);

121
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
122
123
124

  streamDefVlist(streamID2, vlistID2);

125
  int maxrecs = vlistNrecs(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
  int gridsizemax = vlistGridsizeMax(vlistID1);
130

131
  field_type field;
132
  field_init(&field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
133
  field.ptr = (double*) Malloc(gridsizemax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134

135
136
  int tsID = 0;
  int otsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
138
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
139
140
      int vdate = taxisInqVdate(taxisID1);
      int vtime = taxisInqVtime(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
      cdiDecodeDate(vdate, &year, &month, &day);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142

143
      int seas = month_to_season(month);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144

Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
      set_date(vdate, vtime, &datetime[seas]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146
147
148

      if ( vars1[seas] == NULL )
	{
149
150
	  vars1[seas] = field_malloc(vlistID1, FIELD_PTR);
	  samp1[seas] = field_malloc(vlistID1, FIELD_NONE);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
	  if ( lvarstd || lrange )
152
	    vars2[seas] = field_malloc(vlistID1, FIELD_PTR);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
	}

155
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156
157
158
	{
	  streamInqRecord(streamID1, &varID, &levelID);

159
160
	  if ( tsID == 0 )
	    {
161
162
              recinfo[recID].varID   = varID;
              recinfo[recID].levelID = levelID;
163
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164

Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
          field_type *psamp1 = &samp1[seas][varID][levelID];
166
167
          field_type *pvars1 = &vars1[seas][varID][levelID];
          field_type *pvars2 = vars2[seas] ? &vars2[seas][varID][levelID] : NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
          int nsets = seas_nsets[seas];
169

Uwe Schulzweida's avatar
Uwe Schulzweida committed
170
	  int gridsize = pvars1->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
171

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
	      if ( nmiss > 0 || psamp1->ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
186
		  if ( psamp1->ptr == NULL )
		    psamp1->ptr = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187

188
		  for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
                    psamp1->ptr[i] = !DBL_IS_EQUAL(pvars1->ptr[i], pvars1->missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
191
192
193
	    }
	  else
	    {
194
195
	      streamReadRecord(streamID1, field.ptr, &nmiss);
              field.nmiss   = (size_t)nmiss;
196
197
	      field.grid    = pvars1->grid;
	      field.missval = pvars1->missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198

Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
	      if ( field.nmiss > 0 || psamp1->ptr )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
		{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
		  if ( psamp1->ptr == NULL )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
		    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
		      psamp1->ptr = (double*) Malloc(gridsize*sizeof(double));
204
		      for ( int i = 0; i < gridsize; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
			psamp1->ptr[i] = nsets;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
207
		    }
		  
208
209
		  for ( int i = 0; i < gridsize; i++ )
		    if ( !DBL_IS_EQUAL(field.ptr[i], pvars1->missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
210
		      psamp1->ptr[i]++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
212
		}

213
	      if ( lvarstd )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
214
		{
215
216
		  farsumq(pvars2, field);
		  farsum(pvars1, field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
		}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
219
220
221
222
              else if ( lrange )
                {
                  farmin(pvars2, field);
                  farmax(pvars1, field);
                }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223
224
	      else
		{
225
		  farfun(pvars1, field, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
227
228
229
		}
	    }
	}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
230
      if ( seas_nsets[seas] == 0 && lvarstd )
231
232
233
234
235
236
237
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
            field_type *pvars1 = &vars1[seas][varID][levelID];
            field_type *pvars2 = &vars2[seas][varID][levelID];

Uwe Schulzweida's avatar
Uwe Schulzweida committed
238
	    if ( vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
239
240

            farmoq(pvars2, *pvars1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
242
	  }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
      seas_nsets[seas]++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
244
245
246
      tsID++;
    }

247
  for ( int seas = 0; seas < NSEAS; seas++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
    if ( seas_nsets[seas] )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
249
      {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
        int nsets = seas_nsets[seas];
        for ( int recID = 0; recID < maxrecs; recID++ )
          {
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
            field_type *psamp1 = &samp1[seas][varID][levelID];
            field_type *pvars1 = &vars1[seas][varID][levelID];
            field_type *pvars2 = vars2[seas] ? &vars2[seas][varID][levelID] : NULL;

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

            if ( lmean )
              {
                if ( psamp1->ptr ) fardiv(pvars1, *psamp1);
                else               farcdiv(pvars1, (double)nsets);
              }
            else if ( lvarstd )
              {
                if ( psamp1->ptr )
                  {
                    if ( lstd ) farstd(pvars1, *pvars2, *psamp1, divisor);
                    else        farvar(pvars1, *pvars2, *psamp1, divisor);
                  }
                else
                  {
                    if ( lstd ) farcstd(pvars1, *pvars2, nsets, divisor);
                    else        farcvar(pvars1, *pvars2, nsets, divisor);
                  }
              }
            else if ( lrange )
              {
                farsub(pvars1, *pvars2);
              }
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284

285
286
	taxisDefVdate(taxisID2, datetime[seas].vdate);
	taxisDefVtime(taxisID2, datetime[seas].vtime);
287
	streamDefTimestep(streamID2, otsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288

289
	for ( int recID = 0; recID < maxrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
	  {
291
292
293
            int varID   = recinfo[recID].varID;
            int levelID = recinfo[recID].levelID;
            field_type *pvars1 = &vars1[seas][varID][levelID];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294

Uwe Schulzweida's avatar
Uwe Schulzweida committed
295
	    if ( otsID && vlistInqVarTsteptype(vlistID1, varID) == TSTEP_CONSTANT ) continue;
296
297

	    streamDefRecord(streamID2, varID, levelID);
298
	    streamWriteRecord(streamID2, pvars1->ptr, (int)pvars1->nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
299
	  }
300
301

	otsID++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
303
      }

304
  for ( int seas = 0; seas < NSEAS; seas++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
305
306
307
    {
      if ( vars1[seas] != NULL )
	{
308
309
	  field_free(vars1[seas], vlistID1);
	  field_free(samp1[seas], vlistID1);
310
	  if ( lvarstd ) Free(vars2[seas]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
311
312
313
	}
    }

314
  if ( field.ptr ) Free(field.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315

316
  Free(recinfo);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
317
318
319
320
321
322

  streamClose(streamID2);
  streamClose(streamID1);

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
323
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324
}