Change.cc 10.1 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:

21
      Change     chcode          Change code number
22
      Change     chtabnum        Change GRIB1 parameter table number
23
      Change     chparam         Change parameter identifier
24
      Change     chname          Change variable name
25
26
27
      Change     chlevel         Change level
      Change     chlevelc        Change level of one code
      Change     chlevelv        Change level of one variable
28
      Change     chltype         Change GRIB level type 
29
30
*/

Ralf Mueller's avatar
Ralf Mueller committed
31
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
33
34
35
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"

36
int stringToParam(const char *paramstr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
38
39
40
41
42

#define  MAXARG     16384

void *Change(void *argument)
{
  int nrecs, nvars;
43
44
  int varID = 0, levelID;
  int chints[MAXARG];
45
  char *chnames[MAXARG];
46
  char varname[CDI_MAX_NAME];
47
  char *chname = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
48
  int chcode = 0;
49
  int param;
50
  int code, tabnum, i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
52
  int nmiss;
  int nfound;
53
  int nzaxis, zaxisID1, zaxisID2, k, nlevs, index; 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
  double chlevels[MAXARG];
55
  int  chltypes[MAXARG];              
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
  double *levels = NULL;
57
  double *newlevels = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
59
60

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
61
  // clang-format off
62
63
64
65
66
67
68
69
70
  int CHCODE   = cdoOperatorAdd("chcode",   0, 0, "pairs of old and new code numbers");
  int CHTABNUM = cdoOperatorAdd("chtabnum", 0, 0, "pairs of old and new GRIB1 table numbers");
  int CHPARAM  = cdoOperatorAdd("chparam",  0, 0, "pairs of old and new parameter identifiers");
  int CHNAME   = cdoOperatorAdd("chname",   0, 0, "pairs of old and new variable names");
  int CHUNIT   = cdoOperatorAdd("chunit",   0, 0, "pairs of old and new variable units");
  int CHLEVEL  = cdoOperatorAdd("chlevel",  0, 0, "pairs of old and new levels");
  int CHLEVELC = cdoOperatorAdd("chlevelc", 0, 0, "code number, old and new level");
  int CHLEVELV = cdoOperatorAdd("chlevelv", 0, 0, "variable name, old and new level");
  int CHLTYPE  = cdoOperatorAdd("chltype",  0, 0, "pairs of old and new type");          
Uwe Schulzweida's avatar
Uwe Schulzweida committed
71
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72

73
  int operatorID = cdoOperatorID();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
76

  operatorInputArg(cdoOperatorEnter(operatorID));

77
  int nch = operatorArgc();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78

79
  if ( operatorID == CHCODE || operatorID == CHTABNUM )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
83
	chints[i] = parameter2int(operatorArgv()[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
    }
85
  else if ( operatorID == CHPARAM || operatorID == CHNAME || operatorID == CHUNIT )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
87
88
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
89
	chnames[i] = operatorArgv()[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
90
91
92
93
94
    }
  else if ( operatorID == CHLEVEL )
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
95
	chlevels[i] = parameter2double(operatorArgv()[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
97
98
99
100
    }
  else if ( operatorID == CHLEVELC )
    {
      operatorCheckArgc(3);
      
101
      chcode = parameter2int(operatorArgv()[0]);
102
103
      chlevels[0] = parameter2double(operatorArgv()[1]);
      chlevels[1] = parameter2double(operatorArgv()[2]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
105
106
107
108
    }
  else if ( operatorID == CHLEVELV )
    {
      operatorCheckArgc(3);
      
109
      chname = operatorArgv()[0];
110
111
      chlevels[0] = parameter2double(operatorArgv()[1]);
      chlevels[1] = parameter2double(operatorArgv()[2]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
    }
113
114
115
116
  else if ( operatorID == CHLTYPE )                  
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
117
	chltypes[i] = parameter2int(operatorArgv()[i]);
118
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
119

120
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121

122
  int vlistID1 = pstreamInqVlist(streamID1);
123
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124

125
126
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127
128
129
130
131
132
133
134
135
  vlistDefTaxis(vlistID2, taxisID2);

  if ( operatorID == CHCODE )
    {
      nvars = vlistNvars(vlistID2);
      for ( varID = 0; varID < nvars; varID++ )
	{
	  code = vlistInqVarCode(vlistID2, varID);
	  for ( i = 0; i < nch; i += 2 )
136
137
	    if ( code == chints[i] )
	      vlistDefVarCode(vlistID2, varID, chints[i+1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
138
139
	}
    }
140
141
142
143
144
145
146
147
  else if ( operatorID == CHTABNUM )
    {
      int tableID;
      nvars = vlistNvars(vlistID2);
      for ( varID = 0; varID < nvars; varID++ )
	{
	  tabnum = tableInqNum(vlistInqVarTable(vlistID2, varID));
	  for ( i = 0; i < nch; i += 2 )
148
	    if ( tabnum == chints[i] )
149
	      {
150
		tableID = tableDef(-1, chints[i+1], NULL);
151
152
153
154
		vlistDefVarTable(vlistID2, varID, tableID);
	      }
	}
    }
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
  else if ( operatorID == CHPARAM )
    {
      nvars = vlistNvars(vlistID2);
      for ( varID = 0; varID < nvars; varID++ )
	{
	  param = vlistInqVarParam(vlistID2, varID);
	  if ( cdoVerbose )
	    {
	      int pnum, pcat, pdis;
	      cdiDecodeParam(param, &pnum, &pcat, &pdis);
	      cdoPrint("pnum, pcat, pdis: %d.%d.%d", pnum, pcat, pdis);
	    }
	  for ( i = 0; i < nch; i += 2 )
	    if ( param == stringToParam(chnames[i]) )
	      vlistDefVarParam(vlistID2, varID, stringToParam(chnames[i+1]));
	}
    }
172
  else if ( operatorID == CHNAME )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
173
174
175
176
177
178
    {
      nvars = vlistNvars(vlistID2);
      for ( varID = 0; varID < nvars; varID++ )
	{
	  vlistInqVarName(vlistID2, varID, varname);
	  for ( i = 0; i < nch; i += 2 )
179
180
	    if ( strcmp(varname, chnames[i]) == 0 )
	      vlistDefVarName(vlistID2, varID, chnames[i+1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
182
	}
    }
183
184
185
186
187
188
189
190
191
192
193
194
  else if ( operatorID == CHUNIT )
    {
      nvars = vlistNvars(vlistID2);
      for ( varID = 0; varID < nvars; varID++ )
	{

	  vlistInqVarUnits(vlistID2, varID, varname);
	  for ( i = 0; i < nch; i += 2 )
	    if ( strcmp(varname, chnames[i]) == 0 )
	      vlistDefVarUnits(vlistID2, varID, chnames[i+1]);
	}
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
195
196
197
198
199
200
  else if ( operatorID == CHLEVEL )
    {
      nzaxis = vlistNzaxis(vlistID2);
      for ( index = 0; index < nzaxis; index++ )
	{
	  zaxisID1 = vlistZaxis(vlistID2, index);
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
          if ( zaxisInqLevels(zaxisID1, NULL) )
            {
              nlevs = zaxisInqSize(zaxisID1);
              levels = (double*) Malloc(nlevs*sizeof(double));
              newlevels = (double*) Malloc(nlevs*sizeof(double));
              zaxisInqLevels(zaxisID1, levels);

              for ( k = 0; k < nlevs; k++ ) newlevels[k] = levels[k];

              nfound = 0;
              for ( i = 0; i < nch; i += 2 )
                for ( k = 0; k < nlevs; k++ )
                  if ( fabs(levels[k] - chlevels[i]) < 0.0001 ) nfound++;

              if ( nfound )
                {
                  zaxisID2 = zaxisDuplicate(zaxisID1);
                  for ( i = 0; i < nch; i += 2 )
                    for ( k = 0; k < nlevs; k++ )
                      if ( fabs(levels[k] - chlevels[i]) < 0.001 )
                        newlevels[k] = chlevels[i+1];

                  zaxisDefLevels(zaxisID2, newlevels);
                  vlistChangeZaxis(vlistID2, zaxisID1, zaxisID2);
                }

              Free(levels);
              Free(newlevels);
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
    }
  else if ( operatorID == CHLEVELC || operatorID == CHLEVELV )
    {
      nvars = vlistNvars(vlistID2);
      if ( operatorID == CHLEVELC )
	{
	  for ( varID = 0; varID < nvars; varID++ )
	    {
	      code = vlistInqVarCode(vlistID2, varID);
	      if ( code == chcode ) break;
	    }
	  if ( varID == nvars ) cdoAbort("Code %d not found!", chcode);
	}
      else
	{
	  for ( varID = 0; varID < nvars; varID++ )
	    {
	      vlistInqVarName(vlistID2, varID, varname);
249
	      if ( strcmp(varname, chname) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
	    }
251
	  if ( varID == nvars ) cdoAbort("Variable name %s not found!", chname);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
253
254
	}

      zaxisID1 = vlistInqVarZaxis(vlistID2, varID);
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
      if ( zaxisInqLevels(zaxisID1, NULL) )
        {
          nlevs = zaxisInqSize(zaxisID1);
          levels = (double*) Malloc(nlevs*sizeof(double));
          zaxisInqLevels(zaxisID1, levels);
          nfound = 0;
          for ( k = 0; k < nlevs; k++ )
            if ( fabs(levels[k] - chlevels[0]) < 0.0001 ) nfound++;

          if ( nfound )
            {
              zaxisID2 = zaxisDuplicate(zaxisID1);
              for ( k = 0; k < nlevs; k++ )
                if ( fabs(levels[k] - chlevels[0]) < 0.001 )
                  levels[k] = chlevels[1];

              zaxisDefLevels(zaxisID2, levels);
              vlistChangeVarZaxis(vlistID2, varID, zaxisID2);
            }
          else
            cdoAbort("Level %g not found!", chlevels[0]);

          Free(levels);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279
    }
280
281
  else if ( operatorID == CHLTYPE )                
    {
282
      int ltype, ltype1, ltype2;
283
284
285
286
287
288
289

      nzaxis = vlistNzaxis(vlistID2);
      for ( index = 0; index < nzaxis; index++ )
	{
	  zaxisID1 = vlistZaxis(vlistID2, index);
	  zaxisID2 = zaxisDuplicate(zaxisID1);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
	  ltype = zaxisInqLtype(zaxisID1);
291
292
293
294
295
296

	  for ( i = 0; i < nch; i += 2 )
	    {
	      ltype1 = chltypes[i];
	      ltype2 = chltypes[i+1];

297
	      if ( ltype1 == ltype )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298
		{
299
300
		  zaxisChangeType(zaxisID2, ZAXIS_GENERIC);
		  zaxisDefLtype(zaxisID2, ltype2);
301
		  vlistChangeZaxis(vlistID2, zaxisID1, zaxisID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302
		}
303
304
305
	    }
	}
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306

307
308
  int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
  pstreamDefVlist(streamID2, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
309

310
311
  int gridsize = vlistGridsizeMax(vlistID2);
  double *array = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
312

313
  int tsID1 = 0;
314
  while ( (nrecs = pstreamInqTimestep(streamID1, tsID1)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
315
316
317
    {
      taxisCopyTimestep(taxisID2, taxisID1);

318
      pstreamDefTimestep(streamID2, tsID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
	       
320
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
321
	{
322
323
	  pstreamInqRecord(streamID1, &varID, &levelID);
	  pstreamDefRecord(streamID2,  varID,  levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
324
	  
325
326
	  pstreamReadRecord(streamID1, array, &nmiss);
	  pstreamWriteRecord(streamID2, array, nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
327
328
329
330
	}
      tsID1++;
    }

331
332
  pstreamClose(streamID1);
  pstreamClose(streamID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
333

334
  if ( array ) Free(array);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
336
337

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
338
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
339
}