remap_bilinear_scrip.c 13.7 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
int find_ij_weights(double plon, double plat, double* restrict src_lats, double* restrict src_lons, double *ig, double *jg)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
16
{
17
  int lfound = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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
  long iter;                     /*  iteration counters   */
  double iguess, jguess;         /*  current guess for bilinear coordinate  */
  double deli, delj;             /*  corrections to iw,jw                   */
  double dth1, dth2, dth3;       /*  some latitude  differences             */
  double dph1, dph2, dph3;       /*  some longitude differences             */
  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  */

  dth1 = src_lats[1] - src_lats[0];
  dth2 = src_lats[3] - src_lats[0];
  dth3 = src_lats[2] - src_lats[1] - dth2;

  dph1 = src_lons[1] - src_lons[0];
  dph2 = src_lons[3] - src_lons[0];
  dph3 = src_lons[2] - src_lons[1];

  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;

  iguess = HALF;
  jguess = HALF;

  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
81
82
83
84
85
  if ( iter < remap_max_iter ) lfound = 1;

  return (lfound);
}

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

Uwe Schulzweida's avatar
Uwe Schulzweida committed
95
int num_src_points(const int* restrict mask, const int src_add[4], double src_lats[4])
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
{
  int icount = 0;

  for ( int n = 0; n < 4; ++n )
    {
      if ( mask[src_add[n]] )
	icount++;
      else
	src_lats[n] = 0.;
    }

  return (icount);
}

static
Uwe Schulzweida's avatar
Uwe Schulzweida committed
111
void renormalize_weights(const double src_lats[4], double wgts[4])
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
{
  int n;
  double sum_wgts = 0.0; /* sum of weights for normalization */
  /* 2012-05-08 Uwe Schulzweida: using absolute value of src_lats (bug fix) */
  for ( n = 0; n < 4; ++n ) sum_wgts += fabs(src_lats[n]);
  for ( n = 0; n < 4; ++n ) wgts[n] = fabs(src_lats[n])/sum_wgts;
}

static
void bilinear_warning(double plon, double plat, double iw, double jw, int* src_add, double* src_lons, double* src_lats, remapgrid_t* src_grid)
{
  static int lwarn = TRUE;

  if ( cdoVerbose )
    {
      cdoPrint("Point coords: %g %g", plat, plon);
      cdoPrint("Src grid lats: %g %g %g %g", src_lats[0], src_lats[1], src_lats[2], src_lats[3]);
      cdoPrint("Src grid lons: %g %g %g %g", src_lons[0], src_lons[1], src_lons[2], src_lons[3]);
      cdoPrint("Src grid addresses: %d %d %d %d", src_add[0], src_add[1], src_add[2], src_add[3]);
      cdoPrint("Src grid lats: %g %g %g %g",
	       src_grid->cell_center_lat[src_add[0]], src_grid->cell_center_lat[src_add[1]],
	       src_grid->cell_center_lat[src_add[2]], src_grid->cell_center_lat[src_add[3]]);
      cdoPrint("Src grid lons: %g %g %g %g",
	       src_grid->cell_center_lon[src_add[0]], src_grid->cell_center_lon[src_add[1]],
	       src_grid->cell_center_lon[src_add[2]], src_grid->cell_center_lon[src_add[3]]);
      cdoPrint("Current iw,jw : %g %g", iw, jw);
    }

  if ( cdoVerbose || lwarn )
    {
      lwarn = FALSE;
      //  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
Uwe Schulzweida's avatar
Uwe Schulzweida committed
149
void bilinear_remap(double* restrict tgt_point, const double* restrict src_array, const double wgts[4], const int src_add[4])
150
151
152
{
  // *tgt_point = 0.;
  // for ( int n = 0; n < 4; ++n ) *tgt_point += src_array[src_add[n]]*wgts[n];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
153
154
  *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
155
156
157
158
159
160
161
162
163
}

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

  This routine computes the weights for a bilinear interpolation.

  -----------------------------------------------------------------------
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164
void scrip_remap_bilinear_weights(remapgrid_t *src_grid, remapgrid_t *tgt_grid, remapvars_t *rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
166
167
{
  /*   Local variables */
  int  search_result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168
  long tgt_cell_add;             /*  destination addresss                   */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
  int 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      */
  double plat, plon;             /*  lat/lon coords of destination point    */
  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"); 

188
189
190
191
192
193
  long tgt_grid_size = tgt_grid->size;

  weightlinks_t *weightlinks = (weightlinks_t *) malloc(tgt_grid_size*sizeof(weightlinks_t));

  double findex = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195
196
197
  /* Loop over destination grid */

#if defined(_OPENMP)
#pragma omp parallel for default(none) \
198
  shared(ompNumThreads, cdoVerbose, weightlinks, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, rv, findex) \
199
  private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result)    \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200
  schedule(static)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
201
#endif
202
  for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
204
    {
      int lprogress = 1;
205
206
      if ( cdo_omp_get_thread_num() != 0 ) lprogress = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
207
#if defined(_OPENMP)
208
#include "pragma_omp_atomic_update.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
209
210
211
212
#endif
      findex++;
      if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);

213
214
215
      weightlinks[tgt_cell_add].nlinks = 0;	

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

217
218
      plat = tgt_grid->cell_center_lat[tgt_cell_add];
      plon = tgt_grid->cell_center_lon[tgt_cell_add];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233

      /* Find nearest square of grid points on source grid  */
      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);

      /* Check to see if points are land points */
      if ( search_result > 0 )
	{
234
	  for ( int n = 0; n < 4; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
235
236
237
238
239
240
241
242
	    if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
	}

      /* If point found, find local iw,jw coordinates for weights  */
      if ( search_result > 0 )
	{
	  double iw, jw;  /*  current guess for bilinear coordinate  */

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

245
          if ( find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
247
	    {
	      /* Successfully found iw,jw - compute weights */
248
	      set_bilinear_weights(iw, jw, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
249

250
	      store_weightlinks(4, src_add, wgts, tgt_cell_add, weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
252
253
	    }
          else
	    {
254
	      bilinear_warning(plon, plat, iw, jw, src_add, src_lons, src_lats, src_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
256
257
258
259
260
261
262
263
264
265

	      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 )
	{
266
          if ( num_src_points(src_grid->mask, src_add, src_lats) > 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
267
	    {
268
	      renormalize_weights(src_lats, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269

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

272
	      store_weightlinks(4, src_add, wgts, tgt_cell_add, weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
273
274
	    }
        }
275
276
277
278
279
    }

  weightlinks2remaplinks(tgt_grid_size, weightlinks, rv);

  if ( weightlinks ) free(weightlinks);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280
281

  if ( cdoTimer ) timer_stop(timer_remap_bil);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
282
283
} /* scrip_remap_weights_bilinear */

284
285
286
287
/*
  -----------------------------------------------------------------------

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

289
290
  -----------------------------------------------------------------------
*/
291
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
292
293
294
{
  /*   Local variables */
  int  search_result;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
295
  long tgt_cell_add;             /*  destination addresss                   */
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296
297
298
299
300
301
302
303
304
305
306
307
308
309
  int 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      */
  double plat, plon;             /*  lat/lon coords of destination point    */
  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();

310
  long tgt_grid_size = tgt_grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
311
312
313
314
315
316

  /* Compute mappings from source to target grid */

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

317
318
  double findex = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
319
320
321
322
  /* Loop over destination grid */

#if defined(_OPENMP)
#pragma omp parallel for default(none) \
323
  shared(ompNumThreads, cdoVerbose, cdoSilentMode, remap_grid_type, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, findex) \
324
  private(tgt_cell_add, src_add, src_lats, src_lons, wgts, plat, plon, search_result)    \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
325
326
  schedule(static)
#endif
327
  for ( tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
328
329
    {
      int lprogress = 1;
330
      if ( cdo_omp_get_thread_num() != 0 ) lprogress = 0;
331
332
333

      if ( !cdoSilentMode )
	{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
#if defined(_OPENMP)
335
#include "pragma_omp_atomic_update.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
336
#endif
337
338
339
	  findex++;
	  if ( lprogress ) progressStatus(0, 1, findex/tgt_grid_size);
	}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
340

341
      tgt_array[tgt_cell_add] = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
342

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

345
346
347
348
349
350
351
352
353
354
355
356
357
358
      if ( tgt_grid->remap_grid_type == REMAP_GRID_TYPE_REG2D )
        {
          long nx = tgt_grid->dims[0];
	  long iy = tgt_cell_add/nx;
	  long ix = tgt_cell_add - iy*nx;

          plat = tgt_grid->reg2d_center_lat[iy];
          plon = tgt_grid->reg2d_center_lon[ix];
        }
      else
        {
          plat = tgt_grid->cell_center_lat[tgt_cell_add];
          plon = tgt_grid->cell_center_lon[tgt_cell_add];
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373

      /* Find nearest square of grid points on source grid  */
      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);

      /* Check to see if points are land points */
      if ( search_result > 0 )
	{
374
	  for ( int n = 0; n < 4; ++n )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
375
376
377
378
379
380
381
382
	    if ( ! src_grid->mask[src_add[n]] ) search_result = 0;
	}

      /* If point found, find local iw,jw coordinates for weights  */
      if ( search_result > 0 )
	{
	  double iw, jw;  /*  current guess for bilinear coordinate  */

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

385
          if ( find_ij_weights(plon, plat, src_lats, src_lons, &iw, &jw) )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
386
387
	    {
	      /* Successfully found iw,jw - compute weights */
388
	      set_bilinear_weights(iw, jw, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
389

390
	      sort_add_and_wgts(4, src_add, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391

392
	      bilinear_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
393
394
395
	    }
          else
	    {
396
	      bilinear_warning(plon, plat, iw, jw, src_add, src_lons, src_lats, src_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
397
398
399
400
401
402
403
404
405
406
407

	      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 )
	{
408
          if ( num_src_points(src_grid->mask, src_add, src_lats) > 0 )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
	    {
410
	      renormalize_weights(src_lats, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411

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

414
	      sort_add_and_wgts(4, src_add, wgts);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415

416
	      bilinear_remap(&tgt_array[tgt_cell_add], src_array, wgts, src_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
417
418
	    }
        }
419
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
420
421
422

  if ( cdoTimer ) timer_stop(timer_remap_bil);
} /* scrip_remap_bilinear */