Zonstat.c 5.35 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
/*
  This file is part of CDO. CDO is a collection of Operators to
3
  manipulate and analyse Climate model Data.1
Uwe Schulzweida's avatar
Uwe Schulzweida committed
4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
5
  Copyright (C) 2003-2015 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
21
22
/*
   This module contains the following operators:

      Zonstat    zonmin          Zonal minimum
      Zonstat    zonmax          Zonal maximum
Uwe Schulzweida's avatar
Uwe Schulzweida committed
23
      Zonstat    zonrange        Zonal range
Uwe Schulzweida's avatar
Uwe Schulzweida committed
24
25
26
27
      Zonstat    zonsum          Zonal sum
      Zonstat    zonmean         Zonal mean
      Zonstat    zonavg          Zonal average
      Zonstat    zonstd          Zonal standard deviation
28
      Zonstat    zonstd1         Zonal standard deviation [Divisor is (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
      Zonstat    zonvar          Zonal variance
30
      Zonstat    zonvar1         Zonal variance [Divisor is (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
      Zonstat    zonpctl         Zonal percentiles
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32
33
34
*/


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


void *Zonstat(void *argument)
{
44
45
  int gridID1 = -1, gridID2 = -1;
  int zongridID = -1;
46
  int index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
  int recID, nrecs;
48
  int varID, levelID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
53
54
55
56
57
58
  cdoOperatorAdd("zonmin",   func_min,   0, NULL);
  cdoOperatorAdd("zonmax",   func_max,   0, NULL);
  cdoOperatorAdd("zonrange", func_range, 0, NULL);
  cdoOperatorAdd("zonsum",   func_sum,   0, NULL);
  cdoOperatorAdd("zonmean",  func_mean,  0, NULL);
  cdoOperatorAdd("zonavg",   func_avg,   0, NULL);
  cdoOperatorAdd("zonvar",   func_var,   0, NULL);
59
  cdoOperatorAdd("zonvar1",  func_var1,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
  cdoOperatorAdd("zonstd",   func_std,   0, NULL);
61
  cdoOperatorAdd("zonstd1",  func_std1,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62
  cdoOperatorAdd("zonpctl",  func_pctl,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63

64
65
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
66

67
  double pn = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
69
70
  if ( operfunc == func_pctl )
    {
      operatorInputArg("percentile number");
71
72
      pn = parameter2double(operatorArgv()[0]);
      percentile_check_number(pn);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
73
74
    }

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

77
78
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79

80
81
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82
83
  vlistDefTaxis(vlistID2, taxisID2);

84
  int ngrids = vlistNgrids(vlistID1);
85
86
87
88
89
90
91
92
93
94
95
  for ( index = 0; index < ngrids; index++ )
    {
      if ( gridInqXsize(vlistGrid(vlistID1, index)) > 1 ) 
	{
	  if ( gridID1 == -1 ) gridID1 = vlistGrid(vlistID1, index);
	}
      else
	{
	  if ( zongridID == -1 ) zongridID = vlistGrid(vlistID1, index);
	}
    }
96

97
  int ndiffgrids = 0;
98
99
100
101
102
  for ( index = 0; index < ngrids; index++ )
    {
      if ( zongridID != -1 && zongridID == vlistGrid(vlistID1, index) ) continue;
      if ( gridID1 != vlistGrid(vlistID1, index) ) ndiffgrids++;
    }
103
104
  if ( ndiffgrids > 0 ) cdoAbort("Too many different grids!");

105
  if ( gridID1 != -1 )
106
    {
107
108
109
110
111
112
113
114
115
      if ( gridInqType(gridID1) == GRID_LONLAT   ||
	   gridInqType(gridID1) == GRID_GAUSSIAN ||
	   gridInqType(gridID1) == GRID_GENERIC )
	{
	  if ( zongridID != -1 && gridInqYsize(zongridID) == gridInqYsize(gridID1) )
	    gridID2 = zongridID;
	  else
	    gridID2 = gridToZonal(gridID1);
	}
116
      else
117
118
119
	{
	  cdoAbort("Unsupported gridtype: %s", gridNamePtr(gridInqType(gridID1)));
	}
120
121
122
    }
  else
    {
123
124
      gridID2 = zongridID;
      cdoWarning("Input stream contains only zonal data!");
125
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126

127
128
  for ( index = 0; index < ngrids; index++ )
    vlistChangeGridIndex(vlistID2, index, gridID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129

130
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132
133
134

  streamDefVlist(streamID2, vlistID2);

  gridID1 = vlistInqVarGrid(vlistID1, 0);
135
  int nlatmax = gridInqYsize(gridID1); /* max nlat ? */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136

137
138
139
  int lim = vlistGridsizeMax(vlistID1);

  field_t field1, field2;
140
141
  field_init(&field2);
  field_init(&field2);
142
143
  field1.ptr  = (double*) Malloc(lim*sizeof(double));
  field2.ptr  = (double*) Malloc(nlatmax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
145
  field2.grid = gridID2;

146
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
      taxisCopyTimestep(taxisID2, taxisID1);

      streamDefTimestep(streamID2, tsID);

      for ( recID = 0; recID < nrecs; recID++ )
	{
	  streamInqRecord(streamID1, &varID, &levelID);
	  streamReadRecord(streamID1, field1.ptr, &field1.nmiss);

	  field1.grid    = vlistInqVarGrid(vlistID1, varID);
	  field1.missval = vlistInqVarMissval(vlistID1, varID);
	  field2.missval = vlistInqVarMissval(vlistID1, varID);

162
163
164
165
166
167
168
169
170
171
172
173
	  if ( zongridID != -1 && zongridID == field1.grid )
	    {
	      memcpy(field2.ptr, field1.ptr, nlatmax*sizeof(double));
	      field2.nmiss = field1.nmiss;
	    }
	  else
	    {
	      if ( operfunc == func_pctl )
		zonpctl(field1, & field2, pn);
	      else  
		zonfun(field1, &field2, operfunc);
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
176
177
178
179
180
181
182
183

	  streamDefRecord(streamID2, varID,  levelID);
	  streamWriteRecord(streamID2, field2.ptr, field2.nmiss);
	}
      tsID++;
    }

  streamClose(streamID2);
  streamClose(streamID1);

184
185
  if ( field1.ptr ) Free(field1.ptr);
  if ( field2.ptr ) Free(field2.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
187
188

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
}