Gridboxstat.c 22.5 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
18
19
20
21
22
23
24
25
26
  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.
*/

/*
   This module contains the following operators:

      Gridboxstat    gridboxmin          Gridbox minimum
      Gridboxstat    gridboxmax          Gridbox maximum
      Gridboxstat    gridboxsum          Gridbox sum
      Gridboxstat    gridboxmean         Gridbox mean
      Gridboxstat    gridboxavg          Gridbox average
      Gridboxstat    gridboxstd          Gridbox standard deviation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
      Gridboxstat    gridboxstd1         Gridbox standard deviation [Normalize by (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
28
      Gridboxstat    gridboxvar          Gridbox variance
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
      Gridboxstat    gridboxvar1         Gridbox variance [Normalize by (n-1)]
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30
31
32
*/


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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
40
41
42
43

static
int genBoxGrid(int gridID1, int xinc, int yinc)
{
  int i, j, i1;
44
45
  int gridID2 = -1;
  int nlon1 = 0, nlat1 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46
  int gridsize2 = 0, nlon2 = 0, nlat2 = 0;
47
48
49
50
51
  int x1, y1, x2, y2;
  int use_x1, use_y1;
  int corner, add, g2_add; 
  int g1_add;
  int circular = gridIsCircular(gridID1) ;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52
53
  double *grid1_corner_lon = NULL, *grid1_corner_lat = NULL;
  double *grid2_corner_lon = NULL, *grid2_corner_lat = NULL;
54
  double on_up, on_lo, ox_up, ox_lo, an_le, an_ri, ax_le, ax_ri;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
55
  double xvals2_0 = 0;
56
  double area_norm;  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57

58
59
  int gridtype  = gridInqType(gridID1);
  int gridsize1 = gridInqSize(gridID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60

61
  if ( xinc < 1 || yinc < 1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62
    cdoAbort("xinc and yinc must not be smaller than 1!");
63
  
64
65
  if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT ||
       gridtype == GRID_CURVILINEAR || gridtype == GRID_GENERIC )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
66
67
68
69
70
71
72
73
74
75
76
77
78
79
    {
      nlon1 = gridInqXsize(gridID1);
      nlat1 = gridInqYsize(gridID1);

      nlon2 = nlon1/xinc;
      nlat2 = nlat1/yinc;
      if ( nlon1%xinc ) nlon2++;
      if ( nlat1%yinc ) nlat2++;
      gridsize2 = nlon2*nlat2;

      gridID2 = gridCreate(gridtype, gridsize2);
      gridDefXsize(gridID2, nlon2);
      gridDefYsize(gridID2, nlat2);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
81
82
83
  else
    {
      cdoAbort("Unsupported grid: %s", gridNamePtr(gridtype));
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84

85
  if ( xinc > nlon1 || yinc > nlat1 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
    cdoAbort("xinc and/or yinc exceeds gridsize!");
87
88

  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
89
90
  if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT )
    {
91
      bool gridHasBounds = (gridInqXbounds(gridID1, NULL) && gridInqYbounds(gridID1, NULL));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
92

93
94
95
96
      double *xvals1 = (double*) Malloc(nlon1*sizeof(double));
      double *yvals1 = (double*) Malloc(nlat1*sizeof(double));
      double *xvals2 = (double*) Malloc(nlon2*sizeof(double));
      double *yvals2 = (double*) Malloc(nlat2*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
97
98
99
      gridInqXvals(gridID1, xvals1);
      gridInqYvals(gridID1, yvals1);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
100
      if ( gridHasBounds )
101
        {
102
103
104
105
          grid1_corner_lon = (double*) Malloc(2*nlon1*sizeof(double));
          grid1_corner_lat = (double*) Malloc(2*nlat1*sizeof(double));
          grid2_corner_lon = (double*) Malloc(2*nlon2*sizeof(double));
          grid2_corner_lat = (double*) Malloc(2*nlat2*sizeof(double));
106
107
108
          gridInqXbounds(gridID1, grid1_corner_lon);
          gridInqYbounds(gridID1, grid1_corner_lat);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
110
111

      j = 0;
      for ( i = 0; i < nlon1; i += xinc )
112
113
114
115
        {
          i1 = i+(xinc-1);
          if ( i1 >= nlon1-1 ) i1 = nlon1-1; 
          xvals2[j] = xvals1[i] + (xvals1[i1] - xvals1[i])/2.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
          if ( gridHasBounds )
117
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118
              grid2_corner_lon[2*j  ] = grid1_corner_lon[2*i];
119
120
              grid2_corner_lon[2*j+1] = grid1_corner_lon[2*i1+1];
            }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
121
          j++;        
122
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123
124
      j = 0;
      for ( i = 0; i < nlat1; i += yinc )
125
126
127
128
        {
          i1 = i+(yinc-1);
          if ( i1 >= nlat1-1 ) i1 = nlat1-1; 
          yvals2[j] = yvals1[i] + (yvals1[i1] - yvals1[i])/2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
          if ( gridHasBounds )
130
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
              grid2_corner_lat[2*j]   = grid1_corner_lat[2*i];
132
133
134
135
136
              grid2_corner_lat[2*j+1] = grid1_corner_lat[2*i1+1];
            }
          j++;
        }
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
137
138
      gridDefXvals(gridID2, xvals2);
      gridDefYvals(gridID2, yvals2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139

140
141
      Free(xvals2);
      Free(yvals2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
142

143
144
      Free(xvals1);
      Free(yvals1);
145
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
146
      if ( gridHasBounds )
147
148
149
150
        {
          gridDefNvertex(gridID2, 2);
          gridDefXbounds(gridID2, grid2_corner_lon);
          gridDefYbounds(gridID2, grid2_corner_lat);
151
152
          Free(grid2_corner_lon);
          Free(grid2_corner_lat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153

154
155
          Free(grid1_corner_lon);
          Free(grid1_corner_lat);
156
        }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
157
    } /* if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_LONLAT ) */
158
159
160
  else if ( gridtype == GRID_GENERIC )
    {
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
162
  else if ( gridtype == GRID_CURVILINEAR )
    {
163
      bool gridHasBounds = (gridInqXbounds(gridID1, NULL) && gridInqYbounds(gridID1, NULL));
Cedrick Ansorge's avatar
Cedrick Ansorge committed
164
      
165
166
167
168
      double *xvals1 = (double*) Malloc(nlon1*nlat1*sizeof(double));
      double *yvals1 = (double*) Malloc(nlon1*nlat1*sizeof(double));
      double *xvals2 = (double*) Malloc(nlon2*nlat2*sizeof(double));
      double *yvals2 = (double*) Malloc(nlon2*nlat2*sizeof(double));
169
170
      gridInqXvals(gridID1, xvals1);
      gridInqYvals(gridID1, yvals1);
171
172
173

      /* Convert lat/lon units if required */
      {
174
	char units[CDI_MAX_NAME];
175
	gridInqXunits(gridID1, units);
176
	grid_to_degree(units, nlon1*nlat1, xvals1, "grid center lon");
177
	gridInqYunits(gridID1, units);
178
	grid_to_degree(units, nlon1*nlat1, yvals1, "grid center lat");
179
      }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
      
Cedrick Ansorge's avatar
Cedrick Ansorge committed
181
182
      if ( gridHasBounds )
        {
183
184
185
186
          grid1_corner_lon = (double*) Malloc(4*nlon1*nlat1*sizeof(double));
          grid1_corner_lat = (double*) Malloc(4*nlon1*nlat1*sizeof(double));
          grid2_corner_lon = (double*) Malloc(4*nlon2*nlat2*sizeof(double));
          grid2_corner_lat = (double*) Malloc(4*nlon2*nlat2*sizeof(double));
Cedrick Ansorge's avatar
Cedrick Ansorge committed
187
188
          gridInqXbounds(gridID1, grid1_corner_lon);
          gridInqYbounds(gridID1, grid1_corner_lat);
189
190
191

	  /* Convert lat/lon units if required */
	  {
192
	    char units[CDI_MAX_NAME];
193
	    gridInqXunits(gridID1, units);
194
	    grid_to_degree(units, 4*nlon1*nlat1, grid1_corner_lon, "grid corner lon");
195
	    gridInqYunits(gridID1, units);
196
	    grid_to_degree(units, 4*nlon1*nlat1, grid1_corner_lat, "grid corner lat");
197
	  }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
198
        }
199
      
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
      /* Process grid2 bounds */
201
202
203
204
205
      area_norm = xinc*yinc;
      for ( y2 = 0; y2 < nlat2; y2++ )
        {
          for ( x2 = 0; x2 < nlon2; x2++ )
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206
              g2_add = (y2*nlon2+x2);
207
208
              on_up = on_lo = 360.; ox_up = ox_lo = -360.;
              an_ri = an_le =  90.; ax_ri = ax_le =  -90.; 
209
210
211
              
              for ( y1 = y2*yinc; y1 < yinc*(y2+1); y1++ )
                {
212
                  if ( y1 >= nlat1 ) use_y1 = nlat1-1;
213
214
                  else use_y1 = y1;
                  for ( x1 = x2*xinc; x1 < xinc*(x2+1) ; x1++ )
215
                    {                    
216
                      use_x1 = x1;
217
218
219
220
221
222
223
224
225
226
227
228
229
                      if ( x1 >= nlon1  ) 
                        {
                          if ( circular && use_y1 == y1 ) 
                            use_y1 -= 1;
                          else            
                            use_x1  = nlon1-1;
                        }
                      
                      g1_add= (use_y1*nlon1)+use_x1;                      

                      if ( y1 == y2*yinc && x1 == x2*xinc )
                        {
                          xvals2_0 = xvals1[g1_add];
Cedrick Ansorge's avatar
Cedrick Ansorge committed
230
                          xvals2[g2_add] = xvals1[g1_add]/area_norm;                        
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
                          yvals2[g2_add] = yvals1[g1_add]/area_norm;
232
                        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
                      else if ( fabs(xvals1[g1_add] - xvals2_0) > 270. )
234
                        {
235
236
237
238
                          if ( (xvals1[g1_add] - xvals2_0) > 270. )
                            xvals2[g2_add] += (xvals1[g1_add]-360.)/area_norm;
                          else if ( ( xvals1[g1_add] - xvals2_0 ) < -270. )
                            xvals2[g2_add] += (xvals1[g1_add]+360.)/area_norm;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
239
			  yvals2[g2_add] += yvals1[g1_add]/area_norm;  
240
                        }
241
                      else                      
Cedrick Ansorge's avatar
Cedrick Ansorge committed
242
                        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
244
                          xvals2[g2_add] += xvals1[g1_add]/area_norm;
                          yvals2[g2_add] += yvals1[g1_add]/area_norm;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
245
246
                        }
                     
Cedrick Ansorge's avatar
Cedrick Ansorge committed
247
                      if ( gridHasBounds )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
248
249
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
284
                        {
			  int c_flag[4], corner2, g1_add2, g1_add3;
			  double lon, lat, lon2, lat2, lon3, lat3;
			  /* upper left cell */
			  if ( y1 == y2*yinc && x1 == x2*xinc)
			    {                                 
			      c_flag[0] = c_flag[1] = c_flag[2] = c_flag[3] = 0;
			      for ( corner = 0; corner < 4; corner++ )
				{                                      
				  add = 4*g1_add + corner;
				  lon = grid1_corner_lon[add];
				  lat = grid1_corner_lat[add]; 
				  g1_add2 = g1_add+1;
				  if ( g1_add+nlon1 > gridsize1 ) 
				    {
				      cdoWarning("Can't find cell below upper left");
				      continue; 
				    }
				  g1_add3 = g1_add+nlon1;
				  for ( corner2 = 0; corner2 < 4; corner2++ )
				    {                                          
				      lon2 = grid1_corner_lon[4*g1_add2+corner2];
				      lat2 = grid1_corner_lat[4*g1_add2+corner2];
				      lon3 = grid1_corner_lon[4*g1_add3+corner2];
				      lat3 = grid1_corner_lat[4*g1_add3+corner2];
				      if ((IS_EQUAL(lon2, lon) && IS_EQUAL(lat2, lat))  ||
					  (IS_EQUAL(lon3, lon) && IS_EQUAL(lat3, lat)) )
					c_flag[corner] = 1;
				    }
				}
			      for ( corner = 0; corner<4; corner++ )
				if ( !c_flag[corner] ) break;                                 
			      on_up = grid1_corner_lon[4*g1_add + corner];
			      ax_le = grid1_corner_lat[4*g1_add + corner];                                  
			      if ( c_flag[0] + c_flag[1] + c_flag[2] + c_flag[3] < 3 )
				cdoWarning("found two matching corners!");                                  
			    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
285
                              
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
			  /* upper right cell */
			  if ( ( y1 == y2*yinc ) && ( x1 == (x2+1)*xinc - 1 ) )
			    {
			      c_flag[0] = c_flag[1] = c_flag[2] = c_flag[3] = 0;
			      for ( corner = 0; corner < 4; corner++ )
				{                                      
				  add = 4*g1_add + corner;
				  lon = grid1_corner_lon[add];
				  lat = grid1_corner_lat[add]; 
				  g1_add2 = g1_add-1;                                      
				  if ( g1_add+nlon1 > gridsize1 ) 
				    {
				      cdoWarning("Can't find cell below upper left");
				      continue; 
				    }
				  g1_add3 = g1_add+nlon1;
				  for ( corner2 = 0; corner2 < 4; corner2++ )
				    {                                          
				      lon2 = grid1_corner_lon[4*g1_add2+corner2];
				      lat2 = grid1_corner_lat[4*g1_add2+corner2];
				      lon3 = grid1_corner_lon[4*g1_add3+corner2];
				      lat3 = grid1_corner_lat[4*g1_add3+corner2];
				      if ((IS_EQUAL(lon2, lon) && IS_EQUAL(lat2, lat))  ||
					  (IS_EQUAL(lon3, lon) && IS_EQUAL(lat3, lat)) )
					c_flag[corner] = 1;
				    } 
				}                                 
			      for ( corner = 0; corner<4; corner++ )
				if ( !c_flag[corner] ) break;                                 
			      ox_up = grid1_corner_lon[4*g1_add + corner];
			      ax_ri = grid1_corner_lat[4*g1_add + corner];                                  
			      if ( c_flag[0] + c_flag[1] + c_flag[2] + c_flag[3] < 3 )
				cdoWarning("found two matching corners!");                                     
			    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
320
321
                            
                              
Uwe Schulzweida's avatar
Uwe Schulzweida committed
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
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
			  /* lower right cell */
			  if ( ( y1 == (y2+1)*yinc -1 ) && (x1 == (x2+1)*xinc -1) )
			    {                                  
			      c_flag[0] = c_flag[1] = c_flag[2] = c_flag[3] = 0;
			      for ( corner = 0; corner < 4; corner++ )
				{                                      
				  add = 4*g1_add + corner;
				  lon = grid1_corner_lon[add];
				  lat = grid1_corner_lat[add]; 
				  g1_add2 = g1_add-1;
				  if ( g1_add-nlon1 < 0 ) 
				    {
				      cdoWarning("Can't find cell above lower right left");
				      continue; 
				    }
				  g1_add3 = g1_add-nlon1;
				  for ( corner2 = 0; corner2 < 4; corner2++ )
				    {                                          
				      lon2 = grid1_corner_lon[4*g1_add2+corner2];
				      lat2 = grid1_corner_lat[4*g1_add2+corner2];
				      lon3 = grid1_corner_lon[4*g1_add3+corner2];
				      lat3 = grid1_corner_lat[4*g1_add3+corner2];
				      if ((IS_EQUAL(lon2, lon) && IS_EQUAL(lat2, lat))  ||
					  (IS_EQUAL(lon3, lon) && IS_EQUAL(lat3, lat)) )
					c_flag[corner] = 1;
				    } 
				}                                  
			      for ( corner = 0; corner<4; corner++ )
				if ( !c_flag[corner] ) break;                                 
			      ox_lo = grid1_corner_lon[4*g1_add + corner];
			      an_ri = grid1_corner_lat[4*g1_add + corner];                                  
			      if ( c_flag[0] + c_flag[1] + c_flag[2] + c_flag[3] < 3 )
				cdoWarning("found two matching corners!");        
			    }
			  
			  /* lower left cell */
			  if ( ( y1 == (y2+1)*yinc -1 ) && ( x1 == x2*xinc ) )
			    {    
			      c_flag[0] = c_flag[1] = c_flag[2] = c_flag[3] = 0;
			      for ( corner = 0; corner < 4; corner++ )
				{                                      
				  add = 4*g1_add + corner;
				  lon = grid1_corner_lon[add];
				  lat = grid1_corner_lat[add]; 
				  g1_add2 = g1_add+1;
				  if ( g1_add-nlon1 < 0 ) 
				    {
				      cdoWarning("Can't find cell above lower right left");
				      continue; 
				    }
				  g1_add3 = g1_add-nlon1;
				  for ( corner2 = 0; corner2 < 4; corner2++ )
				    {                                          
				      lon2 = grid1_corner_lon[4*g1_add2+corner2];
				      lat2 = grid1_corner_lat[4*g1_add2+corner2];
				      lon3 = grid1_corner_lon[4*g1_add3+corner2];
				      lat3 = grid1_corner_lat[4*g1_add3+corner2];
				      if ((IS_EQUAL(lon2, lon) && IS_EQUAL(lat2, lat))  ||
					  (IS_EQUAL(lon3, lon) && IS_EQUAL(lat3, lat)) )
					c_flag[corner] = 1;
				    }
				}                                 
			      for ( corner = 0; corner<4; corner++ )
				if ( !c_flag[corner] ) break;                                 
			      on_lo = grid1_corner_lon[4*g1_add + corner];
			      an_le = grid1_corner_lat[4*g1_add + corner];                                  
			      if ( c_flag[0] + c_flag[1] + c_flag[2] + c_flag[3] < 3 )
				cdoWarning("found two matching corners!");                                      
			    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
                        }/* if ( gridHasBounds) */                  
                    }/* for ( x1 = x2*xinc; x1 < xinc*(x2+1) ; x1++ ) */
                }/* for ( y1 = y2*yinc; y1 < yinc*(y2+1); y1++ ) */
 
              
              if ( gridHasBounds )
                {
                  /* upper left corner */
                  grid2_corner_lon[4*g2_add+3] = on_up;
                  grid2_corner_lat[4*g2_add+3] = ax_le;
                  /* upper right corner */
                  grid2_corner_lon[4*g2_add+2] = ox_up;
                  grid2_corner_lat[4*g2_add+2] = ax_ri;
                  /* lower right corner */
                  grid2_corner_lon[4*g2_add+1] = ox_lo;
                  grid2_corner_lat[4*g2_add+1] = an_ri;
                  /* lower left corner */
                  grid2_corner_lon[4*g2_add+0] = on_lo;
                  grid2_corner_lat[4*g2_add+0] = an_le;
                }
411
              
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
413
	      //  while ( xvals2[g2_add] >  180. ) xvals2[g2_add] -= 360.;
	      //  while ( xvals2[g2_add] < -180. ) xvals2[g2_add] += 360.;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
414
415
            } /* for ( x2 = 0; x2 < nlon2; x2++ ) */
        } /* for ( y2 = 0; y2 < nlat2; y2++ ) */
416
      
417
418
      gridDefXvals(gridID2, xvals2);
      gridDefYvals(gridID2, yvals2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
419

420
421
      Free(xvals2);
      Free(yvals2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
422

423
424
      Free(xvals1);
      Free(yvals1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
425
426

      if ( gridHasBounds )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
427
428
429
430
431
        {
          gridDefNvertex(gridID2, 4);
          gridDefXbounds(gridID2, grid2_corner_lon);
          gridDefYbounds(gridID2, grid2_corner_lat);
          
432
433
          Free(grid2_corner_lon);
          Free(grid2_corner_lat);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
434
          
435
436
          Free(grid1_corner_lon);
          Free(grid1_corner_lat);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
437
438
        }
    } /* else if ( gridtype == GRID_CURVILINEAR ) */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
439
440
441
442
  else
    {
      cdoAbort("Unsupported grid: %s", gridNamePtr(gridtype));
    }
443
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444
445
446
447
  return gridID2;
}

static
448
void gridboxstat(field_type *field1, field_type *field2, int xinc, int yinc, int statfunc)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
449
{
450
  bool useWeight = (field1->weight != NULL);
451
452
  /*
  double findex = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453

454
455
  progressInit();
  */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456

457
  int gridsize = xinc*yinc;
458
  field_type *field = (field_type*) Malloc(ompNumThreads*sizeof(field_type));
459
  for ( int i = 0; i < ompNumThreads; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
460
461
    {
      field[i].size    = gridsize;
462
      field[i].ptr     = (double*) Malloc(gridsize*sizeof(double));
463
      field[i].weight  = (useWeight) ? (double*) Malloc(gridsize*sizeof(double)) : NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464
465
466
      field[i].missval = field1->missval;
      field[i].nmiss   = 0;
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
  
468
469
  int gridID1 = field1->grid;
  int gridID2 = field2->grid;
470
471
  double *array1 = field1->ptr;
  double *array2 = field2->ptr;
472
  double missval = field1->missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473

474
475
  int nlon1 = gridInqXsize(gridID1);
  int nlat1 = gridInqYsize(gridID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476

477
478
  int nlon2 = gridInqXsize(gridID2);
  int nlat2 = gridInqYsize(gridID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
479

Uwe Schulzweida's avatar
Uwe Schulzweida committed
480
#if defined(_OPENMP)
481
#pragma omp parallel for default(shared)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482
#endif
483
  for ( int ig = 0; ig < nlat2*nlon2; ++ig )
484
    {
485
486
      int ompthID = cdo_omp_get_thread_num();

487
488
489
490
491
492
      /*
      int lprogress = 1;
#if defined(_OPENMP)
      if ( ompthID != 0 ) lprogress = 0;
#endif
#if defined(_OPENMP)
493
#include "pragma_omp_atomic_update.h"
494
495
496
497
#endif
      findex++;
      if ( lprogress ) progressStatus(0, 1, findex/nlat2*nlon2);
      */
498
499
      int ilat = ig/nlon2;
      int ilon = ig - ilat*nlon2;
500

501
      int isize = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
502
      field[ompthID].nmiss = 0;
503
      for ( int j = 0; j < yinc; ++j )
504
	{
505
	  int jj = ilat*yinc+j;
506
	  if ( jj >= nlat1 ) break;
507
	  for ( int i = 0; i < xinc; ++i )
508
	    {
509
510
	      int ii = ilon*xinc+i;
	      int index = jj*nlon1 + ii;
511
	      if ( ii >= nlon1 ) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
512
513
514
	      field[ompthID].ptr[isize] = array1[index];
	      if ( useWeight ) field[ompthID].weight[isize] = field1->weight[index];
	      if ( DBL_IS_EQUAL(field[ompthID].ptr[isize], field[ompthID].missval) ) field[ompthID].nmiss++;
515
516
517
	      isize++;
	    }
	}
Cedrick Ansorge's avatar
Cedrick Ansorge committed
518
        
Uwe Schulzweida's avatar
Uwe Schulzweida committed
519
520
      field[ompthID].size = isize;
      field2->ptr[ig] = fldfun(field[ompthID], statfunc);
521
    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
522
  
523
524
  int nmiss = 0;
  for ( int i = 0; i < nlat2*nlon2; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
525
    if ( DBL_IS_EQUAL(array2[i], missval) ) nmiss++;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
526
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
527
  field2->nmiss = nmiss;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
528
  
529
  for ( int i = 0; i < ompNumThreads; i++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
530
    {
531
532
      if ( field[i].ptr    ) Free(field[i].ptr);
      if ( field[i].weight ) Free(field[i].weight);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
533
534
    }

535
  if ( field ) Free(field);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
536
537
538
539
540
}


void *Gridboxstat(void *argument)
{
541
  int lastgrid = -1;
542
  int nrecs;
543
  int varID, levelID;
544
  bool wstatus = false;
545
  char varname[CDI_MAX_NAME];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
546
547
548
549
550

  cdoInitialize(argument);

  operatorInputArg("xinc, yinc");
  operatorCheckArgc(2);
551
552
  int xinc = parameter2int(operatorArgv()[0]);
  int yinc = parameter2int(operatorArgv()[1]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
553

554
555
556
557
558
  cdoOperatorAdd("gridboxmin",  func_min,   0, NULL);
  cdoOperatorAdd("gridboxmax",  func_max,   0, NULL);
  cdoOperatorAdd("gridboxsum",  func_sum,   0, NULL);
  cdoOperatorAdd("gridboxmean", func_meanw, 1, NULL);
  cdoOperatorAdd("gridboxavg",  func_avgw,  1, NULL);
559
560
561
562
  cdoOperatorAdd("gridboxvar",  func_varw,  1, NULL);
  cdoOperatorAdd("gridboxvar1", func_var1w, 1, NULL);
  cdoOperatorAdd("gridboxstd",  func_stdw,  1, NULL);
  cdoOperatorAdd("gridboxstd1", func_std1w, 1, NULL);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563

564
565
  int operatorID = cdoOperatorID();
  int operfunc = cdoOperatorF1(operatorID);
566
  bool needWeights = cdoOperatorF2(operatorID) != 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
567

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

570
571
  int vlistID1 = streamInqVlist(streamID1);
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
572

573
574
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
575
576
  vlistDefTaxis(vlistID2, taxisID2);

577
  int ngrids = vlistNgrids(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
578
579
  if ( ngrids > 1 )  cdoAbort("Too many different grids!");

580
  int gridID1 = vlistGrid(vlistID1, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
581
582
583
  if ( gridInqType(gridID1) == GRID_GAUSSIAN_REDUCED )
    cdoAbort("Gaussian reduced grid found. Use option -R to convert it to a regular grid!");

584
  int gridID2 = genBoxGrid(gridID1, xinc, yinc);
585
  for ( int index = 0; index < ngrids; index++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
586
587
    vlistChangeGridIndex(vlistID2, index, gridID2);

588
  int streamID2 = streamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
589
590
591

  streamDefVlist(streamID2, vlistID2);

592
  field_type field1, field2;
593
594
595
  field_init(&field1);
  field_init(&field2);

596
  int gridsize1 = gridInqSize(gridID1);
597
  field1.ptr    = (double*) Malloc(gridsize1*sizeof(double));
598
  field1.weight = needWeights ? (double*) Malloc(gridsize1*sizeof(double)) : NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
599

600
  int gridsize2 = gridInqSize(gridID2);
601
  field2.ptr    = (double*) Malloc(gridsize2*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
602
603
  field2.weight = NULL;

604
  int tsID = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
605
606
607
608
609
610
  while ( (nrecs = streamInqTimestep(streamID1, tsID)) )
    {
      taxisCopyTimestep(taxisID2, taxisID1);

      streamDefTimestep(streamID2, tsID);

611
      for ( int recID = 0; recID < nrecs; recID++ )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
612
        {
613
          int nmiss;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
614
          streamInqRecord(streamID1, &varID, &levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
615
616
          streamReadRecord(streamID1, field1.ptr, &nmiss);
          field1.nmiss = (size_t) nmiss;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
617
618
619
620
621
622
623
624
625

          field1.grid = vlistInqVarGrid(vlistID1, varID);
          field1.size = gridInqSize(field1.grid);
          field1.missval = vlistInqVarMissval(vlistID1, varID);
          
          field2.grid = gridID2;
          field2.size = gridsize2;
          field2.missval = field1.missval;

626
          if ( needWeights && field1.grid != lastgrid )
Cedrick Ansorge's avatar
Cedrick Ansorge committed
627
            {
628
	      lastgrid = field1.grid;
Cedrick Ansorge's avatar
Cedrick Ansorge committed
629
630
              wstatus = gridWeights(field1.grid, field1.weight);
            }
631
632
633
634
635
          if ( wstatus != 0 && tsID == 0 && levelID == 0 )
	    {
	      vlistInqVarName(vlistID1, varID, varname);
	      cdoWarning("Using constant grid cell area weights for variable %s!", varname);
	    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
636
637
638
639
          
          gridboxstat(&field1, &field2, xinc, yinc, operfunc);
          
          streamDefRecord(streamID2, varID,  levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
640
          streamWriteRecord(streamID2, field2.ptr, (int)field2.nmiss);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
641
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
642
643
      tsID++;
    }
Cedrick Ansorge's avatar
Cedrick Ansorge committed
644
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
645
646
  streamClose(streamID2);
  streamClose(streamID1);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
647
  
648
649
  if ( field1.ptr )    Free(field1.ptr);
  if ( field1.weight ) Free(field1.weight);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
650
  
651
652
  if ( field2.ptr )    Free(field2.ptr);
  if ( field2.weight ) Free(field2.weight);
Cedrick Ansorge's avatar
Cedrick Ansorge committed
653
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
654
  cdoFinish();
Cedrick Ansorge's avatar
Cedrick Ansorge committed
655
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
656
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
657
}