remap_bilinear.cc 11.1 KB
Newer Older
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-2020 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
6
7
8
9
10
11
12
13
14
15
16
  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.
*/
Oliver Heidmann's avatar
Oliver Heidmann committed
17

18
#include "process_int.h"
19
#include "cdo_wtime.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
#include "remap.h"
21
#include "remap_store_link.h"
22
#include "cdo_options.h"
23
#include "progress.h"
24
#include "cimdOmp.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25

Uwe Schulzweida's avatar
Uwe Schulzweida committed
26
// bilinear interpolation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27

28
29
30
31
32
33
34
static inline void
limit_dphi_bounds(double &dphi)
{
  if (dphi > 3.0 * PIH) dphi -= PI2;
  if (dphi < -3.0 * PIH) dphi += PI2;
}

35
bool
36
remapFindWeights(const LonLatPoint &llpoint, const double (&src_lons)[4], const double (&src_lats)[4], double *ig, double *jg)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
37
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
  constexpr double converge = 1.e-10;  // Convergence criterion
Uwe Schulzweida's avatar
Uwe Schulzweida committed
39
40
  extern long remap_max_iter;

41
  // Iterate to find iw,jw for bilinear approximation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
42

43
  // some latitude  differences
44
45
46
  const auto dth1 = src_lats[1] - src_lats[0];
  const auto dth2 = src_lats[3] - src_lats[0];
  const auto dth3 = src_lats[2] - src_lats[1] - dth2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
47

48
  // some longitude differences
49
50
51
  auto dph1 = src_lons[1] - src_lons[0];
  auto dph2 = src_lons[3] - src_lons[0];
  auto dph3 = src_lons[2] - src_lons[1];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
52

53
54
55
  limit_dphi_bounds(dph1);
  limit_dphi_bounds(dph2);
  limit_dphi_bounds(dph3);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
56
57
58

  dph3 = dph3 - dph2;

59
  // current guess for bilinear coordinate
Uwe Schulzweida's avatar
Uwe Schulzweida committed
60
61
  double iguess = 0.5;
  double jguess = 0.5;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
62

63
  long iter = 0;  // iteration counters
64
  for (iter = 0; iter < remap_max_iter; ++iter)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
65
    {
66
67
      const auto dthp = llpoint.lat - src_lats[0] - dth1 * iguess - dth2 * jguess - dth3 * iguess * jguess;
      auto dphp = llpoint.lon - src_lons[0];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68

69
      limit_dphi_bounds(dphp);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70

71
      dphp = dphp - dph1 * iguess - dph2 * jguess - dph3 * iguess * jguess;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
72

73
74
75
76
      const auto mat1 = dth1 + dth3 * jguess;
      const auto mat2 = dth2 + dth3 * iguess;
      const auto mat3 = dph1 + dph3 * jguess;
      const auto mat4 = dph2 + dph3 * iguess;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
77

78
      const auto determinant = mat1 * mat4 - mat2 * mat3;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
79

80
81
      const auto deli = (dthp * mat4 - dphp * mat2) / determinant;
      const auto delj = (dphp * mat1 - dthp * mat3) / determinant;
82

83
      if (std::fabs(deli) < converge && std::fabs(delj) < converge) break;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
84
85
86
87
88
89
90
91

      iguess += deli;
      jguess += delj;
    }

  *ig = iguess;
  *jg = jguess;

92
  return (iter < remap_max_iter);
93
94
}

95
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
96
bilinearSetWeights(double iw, double jw, double (&weights)[4])
97
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
  // clang-format off
Uwe Schulzweida's avatar
Uwe Schulzweida committed
99
100
101
102
  weights[0] = (1.-iw) * (1.-jw);
  weights[1] =     iw  * (1.-jw);
  weights[2] =     iw  *     jw;
  weights[3] = (1.-iw) *     jw;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
  // clang-format on
104
105
}

106
107
template <typename T>
int
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
num_src_points(const Varray<T> &mask, const size_t (&src_add)[4], double (&src_lats)[4])
109
{
110
  int icount = 0;
111

112
  for (int n = 0; n < 4; ++n)
113
    {
114
      if (mask[src_add[n]] != 0)
115
116
117
118
119
120
121
122
        icount++;
      else
        src_lats[n] = 0.;
    }

  return icount;
}

123
// Explicit instantiation
Uwe Schulzweida's avatar
Uwe Schulzweida committed
124
125
template int num_src_points(const Varray<uint8_t> &mask, const size_t (&src_add)[4], double (&src_lats)[4]);
template int num_src_points(const Varray<int> &mask, const size_t (&src_add)[4], double (&src_lats)[4]);
126

127
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
renormalizeWeights(const double (&src_lats)[4], double (&weights)[4])
129
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
130
131
132
  double sum_weights = 0.0;  // sum of weights for normalization
  for (unsigned n = 0; n < 4; ++n) sum_weights += std::fabs(src_lats[n]);
  for (unsigned n = 0; n < 4; ++n) weights[n] = std::fabs(src_lats[n]) / sum_weights;
133
134
}

135
static void
136
bilinearWarning()
137
{
138
  static bool lwarn = true;
139

140
  if (Options::cdoVerbose || lwarn)
141
    {
142
      lwarn = false;
143
      cdoWarning("Bilinear interpolation failed for some grid points - used a distance-weighted average instead!");
144
145
146
    }
}

147
static inline double
Uwe Schulzweida's avatar
Uwe Schulzweida committed
148
bilinearRemap(const double *restrict src_array, const double (&weights)[4], const size_t (&src_add)[4])
149
150
{
  // *tgt_point = 0.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
152
153
  // for (unsigned n = 0; n < 4; ++n) *tgt_point += src_array[src_add[n]]*weights[n];
  return src_array[src_add[0]] * weights[0] + src_array[src_add[1]] * weights[1] + src_array[src_add[2]] * weights[2]
         + src_array[src_add[3]] * weights[3];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154
155
}

156
// This routine computes the weights for a bilinear interpolation.
157
void
158
remapBilinearWeights(RemapSearch &rsearch, RemapVars &rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
{
160
161
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;
162

163
  if (Options::cdoVerbose) cdoPrint("Called %s()", __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
164

Uwe Schulzweida's avatar
Uwe Schulzweida committed
165
166
  if (src_grid->rank != 2) cdoAbort("Can't do bilinear interpolation when source grid rank != 2");

167
  auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
168

169
  progress::init();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170

171
  // Compute mappings from source to target grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
172

173
  auto tgt_grid_size = tgt_grid->size;
174

175
  std::vector<WeightLinks> weightLinks(tgt_grid_size);
176
  weightLinksAlloc(4, tgt_grid_size, weightLinks);
177

178
  auto findex = 0.0;
179

Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
  // Loop over destination grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181

182
183
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) shared(findex, rsearch, weightLinks, tgt_grid_size, src_grid, tgt_grid, rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
184
#endif
185
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
    {
187
188
189
#ifdef _OPENMP
#pragma omp atomic
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
190
      findex++;
191
      if (cdo_omp_get_thread_num() == 0) progress::update(0, 1, findex / tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
192

193
      weightLinks[tgt_cell_add].nlinks = 0;
194

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

197
      const auto llpoint = remapgrid_get_lonlat(tgt_grid, tgt_cell_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198

199
200
201
      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
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
      double weights[4];   //  bilinear weights for four corners
203

204
      // Find nearest square of grid points on source grid
205
      auto search_result = remapSearchSquare(rsearch, llpoint, src_add, src_lats, src_lons);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
206

207
      // Check to see if points are mask points
208
209
210
211
212
      if (search_result > 0)
        {
          for (unsigned n = 0; n < 4; ++n)
            if (!src_grid->mask[src_add[n]]) search_result = 0;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
213

214
      // If point found, find local iw,jw coordinates for weights
215
216
      if (search_result > 0)
        {
217
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218

219
          double iw = 0.0, jw = 0.0;  // current guess for bilinear coordinate
220
          if (remapFindWeights(llpoint, src_lons, src_lats, &iw, &jw))
221
222
            {
              // Successfully found iw,jw - compute weights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
223
224
              bilinearSetWeights(iw, jw, weights);
              store_weightlinks(0, 4, src_add, weights, tgt_cell_add, weightLinks);
225
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
226
          else
227
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
228
              bilinearWarning();
229
230
231
              search_result = -1;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
233

      /*
234
235
        Search for bilinear failed - use a distance-weighted average instead
        (this is typically near the pole) Distance was stored in src_lats!
Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
      */
237
238
      if (search_result < 0)
        {
239
          if (num_src_points(src_grid->mask, src_add, src_lats) > 0)
240
241
            {
              tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
243
              renormalizeWeights(src_lats, weights);
              store_weightlinks(0, 4, src_add, weights, tgt_cell_add, weightLinks);
244
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
245
        }
246
247
    }

248
  progress::update(0, 1, 1);
249

250
  weightLinksToRemapLinks(0, tgt_grid_size, weightLinks, rv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251

252
  if (Options::cdoVerbose) cdoPrint("Square search nearest: %.2f seconds", cdo_get_wtime() - start);
253
}  // remapBilinearWeights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
254

255
// This routine computes and apply the weights for a bilinear interpolation.
256
void
257
remapBilinear(RemapSearch &rsearch, const double *restrict src_array, double *restrict tgt_array, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
{
259
260
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;
261

262
  if (Options::cdoVerbose) cdoPrint("Called %s()", __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
263

Uwe Schulzweida's avatar
Uwe Schulzweida committed
264
265
  if (src_grid->rank != 2) cdoAbort("Can't do bilinear interpolation when source grid rank != 2");

266
  auto start = Options::cdoVerbose ? cdo_get_wtime() : 0.0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
267

268
  progress::init();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
269

270
271
  auto tgt_grid_size = tgt_grid->size;
  auto src_grid_size = src_grid->size;
272

Uwe Schulzweida's avatar
Uwe Schulzweida committed
273
  Varray<uint8_t> src_grid_mask(src_grid_size);
274
275
276
#ifdef _OPENMP
#pragma omp parallel for default(none) schedule(static) shared(src_grid_size, src_array, src_grid_mask, missval)
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277
  for (size_t i = 0; i < src_grid_size; ++i) src_grid_mask[i] = !DBL_IS_EQUAL(src_array[i], missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
278

279
  // Compute mappings from source to target grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
280

281
  auto findex = 0.0;
282

Uwe Schulzweida's avatar
Uwe Schulzweida committed
283
  // Loop over destination grid
Uwe Schulzweida's avatar
Uwe Schulzweida committed
284

285
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
286
#pragma omp parallel for default(none) schedule(static) \
Uwe Schulzweida's avatar
Uwe Schulzweida committed
287
    shared(findex, rsearch, tgt_grid_size, src_grid, tgt_grid, src_array, tgt_array, missval, src_grid_mask)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
288
#endif
289
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
290
    {
291
292
293
#ifdef _OPENMP
#pragma omp atomic
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
294
      findex++;
295
      if (cdo_omp_get_thread_num() == 0) progress::update(0, 1, findex / tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
296

297
      tgt_array[tgt_cell_add] = missval;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
298

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

301
      const auto llpoint = remapgrid_get_lonlat(tgt_grid, tgt_cell_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
302

303
304
305
      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
Uwe Schulzweida's avatar
Uwe Schulzweida committed
306
      double weights[4];   //  bilinear weights for four corners
307

308
      // Find nearest square of grid points on source grid
309
      auto search_result = remapSearchSquare(rsearch, llpoint, src_add, src_lats, src_lons);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
310

311
      // Check to see if points are mask points
312
313
314
      if (search_result > 0)
        {
          for (unsigned n = 0; n < 4; ++n)
315
            if (src_grid_mask[src_add[n]] == 0) search_result = 0;
316
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
317

318
      // If point found, find local iw,jw coordinates for weights
319
320
      if (search_result > 0)
        {
321
          tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
322

323
          double iw = 0.0, jw = 0.0;  // current guess for bilinear coordinate
324
          if (remapFindWeights(llpoint, src_lons, src_lats, &iw, &jw))
325
326
            {
              // Successfully found iw,jw - compute weights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
327
              bilinearSetWeights(iw, jw, weights);
328
              sort_weights_n4(src_add, weights);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
329
              tgt_array[tgt_cell_add] = bilinearRemap(src_array, weights, src_add);
330
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
331
          else
332
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
333
              bilinearWarning();
334
335
336
              search_result = -1;
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
337
338

      /*
339
340
        Search for bilinear failed - use a distance-weighted average instead
        (this is typically near the pole) Distance was stored in src_lats!
Uwe Schulzweida's avatar
Uwe Schulzweida committed
341
      */
342
343
      if (search_result < 0)
        {
344
          if (num_src_points(src_grid_mask, src_add, src_lats) > 0)
345
346
            {
              tgt_grid->cell_frac[tgt_cell_add] = 1.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
              renormalizeWeights(src_lats, weights);
348
              sort_weights_n4(src_add, weights);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
349
              tgt_array[tgt_cell_add] = bilinearRemap(src_array, weights, src_add);
350
            }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
351
        }
352
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353

354
  progress::update(0, 1, 1);
355

356
  if (Options::cdoVerbose) cdoPrint("Square search nearest: %.2f seconds", cdo_get_wtime() - start);
357
}  // remapBilinear