Change.c 9.95 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-2014 Uwe Schulzweida, Uwe.Schulzweida@zmaw.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

#define  MAXARG     16384

void *Change(void *argument)
{
42
  int CHCODE, CHTABNUM, CHPARAM, CHNAME, CHUNIT, CHLEVEL, CHLEVELC, CHLEVELV, CHLTYPE;  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
45
46
47
48
  int operatorID;
  int streamID1, streamID2 = CDI_UNDEFID;
  int nrecs, nvars;
  int tsID1, recID, varID = 0, levelID;
  int vlistID1, vlistID2;
  int taxisID1, taxisID2;
49
  int chints[MAXARG], nch = 0;
50
  char *chnames[MAXARG];
51
  char varname[CDI_MAX_NAME];
52
  char *chname = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
  int chcode = 0;
54
  int param;
55
  int code, tabnum, i;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
57
58
  int nmiss;
  int gridsize;
  int nfound;
59
  int nzaxis, zaxisID1, zaxisID2, k, nlevs, index; 
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
  double chlevels[MAXARG];
61
  int  chltypes[MAXARG];              
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62
  double *levels = NULL;
63
  double *newlevels = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
64
65
66
67
  double *array = NULL;

  cdoInitialize(argument);

68
  CHCODE   = cdoOperatorAdd("chcode",   0, 0, "pairs of old and new code numbers");
69
  CHTABNUM = cdoOperatorAdd("chtabnum", 0, 0, "pairs of old and new GRIB1 table numbers");
70
  CHPARAM  = cdoOperatorAdd("chparam",  0, 0, "pairs of old and new parameter identifiers");
71
  CHNAME   = cdoOperatorAdd("chname",   0, 0, "pairs of old and new variable names");
72
  CHUNIT   = cdoOperatorAdd("chunit",   0, 0, "pairs of old and new variable units");
73
  CHLEVEL  = cdoOperatorAdd("chlevel",  0, 0, "pairs of old and new levels");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
74
75
  CHLEVELC = cdoOperatorAdd("chlevelc", 0, 0, "code number, old and new level");
  CHLEVELV = cdoOperatorAdd("chlevelv", 0, 0, "variable name, old and new level");
76
  CHLTYPE  = cdoOperatorAdd("chltype",  0, 0, "pairs of old and new type");          
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77
78
79
80
81
82
83

  operatorID = cdoOperatorID();

  operatorInputArg(cdoOperatorEnter(operatorID));

  nch = operatorArgc();

84
  if ( operatorID == CHCODE || operatorID == CHTABNUM )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
86
87
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
88
	chints[i] = atoi(operatorArgv()[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
    }
90
  else if ( operatorID == CHPARAM || operatorID == CHNAME || operatorID == CHUNIT )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
93
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
94
	chnames[i] = operatorArgv()[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
    }
  else if ( operatorID == CHLEVEL )
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
	chlevels[i] = atof(operatorArgv()[i]);
    }
  else if ( operatorID == CHLEVELC )
    {
      operatorCheckArgc(3);
      
      chcode = atoi(operatorArgv()[0]);
      chlevels[0] = atof(operatorArgv()[1]);
      chlevels[1] = atof(operatorArgv()[2]);
    }
  else if ( operatorID == CHLEVELV )
    {
      operatorCheckArgc(3);
      
114
      chname = operatorArgv()[0];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
116
117
      chlevels[0] = atof(operatorArgv()[1]);
      chlevels[1] = atof(operatorArgv()[2]);
    }
118
119
120
121
122
123
  else if ( operatorID == CHLTYPE )                  
    {
      if ( nch%2 ) cdoAbort("Odd number of input arguments!");
      for ( i = 0; i < nch; i++ )
	chltypes[i] = atoi(operatorArgv()[i]);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140

  streamID1 = streamOpenRead(cdoStreamName(0));

  vlistID1 = streamInqVlist(streamID1);
  vlistID2 = vlistDuplicate(vlistID1);

  taxisID1 = vlistInqTaxis(vlistID1);
  taxisID2 = taxisDuplicate(taxisID1);
  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 )
141
142
	    if ( code == chints[i] )
	      vlistDefVarCode(vlistID2, varID, chints[i+1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
144
	}
    }
145
146
147
148
149
150
151
152
  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 )
153
	    if ( tabnum == chints[i] )
154
	      {
155
		tableID = tableDef(-1, chints[i+1], NULL);
156
157
158
159
		vlistDefVarTable(vlistID2, varID, tableID);
	      }
	}
    }
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
  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]));
	}
    }
177
  else if ( operatorID == CHNAME )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
179
180
181
182
183
    {
      nvars = vlistNvars(vlistID2);
      for ( varID = 0; varID < nvars; varID++ )
	{
	  vlistInqVarName(vlistID2, varID, varname);
	  for ( i = 0; i < nch; i += 2 )
184
185
	    if ( strcmp(varname, chnames[i]) == 0 )
	      vlistDefVarName(vlistID2, varID, chnames[i+1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
187
	}
    }
188
189
190
191
192
193
194
195
196
197
198
199
  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
200
201
202
203
204
205
206
  else if ( operatorID == CHLEVEL )
    {
      nzaxis = vlistNzaxis(vlistID2);
      for ( index = 0; index < nzaxis; index++ )
	{
	  zaxisID1 = vlistZaxis(vlistID2, index);
	  nlevs = zaxisInqSize(zaxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
208
	  levels = (double*) malloc(nlevs*sizeof(double));
	  newlevels = (double*) malloc(nlevs*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209
	  zaxisInqLevels(zaxisID1, levels);
210
211
212

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
214
215
216
217
218
219
220
221
222
223
	  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 )
224
		    newlevels[k] = chlevels[i+1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
225

226
	      zaxisDefLevels(zaxisID2, newlevels);
227
	      vlistChangeZaxis(vlistID2, zaxisID1, zaxisID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
229
230
	    }

	  free(levels);
231
	  free(newlevels);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
	}
    }
  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);
251
	      if ( strcmp(varname, chname) == 0 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
	    }
253
	  if ( varID == nvars ) cdoAbort("Variable name %s not found!", chname);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
254
255
256
257
	}

      zaxisID1 = vlistInqVarZaxis(vlistID2, varID);
      nlevs = zaxisInqSize(zaxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
      levels = (double*) malloc(nlevs*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
      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);
    }
279
280
  else if ( operatorID == CHLTYPE )                
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
281
      int zaxistype, zaxistype1, zaxistype2, ltype, ltype1, ltype2, sameltype;
282
283
284
285
286
287
288

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
290
	  zaxistype = zaxisInqType(zaxisID1);
	  ltype = zaxisInqLtype(zaxisID1);
291
292
293

	  for ( i = 0; i < nch; i += 2 )
	    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
	      sameltype = FALSE;
295
296
297
298
299
300
	      ltype1 = chltypes[i];
	      ltype2 = chltypes[i+1];

	      zaxistype1 = ltype2ztype(ltype1);
	      zaxistype2 = ltype2ztype(ltype2);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
301
302
303
304
305
306
307
308
	      if ( zaxistype1 == zaxistype ) sameltype = TRUE;

	      if ( !(zaxistype1 == ZAXIS_GENERIC && ltype1 == ltype) ) sameltype = FALSE;

	      if ( sameltype)
		{
		  zaxisChangeType(zaxisID2, zaxistype2);
		  if ( zaxistype == ZAXIS_GENERIC ) zaxisDefLtype(zaxisID2, ltype2);
309
		  vlistChangeZaxis(vlistID2, zaxisID1, zaxisID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310
		}
311
312
313
	    }
	}
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
314
315
316
317
318
319

  streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());

  streamDefVlist(streamID2, vlistID2);

  gridsize = vlistGridsizeMax(vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320
  array = (double*) malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348

  tsID1 = 0;
  while ( (nrecs = streamInqTimestep(streamID1, tsID1)) )
    {
      taxisCopyTimestep(taxisID2, taxisID1);

      streamDefTimestep(streamID2, tsID1);
	       
      for ( recID = 0; recID < nrecs; recID++ )
	{
	  streamInqRecord(streamID1, &varID, &levelID);
	  streamDefRecord(streamID2,  varID,  levelID);
	  
	  streamReadRecord(streamID1, array, &nmiss);
	  streamWriteRecord(streamID2, array, nmiss);
	}
      tsID1++;
    }

  streamClose(streamID1);
  streamClose(streamID2);

  if ( array ) free(array);

  cdoFinish();

  return (0);
}