remap_bilinear_scrip.cc 16 KB
Newer Older
Uwe Schulzweida's avatar
Uwe Schulzweida committed
1
2
3
4
#include "cdo.h"
#include "cdo_int.h"
#include "grid.h"
#include "remap.h"
5
#include "remap_store_link.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
6
7
8
9
10
11
12
13
14


/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */
/*                                                                         */
/*      BILINEAR INTERPOLATION                                             */
/*                                                                         */
/* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! */


15
bool find_ij_weights(double plon, double plat, double *restrict src_lons, double *restrict src_lats, double *ig, double *jg)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16
{
17
  bool lfound = false;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
19
20
21
22
23
24
25
26
27
  long iter;                     /*  iteration counters   */
  double deli, delj;             /*  corrections to iw,jw                   */
  double dthp, dphp;             /*  difference between point and sw corner */
  double mat1, mat2, mat3, mat4; /*  matrix elements                        */
  double determinant;            /*  matrix determinant                     */
  double converge = 1.e-10;      /* Convergence criterion                   */
  extern long remap_max_iter;

  /* Iterate to find iw,jw for bilinear approximation  */

28
29
30
31
  // some latitude  differences 
  double dth1 = src_lats[1] - src_lats[0];
  double dth2 = src_lats[3] - src_lats[0];
  double dth3 = src_lats[2] - src_lats[1] - dth2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
32

33
34
35
36
  // some longitude differences
  double dph1 = src_lons[1] - src_lons[0];
  double dph2 = src_lons[3] - src_lons[0];
  double dph3 = src_lons[2] - src_lons[1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
38
39
40
41
42
43
44
45
46

  if ( dph1 >  THREE*PIH ) dph1 -= PI2;
  if ( dph2 >  THREE*PIH ) dph2 -= PI2;
  if ( dph3 >  THREE*PIH ) dph3 -= PI2;
  if ( dph1 < -THREE*PIH ) dph1 += PI2;
  if ( dph2 < -THREE*PIH ) dph2 += PI2;
  if ( dph3 < -THREE*PIH ) dph3 += PI2;

  dph3 = dph3 - dph2;

47
48
49
  // current guess for bilinear coordinate
  double iguess = HALF;
  double jguess = HALF;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

  for ( iter = 0; iter < remap_max_iter; ++iter )
    {
      dthp = plat - src_lats[0] - dth1*iguess - dth2*jguess - dth3*iguess*jguess;
      dphp = plon - src_lons[0];
      
      if ( dphp >  THREE*PIH ) dphp -= PI2;
      if ( dphp < -THREE*PIH ) dphp += PI2;

      dphp = dphp - dph1*iguess - dph2*jguess - dph3*iguess*jguess;

      mat1 = dth1 + dth3*jguess;
      mat2 = dth2 + dth3*iguess;
      mat3 = dph1 + dph3*jguess;
      mat4 = dph2 + dph3*iguess;

      determinant = mat1*mat4 - mat2*mat3;

      deli = (dthp*mat4 - dphp*mat2)/determinant;
      delj = (dphp*mat1 - dthp*mat3)/determinant;

      if ( fabs(deli) < converge && fabs(delj) < converge ) break;

      iguess += deli;
      jguess += delj;
    }

  *ig = iguess;
  *jg = jguess;

80
  if ( iter < remap_max_iter ) lfound = true;
81

82
  return lfound;
83
84
85
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
void set_bilinear_weights(double iw, double jw, double wgts[4])
87
88
89
90
91
92
93
{
  wgts[0] = (1.-iw) * (1.-jw);
  wgts[1] =     iw  * (1.-jw);
  wgts[2] =     iw  *     jw;
  wgts[3] = (1.-iw) *     jw;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
94

95
unsigned num_src_points(const int* restrict mask, const size_t src_add[4], double src_lats[4])
96
{
97
  unsigned icount = 0;
98

99
  for ( unsigned n = 0; n < 4; ++n )
100
101
102
103
104
105
106
    {
      if ( mask[src_add[n]] )
	icount++;
      else
	src_lats[n] = 0.;
    }

107
  return icount;
108
109
110
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
void renormalize_weights(const double src_lats[4], double wgts[4])
112
113
114
{
  double sum_wgts = 0.0; /* sum of weights for normalization */
  /* 2012-05-08 Uwe Schulzweida: using absolute value of src_lats (bug fix) */
115
116
  for ( unsigned n = 0; n < 4; ++n ) sum_wgts += fabs(src_lats[n]);
  for ( unsigned n = 0; n < 4; ++n ) wgts[n] = fabs(src_lats[n])/sum_wgts;
117
118
119
}

static
120
void bilinear_warning(double plon, double plat, double iw, double jw, size_t* src_add, double* src_lons, double* src_lats, remapgrid_t* src_grid)
121
{
122
  static bool lwarn = true;
123
124
125

  if ( cdoVerbose )
    {
126
127
128
      cdoPrint("Point coords: %g %g", plat*RAD2DEG, plon*RAD2DEG);
      cdoPrint("Src grid lats: %g %g %g %g", src_lats[0]*RAD2DEG, src_lats[1]*RAD2DEG, src_lats[2]*RAD2DEG, src_lats[3]*RAD2DEG);
      cdoPrint("Src grid lons: %g %g %g %g", src_lons[0]*RAD2DEG, src_lons[1]*RAD2DEG, src_lons[2]*RAD2DEG, src_lons[3]*RAD2DEG);
129
      cdoPrint("Src grid addresses: %zu %zu %zu %zu", src_add[0], src_add[1], src_add[2], src_add[3]);
130
      cdoPrint("Src grid lats: %g %g %g %g",
131
132
	       src_grid->cell_center_lat[src_add[0]]*RAD2DEG, src_grid->cell_center_lat[src_add[1]]*RAD2DEG,
	       src_grid->cell_center_lat[src_add[2]]*RAD2DEG, src_grid->cell_center_lat[src_add[3]]*RAD2DEG);
133
      cdoPrint("Src grid lons: %g %g %g %g",
134
135
	       src_grid->cell_center_lon[src_add[0]]*RAD2DEG, src_grid->cell_center_lon[src_add[1]]*RAD2DEG,
	       src_grid->cell_center_lon[src_add[2]]*RAD2DEG, src_grid->cell_center_lon[src_add[3]]*RAD2DEG);
136
137
138
139
140
      cdoPrint("Current iw,jw : %g %g", iw, jw);
    }

  if ( cdoVerbose || lwarn )
    {
141
      lwarn = false;
142
143
144
145
146
147
      //  cdoWarning("Iteration for iw,jw exceed max iteration count of %d!", remap_max_iter);
      cdoWarning("Bilinear interpolation failed for some grid points - used a distance-weighted average instead!");
    }
}

static
148
void bilinear_remap(double* restrict tgt_point, const double *restrict src_array, const double wgts[4], const size_t src_add[4])
149
150
{
  // *tgt_point = 0.;
151
  // for ( unsigned n = 0; n < 4; ++n ) *tgt_point += src_array[src_add[n]]*wgts[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
152
153
  *tgt_point = src_array[src_add[0]]*wgts[0] + src_array[src_add[1]]*wgts[1]
             + src_array[src_add[2]]*wgts[2] + src_array[src_add[3]]*wgts[3];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
155
156
157
158
159
160
161
162
}

/*
  -----------------------------------------------------------------------

  This routine computes the weights for a bilinear interpolation.

  -----------------------------------------------------------------------
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
163
void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
{
  extern int timer_remap_bil;
  int remap_grid_type = src_grid->remap_grid_type;

  if ( cdoVerbose ) cdoPrint("Called %s()", __func__);

  if ( cdoTimer ) timer_start(timer_remap_bil);

  progressInit();

  /* Compute mappings from source to target grid */

  if ( src_grid->rank != 2 )
    cdoAbort("Can not do bilinear interpolation when source grid rank != 2"); 

179
  size_t tgt_grid_size = tgt_grid->size;
180

181
  weightlinks_t *weightlinks = (weightlinks_t *) Malloc(tgt_grid_size*sizeof(weightlinks_t));
182
  weightlinks[0].addweights = (addweight_t *) Malloc(4*tgt_grid_size*sizeof(addweight_t));
183
  for ( size_t tgt_cell_add = 1; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
184
    weightlinks[tgt_cell_add].addweights = weightlinks[0].addweights + 4*tgt_cell_add;
185
186
187

  double findex = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
189
  /* Loop over destination grid */

190
191
192
#ifdef  HAVE_OPENMP4
#pragma omp parallel for default(none)  schedule(static)  reduction(+:findex) \
  shared(weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
#endif
194
  for ( size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
195
196
    {
      findex++;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
197
      if ( cdo_omp_get_thread_num() == 0 ) progressStatus(0, 1, findex/tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198

199
200
201
      weightlinks[tgt_cell_add].nlinks = 0;	

      if ( ! tgt_grid->mask[tgt_cell_add] ) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202

203
204
      double plon = 0, plat = 0;
      remapgrid_get_lonlat(tgt_grid, tgt_cell_add, &plon, &plat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205

206
207
208
209
210
      size_t src_add[4];    //  address for the four source points
      double src_lats[4];   //  latitudes  of four bilinear corners
      double src_lons[4];   //  longitudes of four bilinear corners
      double wgts[4];       //  bilinear weights for four corners

211
      // Find nearest square of grid points on source grid
212
      int search_result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213
214
215
216
217
218
219
220
221
222
      if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
	search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons, 
					  plat, plon, src_grid->dims,
					  src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
      else
	search_result = grid_search(src_grid, src_add, src_lats, src_lons, 
				    plat, plon, src_grid->dims,
				    src_grid->cell_center_lat, src_grid->cell_center_lon,
				    src_grid->cell_bound_box, src_grid->bin_addr);

223
      // Check to see if points are mask points
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
225
      if ( search_result > 0 )
	{
226
	  for ( unsigned n = 0; n < 4; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
228
229
	    if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
	}

230
      // If point found, find local iw,jw coordinates for weights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
232
      if ( search_result > 0 )
	{
233
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234

235
	  double iw, jw;  // current guess for bilinear coordinate 
236
          if ( find_ij_weights(plon, plat, src_lons, src_lats, &iw, &jw) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
237
	    {
238
	      // Successfully found iw,jw - compute weights
239
	      set_bilinear_weights(iw, jw, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
240

241
	      store_weightlinks(0, 4, src_add, wgts, tgt_cell_add, weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
243
244
	    }
          else
	    {
245
	      bilinear_warning(plon, plat, iw, jw, src_add, src_lons, src_lats, src_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
247
248
249
250
251
252
253
254
255
256

	      search_result = -1;
	    }
	}

      /*
	Search for bilinear failed - use a distance-weighted average instead (this is typically near the pole)
	Distance was stored in src_lats!
      */
      if ( search_result < 0 )
	{
257
          if ( num_src_points(src_grid->mask, src_add, src_lats) > 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
	    {
259
	      renormalize_weights(src_lats, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
260

261
	      tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
262

263
	      store_weightlinks(0, 4, src_add, wgts, tgt_cell_add, weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
265
	    }
        }
266
267
    }

268
  weightlinks2remaplinks(0, tgt_grid_size, weightlinks, rv);
269

270
  if ( weightlinks ) Free(weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
271
272

  if ( cdoTimer ) timer_stop(timer_remap_bil);
273
} // scrip_remap_weights_bilinear
Uwe Schulzweida's avatar
Uwe Schulzweida committed
274

275
276
277
278
/*
  -----------------------------------------------------------------------

  This routine computes and apply the weights for a bilinear interpolation.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
279

280
281
  -----------------------------------------------------------------------
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
282
283
284
285
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
320
321
322
323
324
325
326
327
328
329
330
331

//#define TEST_KDTREE
#ifdef  TEST_KDTREE
#include "grid_search.h"
int grid_search_test(struct gridsearch *gs, size_t *restrict src_add, double *restrict src_lats, 
                     double *restrict src_lons,  double plat, double plon, const size_t *restrict src_grid_dims,
                     double *restrict src_center_lat, double *restrict src_center_lon)
{
  /*
    Output variables:

    int    src_add[4]              ! address of each corner point enclosing P
    double src_lats[4]             ! latitudes  of the four corner points
    double src_lons[4]             ! longitudes of the four corner points

    Input variables:

    double plat                    ! latitude  of the search point
    double plon                    ! longitude of the search point

    int src_grid_dims[2]           ! size of each src grid dimension

  */
  bool is_cyclic = true;
  int search_result = 0;

  for ( unsigned n = 0; n < 4; ++n ) src_add[n] = 0;
 
  /* Now perform a more detailed search */

  size_t nx = src_grid_dims[0];
  size_t ny = src_grid_dims[1];

  double search_radius = gs->search_radius;
  const double range0 = SQR(search_radius);
  double range = range0;
  size_t add = gridsearch_nearest(gs, plon, plat, &range);
  // printf("plon, plat, add, range %g %g %g %g %zu %g\n", plon*RAD2DEG, plat*RAD2DEG,
  //     src_center_lon[add]*RAD2DEG, src_center_lat[add]*RAD2DEG,add, range);
  if ( add != GS_NOT_FOUND )
    {
      size_t idx[4];
      for ( unsigned k = 0; k < 4; ++k )
        {
          /* Determine neighbor addresses */
          size_t j = add/nx;
          size_t i = add - j*nx;
          if ( k == 1 || k == 3 )  i = (i > 0) ? i - 1 : (is_cyclic) ? nx-1 : 0;
          if ( k == 2 || k == 3 )  j = (j > 0) ? j - 1 : 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
332
333
          if ( point_in_quad(is_cyclic, nx, ny, i, j, src_add, src_lons, src_lats,
                             plon, plat, src_center_lon, src_center_lat) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
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
	    {
	      search_result = 1;
	      return search_result;
	    }
	  /* Otherwise move on to next cell */
        }
      /*
        If no cell found, point is likely either in a box that straddles either pole or is outside 
        the grid. Fall back to a distance-weighted average of the four closest points.
        Go ahead and compute weights here, but store in src_lats and return -add to prevent the 
        parent routine from computing bilinear weights.
      */
      // if ( !src_grid->lextrapolate ) return search_result;

      /*
        printf("Could not find location for %g %g\n", plat*RAD2DEG, plon*RAD2DEG);
        printf("Using nearest-neighbor for this point\n");
      */
      search_result = add;
    }

  return search_result;
}  /* grid_search_test */
#endif

359
void scrip_remap_bilinear(remapgrid_t *src_grid, remapgrid_t *tgt_grid, const double *restrict src_array, double *restrict tgt_array, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
360
361
362
363
364
365
366
367
368
369
{
  extern int timer_remap_bil;
  int remap_grid_type = src_grid->remap_grid_type;

  if ( cdoVerbose ) cdoPrint("Called %s()", __func__);

  if ( cdoTimer ) timer_start(timer_remap_bil);

  progressInit();

Uwe Schulzweida's avatar
Uwe Schulzweida committed
370
#ifdef  TEST_KDTREE
371
372
  bool xIsCyclic = false;
  size_t dims[2] = {src_grid->size, 0};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
373
  struct gridsearch *gs = NULL;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
374
  if ( remap_grid_type != REMAP_GRID_TYPE_REG2D )
375
    gs = gridsearch_create_nn(xIsCyclic, dims, src_grid->size, src_grid->cell_center_lon, src_grid->cell_center_lat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
376
377
#endif

378
  size_t tgt_grid_size = tgt_grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
379
380
381
382
383
384

  /* Compute mappings from source to target grid */

  if ( src_grid->rank != 2 )
    cdoAbort("Can not do bilinear interpolation when source grid rank != 2"); 

385
386
  double findex = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
387
388
  /* Loop over destination grid */

389
390
391
#ifdef  HAVE_OPENMP4
#pragma omp parallel for default(none)  schedule(static)  reduction(+:findex) \
  shared(cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
392
#endif
393
  for ( size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
394
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
396
      findex++;
      if ( cdo_omp_get_thread_num() == 0 ) progressStatus(0, 1, findex/tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397

398
      tgt_array[tgt_cell_add] = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
399

400
      if ( ! tgt_grid->mask[tgt_cell_add] ) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
401

402
403
      double plon = 0, plat = 0;
      remapgrid_get_lonlat(tgt_grid, tgt_cell_add, &plon, &plat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
404

405
406
407
408
409
      size_t src_add[4];    //  address for the four source points
      double src_lats[4];   //  latitudes  of four bilinear corners
      double src_lons[4];   //  longitudes of four bilinear corners
      double wgts[4];       //  bilinear weights for four corners

410
      // Find nearest square of grid points on source grid
411
      int search_result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
412
413
414
415
416
      if ( remap_grid_type == REMAP_GRID_TYPE_REG2D )
	search_result = grid_search_reg2d(src_grid, src_add, src_lats, src_lons, 
					  plat, plon, src_grid->dims,
					  src_grid->reg2d_center_lat, src_grid->reg2d_center_lon);
      else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
418
419
420
#ifdef  TEST_KDTREE
        search_result = grid_search_test(gs, src_add, src_lats, src_lons, plat, plon, src_grid->dims,
                                         src_grid->cell_center_lat, src_grid->cell_center_lon);
#else
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
422
423
424
	search_result = grid_search(src_grid, src_add, src_lats, src_lons, 
				    plat, plon, src_grid->dims,
				    src_grid->cell_center_lat, src_grid->cell_center_lon,
				    src_grid->cell_bound_box, src_grid->bin_addr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
425
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
426

427
      // Check to see if points are mask points
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428
429
      if ( search_result > 0 )
	{
430
	  for ( unsigned n = 0; n < 4; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
431
432
433
	    if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
	}

434
      // If point found, find local iw,jw coordinates for weights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
435
436
      if ( search_result > 0 )
	{
437
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
438

439
	  double iw, jw;  // current guess for bilinear coordinate
440
          if ( find_ij_weights(plon, plat, src_lons, src_lats, &iw, &jw) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
441
	    {
442
	      // Successfully found iw,jw - compute weights
443
	      set_bilinear_weights(iw, jw, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444

445
	      sort_add_and_wgts(4, src_add, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
446

447
	      bilinear_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
448
449
450
	    }
          else
	    {
451
	      bilinear_warning(plon, plat, iw, jw, src_add, src_lons, src_lats, src_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
452
453
454
455
456
457
458
459
460
461
462

	      search_result = -1;
	    }
	}

      /*
	Search for bilinear failed - use a distance-weighted average instead (this is typically near the pole)
	Distance was stored in src_lats!
      */
      if ( search_result < 0 )
	{
463
          if ( num_src_points(src_grid->mask, src_add, src_lats) > 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
464
	    {
465
	      renormalize_weights(src_lats, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
466

467
	      tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
468

469
	      sort_add_and_wgts(4, src_add, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
470

471
	      bilinear_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
472
473
	    }
        }
474
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
475

Uwe Schulzweida's avatar
Uwe Schulzweida committed
476
477
478
479
#ifdef  TEST_KDTREE
  if ( gs ) gridsearch_delete(gs);
#endif

Uwe Schulzweida's avatar
Uwe Schulzweida committed
480
  if ( cdoTimer ) timer_stop(timer_remap_bil);
481
} // scrip_remap_bilinear