Fillmiss.cc 15.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
21
/*
   This module contains the following operators:

*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
22
#include <time.h> // clock()
23

24
#include <cdi.h>
Ralf Mueller's avatar
Ralf Mueller committed
25
#include <cdi.h>
Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
27
28
#include "cdo.h"
#include "cdo_int.h"
#include "pstream.h"
29
#include "grid.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
30

Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
#include "grid_search.h"
32
33
34
#ifdef __cplusplus
extern "C" {
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
35
#include "clipping/geometry.h"
36
37
38
#ifdef __cplusplus
}
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39

Uwe Schulzweida's avatar
Uwe Schulzweida committed
40

41
void fillmiss(field_type *field1, field_type *field2, int nfill)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
43
44
  int nx, ny, i, j;
  int nmiss2 = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
45
46
47
  int kr, ku, kl, ko;
  int ir, iu, il, io;
  int kh, kv, k1, k2, kk;
48
  int globgrid = FALSE,gridtype;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
49
50
51
52
  double s1, s2;
  double xr, xu, xl, xo;
  double **matrix1, **matrix2;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
53
54
55
56
57
  int gridID = field1->grid;
  int nmiss1 = field1->nmiss;
  double missval  = field1->missval;
  double *array1  = field1->ptr;
  double *array2  = field2->ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58

59
60
61
62
63
64
65
  nx       = gridInqXsize(gridID);
  ny       = gridInqYsize(gridID);
  globgrid = gridIsCircular(gridID);

  gridtype = gridInqType(gridID);
  if ( !(gridtype == GRID_LONLAT || gridtype == GRID_GAUSSIAN ) )
    cdoAbort("Unsupported grid type: %s!", gridNamePtr(gridtype));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
66

67
68
  matrix1 = (double **) Malloc(ny * sizeof(double *));
  matrix2 = (double **) Malloc(ny * sizeof(double *));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

  for ( j = 0; j < ny; j++ )
    {
      matrix1[j] = array1 + j*nx;
      matrix2[j] = array2 + j*nx;
    }

  for ( j = 0; j < ny; j++ )
    for ( i = 0; i < nx; i++ )
      {
	if ( DBL_IS_EQUAL(matrix1[j][i], missval) )
	  {
	    nmiss2++;

	    kr = ku = kl = ko = 0;
	    xr = xu = xl = xo = 0.;

	    for ( ir = i + 1; ir < nx; ir++ )
	      if ( !DBL_IS_EQUAL(matrix1[j][ir], missval) )
		{ kr = ir - i; xr = matrix1[j][ir]; break; }

	    if ( globgrid && ir == nx )
	      {
		for ( ir = 0; ir < i; ir++ )
		  if ( !DBL_IS_EQUAL(matrix1[j][ir], missval) )
		    { kr = nx + ir - i; xr = matrix1[j][ir]; break; }
	      }

	    for ( il = i-1; il >= 0; il-- )
	      if ( !DBL_IS_EQUAL(matrix1[j][il], missval) )
		{ kl = i - il; xl = matrix1[j][il]; break; }

	    if ( globgrid && il == -1 )
	      {
		for ( il = nx-1; il > i; il-- )
		  if ( !DBL_IS_EQUAL(matrix1[j][il], missval) )
		    { kl = nx + i - il; xl = matrix1[j][il]; break; }
	      }

	    for ( iu = j + 1; iu < ny; iu++ )
	      if ( !DBL_IS_EQUAL(matrix1[iu][i], missval) )
		{ ku = iu - j; xu = matrix1[iu][i]; break; }
	    
	    for ( io = j - 1; io >= 0; io-- )
	      if ( !DBL_IS_EQUAL(matrix1[io][i], missval) )
		{ ko = j - io; xo = matrix1[io][i]; break; }
	    
116
	    /*  printf("%d %d %d %d %d %d %g %g %g %g\n", j,i,kr,kl,ku,ko,xr,xl,xu,xo);*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
117
118
119
120
121
122
123
124
125
126
127
128
129
130

	    kh = kl + kr;
	    kv = ko + ku;
	    if      ( kh == 0 ) { s1 = 0.; k1 = 0; }
	    else if ( kl == 0 ) { s1 = xr; k1 = 1; }
	    else if ( kr == 0 ) { s1 = xl; k1 = 1; }
	    else { s1 = xr*kl/kh + xl*kr/kh; k1 = 2; }

	    if      ( kv == 0 ) { s2 = 0.; k2 = 0; }
	    else if ( ku == 0 ) { s2 = xo; k2 = 1; }
	    else if ( ko == 0 ) { s2 = xu; k2 = 1; }
	    else { s2 = xu*ko/kv + xo*ku/kv; k2 = 2; }

	    kk = k1 + k2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
132
133
134
135
136
137
138
139
	    if ( kk >= nfill )
	      {
		if      ( kk == 0 ) cdoAbort("no point found!");
		else if ( k1 == 0 ) matrix2[j][i] = s2;
		else if ( k2 == 0 ) matrix2[j][i] = s1;
		else  matrix2[j][i] = s1*k2/kk + s2*k1/kk;
	      }
	    else
	      matrix2[j][i] = matrix1[j][i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
140
141
142
143
144
145
146
147
148
149
150

	    /* matrix1[j][i] = matrix2[j][i]; */
	  }
	else
	  {
	    matrix2[j][i] = matrix1[j][i];
	  }
      }

  if ( nmiss1 != nmiss2 ) cdoAbort("found only %d of %d missing values!", nmiss2, nmiss1);

151
152
  Free(matrix2);
  Free(matrix1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
155

156
void fillmiss_one_step(field_type *field1, field_type *field2, int maxfill)
157
158
{
  int gridID, nx, ny, i, j;
159
  int nmiss2 = 0;
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
  int kr, ku, kl, ko;
  int ir, iu, il, io;
  int kh, kv, k1, k2, kk;
  double s1, s2;
  double xr, xu, xl, xo;
  double missval;
  double *array1, *array2;
  double **matrix1, **matrix2;

  gridID  = field1->grid;
  missval = field1->missval;
  array1  = field1->ptr;
  array2  = field2->ptr;

  nx  = gridInqXsize(gridID);
  ny  = gridInqYsize(gridID);

177
178
  matrix1 = (double **) Malloc(ny * sizeof(double *));
  matrix2 = (double **) Malloc(ny * sizeof(double *));
179

180
  for ( j = 0; j < ny; j++ ) { matrix1[j] = array1 + j*nx; matrix2[j] = array2 + j*nx; }
181

182
  for (int fill_iterations=0; fill_iterations < maxfill; fill_iterations++) {
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
  for ( j = 0; j < ny; j++ )
    for ( i = 0; i < nx; i++ )
      {
        if ( DBL_IS_EQUAL(matrix1[j][i], missval) )
          {
            nmiss2++;

            kr = ku = kl = ko = 0;
            xr = xu = xl = xo = 0.;

            for ( ir = i + 1; ir < nx; ir++ )
              if ( !DBL_IS_EQUAL(matrix1[j][ir], missval) )
                { kr = ir - i; xr = matrix1[j][ir]; break; }

            for ( il = i-1; il >= 0; il-- )
              if ( !DBL_IS_EQUAL(matrix1[j][il], missval) )
                { kl = i - il; xl = matrix1[j][il]; break; }


            for ( iu = j + 1; iu < ny; iu++ )
              if ( !DBL_IS_EQUAL(matrix1[iu][i], missval) )
                { ku = iu - j; xu = matrix1[iu][i]; break; }

            for ( io = j - 1; io >= 0; io-- )
              if ( !DBL_IS_EQUAL(matrix1[io][i], missval) )
                { ko = j - io; xo = matrix1[io][i]; break; }


            kh = kl + kr;
            kv = ko + ku;
            if      ( kh == 0 ) { s1 = 0.; k1 = 0; }
214
215
            else if ( kl == 0 ) { s1 = xr; k1 = kr; }
            else if ( kr == 0 ) { s1 = xl; k1 = kl; }
216
217
218
219
220
            else
              {
                if ( kl < kr )
                {
                  s1 = xl;
221
                  k1 = kl;
222
223
224
225
                }
              else
                {
                  s1 = xr;
226
                  k1 = kr;
227
228
229
230
                }
              }

            if      ( kv == 0 ) { s2 = 0.; k2 = 0; }
231
232
            else if ( ku == 0 ) { s2 = xo; k2 = ko; }
            else if ( ko == 0 ) { s2 = xu; k2 = ku; }
233
234
235
236
            else
              {
                if ( ku < ko )
                  {
237
238
                    s2 = xu;
                    k2 = ku;
239
240
241
                  }
                else
                  {
242
243
                    s2 = xo;
                    k2 = ko;
244
245
246
247
                  }
              }

            kk = k1 + k2;
248
            if      ( kk == 0 ) matrix2[j][i] = matrix1[j][i];
249
250
251
252
            else if ( k1 == 0 ) matrix2[j][i] = s2;
            else if ( k2 == 0 ) matrix2[j][i] = s1;
            else
              {
253
                if ( k1 <= k2 )
254
255
256
257
258
259
260
261
262
263
                {
                  matrix2[j][i] = s1;
                }
                else
                {
                  matrix2[j][i] = s2;
                }

              }

264
            //printf("%d %d %2d %2d %2d %2d %2g %2g %2g %2g %2g %2g %2g\n", j,i,kr,kl,ku,ko,xr,xl,xu,xo,s1,s2,matrix2[j][i]);
265
266
267
268
269
270
271
            /* matrix1[j][i] = matrix2[j][i]; */
          }
        else
          {
            matrix2[j][i] = matrix1[j][i];
          }
      }
272
273
  for ( j = 0; j < ny; j++ ) for ( i = 0; i < nx; i++ ) matrix1[j][i] = matrix2[j][i];
  }
274

275
276
  Free(matrix2);
  Free(matrix1);
277
278
}

279

280
int grid_search_nbr(struct gridsearch *gs, int num_neighbors, int *restrict nbr_add, double *restrict nbr_dist, double plon, double plat);
281
282
double nbr_compute_weights(unsigned num_neighbors, const int *restrict src_grid_mask, bool *restrict nbr_mask, const int *restrict nbr_add, double *restrict nbr_dist);
unsigned nbr_normalize_weights(unsigned num_neighbors, double dist_tot, const bool *restrict nbr_mask, int *restrict nbr_add, double *restrict nbr_dist);
283
284

static
285
void setmisstodis(field_type *field1, field_type *field2, int num_neighbors)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
287
{
  int gridID = field1->grid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
  int gridID0 = gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
289
290
291
  double missval = field1->missval;
  double *array1 = field1->ptr;
  double *array2 = field2->ptr;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
292
293
294

  unsigned gridsize = gridInqSize(gridID);

295
296
297
  unsigned nmiss = field1->nmiss;
  unsigned nvals = gridsize - nmiss;

298
299
  double *xvals = (double*) Malloc(gridsize*sizeof(double));
  double *yvals = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315

  if ( gridInqType(gridID) == GRID_GME ) gridID = gridToUnstructured(gridID, 0);

  if ( gridInqType(gridID) != GRID_UNSTRUCTURED && gridInqType(gridID) != GRID_CURVILINEAR )
    gridID = gridToCurvilinear(gridID, 0);

  gridInqXvals(gridID, xvals);
  gridInqYvals(gridID, yvals);

  /* Convert lat/lon units if required */
  char units[CDI_MAX_NAME];
  gridInqXunits(gridID, units);
  grid_to_radian(units, gridsize, xvals, "grid center lon");
  gridInqYunits(gridID, units);
  grid_to_radian(units, gridsize, yvals, "grid center lat");

316
317
  unsigned *mindex = NULL;
  unsigned *vindex = NULL;
318
319
  double *lons = NULL;
  double *lats = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
320

321
  if ( nmiss ) mindex = (unsigned *) Calloc(1, nmiss*sizeof(unsigned));
322
323
  if ( nvals )
    {
324
325
326
      vindex = (unsigned *) Calloc(1, nvals*sizeof(unsigned));
      lons   = (double *) Malloc(nvals*sizeof(double));
      lats   = (double *) Malloc(nvals*sizeof(double));
327
328
    }
  
Uwe Schulzweida's avatar
Uwe Schulzweida committed
329
  unsigned nv = 0, nm = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
330
331
  for ( unsigned i = 0; i < gridsize; ++i ) 
    {
332
      array2[i] = array1[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
333
      if ( DBL_IS_EQUAL(array1[i], missval) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
335
336
337
338
339
          mindex[nm] = i;
          nm++;
        }
      else
        {
340
341
          lons[nv] = xvals[i];
          lats[nv] = yvals[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342
343
          vindex[nv] = i;
          nv++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
344
345
        }
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
346

Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
  if ( nv != nvals ) cdoAbort("Internal problem, number of valid values differ!");
348
  
349
  bool nbr_mask[num_neighbors];   /* mask at nearest neighbors                   */
350
351
  int nbr_add[num_neighbors];     /* source address at nearest neighbors         */
  double nbr_dist[num_neighbors]; /* angular distance four nearest neighbors     */
352
353
354
355

  clock_t start, finish;
  start = clock();

356
357
  struct gridsearch *gs = NULL;

358
359
360
361
362
363
364
  if ( nmiss )
    {
      if ( num_neighbors == 1 )
        gs = gridsearch_create_nn(nvals, lons, lats);
      else
        gs = gridsearch_create(nvals, lons, lats);
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
365
366
367
  
  finish = clock();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
368
  if ( cdoVerbose ) printf("gridsearch created: %.2f seconds\n", ((double)(finish-start))/CLOCKS_PER_SEC);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
369

370
371
  progressInit();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
372
  start = clock();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373

374
  double findex = 0;
375

376
377
#pragma omp parallel for default(none) shared(findex, mindex, vindex, array1, array2, xvals, yvals, gs, nmiss, num_neighbors) \
                                      private(nbr_mask, nbr_add, nbr_dist)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
378
  for ( unsigned i = 0; i < nmiss; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
379
    {
380
381
382
#if defined(_OPENMP)
#include "pragma_omp_atomic_update.h"
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
383
384
      findex++;
      if ( cdo_omp_get_thread_num() == 0 ) progressStatus(0, 1, findex/nmiss);
385

386
      grid_search_nbr(gs, num_neighbors, nbr_add, nbr_dist, xvals[mindex[i]], yvals[mindex[i]]);
387

388
389
      /* Compute weights based on inverse distance if mask is false, eliminate those points */
      double dist_tot = nbr_compute_weights(num_neighbors, NULL, nbr_mask, nbr_add, nbr_dist);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
390

391
392
393
394
395
396
397
398
      /* Normalize weights and store the link */
      unsigned nadds = nbr_normalize_weights(num_neighbors, dist_tot, nbr_mask, nbr_add, nbr_dist);
      if ( nadds )
        {
          double result = 0;
          for ( unsigned n = 0; n < nadds; ++n ) result += array1[vindex[nbr_add[n]]]*nbr_dist[n];
          array2[mindex[i]] = result;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
399
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
400

401
402
  if ( mindex ) Free(mindex);
  if ( vindex ) Free(vindex);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
403

Uwe Schulzweida's avatar
Uwe Schulzweida committed
404
405
  finish = clock();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
406
  if ( cdoVerbose ) printf("gridsearch nearest: %.2f seconds\n", ((double)(finish-start))/CLOCKS_PER_SEC);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
407

408
  if ( gs ) gridsearch_delete(gs);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
409

Uwe Schulzweida's avatar
Uwe Schulzweida committed
410
411
  if ( gridID0 != gridID ) gridDestroy(gridID);

412
413
414
415
  if ( lons ) Free(lons);
  if ( lats ) Free(lats);
  Free(xvals);
  Free(yvals);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416
417
418
}


Uwe Schulzweida's avatar
Uwe Schulzweida committed
419
420
void *Fillmiss(void *argument)
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
  int nmiss;
422
  int nrecs, varID, levelID;
423
  void (*fill_method) (field_type *fin , field_type *fout , int) = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
424
425
426

  cdoInitialize(argument);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
427
  // clang-format off
428
429
430
431
  int FILLMISS        = cdoOperatorAdd("fillmiss"   ,   0, 0, "nfill");
  int FILLMISSONESTEP = cdoOperatorAdd("fillmiss2"  ,   0, 0, "nfill");
  int SETMISSTONN     = cdoOperatorAdd("setmisstonn" ,  0, 0, "");
  int SETMISSTODIS    = cdoOperatorAdd("setmisstodis" , 0, 0, "number of neighbors");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
432
  // clang-format on
Uwe Schulzweida's avatar
Uwe Schulzweida committed
433
434
435

  int operatorID      = cdoOperatorID();

436
  int nfill = 1;
437
438
439
440
  if ( operatorID == FILLMISS )
     {
       fill_method = &fillmiss;
     }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
  else if ( operatorID == FILLMISSONESTEP )
442
443
444
     {
       fill_method = &fillmiss_one_step;
     }
445
  else if ( operatorID == SETMISSTONN )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
446
     {
447
448
449
450
451
452
       fill_method = &setmisstodis;
     }
  else if ( operatorID == SETMISSTODIS )
     {
       nfill = 4;
       fill_method = &setmisstodis;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
453
     }
454
455

  /* Argument handling */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
456
457
458
459
460
461
  {
    int oargc = operatorArgc();
    char **oargv = operatorArgv();

    if ( oargc == 1 )
      {
462
463
        nfill = parameter2int(oargv[0]);
        if ( operatorID == FILLMISS ) 
464
465
466
          {
            if ( nfill < 1 || nfill > 4 ) cdoAbort("nfill out of range!");
          }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
467
468
469
470
471
      }
    else if ( oargc > 1 )
      cdoAbort("Too many arguments!");
  }

472
  int streamID1 = pstreamOpenRead(cdoStreamName(0));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
473

474
  int vlistID1 = pstreamInqVlist(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475
  int vlistID2 = vlistDuplicate(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
476

Uwe Schulzweida's avatar
Uwe Schulzweida committed
477
478
  int taxisID1 = vlistInqTaxis(vlistID1);
  int taxisID2 = taxisDuplicate(taxisID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
479
480
  vlistDefTaxis(vlistID2, taxisID2);

481
  int streamID2 = pstreamOpenWrite(cdoStreamName(1), cdoFiletype());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
482

483
  pstreamDefVlist(streamID2, vlistID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
484

Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
  int gridsize = vlistGridsizeMax(vlistID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486

487
  field_type field1, field2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
488
489
  field_init(&field1);
  field_init(&field2);
490
491
  field1.ptr = (double*) Malloc(gridsize*sizeof(double));
  field2.ptr = (double*) Malloc(gridsize*sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
492

Uwe Schulzweida's avatar
Uwe Schulzweida committed
493
  int tsID = 0;
494
  while ( (nrecs = pstreamInqTimestep(streamID1, tsID)) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
495
496
497
    {
      taxisCopyTimestep(taxisID2, taxisID1);

498
      pstreamDefTimestep(streamID2, tsID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
499
	       
500
      for ( int recID = 0; recID < nrecs; recID++ )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
501
	{
502
503
	  pstreamInqRecord(streamID1, &varID, &levelID);
	  pstreamReadRecord(streamID1, field1.ptr, &nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
504
          field1.nmiss = (size_t) nmiss;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
505

506
	  pstreamDefRecord(streamID2, varID, levelID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
507
508
509

          if ( field1.nmiss == 0 )
            {
510
              pstreamWriteRecord(streamID2, field1.ptr, 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
511
512
513
            }
          else
            {
514
              int gridID = vlistInqVarGrid(vlistID1, varID);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
515

516
517
              if ( (operatorID == FILLMISS || operatorID == FILLMISSONESTEP) && 
                   (gridInqType(gridID) == GRID_GME || gridInqType(gridID) == GRID_UNSTRUCTURED) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
518
519
520
521
522
523
524
525
526
527
528
                cdoAbort("%s data unsupported!", gridNamePtr(gridInqType(gridID)) );
                
              field1.grid    = gridID;
              field1.missval = vlistInqVarMissval(vlistID1, varID);

              field2.grid    = field1.grid;
              field2.nmiss   = 0;
              field2.missval = field1.missval;

              fill_method(&field1, &field2, nfill);

529
530
531
              int gridsize = gridInqSize(field2.grid);
              int nmiss = 0;
              for ( int i = 0; i < gridsize; ++i )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
532
533
                if ( DBL_IS_EQUAL(field2.ptr[i], field2.missval) ) nmiss++;
              
534
              pstreamWriteRecord(streamID2, field2.ptr, nmiss);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
535
536
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
537
538
539
      tsID++;
    }

540
541
  pstreamClose(streamID2);
  pstreamClose(streamID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542

543
544
  if ( field2.ptr ) Free(field2.ptr);
  if ( field1.ptr ) Free(field1.ptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
545
546
547

  cdoFinish();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
548
  return 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
}