remap.h 7.36 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-2019 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.
*/
17
18
#ifndef REMAP_H
#define REMAP_H
Uwe Schulzweida's avatar
Uwe Schulzweida committed
19

20
#include <stdint.h>
21
#include <cmath>
22
#include "remap_vars.h"
23
#include "remap_grid_cell_search.h"
24
#include "grid_point_search.h"
Uwe Schulzweida's avatar
Uwe Schulzweida committed
25

26
#ifndef M_PI
Uwe Schulzweida's avatar
Uwe Schulzweida committed
27
#define M_PI 3.14159265358979323846264338327950288
28
29
#endif

30
31
32
33
template <typename T>
inline void
vectorFree(std::vector<T> &v)
{
34
35
  v.clear();
  v.shrink_to_fit();
36
37
}

Uwe Schulzweida's avatar
Uwe Schulzweida committed
38
39
40
constexpr double PI = M_PI;
constexpr double PI2 = (2.0 * PI);
constexpr double PIH = (0.5 * PI);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
41
42
43
constexpr float PI_f = PI;
constexpr float PI2_f = PI2;
constexpr float PIH_f = PIH;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
44

45
#define TINY 1.e-14
Uwe Schulzweida's avatar
Uwe Schulzweida committed
46

Uwe Schulzweida's avatar
Uwe Schulzweida committed
47
48
49
50
51
52
53
enum struct RemapGridType
{
  Undefined,
  Reg2D,
  Curve2D,
  Unstruct
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
54

55
56
#define REMAP_GRID_BASIS_SRC 1
#define REMAP_GRID_BASIS_TGT 2
Uwe Schulzweida's avatar
Uwe Schulzweida committed
57

58
59
60
61
62
63
64
enum struct SubmapType
{
  NONE,
  LAF,
  SUM
};

65
struct RemapGrid
66
67
{
  int gridID;
68
  int tmpgridID;
69
  RemapGridType type;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
70
71
72
  int rank;                 // rank of the grid
  size_t size;              // total points on the grid
  size_t num_cell_corners;  // number of corners for each grid cell
73
74

  bool lneed_cell_corners;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
75
  bool luse_cell_corners;  // use corners for bounding boxes
76
77
78
79
80

  bool lextrapolate;
  bool non_global;
  bool is_cyclic;

Uwe Schulzweida's avatar
Uwe Schulzweida committed
81
  size_t dims[2];  // size of grid dimension
82

Uwe Schulzweida's avatar
Uwe Schulzweida committed
83
84
  int nvgp;               // size of vgpm
  std::vector<int> vgpm;  // flag which cells are valid
85

Uwe Schulzweida's avatar
Uwe Schulzweida committed
86
  std::vector<int> mask;  // flag which cells participate
87

Uwe Schulzweida's avatar
Uwe Schulzweida committed
88
89
90
91
  double *reg2d_center_lon;  // reg2d lon/lat coordinates for
  double *reg2d_center_lat;  // each grid center in radians
  double *reg2d_corner_lon;  // reg2d lon/lat coordinates for
  double *reg2d_corner_lat;  // each grid corner in radians
92

93
94
95
96
  double *cell_center_lon;  // lon/lat coordinates for
  double *cell_center_lat;  // each grid center in radians
  double *cell_corner_lon;  // lon/lat coordinates for
  double *cell_corner_lat;  // each grid corner in radians
97

Uwe Schulzweida's avatar
Uwe Schulzweida committed
98
99
  std::vector<double> cell_area;  // total area of each grid cell
  std::vector<double> cell_frac;  // fractional area of grid cells participating in remapping
100
};
101

Uwe Schulzweida's avatar
Uwe Schulzweida committed
102
103
struct GridSearchBins
{
Uwe Schulzweida's avatar
Uwe Schulzweida committed
104
105
106
107
108
  unsigned nbins;                     // num of bins for restricted search
  size_t ncells;                      // total number of grid cells (cell_bound_box)
  std::vector<size_t> bin_addr;       // min,max adds for grid cells in this lat bin
  std::vector<float> bin_lats;        // min,max latitude for each search bin
  std::vector<float> cell_bound_box;  // lon/lat bounding box for use
Uwe Schulzweida's avatar
Uwe Schulzweida committed
109
110
};

111
struct RemapSearch
Uwe Schulzweida's avatar
Uwe Schulzweida committed
112
{
113
114
115
116
117
118
  RemapGrid *srcGrid;
  RemapGrid *tgtGrid;

  GridSearchBins srcBins;
  GridSearchBins tgtBins;

119
  GridPointSearch *gps;
120
  GridCellSearch *gcs;
Uwe Schulzweida's avatar
Uwe Schulzweida committed
121
};
122

123
struct remapType
124
125
126
127
128
{
  int nused;
  int gridID;
  size_t gridsize;
  size_t nmiss;
129
130
  RemapGrid src_grid;
  RemapGrid tgt_grid;
131
  RemapVars vars;
132
  RemapSearch search;
133
};
Uwe Schulzweida's avatar
Uwe Schulzweida committed
134

135
136
137
138
#define REMAP_WRITE_REMAP 2
#define REMAP_MAX_ITER 3
#define REMAP_NUM_SRCH_BINS 4
#define REMAP_GENWEIGHTS 5
139

140
void remap_set_threshhold(double threshhold);
141
void remap_set_int(int remapvar, int value);
142

Uwe Schulzweida's avatar
Uwe Schulzweida committed
143
void remapInitGrids(RemapMethod mapType, bool lextrapolate, int gridID1, RemapGrid &src_grid, int gridID2, RemapGrid &tgt_grid);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
144

145
146
void remapGridInit(RemapGrid &grid);
void remapGridFree(RemapGrid &grid);
147
148
void remapGridAlloc(RemapMethod mapType, RemapGrid &grid);
void remapSearchInit(RemapMethod mapType, RemapSearch &search, RemapGrid &src_grid, RemapGrid &tgt_grid);
149
void remapSearchFree(RemapSearch &search);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
150

Uwe Schulzweida's avatar
Uwe Schulzweida committed
151
void remapSearchPoints(RemapSearch &rsearch, double plon, double plat, knnWeightsType &knnWeights);
152
int remapSearchSquare(RemapSearch &rsearch, double plon, double plat, size_t *src_add, double *src_lats, double *src_lons);
153
size_t remapSearchCells(RemapSearch &rsearch, bool isReg2dCell, grid_cell &gridCell, size_t *srchAddr);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
154

155
156
void remapBilinearWeights(RemapSearch &rsearch, RemapVars &rv);
void remapBicubicWeights(RemapSearch &rsearch, RemapVars &rv);
157
void remapDistwgtWeights(size_t numNeighbors, RemapSearch &rsearch, RemapVars &rv);
158
void remapConservWeights(RemapSearch &rsearch, RemapVars &rv);
159
void remapConservWeightsScrip(RemapSearch &rsearch, RemapVars &rv);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
160

Uwe Schulzweida's avatar
Uwe Schulzweida committed
161
162
163
164
void remapBilinear(RemapSearch &rsearch, const double *src_array, double *tgt_array, double missval);
void remapBicubic(RemapSearch &rsearch, const double *src_array, double *tgt_array, double missval);
void remapDistwgt(size_t numNeighbors, RemapSearch &rsearch, const double *src_array, double *tgt_array, double missval);
void remapConserv(RemapSearch &rsearch, const double *src_array, double *tgt_array, double missval);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
165

Uwe Schulzweida's avatar
Uwe Schulzweida committed
166
167
void remapStat(int remapOrder, RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv, const double *array1, const double *array2,
               double missval);
168
void remapGradients(RemapGrid &grid, const std::vector<bool> &mask, const double *array, RemapGradients &gradients);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
169
void remapGradients(RemapGrid &grid, const double *array, RemapGradients &gradients);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
170

Uwe Schulzweida's avatar
Uwe Schulzweida committed
171
172
void sort_add(size_t num_links, size_t num_wts, size_t *add1, size_t *add2, double *weights);
void sort_iter(size_t num_links, size_t num_wts, size_t *add1, size_t *add2, double *weights, int parent);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
173

Uwe Schulzweida's avatar
Uwe Schulzweida committed
174
175
void remapWriteDataScrip(const char *interp_file, RemapMethod mapType, SubmapType submapType, int numNeighbors, int remapOrder,
                         RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv);
176
void remapReadDataScrip(const char *interp_file, int gridID1, int gridID2, RemapMethod *mapType, SubmapType *submapType,
177
                        int *numNeighbors, int &remapOrder, RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv);
178

179
void calcLatBins(GridSearchBins &searchBins);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
180
181
size_t get_srch_cells(size_t tgt_cell_addr, GridSearchBins &tgtBins, GridSearchBins &srcBins, float *tgt_cell_bound_box,
                      size_t *srch_add);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
182

Uwe Schulzweida's avatar
Uwe Schulzweida committed
183
184
int gridSearchSquareReg2dNN(size_t nx, size_t ny, size_t *nbr_add, double *nbr_dist, double plat, double plon,
                            const double *src_center_lat, const double *src_center_lon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
185

Uwe Schulzweida's avatar
Uwe Schulzweida committed
186
int gridSearchSquareReg2d(RemapGrid *src_grid, size_t *src_add, double *src_lats, double *src_lons, double plat, double plon);
Uwe Schulzweida's avatar
Uwe Schulzweida committed
187

188
bool pointInQuad(bool isCyclic, size_t nx, size_t ny, size_t i, size_t j, size_t adds[4], double lons[4], double lats[4],
Uwe Schulzweida's avatar
Uwe Schulzweida committed
189
                 double plon, double plat, const double *centerLon, const double *centerLat);
190

Uwe Schulzweida's avatar
Uwe Schulzweida committed
191
192
int gridSearchSquareCurv2dScrip(RemapGrid *src_grid, size_t *src_add, double *src_lats, double *src_lons, double plat, double plon,
                                GridSearchBins &srcBins);
193

Uwe Schulzweida's avatar
Uwe Schulzweida committed
194
195
bool remapFindWeights(double plon, double plat, double *src_lons, double *src_lats, double *ig, double *jg);
int rect_grid_search(size_t *ii, size_t *jj, double x, double y, size_t nxm, size_t nym, const double *xm, const double *ym);
196

197
void remapgrid_get_lonlat(RemapGrid *grid, size_t cell_add, double *plon, double *plat);
198

Uwe Schulzweida's avatar
Uwe Schulzweida committed
199
void remapCheckArea(size_t grid_size, double *cell_area, const char *name);
200
201

#endif /* REMAP_H */