remaplib.cc 33.2 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.
*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
17
18
19
/*
  This is a C library of the Fortran SCRIP version 1.4

Uwe Schulzweida's avatar
Uwe Schulzweida committed
20
  ===>>> Please send bug reports to <https://mpimet.mpg.de/cdo> <<<===
Uwe Schulzweida's avatar
Uwe Schulzweida committed
21
22
23
24
25
26
27

  Spherical Coordinate Remapping and Interpolation Package (SCRIP)
  ================================================================

  SCRIP is a software package which computes addresses and weights for
  remapping and interpolating fields between grids in spherical coordinates.
  It was written originally for remapping fields to other grids in a coupled
28
  climate model, but is sufficiently general that it can be used in other
Uwe Schulzweida's avatar
Uwe Schulzweida committed
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
  applications as well. The package should work for any grid on the surface
  of a sphere. SCRIP currently supports four remapping options:

  Conservative remapping
  ----------------------
  First- and second-order conservative remapping as described in
  Jones (1999, Monthly Weather Review, 127, 2204-2210).

  Bilinear interpolation
  ----------------------
  Slightly generalized to use a local bilinear approximation
  (only logically-rectangular grids).

  Bicubic interpolation
  ----------------------
  Similarly generalized (only logically-rectangular grids).

  Distance-weighted averaging
  ---------------------------
48
49
  Distance-weighted average of a user-specified number of nearest neighbor
  values.
Uwe Schulzweida's avatar
Uwe Schulzweida committed
50
51
52
53
54
55
56

  Documentation
  =============

  http://climate.lanl.gov/Software/SCRIP/SCRIPusers.pdf

*/
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57
/*
58
  2013-11-08 Uwe Schulzweida: split remapgrid class to src_grid and tgt_grid
59
60
61
62
  2012-01-16 Uwe Schulzweida: alloc grid2_bound_box only for conservative
  remapping 2011-01-07 Uwe Schulzweida: Changed remap weights from 2D to 1D
  array 2009-05-25 Uwe Schulzweida: Changed restrict data type from double to
  int 2009-01-11 Uwe Schulzweida: OpenMP parallelization
Uwe Schulzweida's avatar
Uwe Schulzweida committed
63
 */
Cedrick Ansorge's avatar
Cedrick Ansorge committed
64

65
66
#ifdef HAVE_CONFIG_H
#include "config.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
67
68
#endif

Ralf Mueller's avatar
Ralf Mueller committed
69
#include <cdi.h>
70

71
#include "cdo_options.h"
72
#include "cimdOmp.h"
73
#include "dmemory.h"
74
#include "process_int.h"
75
#include "cdo_wtime.h"
76
#include <mpim_grid.h>
77
#include "gridreference.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
78
79
#include "remap.h"

Uwe Schulzweida's avatar
Uwe Schulzweida committed
80
#define IS_REG2D_GRID(gridID) (gridInqType(gridID) == GRID_LONLAT || gridInqType(gridID) == GRID_GAUSSIAN)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
81

82
83
84
85
static bool remap_gen_weights = true;
static bool remap_write_remap = false;
static int remap_num_srch_bins = 180;
#define DEFAULT_MAX_ITER 100
Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
size_t remap_max_iter = DEFAULT_MAX_ITER;  // Max iteration count for i, j iteration
87

88
89
void
remap_set_int(int remapvar, int value)
90
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
91
92
93
94
95
96
97
  // clang-format off
  if      (remapvar == REMAP_WRITE_REMAP)   remap_write_remap = value > 0;
  else if (remapvar == REMAP_MAX_ITER)      remap_max_iter = value;
  else if (remapvar == REMAP_NUM_SRCH_BINS) remap_num_srch_bins = value;
  else if (remapvar == REMAP_GENWEIGHTS)    remap_gen_weights = value > 0;
  else    cdoAbort("Unsupported remap variable (%d)!", remapvar);
  // clang-format on
98
99
}

100
/*****************************************************************************/
101
static void
102
103
boundboxFromCorners(size_t size, size_t nc, const double *restrict corner_lon, const double *restrict corner_lat,
                    float *restrict bound_box)
104
{
105
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
106
#pragma omp parallel for default(none) shared(bound_box, corner_lat, corner_lon, nc, size)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
107
#endif
108
  for (size_t i = 0; i < size; ++i)
109
    {
110
111
      size_t i4 = i << 2;  // *4
      size_t inc = i * nc;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
113
      float clat = corner_lat[inc];
      float clon = corner_lon[inc];
114
      bound_box[i4 + 0] = clat;
115
116
117
118
119
      bound_box[i4 + 1] = clat;
      bound_box[i4 + 2] = clon;
      bound_box[i4 + 3] = clon;
      for (size_t j = 1; j < nc; ++j)
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
120
121
          clat = corner_lat[inc + j];
          clon = corner_lon[inc + j];
122
          if (clat < bound_box[i4 + 0]) bound_box[i4 + 0] = clat;
123
124
125
126
          if (clat > bound_box[i4 + 1]) bound_box[i4 + 1] = clat;
          if (clon < bound_box[i4 + 2]) bound_box[i4 + 2] = clon;
          if (clon > bound_box[i4 + 3]) bound_box[i4 + 3] = clon;
        }
127
128
129
    }
}

130
static void
131
132
boundboxFromCenter(bool lonIsCyclic, size_t size, size_t nx, size_t ny, const double *restrict center_lon,
                   const double *restrict center_lat, float *restrict bound_box)
133
{
134
#ifdef _OPENMP
135
#pragma omp parallel for default(none)  shared(lonIsCyclic, size, nx, ny, center_lon, center_lat, bound_box)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
136
#endif
137
  for (size_t n = 0; n < size; n++)
138
    {
139
140
141
      size_t idx[4];
      float tmp_lats[4], tmp_lons[4];

142
      size_t n4 = n << 2;
143

Uwe Schulzweida's avatar
Uwe Schulzweida committed
144
      // Find N,S and NE points to this grid point
145

146
147
      size_t j = n / nx;
      size_t i = n - j * nx;
148

149
150
      size_t ip1 = (i < (nx - 1)) ? i + 1 : lonIsCyclic ? 0 : i;
      size_t jp1 = (j < (ny - 1)) ? j + 1 : j;
151

152
153
154
155
      idx[0] = n;
      idx[1] = j * nx + ip1;    // east
      idx[2] = jp1 * nx + ip1;  // north-east
      idx[3] = jp1 * nx + i;    // north
156

Uwe Schulzweida's avatar
Uwe Schulzweida committed
157
      // Find N,S and NE lat/lon coords and check bounding box
158

Uwe Schulzweida's avatar
Uwe Schulzweida committed
159
160
      for (unsigned k = 0; k < 4; ++k) tmp_lons[k] = center_lon[idx[k]];
      for (unsigned k = 0; k < 4; ++k) tmp_lats[k] = center_lat[idx[k]];
161

162
      bound_box[n4 + 0] = tmp_lats[0];
163
164
165
166
      bound_box[n4 + 1] = tmp_lats[0];
      bound_box[n4 + 2] = tmp_lons[0];
      bound_box[n4 + 3] = tmp_lons[0];

167
      for (unsigned k = 1; k < 4; k++)
168
        {
169
          if (tmp_lats[k] < bound_box[n4 + 0]) bound_box[n4 + 0] = tmp_lats[k];
170
171
172
173
          if (tmp_lats[k] > bound_box[n4 + 1]) bound_box[n4 + 1] = tmp_lats[k];
          if (tmp_lons[k] < bound_box[n4 + 2]) bound_box[n4 + 2] = tmp_lons[k];
          if (tmp_lons[k] > bound_box[n4 + 3]) bound_box[n4 + 3] = tmp_lons[k];
        }
174
175
176
    }
}

177
void
178
remapgrid_get_lonlat(RemapGrid *grid, size_t cell_add, double *plon, double *plat)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
179
{
180
  if (grid->type == RemapGridType::Reg2D)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
181
    {
182
      size_t nx = grid->dims[0];
183
184
      size_t iy = cell_add / nx;
      size_t ix = cell_add - iy * nx;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185
186
      *plat = grid->reg2d_center_lat[iy];
      *plon = grid->reg2d_center_lon[ix];
187
      if (*plon < 0) *plon += PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
188
189
190
191
192
193
194
195
    }
  else
    {
      *plat = grid->cell_center_lat[cell_add];
      *plon = grid->cell_center_lon[cell_add];
    }
}

196
197
void
check_lon_range(size_t nlons, double *lons)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
198
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
  assert(lons != nullptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
200

201
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
202
#pragma omp parallel for default(none) shared(nlons, lons)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
203
#endif
204
  for (size_t n = 0; n < nlons; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
205
    {
206
      // remove missing values
207
208
      if (lons[n] < -PI2) lons[n] = 0;
      if (lons[n] > 2 * PI2) lons[n] = PI2;
209

210
      if (lons[n] > PI2) lons[n] -= PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
211
      if (lons[n] < 0.0) lons[n] += PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
212
213
214
    }
}

215
216
void
check_lat_range(size_t nlats, double *lats)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
217
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
218
  assert(lats != nullptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
219

220
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
221
#pragma omp parallel for default(none) shared(nlats, lats)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
222
#endif
223
  for (size_t n = 0; n < nlats; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
224
    {
225
226
      if (lats[n] > PIH) lats[n] = PIH;
      if (lats[n] < -PIH) lats[n] = -PIH;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
227
228
229
    }
}

230
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
231
check_lon_boundbox_range(size_t nlons, float *bound_box)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
232
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
233
  assert(bound_box != nullptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
234

235
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
236
#pragma omp parallel for default(none) shared(nlons, bound_box)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
237
#endif
238
  for (size_t n = 0; n < nlons; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
239
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
240
      size_t n4 = n << 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
241
      if (fabsf(bound_box[n4 + 3] - bound_box[n4 + 2]) > PI_f)
242
        {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
243
244
          bound_box[n4 + 2] = 0.0f;
          bound_box[n4 + 3] = PI2_f;
245
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
246
247
248
    }
}

249
static void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
250
check_lat_boundbox_range(size_t nlats, float *restrict bound_box, double *restrict lats)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
251
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
252
  assert(bound_box != nullptr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
253

254
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
255
#pragma omp parallel for default(none) shared(nlats, bound_box, lats)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
256
#endif
257
  for (size_t n = 0; n < nlats; ++n)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
258
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
259
260
261
      size_t n4 = n << 2;
      if ((float) lats[n] < bound_box[n4]) bound_box[n4] = -PIH_f;
      if ((float) lats[n] > bound_box[n4 + 1]) bound_box[n4 + 1] = PIH_f;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
262
263
264
    }
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
265
/*****************************************************************************/
266

267
268
static void
grid_check_lat_borders_rad(size_t n, double *ybounds)
269
{
270
271
  constexpr double YLIM = 88 * DEG2RAD;
  if (ybounds[0] > ybounds[n - 1])
272
    {
273
274
      if (ybounds[0] > YLIM) ybounds[0] = PIH;
      if (ybounds[n - 1] < -YLIM) ybounds[n - 1] = -PIH;
275
276
277
    }
  else
    {
278
279
      if (ybounds[0] < -YLIM) ybounds[0] = -PIH;
      if (ybounds[n - 1] > YLIM) ybounds[n - 1] = PIH;
280
281
    }
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
282

283
static void
284
remap_define_reg2d(int gridID, RemapGrid &grid)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
285
{
286
287
  size_t nx = grid.dims[0];
  size_t ny = grid.dims[1];
288
289
  size_t nxp1 = nx + 1;
  size_t nyp1 = ny + 1;
290

291
  size_t nxm = nx;
292
  if (grid.is_cyclic) nxm++;
293

294
  if (grid.size != nx * ny) cdoAbort("Internal error, wrong dimensions!");
295

296
297
  grid.reg2d_center_lon = (double *) Malloc(nxm * sizeof(double));
  grid.reg2d_center_lat = (double *) Malloc(ny * sizeof(double));
298

299
300
301
302
  grid.reg2d_center_lon[0] = 0;
  grid.reg2d_center_lat[0] = 0;
  gridInqXvals(gridID, grid.reg2d_center_lon);
  gridInqYvals(gridID, grid.reg2d_center_lat);
303

304
  // Convert lat/lon units if required
305
306
  cdo_grid_to_radian(gridID, CDI_XAXIS, nx, grid.reg2d_center_lon, "grid reg2d center lon");
  cdo_grid_to_radian(gridID, CDI_YAXIS, ny, grid.reg2d_center_lat, "grid reg2d center lat");
307

308
  if (grid.reg2d_center_lon[nx - 1] < grid.reg2d_center_lon[0])
309
    for (size_t i = 1; i < nx; ++i)
310
      if (grid.reg2d_center_lon[i] < grid.reg2d_center_lon[i - 1]) grid.reg2d_center_lon[i] += PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
311

312
  if (grid.is_cyclic) grid.reg2d_center_lon[nx] = grid.reg2d_center_lon[0] + PI2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
313

314
315
  grid.reg2d_corner_lon = (double *) Malloc(nxp1 * sizeof(double));
  grid.reg2d_corner_lat = (double *) Malloc(nyp1 * sizeof(double));
Uwe Schulzweida's avatar
Uwe Schulzweida committed
316

317
318
319
320
321
322
  if (gridInqXbounds(gridID, nullptr))
    {
      std::vector<double> xbounds(2 * nx);
      gridInqXbounds(gridID, xbounds.data());
      grid.reg2d_corner_lon[0] = xbounds[1];
      for (size_t i = 0; i < nx; ++i) grid.reg2d_corner_lon[i + 1] = xbounds[2 * i];
323
      cdo_grid_to_radian(gridID, CDI_XAXIS, nx + 1, grid.reg2d_corner_lon, "grid reg2d corner lon");
324
325
326
327
328
329
330
331
332
333
334
335
    }
  else
    {
      grid_gen_corners(nx, grid.reg2d_center_lon, grid.reg2d_corner_lon);
    }

  if (gridInqYbounds(gridID, nullptr))
    {
      std::vector<double> ybounds(2 * ny);
      gridInqYbounds(gridID, ybounds.data());
      grid.reg2d_corner_lat[0] = ybounds[1];
      for (size_t j = 0; j < ny; ++j) grid.reg2d_corner_lat[j + 1] = ybounds[2 * j];
336
      cdo_grid_to_radian(gridID, CDI_YAXIS, ny + 1, grid.reg2d_corner_lat, "grid reg2d corner lat");
337
338
339
340
341
342
    }
  else
    {
      grid_gen_corners(ny, grid.reg2d_center_lat, grid.reg2d_corner_lat);
      grid_check_lat_borders_rad(ny + 1, grid.reg2d_corner_lat);
    }
343
}
344

345
static void
346
remapDefineGrid(RemapMethod mapType, int gridID, RemapGrid &grid, const char *txt)
347
{
348
  bool lgrid_destroy = false;
349
  int gridID_gme = -1;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
350

351
352
  auto gridtype = gridInqType(grid.gridID);
  if (gridtype != GRID_UNSTRUCTURED && gridtype != GRID_CURVILINEAR)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
353
    {
354
      if (gridtype == GRID_GME)
355
        {
356
357
          gridID_gme = gridToUnstructured(grid.gridID, 1);
          grid.nvgp = gridInqSize(gridID_gme);
358
359
          gridID = gridDuplicate(gridID_gme);
          gridCompress(gridID);
360
          grid.luse_cell_corners = true;
361
        }
362
363
364
365
366
      else if (gridtype == GRID_GAUSSIAN_REDUCED)
        {
          lgrid_destroy = true;
          gridID = gridToUnstructured(grid.gridID, 1);
        }
367
      else if (remap_write_remap || grid.type != RemapGridType::Reg2D)
368
369
        {
          lgrid_destroy = true;
370
          gridID = gridToCurvilinear(grid.gridID, 1);
371
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
372
373
    }

374
  auto gridsize = grid.size = gridInqSize(gridID);
375

376
377
378
  grid.dims[0] = gridInqXsize(gridID);
  grid.dims[1] = gridInqYsize(gridID);
  if (gridInqType(grid.gridID) != GRID_UNSTRUCTURED)
379
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
380
      if (grid.dims[0] == 0) cdoAbort("%s grid without longitude coordinates!", gridNamePtr(gridInqType(grid.gridID)));
381
      if (grid.dims[1] == 0) cdoAbort("%s grid without latitude coordinates!", gridNamePtr(gridInqType(grid.gridID)));
382
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
383

384
  grid.is_cyclic = (gridIsCircular(gridID) > 0);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
385

386
  grid.rank = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? 1 : 2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
387

388
  grid.num_cell_corners = (gridInqType(gridID) == GRID_UNSTRUCTURED) ? gridInqNvertex(gridID) : 4;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
389

390
  remapGridAlloc(mapType, grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
391

392
  // Initialize logical mask
Uwe Schulzweida's avatar
Uwe Schulzweida committed
393

394
#ifdef _OPENMP
Uwe Schulzweida's avatar
Uwe Schulzweida committed
395
396
#pragma omp parallel for default(none) shared(gridsize, grid)
#endif
397
  for (size_t i = 0; i < gridsize; ++i) grid.mask[i] = true;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
398

Uwe Schulzweida's avatar
Uwe Schulzweida committed
399
  if (gridInqMask(gridID, nullptr))
400
    {
401
402
      std::vector<int> mask(gridsize);
      gridInqMask(gridID, &mask[0]);
403
      for (size_t i = 0; i < gridsize; ++i)
404
        if (mask[i] == 0) grid.mask[i] = false;
405
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
406

407
  if (!remap_write_remap && grid.type == RemapGridType::Reg2D) return;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
408

Uwe Schulzweida's avatar
Uwe Schulzweida committed
409
  if (!(gridInqXvals(gridID, nullptr) && gridInqYvals(gridID, nullptr))) cdoAbort("%s grid cell center coordinates missing!", txt);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
410

411
412
  gridInqXvals(gridID, grid.cell_center_lon);
  gridInqYvals(gridID, grid.cell_center_lat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
413

414
  if (grid.lneed_cell_corners)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
415
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
416
      if (gridInqXbounds(gridID, nullptr) && gridInqYbounds(gridID, nullptr))
417
        {
418
419
          gridInqXbounds(gridID, grid.cell_corner_lon);
          gridInqYbounds(gridID, grid.cell_corner_lat);
420
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
421
      else
422
423
424
        {
          cdoAbort("%s grid cell corner coordinates missing!", txt);
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
425
426
    }

427
  if (gridInqType(grid.gridID) == GRID_GME) gridInqMaskGME(gridID_gme, &grid.vgpm[0]);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
428

429
  // Convert lat/lon units if required
430
431
  cdo_grid_to_radian(gridID, CDI_XAXIS, grid.size, grid.cell_center_lon, "grid center lon");
  cdo_grid_to_radian(gridID, CDI_YAXIS, grid.size, grid.cell_center_lat, "grid center lat");
432
  if (grid.num_cell_corners && grid.lneed_cell_corners)
433
    {
434
435
      cdo_grid_to_radian(gridID, CDI_XAXIS, grid.num_cell_corners * grid.size, grid.cell_corner_lon, "grid corner lon");
      cdo_grid_to_radian(gridID, CDI_YAXIS, grid.num_cell_corners * grid.size, grid.cell_corner_lat, "grid corner lat");
436
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
437

438
  if (lgrid_destroy) gridDestroy(gridID);
439

440
  // Convert longitudes to 0,2pi interval
441
  check_lon_range(grid.size, grid.cell_center_lon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
442
  if (grid.num_cell_corners && grid.lneed_cell_corners) check_lon_range(grid.num_cell_corners * grid.size, grid.cell_corner_lon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
443

444
  //  Make sure input latitude range is within the machine values for +/- pi/2
445
  check_lat_range(grid.size, grid.cell_center_lat);
446

Uwe Schulzweida's avatar
Uwe Schulzweida committed
447
  if (grid.num_cell_corners && grid.lneed_cell_corners) check_lat_range(grid.num_cell_corners * grid.size, grid.cell_corner_lat);
448
449
450
}

/*  Compute bounding boxes for restricting future grid searches */
451
static void
452
cell_bounding_boxes(RemapGrid &grid, float *cell_bound_box, int remap_grid_basis)
453
{
454
  if (grid.luse_cell_corners)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
455
    {
456
      if (grid.lneed_cell_corners)
457
        {
458
          if (Options::cdoVerbose) cdoPrint("Grid: boundboxFromCorners");
459
          boundboxFromCorners(grid.size, grid.num_cell_corners, grid.cell_corner_lon, grid.cell_corner_lat, cell_bound_box);
460
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
461
      else  // full grid search
462
        {
463
          if (Options::cdoVerbose) cdoPrint("Grid: bounds missing -> full grid search!");
464

465
          size_t gridsize = grid.size;
466
467
          for (size_t i = 0; i < gridsize; ++i)
            {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
468
469
470
471
              cell_bound_box[i * 4] = -PIH_f;
              cell_bound_box[i * 4 + 1] = PIH_f;
              cell_bound_box[i * 4 + 2] = 0.0f;
              cell_bound_box[i * 4 + 3] = PI2_f;
472
473
            }
        }
474
    }
475
  else if (remap_grid_basis == REMAP_GRID_BASIS_SRC)
476
    {
477
      if (Options::cdoVerbose) cdoPrint("Grid: boundboxFromCenter");
478
      if (grid.rank != 2) cdoAbort("Internal problem, grid rank = %d!", grid.rank);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
479

480
481
      size_t nx = grid.dims[0];
      size_t ny = grid.dims[1];
482
      boundboxFromCenter(grid.is_cyclic, grid.size, nx, ny, grid.cell_center_lon, grid.cell_center_lat, cell_bound_box);
483
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
484

Uwe Schulzweida's avatar
Uwe Schulzweida committed
485
  if (remap_grid_basis == REMAP_GRID_BASIS_SRC || grid.lneed_cell_corners) check_lon_boundbox_range(grid.size, cell_bound_box);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
486

487
  // Try to check for cells that overlap poles
488
  if (remap_grid_basis == REMAP_GRID_BASIS_SRC || grid.lneed_cell_corners)
489
490
491
    check_lat_boundbox_range(grid.size, cell_bound_box, grid.cell_center_lat);
}

492
493
494
495
496
497
498
void
remapGridAlloc(RemapMethod mapType, RemapGrid &grid)
{
  if (grid.nvgp) grid.vgpm.resize(grid.nvgp);

  grid.mask.resize(grid.size);

499
  if (remap_write_remap || grid.type != RemapGridType::Reg2D)
500
501
502
503
504
    {
      grid.cell_center_lon = (double *) Malloc(grid.size * sizeof(double));
      grid.cell_center_lat = (double *) Malloc(grid.size * sizeof(double));
    }

505
  if (mapType == RemapMethod::CONSERV_SCRIP || mapType == RemapMethod::CONSERV)
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
    {
      grid.cell_area.resize(grid.size, 0.0);
    }

  grid.cell_frac.resize(grid.size, 0.0);

  if (grid.lneed_cell_corners)
    {
      if (grid.num_cell_corners > 0)
        {
          size_t nalloc = grid.num_cell_corners * grid.size;
          grid.cell_corner_lon = (double *) Calloc(nalloc, sizeof(double));
          grid.cell_corner_lat = (double *) Calloc(nalloc, sizeof(double));
        }
    }
}

523
void
524
remapGridInit(RemapGrid &grid)
525
{
526
  grid.tmpgridID = -1;
527
  grid.type = RemapGridType::Undefined;
528
529
530
531
532
533
534

  grid.num_cell_corners = 0;
  grid.luse_cell_corners = false;
  grid.lneed_cell_corners = false;

  grid.nvgp = 0;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
535
536
537
538
  grid.reg2d_center_lon = nullptr;
  grid.reg2d_center_lat = nullptr;
  grid.reg2d_corner_lon = nullptr;
  grid.reg2d_corner_lat = nullptr;
539

Uwe Schulzweida's avatar
Uwe Schulzweida committed
540
541
542
543
  grid.cell_center_lon = nullptr;
  grid.cell_center_lat = nullptr;
  grid.cell_corner_lon = nullptr;
  grid.cell_corner_lat = nullptr;
544
545
546
}

void
547
remapGridFree(RemapGrid &grid)
548
{
549
550
  varrayFree(grid.vgpm);
  varrayFree(grid.mask);
551
552
553
554
555
556
557
558
559
560
561

  if (grid.reg2d_center_lat) Free(grid.reg2d_center_lat);
  if (grid.reg2d_center_lon) Free(grid.reg2d_center_lon);
  if (grid.reg2d_corner_lat) Free(grid.reg2d_corner_lat);
  if (grid.reg2d_corner_lon) Free(grid.reg2d_corner_lon);

  if (grid.cell_center_lat) Free(grid.cell_center_lat);
  if (grid.cell_center_lon) Free(grid.cell_center_lon);
  if (grid.cell_corner_lat) Free(grid.cell_corner_lat);
  if (grid.cell_corner_lon) Free(grid.cell_corner_lon);

562
563
  varrayFree(grid.cell_area);
  varrayFree(grid.cell_frac);
564
565

  if (grid.tmpgridID != -1) gridDestroy(grid.tmpgridID);
566
567
568
}

void
569
remapSearchInit(RemapMethod mapType, RemapSearch &search, RemapGrid &src_grid, RemapGrid &tgt_grid)
570
{
571
  extern PointSearchMethod pointSearchMethod;
572
  extern CellSearchMethod cellSearchMethod;
573

574
575
576
577
578
  search.srcGrid = &src_grid;
  search.tgtGrid = &tgt_grid;

  search.srcBins.ncells = src_grid.size;
  search.tgtBins.ncells = tgt_grid.size;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
579

580
581
  search.srcBins.nbins = remap_num_srch_bins;
  search.tgtBins.nbins = remap_num_srch_bins;
582

583
  bool usePointsearch = mapType == RemapMethod::DISTWGT;
584
  if (src_grid.type != RemapGridType::Reg2D && pointSearchMethod != PointSearchMethod::latbins)
585
    {
586
587
      usePointsearch |= mapType == RemapMethod::BILINEAR;
      usePointsearch |= mapType == RemapMethod::BICUBIC;
588
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
589

590
  bool useCellsearch = (mapType == RemapMethod::CONSERV)
591
                       && (cellSearchMethod == CellSearchMethod::spherepart || src_grid.type == RemapGridType::Reg2D);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
592

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

595
  if (usePointsearch)
596
    {
597
      bool xIsCyclic = src_grid.is_cyclic;
598
      if (src_grid.type == RemapGridType::Reg2D)
599
        gridPointSearchCreateReg2d(search.gps, xIsCyclic, src_grid.dims, src_grid.reg2d_center_lon, src_grid.reg2d_center_lat);
600
      else
601
        gridPointSearchCreate(search.gps, xIsCyclic, src_grid.dims, src_grid.size, src_grid.cell_center_lon, src_grid.cell_center_lat);
602

603
      if (src_grid.lextrapolate) gridPointSearchExtrapolate(search.gps);
604

605
      if (Options::cdoVerbose) cdoPrint("Point search created: %.2f seconds", cdo_get_wtime() - start);
606
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
607
  else if (useCellsearch)
608
    {
609
      if (src_grid.type == RemapGridType::Reg2D)
610
        gridCellSearchCreateReg2d(search.gcs, src_grid.dims, src_grid.reg2d_corner_lon, src_grid.reg2d_corner_lat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
611
      else
612
        gridCellSearchCreate(search.gcs, src_grid.size, src_grid.num_cell_corners, src_grid.cell_corner_lon, src_grid.cell_corner_lat);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
613

614
      if (Options::cdoVerbose) cdoPrint("Cell search created: %.2f seconds", cdo_get_wtime() - start);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
615
    }
616
  else if (!(src_grid.type == RemapGridType::Reg2D || tgt_grid.type == RemapGridType::Reg2D))
617
    {
618
619
620
621
622
623
624
      search.srcBins.cell_bound_box.resize(4 * src_grid.size);
      if (tgt_grid.luse_cell_corners) search.tgtBins.cell_bound_box.resize(4 * tgt_grid.size);

      cell_bounding_boxes(src_grid, &search.srcBins.cell_bound_box[0], REMAP_GRID_BASIS_SRC);
      cell_bounding_boxes(tgt_grid, &search.tgtBins.cell_bound_box[0], REMAP_GRID_BASIS_TGT);
      // Set up and assign address ranges to search bins in order to further restrict later searches
      calcLatBins(search.srcBins);
625
      if (mapType == RemapMethod::CONSERV_SCRIP || mapType == RemapMethod::CONSERV)
626
        {
627
          calcLatBins(search.tgtBins);
628
629
630
          varrayFree(search.tgtBins.bin_lats);
          varrayFree(search.srcBins.bin_lats);
          if (mapType == RemapMethod::CONSERV) varrayFree(search.tgtBins.cell_bound_box);
631
        }
632

633
      if (Options::cdoVerbose) cdoPrint("Latitude bins created: %.2f seconds", cdo_get_wtime() - start);
634
635
636
637
638
639
    }
}

void
remapSearchFree(RemapSearch &search)
{
640
641
642
  varrayFree(search.srcBins.bin_addr);
  varrayFree(search.srcBins.bin_lats);
  varrayFree(search.srcBins.cell_bound_box);
643

644
645
646
  varrayFree(search.tgtBins.bin_addr);
  varrayFree(search.tgtBins.bin_lats);
  varrayFree(search.tgtBins.cell_bound_box);
647

648
  gridPointSearchDelete(search.gps);
649
  gridCellSearchDelete(search.gcs);
650
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
651

652
void
Uwe Schulzweida's avatar
Uwe Schulzweida committed
653
remapInitGrids(RemapMethod mapType, bool lextrapolate, int gridID1, RemapGrid &src_grid, int gridID2, RemapGrid &tgt_grid)
654
655
{
  int reg2d_src_gridID = gridID1;
656
  int reg2d_tgt_gridID = gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
657

658
659
  remapGridInit(src_grid);
  remapGridInit(tgt_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
660

661
  if (mapType == RemapMethod::BILINEAR || mapType == RemapMethod::BICUBIC || mapType == RemapMethod::DISTWGT
662
      || mapType == RemapMethod::CONSERV)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
663
    {
664
665
      if (IS_REG2D_GRID(gridID1)) src_grid.type = RemapGridType::Reg2D;
      // src_grid.type = 0;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
666
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
667

Uwe Schulzweida's avatar
Uwe Schulzweida committed
668
  //#ifdef YAC_CELL_SEARCH
669
  // if (IS_REG2D_GRID(gridID2) && mapType == RemapMethod::CONSERV) tgt_grid.type = RemapGridType::Reg2D;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
670
  //#else
671
  if (src_grid.type == RemapGridType::Reg2D)
672
    {
673
      if (IS_REG2D_GRID(gridID2) && mapType == RemapMethod::CONSERV) tgt_grid.type = RemapGridType::Reg2D;
674
      // else src_grid.type = -1;
675
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
676
  //#endif
677

678
  if (!remap_gen_weights && IS_REG2D_GRID(gridID2) && tgt_grid.type != RemapGridType::Reg2D)
679
    {
680
      if (mapType == RemapMethod::DISTWGT) tgt_grid.type = RemapGridType::Reg2D;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
681
      if (mapType == RemapMethod::BILINEAR && src_grid.type == RemapGridType::Reg2D) tgt_grid.type = RemapGridType::Reg2D;
682
683
    }

684
  src_grid.lextrapolate = lextrapolate;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
685

686
  if (mapType == RemapMethod::CONSERV_SCRIP || mapType == RemapMethod::CONSERV)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
687
    {
688
      if (src_grid.type != RemapGridType::Reg2D)
689
        {
690
691
          src_grid.luse_cell_corners = true;
          src_grid.lneed_cell_corners = true;
692
693
        }

694
      if (tgt_grid.type != RemapGridType::Reg2D)
695
        {
696
697
          tgt_grid.luse_cell_corners = true;
          tgt_grid.lneed_cell_corners = true;
698
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
699
700
    }

701
702
  src_grid.gridID = gridID1;
  tgt_grid.gridID = gridID2;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
703

704
  if (gridInqType(gridID1) == GRID_UNSTRUCTURED)
Uwe Schulzweida's avatar
Uwe Schulzweida committed
705
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
706
      if (gridInqYvals(gridID1, nullptr) == 0 || gridInqXvals(gridID1, nullptr) == 0)
707
        {
708
709
710
          int number = 0;
          cdiInqKeyInt(gridID1, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, &number);
          if (number > 0)
711
            {
712
              src_grid.gridID = gridID1 = referenceToGrid(gridID1);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
713
              if (gridID1 == -1) cdoAbort("Reference to source grid not found!");
714
715
            }
        }
716
    }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
717

718
  if (gridInqType(gridID2) == GRID_UNSTRUCTURED)
719
    {
Uwe Schulzweida's avatar
Uwe Schulzweida committed
720
      if (gridInqYvals(gridID2, nullptr) == 0 || gridInqXvals(gridID2, nullptr) == 0)
721
        {
722
723
724
          int number = 0;
          cdiInqKeyInt(gridID2, CDI_GLOBAL, CDI_KEY_NUMBEROFGRIDUSED, &number);
          if (number > 0)
725
            {
726
              tgt_grid.gridID = gridID2 = referenceToGrid(gridID2);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
727
              if (gridID2 == -1) cdoAbort("Reference to target grid not found!");
728
729
            }
        }
Uwe Schulzweida's avatar
Uwe Schulzweida committed
730
731
    }

732
  int sgridID = src_grid.gridID;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
733
734
735
736
737
738
  if (gridInqSize(sgridID) > 1
      && ((gridInqType(sgridID) == GRID_PROJECTION && gridInqProjType(sgridID) == CDI_PROJ_LCC)
          || (gridInqType(sgridID) == GRID_PROJECTION && gridInqProjType(sgridID) == CDI_PROJ_STERE)
          || (gridInqType(sgridID) == GRID_PROJECTION && gridInqProjType(sgridID) == CDI_PROJ_RLL)
          || (gridInqType(sgridID) == GRID_PROJECTION && gridInqProjType(sgridID) == CDI_PROJ_LAEA)
          || (gridInqType(sgridID) == GRID_PROJECTION && gridInqProjType(sgridID) == CDI_PROJ_SINU)))
739
    {
740
      bool lbounds = true;
741
      src_grid.gridID = gridID1 = gridToCurvilinear(src_grid.gridID, lbounds);
742
      src_grid.tmpgridID = src_grid.gridID;
743
    }
744

745
  // if ( src_grid.type != RemapGridType::Reg2D )
746
747
  remapDefineGrid(mapType, gridID1, src_grid, "Source");
  remapDefineGrid(mapType, gridID2, tgt_grid, "Target");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
748

749
  if (src_grid.type == RemapGridType::Reg2D || tgt_grid.type == RemapGridType::Reg2D)
750
    {
751
752
      if (src_grid.type == RemapGridType::Reg2D) remap_define_reg2d(reg2d_src_gridID, src_grid);
      if (tgt_grid.type == RemapGridType::Reg2D) remap_define_reg2d(reg2d_tgt_gridID, tgt_grid);
753
    }
754
}
Uwe Schulzweida's avatar
Uwe Schulzweida committed
755
756
757

/*****************************************************************************/

758
void
759
760
remapStat(int remapOrder, RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv, const Varray<double> &array1,
          const Varray<double> &array2, double missval)
761
{
762
  if (remapOrder == 2)
763
764
765
    cdoPrint("Second order mapping from grid1 to grid2:");
  else
    cdoPrint("First order mapping from grid1 to grid2:");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
766
  cdoPrint("----------------------------------------------");
Uwe Schulzweida's avatar
Uwe Schulzweida committed
767

768
769
  auto mmm = varrayMinMaxMeanMV(src_grid.size, array1, missval);
  cdoPrint(