Merstat.cc 5.3 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
      Merstat    merrange        Meridional range
22
23
24
25
26
27
      Merstat    mermin          Meridional minimum
      Merstat    mermax          Meridional maximum
      Merstat    mersum          Meridional sum
      Merstat    mermean         Meridional mean
      Merstat    meravg          Meridional average
      Merstat    merstd          Meridional standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
      Merstat    merstd          Meridional standard deviation [Normalize by (n-1)]
29
      Merstat    mervar          Meridional variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
      Merstat    mervar          Meridional variance [Normalize by (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
      Merstat    merpctl         Meridional percentiles
32
33
34
*/


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


void *Merstat(void *argument)
{
44
  int gridID1, gridID2 = -1, lastgrid = -1;
45
  int wstatus = FALSE;
46
  int index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
  int nmiss;
48
  int nrecs;
49
  int varID, levelID;
50
  char varname[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
51
52
53

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
54
55
56
57
58
59
60
61
62
63
64
  cdoOperatorAdd("merrange", func_range, 0, NULL);
  cdoOperatorAdd("mermin",   func_min,   0, NULL);
  cdoOperatorAdd("mermax",   func_max,   0, NULL);
  cdoOperatorAdd("mersum",   func_sum,   0, NULL);
  cdoOperatorAdd("mermean",  func_meanw, 1, NULL);
  cdoOperatorAdd("meravg",   func_avgw,  1, NULL);
  cdoOperatorAdd("mervar",   func_varw,  1, NULL);
  cdoOperatorAdd("mervar1",  func_var1w, 1, NULL);
  cdoOperatorAdd("merstd",   func_stdw,  1, NULL);
  cdoOperatorAdd("merstd1",  func_std1w, 1, NULL);
  cdoOperatorAdd("merpctl",  func_pctl,  0, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65
 
66
67
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
68
  bool needWeights = cdoOperatorF2(operatorID) != 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69

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

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

80
81
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
82

83
84
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
85
86
  vlistDefTaxis(vlistID2, taxisID2);

87
88
  int ngrids = vlistNgrids(vlistID1);
  int ndiffgrids = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
90
91
92
93
94
  for ( index = 1; index < ngrids; index++ )
    if ( vlistGrid(vlistID1, 0) != vlistGrid(vlistID1, index))
      ndiffgrids++;

  if ( ndiffgrids > 0 ) cdoAbort("Too many different grids!");

Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
96
97
  index = 0;
  gridID1 = vlistGrid(vlistID1, index);

98
99
100
101
102
103
104
105
106
107
108
  if ( gridInqType(gridID1) == GRID_LONLAT   ||
       gridInqType(gridID1) == GRID_GAUSSIAN ||
       gridInqType(gridID1) == GRID_GENERIC )
    {
      gridID2 = gridToMeridional(gridID1);
    }
  else
    {
      cdoAbort("Unsupported gridtype: %s", gridNamePtr(gridInqType(gridID1)));
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
110
  vlistChangeGridIndex(vlistID2, index, gridID2);

111
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
113
114
115

  streamDefVlist(streamID2, vlistID2);

  gridID1 = vlistInqVarGrid(vlistID1, 0);
116
117
  int nlonmax = gridInqXsize(gridID1); /* max nlon ? */
  int lim = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118

119
  field_type field1, field2;
120
121
122
  field_init(&field1);
  field_init(&field2);

123
  field1.ptr    = (double*) Malloc(lim*sizeof(double));
124
125
  field1.weight = NULL;
  if ( needWeights )
126
    field1.weight = (double*) Malloc(lim*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
127

128
  field2.ptr  = (double*) Malloc(nlonmax*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
130
  field2.grid = gridID2;

131
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
133
134
135
136
137
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
      taxisCopyTimestep(taxisID2, taxisID1);

      streamDefTimestep(streamID2, tsID);

138
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
	{
	  streamInqRecord(streamID1, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
141
142
	  streamReadRecord(streamID1, field1.ptr, &nmiss);
          field1.nmiss = (size_t) nmiss;
143
	  field1.grid = vlistInqVarGrid(vlistID1, varID);
144
	  if ( needWeights && field1.grid != lastgrid )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
145
146
	    {
	      lastgrid = field1.grid;
147
	      wstatus = gridWeights(field1.grid, field1.weight);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
	    }
149
150
151
152
153
	  if ( wstatus != 0 && tsID == 0 && levelID == 0 )
	    {
	      vlistInqVarName(vlistID1, varID, varname);
	      cdoWarning("Using constant grid cell area weights for variable %s!", varname);
	    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
155
156
	  field1.missval = vlistInqVarMissval(vlistID1, varID);
	  field2.missval = vlistInqVarMissval(vlistID1, varID);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
158
159
160
	  if ( operfunc == func_pctl )
	    merpctl(field1, & field2, pn);
	  else  
	    merfun(field1, &field2, operfunc);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
162

	  streamDefRecord(streamID2, varID,  levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
	  streamWriteRecord(streamID2, field2.ptr, (int)field2.nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
165
166
167
168
169
170
	}
      tsID++;
    }

  streamClose(streamID2);
  streamClose(streamID1);

171
172
173
  if ( field1.ptr )    Free(field1.ptr);
  if ( field1.weight ) Free(field1.weight);
  if ( field2.ptr )    Free(field2.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
176

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
177
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
}