remap_conserv.cc 34.3 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"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
#include "remap_store_link.h"
22
#include "cdo_options.h"
23
#include "progress.h"
24
#include "cimdOmp.h"
25
#include <mpim_grid.h>
26

Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
28
extern "C"
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
30
#include "lib/yac/clipping.h"
#include "lib/yac/area.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
31
}
32

33
struct search_t
34
{
35
  enum yac_edge_type *src_edge_type;
36
  size_t srch_corners;
37
  size_t max_srch_cells;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
39
40
41
  Varray<double> partial_areas;
  Varray<double> partial_weights;
  Varray<grid_cell> src_grid_cells;
  Varray<grid_cell> overlap_buffer;
42
};
43

44
static void
45
cellSearchRealloc(size_t num_srch_cells, search_t &search)
46
{
47
  if (num_srch_cells > search.max_srch_cells)
48
    {
49
50
51
52
      search.partial_areas.resize(num_srch_cells);
      search.partial_weights.resize(num_srch_cells);
      search.overlap_buffer.resize(num_srch_cells);
      search.src_grid_cells.resize(num_srch_cells);
53

54
      for (size_t n = search.max_srch_cells; n < num_srch_cells; ++n)
55
        {
56
57
          search.overlap_buffer[n].array_size = 0;
          search.overlap_buffer[n].num_corners = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
58
59
60
61
          search.overlap_buffer[n].edge_type = nullptr;
          search.overlap_buffer[n].coordinates_x = nullptr;
          search.overlap_buffer[n].coordinates_y = nullptr;
          search.overlap_buffer[n].coordinates_xyz = nullptr;
62
63
64
65
66
67

          search.src_grid_cells[n].array_size = search.srch_corners;
          search.src_grid_cells[n].num_corners = search.srch_corners;
          search.src_grid_cells[n].edge_type = search.src_edge_type;
          search.src_grid_cells[n].coordinates_x = new double[search.srch_corners];
          search.src_grid_cells[n].coordinates_y = new double[search.srch_corners];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
68
          search.src_grid_cells[n].coordinates_xyz = new double[search.srch_corners][3];
69
70
        }

71
      search.max_srch_cells = num_srch_cells;
72
73
74
    }
}

75
static void
76
cellSearchFree(search_t &search)
77
{
78
  const auto max_srch_cells = search.max_srch_cells;
79
  for (size_t n = 0; n < max_srch_cells; n++)
80
    {
81
      if (search.overlap_buffer[n].array_size > 0)
82
        {
83
84
85
86
          free(search.overlap_buffer[n].coordinates_x);
          free(search.overlap_buffer[n].coordinates_y);
          if (search.overlap_buffer[n].coordinates_xyz) free(search.overlap_buffer[n].coordinates_xyz);
          if (search.overlap_buffer[n].edge_type) free(search.overlap_buffer[n].edge_type);
87
        }
88

89
90
91
      delete[] search.src_grid_cells[n].coordinates_x;
      delete[] search.src_grid_cells[n].coordinates_y;
      delete[] search.src_grid_cells[n].coordinates_xyz;
92
    }
93

94
95
96
97
  varrayFree(search.partial_areas);
  varrayFree(search.partial_weights);
  varrayFree(search.overlap_buffer);
  varrayFree(search.src_grid_cells);
98
99
}

100
static void
101
boundbox_from_corners1r(size_t ic, size_t nc, const Varray<double> &corner_lon, const Varray<double> &corner_lat,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
                        float *restrict bound_box)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
103
{
104
  const auto inc = ic * nc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
105

Uwe Schulzweida's avatar
Uwe Schulzweida committed
106
107
  float clat = corner_lat[inc];
  float clon = corner_lon[inc];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
108
109
110
111
112
113

  bound_box[0] = clat;
  bound_box[1] = clat;
  bound_box[2] = clon;
  bound_box[3] = clon;

114
  for (size_t j = 1; j < nc; ++j)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
115
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
116
117
      clat = corner_lat[inc + j];
      clon = corner_lon[inc + j];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
118

119
120
121
122
      if (clat < bound_box[0]) bound_box[0] = clat;
      if (clat > bound_box[1]) bound_box[1] = clat;
      if (clon < bound_box[2]) bound_box[2] = clon;
      if (clon > bound_box[3]) bound_box[3] = clon;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
123
124
    }

Uwe Schulzweida's avatar
Uwe Schulzweida committed
125
  if (fabsf(bound_box[3] - bound_box[2]) > PI_f)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
126
127
    {
      bound_box[2] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
128
      bound_box[3] = PI2_f;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
129
130
    }
  /*
Uwe Schulzweida's avatar
Uwe Schulzweida committed
131
  if ( fabsf(bound_box[3] - bound_box[2]) > PI_f )
Uwe Schulzweida's avatar
Uwe Schulzweida committed
132
    {
133
      if ( bound_box[3] > bound_box[2] && (bound_box[3]-PI2_f) < 0.0f )
134
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
135
136
          float tmp = bound_box[2];
          bound_box[2] = bound_box[3] - PI2_f;
137
138
          bound_box[3] = tmp;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
139
140
141
142
    }
  */
}

143
static inline double
144
gridcell_area(const grid_cell &cell)
145
{
146
  return yac_huiliers_area(cell);
147
148
}

149
static void
150
cdo_compute_overlap_areas(size_t n, search_t &search, const grid_cell &target_cell)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
{
152
  auto overlap_cells = search.overlap_buffer.data();
153

154
  // Do the clipping and get the cell for the overlapping area
155
  yac_cell_clipping(n, search.src_grid_cells.data(), target_cell, overlap_cells);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
156

157
  // Get the partial areas for the overlapping regions
158
  for (size_t i = 0; i < n; i++) search.partial_areas[i] = gridcell_area(overlap_cells[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
160

#ifdef VERBOSE
161
  for (size_t i = 0; i < n; i++) printf("overlap area : %lf\n", search.partial_areas[i]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
162
163
#endif
}
164

165
static constexpr double tol = 1.0e-12;
166

167
168
enum cell_type
{
169
170
171
172
173
174
  UNDEF_CELL,
  LON_LAT_CELL,
  LAT_CELL,
  GREAT_CIRCLE_CELL,
  MIXED_CELL
};
175

176
static void
177
cdo_compute_concave_overlap_areas(size_t N, search_t &search, const grid_cell &target_cell, double target_node_x,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
178
                                  double target_node_y)
179
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
  auto &partial_areas = search.partial_areas;
181
182
  auto overlap_cells = search.overlap_buffer.data();
  auto source_cell = search.src_grid_cells.data();
183

184
185
  double coordinates_x[3] = { -1, -1, -1 };
  double coordinates_y[3] = { -1, -1, -1 };
186
  double coordinates_xyz[3][3] = { { -1, -1, -1 }, { -1, -1, -1 }, { -1, -1, -1 } };
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187
  enum yac_edge_type edge_types[3] = { GREAT_CIRCLE, GREAT_CIRCLE, GREAT_CIRCLE };
188

189
  grid_cell target_partial_cell;
190
  target_partial_cell.array_size = 3;
191
192
  target_partial_cell.coordinates_x = coordinates_x;
  target_partial_cell.coordinates_y = coordinates_y;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
193
  target_partial_cell.coordinates_xyz = coordinates_xyz;
194
195
  target_partial_cell.edge_type = edge_types;
  target_partial_cell.num_corners = 3;
196

197
  // Do the clipping and get the cell for the overlapping area
198

Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
  for (size_t n = 0; n < N; n++) partial_areas[n] = 0.0;
200
201
202
203
204

  // common node point to all partial target cells
  target_partial_cell.coordinates_x[0] = target_node_x;
  target_partial_cell.coordinates_y[0] = target_node_y;

205
206
  auto partial_cell_xyz = target_partial_cell.coordinates_xyz;
  auto cell_xyz = target_cell.coordinates_xyz;
207

208
209
210
211
  gcLLtoXYZ(target_node_x, target_node_y, partial_cell_xyz[0]);


  for (size_t num_corners = 0; num_corners < target_cell.num_corners; ++num_corners)
212
    {
213
214
      const auto corner_a = num_corners;
      const auto corner_b = (num_corners + 1) % target_cell.num_corners;
215
216
217
218
219

      // skip clipping and area calculation for degenerated triangles
      //
      // If this is not sufficient, instead we can try something like:
      //
220
      //     point_list target_list
221
222
      //     init_point_list(&target_list);
      //     generate_point_list(&target_list, target_cell);
223
      //     grid_cell temp_target_cell;
224
225
226
227
228
      //     generate_overlap_cell(target_list, temp_target_cell);
      //     free_point_list(&target_list);
      //
      // and use temp_target_cell for triangulation.
      //
Uwe Schulzweida's avatar
Uwe Schulzweida committed
229
      // Compared to the if statement below the alternative seems to be quite costly.
230

231
232
233
234
235
236
237
238
239
      if (((std::fabs(cell_xyz[corner_a][0] - cell_xyz[corner_b][0]) < tol)
           && (std::fabs(cell_xyz[corner_a][1] - cell_xyz[corner_b][1]) < tol)
           && (std::fabs(cell_xyz[corner_a][2] - cell_xyz[corner_b][2]) < tol))
          || ((std::fabs(cell_xyz[corner_a][0] - partial_cell_xyz[0][0]) < tol)
              && (std::fabs(cell_xyz[corner_a][1] - partial_cell_xyz[0][1]) < tol)
              && (std::fabs(cell_xyz[corner_a][2] - partial_cell_xyz[0][2]) < tol))
          || ((std::fabs(cell_xyz[corner_b][0] - partial_cell_xyz[0][0]) < tol)
              && (std::fabs(cell_xyz[corner_b][1] - partial_cell_xyz[0][1]) < tol)
              && (std::fabs(cell_xyz[corner_b][2] - partial_cell_xyz[0][2]) < tol)))
240
241
        continue;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
242
243
244
245
246
      target_partial_cell.coordinates_x[1] = target_cell.coordinates_x[corner_a];
      target_partial_cell.coordinates_y[1] = target_cell.coordinates_y[corner_a];
      target_partial_cell.coordinates_x[2] = target_cell.coordinates_x[corner_b];
      target_partial_cell.coordinates_y[2] = target_cell.coordinates_y[corner_b];

247
248
249
250
251
252
      partial_cell_xyz[1][0] = cell_xyz[corner_a][0];
      partial_cell_xyz[1][1] = cell_xyz[corner_a][1];
      partial_cell_xyz[1][2] = cell_xyz[corner_a][2];
      partial_cell_xyz[2][0] = cell_xyz[corner_b][0];
      partial_cell_xyz[2][1] = cell_xyz[corner_b][1];
      partial_cell_xyz[2][2] = cell_xyz[corner_b][2];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
253

254
      yac_cell_clipping(N, source_cell, target_partial_cell, overlap_cells);
255

256
      // Get the partial areas for the overlapping regions as sum over the partial target cells.
257
      for (size_t n = 0; n < N; n++)
258
        {
259
          partial_areas[n] += gridcell_area(overlap_cells[n]);
260
261
          // we cannot use pole_area because it is rather inaccurate for great
          // circle edges that are nearly circles of longitude
262
          // partial_areas[n] = pole_area (overlap_cells[n]);
263
        }
264
265
266
    }

#ifdef VERBOSE
Uwe Schulzweida's avatar
Uwe Schulzweida committed
267
  for (size_t n = 0; n < N; n++) printf("overlap area %zu: %lf\n", n, partial_areas[n]);
268
269
270
#endif
}

271
static int
272
get_lonlat_circle_index(RemapGrid *remap_grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
273
274
275
{
  int lonlat_circle_index = -1;

276
  if (remap_grid->num_cell_corners == 4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
277
    {
278
      if (remap_grid->type == RemapGridType::Reg2D)
279
280
281
        {
          lonlat_circle_index = 1;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
282
      else
283
        {
284
285
          const auto &clon = remap_grid->cell_corner_lon;
          const auto &clat = remap_grid->cell_corner_lat;
286
          const auto gridsize = remap_grid->size;
287
          size_t iadd = (gridsize < 100) ? 1 : gridsize / 30 - 1;
288
289
290
291
          size_t num_i = 0, num_eq0 = 0, num_eq1 = 0;

          for (size_t i = 0; i < gridsize; i += iadd)
            {
292
              const auto i4 = i * 4;
293
              num_i++;
294
295
296
              // clang-format off
              if (IS_EQUAL(clon[i4 + 1], clon[i4 + 2]) && IS_EQUAL(clon[i4 + 3], clon[i4 + 0]) &&
                  IS_EQUAL(clat[i4 + 0], clat[i4 + 1]) && IS_EQUAL(clat[i4 + 2], clat[i4 + 3]))
297
298
299
                {
                  num_eq1++;
                }
300
301
              else if (IS_EQUAL(clon[i4 + 0], clon[i4 + 1]) && IS_EQUAL(clon[i4 + 2], clon[i4 + 3]) &&
                       IS_EQUAL(clat[i4 + 1], clat[i4 + 2]) && IS_EQUAL(clat[i4 + 3], clat[i4 + 0]))
302
303
304
                {
                  num_eq0++;
                }
305
              // clang-format on
306
307
308
309
310
            }

          if (num_i == num_eq1) lonlat_circle_index = 1;
          if (num_i == num_eq0) lonlat_circle_index = 0;
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
311
312
    }

313
  // printf("lonlat_circle_index %d\n", lonlat_circle_index);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
314

315
  return lonlat_circle_index;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316
317
}

318
static void
319
remapNormalizeWeights(RemapGrid *tgt_grid, RemapVars &rv)
320
{
321
  // Include centroids in weights and normalize using destination area if requested
322
323
  auto num_links = rv.num_links;
  auto num_wts = rv.num_wts;
324

325
  if (rv.normOpt == NormOpt::DESTAREA)
326
    {
327
      const auto &cell_area = tgt_grid->cell_area;
328
#ifdef _OPENMP
329
#pragma omp parallel for default(none) shared(num_wts, num_links, rv, cell_area)
330
#endif
331
332
      for (size_t n = 0; n < num_links; ++n)
        {
333
          const auto i = rv.tgt_cell_add[n];  // current linear address for target grid cell
Uwe Schulzweida's avatar
Uwe Schulzweida committed
334
          const auto norm_factor = IS_NOT_EQUAL(cell_area[i], 0.0) ? 1.0 / cell_area[i] : 0.0;
335
          rv.wts[n * num_wts] *= norm_factor;
336
        }
337
    }
338
  else if (rv.normOpt == NormOpt::FRACAREA)
339
    {
340
      const auto &cell_frac = tgt_grid->cell_frac;
341
#ifdef _OPENMP
342
#pragma omp parallel for default(none) shared(num_wts, num_links, rv, cell_frac)
343
#endif
344
345
      for (size_t n = 0; n < num_links; ++n)
        {
346
          const auto i = rv.tgt_cell_add[n];  // current linear address for target grid cell
Uwe Schulzweida's avatar
Uwe Schulzweida committed
347
          const auto norm_factor = IS_NOT_EQUAL(cell_frac[i], 0.0) ? 1.0 / cell_frac[i] : 0.0;
348
          rv.wts[n * num_wts] *= norm_factor;
349
        }
350
    }
351
  else if (rv.normOpt == NormOpt::NONE)
352
353
354
355
    {
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
356
357
358
359
360
static void
remapNormalizeWeights(NormOpt normOpt, double cell_area, double cell_frac, size_t num_weights, double *weights)
{
  if (normOpt == NormOpt::DESTAREA)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
361
      const auto norm_factor = IS_NOT_EQUAL(cell_area, 0.0) ? 1.0 / cell_area : 0.0;
362
      for (size_t i = 0; i < num_weights; ++i) weights[i] *= norm_factor;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
363
364
365
    }
  else if (normOpt == NormOpt::FRACAREA)
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
366
      const auto norm_factor = IS_NOT_EQUAL(cell_frac, 0.0) ? 1.0 / cell_frac : 0.0;
367
      for (size_t i = 0; i < num_weights; ++i) weights[i] *= norm_factor;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
368
369
370
371
372
373
    }
  else if (normOpt == NormOpt::NONE)
    {
    }
}

374
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
375
setCellCoordinatesYac(RemapGridType remapGridType, const size_t cellAddr, const size_t numCellCorners, RemapGrid *remap_grid,
376
                      const grid_cell &yac_grid_cell)
377
{
378
379
  auto x = yac_grid_cell.coordinates_x;
  auto y = yac_grid_cell.coordinates_y;
380

381
  if (remapGridType == RemapGridType::Reg2D)
382
    {
383
384
385
      const auto nx = remap_grid->dims[0];
      const auto iy = cellAddr / nx;
      const auto ix = cellAddr - iy * nx;
386
387
388
389
390
391
      const auto reg2d_corner_lon = &remap_grid->reg2d_corner_lon[ix];
      const auto reg2d_corner_lat = &remap_grid->reg2d_corner_lat[iy];
      static const int xi[4] = { 0, 1, 1, 0 };
      static const int yi[4] = { 0, 0, 1, 1 };
      for (int i = 0; i < 4; ++i) x[i] = reg2d_corner_lon[xi[i]];
      for (int i = 0; i < 4; ++i) y[i] = reg2d_corner_lat[yi[i]];
392
393
394
    }
  else
    {
395
396
397
398
      const auto cell_corner_lon = &remap_grid->cell_corner_lon[cellAddr * numCellCorners];
      const auto cell_corner_lat = &remap_grid->cell_corner_lat[cellAddr * numCellCorners];
      for (size_t ic = 0; ic < numCellCorners; ++ic) x[ic] = cell_corner_lon[ic];
      for (size_t ic = 0; ic < numCellCorners; ++ic) y[ic] = cell_corner_lat[ic];
399
    }
400

Uwe Schulzweida's avatar
Uwe Schulzweida committed
401
  double(*restrict xyz)[3] = yac_grid_cell.coordinates_xyz;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
402
  for (size_t ic = 0; ic < numCellCorners; ++ic) gcLLtoXYZ(x[ic], y[ic], xyz[ic]);
403
404
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
405
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
406
setCoordinatesYac(const size_t numCells, RemapGridType remapGridType, const std::vector<size_t> &cellIndices,
Uwe Schulzweida's avatar
Uwe Schulzweida committed
407
                  const size_t numCellCorners, RemapGrid *remap_grid, Varray<grid_cell> &yac_grid_cell)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
408
409
{
  for (size_t n = 0; n < numCells; ++n)
410
    setCellCoordinatesYac(remapGridType, cellIndices[n], numCellCorners, remap_grid, yac_grid_cell[n]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
411
412
}

413
static void
414
reg2d_bound_box(RemapGrid *remap_grid, double *grid_bound_box)
415
{
416
417
  const auto nx = remap_grid->dims[0];
  const auto ny = remap_grid->dims[1];
418
419
  const auto &reg2d_corner_lon = remap_grid->reg2d_corner_lon;
  const auto &reg2d_corner_lat = remap_grid->reg2d_corner_lat;
420
421
422

  grid_bound_box[0] = reg2d_corner_lat[0];
  grid_bound_box[1] = reg2d_corner_lat[ny];
423
  if (grid_bound_box[0] > grid_bound_box[1])
424
425
426
427
428
429
430
431
    {
      grid_bound_box[0] = reg2d_corner_lat[ny];
      grid_bound_box[1] = reg2d_corner_lat[0];
    }
  grid_bound_box[2] = reg2d_corner_lon[0];
  grid_bound_box[3] = reg2d_corner_lon[nx];
}

432
433
434
static void
greatCircleTypeInit(std::vector<enum yac_edge_type> &great_circle_type, size_t src_num_cell_corners, size_t tgt_num_cell_corners)
{
435
  auto max_num_cell_corners = src_num_cell_corners;
436
437
438
439
440
441
  if (tgt_num_cell_corners > max_num_cell_corners) max_num_cell_corners = tgt_num_cell_corners;

  great_circle_type.resize(max_num_cell_corners);
  for (size_t i = 0; i < max_num_cell_corners; ++i) great_circle_type[i] = GREAT_CIRCLE;
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
443
scaleCellFrac(size_t numCells, Varray<double> &cellFrac, const Varray<double> &cellArea)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
444
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
445
446
  for (size_t i = 0; i < numCells; ++i)
    if (IS_NOT_EQUAL(cellArea[i], 0)) cellFrac[i] /= cellArea[i];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
448
}

449
static void
450
yacGridCellInit(grid_cell &cell, size_t numCellCorners, enum yac_edge_type *edgeType)
451
{
452
453
454
455
456
  cell.array_size = numCellCorners;
  cell.num_corners = numCellCorners;
  cell.edge_type = edgeType;
  cell.coordinates_x = new double[numCellCorners];
  cell.coordinates_y = new double[numCellCorners];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
457
  cell.coordinates_xyz = new double[numCellCorners][3];
458
459
}

460
static void
461
yacGridCellFree(const grid_cell &cell)
462
463
464
465
466
467
{
  delete[] cell.coordinates_x;
  delete[] cell.coordinates_y;
  delete[] cell.coordinates_xyz;
}

468
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
469
vec_add_weights(Varray<double> &vec, size_t num_weights, const Varray<double> &weights, const std::vector<size_t> &srch_add)
470
471
472
473
474
475
476
477
478
479
480
481
{
  for (size_t i = 0; i < num_weights; ++i)
    {
      const auto partial_weight = weights[i];
      const auto cell_add = srch_add[i];
#ifdef _OPENMP
#pragma omp atomic
#endif
      vec[cell_add] += partial_weight;
    }
}

482
static size_t
Uwe Schulzweida's avatar
Uwe Schulzweida committed
483
remove_invalid_areas(size_t num_srch_cells, Varray<double> &partial_areas, std::vector<size_t> &srch_add)
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
{
  size_t n = 0;
  for (size_t i = 0; i < num_srch_cells; ++i)
    {
      if (partial_areas[i] > 0)
        {
          partial_areas[n] = partial_areas[i];
          srch_add[n] = srch_add[i];
          n++;
        }
    }

  return n;
}

499
static size_t
Uwe Schulzweida's avatar
Uwe Schulzweida committed
500
remove_invalid_weights(size_t grid_size, size_t num_weights, Varray<double> &partial_weights, std::vector<size_t> &srch_add)
501
502
503
504
{
  size_t n = 0;
  for (size_t i = 0; i < num_weights; ++i)
    {
505
      const auto cell_add = (partial_weights[i] > 0) ? srch_add[i] : grid_size;
506
507
508
509
510
511
512
513
514
515
516
517
      if (cell_add != grid_size)
        {
          partial_weights[n] = partial_weights[i];
          srch_add[n] = cell_add;
          n++;
        }
    }

  return n;
}

static size_t
Uwe Schulzweida's avatar
Uwe Schulzweida committed
518
remove_unmask_weights(const Varray<int> &grid_mask, size_t num_weights, Varray<double> &partial_weights, std::vector<size_t> &srch_add)
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
{
  size_t n = 0;
  for (size_t i = 0; i < num_weights; ++i)
    {
      const auto cell_add = srch_add[i];
      /*
        Store the appropriate addresses and weights.
        Also add contributions to cell areas.
        The source grid mask is the master mask.
      */
      if (grid_mask[cell_add])
        {
          partial_weights[n] = partial_weights[i];
          srch_add[n] = cell_add;
          n++;
        }
    }

  return n;
}

540
void
541
remapConservWeights(RemapSearch &rsearch, RemapVars &rv)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
542
{
543
544
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;
545

546
  bool lcheck = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
547

548
  // Variables necessary if segment manages to hit pole
Uwe Schulzweida's avatar
Uwe Schulzweida committed
549
550
  auto srcGridType = src_grid->type;
  auto tgtGridType = tgt_grid->type;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
551

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

554
  progress::init();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
555

556
  double start = Options::cdoVerbose ? cdo_get_wtime() : 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
557

558
559
  auto src_grid_size = src_grid->size;
  auto tgt_grid_size = tgt_grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
560

561
562
  auto src_num_cell_corners = src_grid->num_cell_corners;
  auto tgt_num_cell_corners = tgt_grid->num_cell_corners;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
563

Uwe Schulzweida's avatar
Uwe Schulzweida committed
564
  enum yac_edge_type lonlat_circle_type[] = { LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE };
565
566
  std::vector<enum yac_edge_type> great_circle_type;
  greatCircleTypeInit(great_circle_type, src_num_cell_corners, tgt_num_cell_corners);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
567

568
569
  auto src_edge_type = great_circle_type.data();
  auto tgt_edge_type = great_circle_type.data();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
570

571
572
  enum cell_type target_cell_type = UNDEF_CELL;

573
  if (src_num_cell_corners == 4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
574
    {
575
      const auto lonlat_circle_index = get_lonlat_circle_index(src_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
576
      if (lonlat_circle_index >= 0) src_edge_type = &lonlat_circle_type[lonlat_circle_index];
Uwe Schulzweida's avatar
Uwe Schulzweida committed
577
578
    }

579
  if (tgt_num_cell_corners == 4)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
580
    {
581
      const auto lonlat_circle_index = get_lonlat_circle_index(tgt_grid);
582
583
584
585
586
      if (lonlat_circle_index >= 0)
        {
          target_cell_type = LON_LAT_CELL;
          tgt_edge_type = &lonlat_circle_type[lonlat_circle_index];
        }
587
588
    }

589
  if (!(tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL))
590
    {
591
      if (tgt_grid->cell_center_lon.empty() || tgt_grid->cell_center_lat.empty())
Uwe Schulzweida's avatar
Uwe Schulzweida committed
592
        cdoAbort("Internal problem (%s): missing target point coordinates!", __func__);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
593
594
    }

595
  std::vector<grid_cell> tgt_grid_cell(Threading::ompNumThreads);
596
  for (int i = 0; i < Threading::ompNumThreads; ++i) yacGridCellInit(tgt_grid_cell[i], tgt_num_cell_corners, tgt_edge_type);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
597

598
  std::vector<search_t> search(Threading::ompNumThreads);
599
  for (int i = 0; i < Threading::ompNumThreads; ++i)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
600
    {
601
602
603
      search[i].srch_corners = src_num_cell_corners;
      search[i].src_edge_type = src_edge_type;
      search[i].max_srch_cells = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
604
605
    }

606
607
  std::vector<std::vector<size_t>> srch_add(Threading::ompNumThreads);
  for (int i = 0; i < Threading::ompNumThreads; ++i) srch_add[i].resize(src_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
608

609
  auto srch_corners = src_num_cell_corners;  // num of corners of srch cells
Uwe Schulzweida's avatar
Uwe Schulzweida committed
610

611
  double src_grid_bound_box[4];
612
  if (srcGridType == RemapGridType::Reg2D) reg2d_bound_box(src_grid, src_grid_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
613

614
  std::vector<WeightLinks> weightLinks(tgt_grid_size);
615

616
  double findex = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
617

618
619
620
621
  size_t num_srch_cells_stat[3];
  num_srch_cells_stat[0] = 0;
  num_srch_cells_stat[1] = 100000;
  num_srch_cells_stat[2] = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
622

Uwe Schulzweida's avatar
Uwe Schulzweida committed
623
  extern CellSearchMethod cellSearchMethod;
624
  bool useCellsearch = cellSearchMethod == CellSearchMethod::spherepart || srcGridType == RemapGridType::Reg2D;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
625

Uwe Schulzweida's avatar
Uwe Schulzweida committed
626
  // Loop over destination grid
627

628
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
629
630
631
#pragma omp parallel for schedule(dynamic) default(none)                                                                    \
    shared(findex, rsearch, srcGridType, tgtGridType, src_grid_bound_box, rv, Options::cdoVerbose, tgt_num_cell_corners,    \
           target_cell_type, weightLinks, srch_corners, src_grid, tgt_grid, tgt_grid_size, src_grid_size, search, srch_add, \
632
           tgt_grid_cell, num_srch_cells_stat, useCellsearch)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
633
#endif
634
  for (size_t tgt_cell_add = 0; tgt_cell_add < tgt_grid_size; ++tgt_cell_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
635
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
636
      const auto ompthID = cdo_omp_get_thread_num();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
637

638
639
640
#ifdef _OPENMP
#pragma omp atomic
#endif
Uwe Schulzweida's avatar
Uwe Schulzweida committed
641
      findex++;
642
      if (ompthID == 0) progress::update(0, 1, findex / tgt_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
643

644
      weightLinks[tgt_cell_add].nlinks = 0;
645

646
647
      if (!tgt_grid->mask[tgt_cell_add]) continue;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
648
      setCellCoordinatesYac(tgtGridType, tgt_cell_add, tgt_num_cell_corners, tgt_grid, tgt_grid_cell[ompthID]);
649

650
      // Get search cells
Uwe Schulzweida's avatar
Uwe Schulzweida committed
651
      size_t num_srch_cells;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
652
      if (useCellsearch)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
653
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
654
655
          // num_srch_cells = remapSearchCells(rsearch, target_cell_type == LON_LAT_CELL, tgt_grid_cell[ompthID],
          // srch_add[ompthID].data());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
656
657
          num_srch_cells
              = remapSearchCells(rsearch, tgtGridType == RemapGridType::Reg2D, tgt_grid_cell[ompthID], srch_add[ompthID].data());
658
659
660
661
662
663
664
        }
      else
        {
          float tgt_cell_bound_box_r[4];
          boundbox_from_corners1r(tgt_cell_add, tgt_num_cell_corners, tgt_grid->cell_corner_lon, tgt_grid->cell_corner_lat,
                                  tgt_cell_bound_box_r);

Uwe Schulzweida's avatar
Uwe Schulzweida committed
665
666
          num_srch_cells
              = get_srch_cells(tgt_cell_add, rsearch.tgtBins, rsearch.srcBins, tgt_cell_bound_box_r, srch_add[ompthID].data());
667
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
668

669
      if (1 && Options::cdoVerbose)
670
671
        {
          num_srch_cells_stat[0] += num_srch_cells;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
672
673
          if (num_srch_cells < num_srch_cells_stat[1]) num_srch_cells_stat[1] = num_srch_cells;
          if (num_srch_cells > num_srch_cells_stat[2]) num_srch_cells_stat[2] = num_srch_cells;
674
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
675

676
      if (0 && Options::cdoVerbose) printf("tgt_cell_add %zu  num_srch_cells %zu\n", tgt_cell_add, num_srch_cells);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
677

678
      if (num_srch_cells == 0) continue;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
679

680
      // Create search arrays
Uwe Schulzweida's avatar
Uwe Schulzweida committed
681

682
      cellSearchRealloc(num_srch_cells, search[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
683

Uwe Schulzweida's avatar
Uwe Schulzweida committed
684
      setCoordinatesYac(num_srch_cells, srcGridType, srch_add[ompthID], srch_corners, src_grid, search[ompthID].src_grid_cells);
685
686
687

      if (tgt_num_cell_corners < 4 || target_cell_type == LON_LAT_CELL)
        {
688
          cdo_compute_overlap_areas(num_srch_cells, search[ompthID], tgt_grid_cell[ompthID]);
689
        }
690
      else
691
        {
692
          cdo_compute_concave_overlap_areas(num_srch_cells, search[ompthID], tgt_grid_cell[ompthID],
Uwe Schulzweida's avatar
Uwe Schulzweida committed
693
                                            tgt_grid->cell_center_lon[tgt_cell_add], tgt_grid->cell_center_lat[tgt_cell_add]);
694
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
695

696
697
698
699
      auto &partial_areas = search[ompthID].partial_areas;
      auto &partial_weights = search[ompthID].partial_weights;

      auto num_weights = remove_invalid_areas(num_srch_cells, partial_areas, srch_add[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
700

701
      const auto tgt_area = gridcell_area(tgt_grid_cell[ompthID]);
702
      tgt_grid->cell_area[tgt_cell_add] = tgt_area;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
703

704
      for (size_t i = 0; i < num_weights; ++i) partial_weights[i] = partial_areas[i] / tgt_area;
705

706
      if (rv.normOpt == NormOpt::FRACAREA) yac_correct_weights(num_weights, partial_weights.data());
Uwe Schulzweida's avatar
Uwe Schulzweida committed
707

708
      for (size_t i = 0; i < num_weights; ++i) partial_weights[i] *= tgt_area;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
709

710
      num_weights = remove_invalid_weights(src_grid_size, num_weights, partial_weights, srch_add[ompthID]);
711

712
      vec_add_weights(src_grid->cell_area, num_weights, partial_weights, srch_add[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
713

714
      num_weights = remove_unmask_weights(src_grid->mask, num_weights, partial_weights, srch_add[ompthID]);
715

716
      vec_add_weights(src_grid->cell_frac, num_weights, partial_weights, srch_add[ompthID]);
717

718
      tgt_grid->cell_frac[tgt_cell_add] = varraySum(num_weights, partial_weights);
719

720
      store_weightlinks(1, num_weights, srch_add[ompthID].data(), partial_weights.data(), tgt_cell_add, weightLinks);
721
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
722

723
  progress::update(0, 1, 1);
724

725
  if (1 && Options::cdoVerbose)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
726
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
728
      cdoPrint("Num search cells min,mean,max :  %zu  %3.1f  %zu", num_srch_cells_stat[1],
               num_srch_cells_stat[0] / (double) tgt_grid_size, num_srch_cells_stat[2]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
729
730
    }

731
  // Finished with all cells: deallocate search arrays
Uwe Schulzweida's avatar
Uwe Schulzweida committed
732
  for (auto ompthID = 0; ompthID < Threading::ompNumThreads; ++ompthID)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
733
    {
734
      cellSearchFree(search[ompthID]);
735
      yacGridCellFree(tgt_grid_cell[ompthID]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
736
    }
737

738
  weightLinksToRemapLinks(1, tgt_grid_size, weightLinks, rv);
739

740
741
  // Normalize weights using destination area if requested
  remapNormalizeWeights(tgt_grid, rv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
742

743
  if (Options::cdoVerbose) cdoPrint("Total number of links = %zu", rv.num_links);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
744

Uwe Schulzweida's avatar
Uwe Schulzweida committed
745
746
  scaleCellFrac(src_grid_size, src_grid->cell_frac, src_grid->cell_area);
  scaleCellFrac(tgt_grid_size, tgt_grid->cell_frac, tgt_grid->cell_area);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
747

748
  // Perform some error checking on final weights
749
  if (lcheck)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
750
    {
751
752
      remapCheckArea(src_grid_size, &src_grid->cell_area[0], "Source");
      remapCheckArea(tgt_grid_size, &tgt_grid->cell_area[0], "Target");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
753

754
      remapVarsCheckWeights(rv);
755
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
756

757
  if (Options::cdoVerbose) cdoPrint("Cells search: %.2f seconds", cdo_get_wtime() - start);
758
}  // remapConservWeights
Uwe Schulzweida's avatar
Uwe Schulzweida committed
759

760
static double
Uwe Schulzweida's avatar
Uwe Schulzweida committed
761
conservRemap(const double *restrict src_array, size_t num_weights, const Varray<double> &weights, const std::vector<size_t> &src_add)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
762
{
763
  double tgt_point = 0.;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
764
  for (size_t i = 0; i < num_weights; ++i) tgt_point += src_array[src_add[i]] * weights[i];
765
766

  return tgt_point;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
767
768
}

769
void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
770
remapConserv(NormOpt normOpt, RemapSearch &rsearch, const double *restrict src_array, double *restrict tgt_array, double missval)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
771
{
772
773
  auto src_grid = rsearch.srcGrid;
  auto tgt_grid = rsearch.tgtGrid;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
774
775
776
777
778
779
780
781
782

  bool lcheck = true;

  // Variables necessary if segment manages to hit pole
  auto srcGridType = src_grid->type;
  auto tgtGridType = tgt_grid->type;

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

783
  progress::init();
Uwe Schulzweida's avatar
Uwe Schulzweida committed
784
785
786

  double start = Options::cdoVerbose ? cdo_get_wtime() : 0;

787
788
  auto src_grid_size = src_grid->size;
  auto tgt_grid_size = tgt_grid->size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
789

Uwe Schulzweida's avatar
Uwe Schulzweida committed
790
  Varray<int> src_grid_mask(src_grid_size);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
791
792
793
#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
794
  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
795

796
797
  auto src_num_cell_corners = src_grid->num_cell_corners;
  auto tgt_num_cell_corners = tgt_grid->num_cell_corners;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
798
799
800
801
802

  enum yac_edge_type lonlat_circle_type[] = { LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE, LAT_CIRCLE, LON_CIRCLE };
  std::vector<enum yac_edge_type> great_circle_type;
  greatCircleTypeInit(great_circle_type, src_num_cell_corners, tgt_num_cell_corners);