1 /* 2 This file is part of CDO. CDO is a collection of Operators to manipulate and analyse Climate model Data. 3 4 Author: Uwe Schulzweida 5 6 */ 7 #ifndef REMAP_H 8 #define REMAP_H 9 10 #include <cstdint> 11 #include <cmath> 12 13 #include "varray.h" 14 #include "remap_vars.h" 15 #include "remap_grid_cell_search.h" 16 #include "grid_point_search.h" 17 18 #ifndef M_PI 19 #define M_PI 3.14159265358979323846264338327950288 20 #endif 21 22 constexpr double PI = M_PI; 23 constexpr double PI2 = (2.0 * PI); 24 constexpr double PIH = (0.5 * PI); 25 constexpr float PI_f = PI; 26 constexpr float PI2_f = PI2; 27 constexpr float PIH_f = PIH; 28 29 constexpr double TINY = 1.e-14; 30 31 enum class RemapGridType 32 { 33 Undefined, 34 Reg2D, 35 Curve2D, 36 Unstruct 37 }; 38 39 #define REMAP_GRID_BASIS_SRC 1 40 #define REMAP_GRID_BASIS_TGT 2 41 42 enum class SubmapType 43 { 44 NONE, 45 LAF, 46 SUM, 47 AVG 48 }; 49 50 struct LonLatPoint 51 { 52 double lon, lat; LonLatPointLonLatPoint53 LonLatPoint() {}; LonLatPointLonLatPoint54 LonLatPoint(double _lon, double _lat) : lon(_lon), lat(_lat) {}; 55 }; 56 57 struct // RemapGrid 58 #ifdef WARN_UNUSED 59 [[gnu::warn_unused]] 60 #endif 61 RemapGrid 62 { 63 int gridID; 64 int tmpgridID; 65 RemapGridType type; 66 int rank; // rank of the grid 67 size_t size; // total points on the grid 68 size_t num_cell_corners; // number of corners for each grid cell 69 70 bool lneed_cell_corners; 71 bool luse_cell_corners; // use corners for bounding boxes 72 73 bool lextrapolate; 74 bool non_global; 75 bool is_cyclic; 76 77 size_t dims[2]; // size of grid dimension 78 79 int nvgp; // size of vgpm 80 Varray<int> vgpm; // flag which cells are valid 81 Varray<short> mask; // flag which cells participate 82 83 Varray<double> reg2d_center_lon; // reg2d lon/lat coordinates for 84 Varray<double> reg2d_center_lat; // each grid center in radians 85 Varray<double> reg2d_corner_lon; // reg2d lon/lat coordinates for 86 Varray<double> reg2d_corner_lat; // each grid corner in radians 87 88 Varray<double> cell_center_lon; // lon/lat coordinates for 89 Varray<double> cell_center_lat; // each grid center in radians 90 Varray<double> cell_corner_lon; // lon/lat coordinates for 91 Varray<double> cell_corner_lat; // each grid corner in radians 92 93 Varray<double> cell_area; // total area of each grid cell 94 Varray<double> cell_frac; // fractional area of grid cells participating in remapping 95 }; 96 97 struct // GridSearchBins 98 #ifdef WARN_UNUSED 99 [[gnu::warn_unused]] 100 #endif 101 GridSearchBins 102 { 103 unsigned nbins; // num of bins for restricted search 104 size_t ncells; // total number of grid cells (cell_bound_box) 105 Varray<size_t> bin_addr; // min,max adds for grid cells in this lat bin 106 Varray<float> bin_lats; // min,max latitude for each search bin 107 Varray<float> cell_bound_box; // lon/lat bounding box for use 108 }; 109 110 struct // RemapSearch 111 #ifdef WARN_UNUSED 112 [[gnu::warn_unused]] 113 #endif 114 RemapSearch 115 { 116 RemapGrid *srcGrid; 117 RemapGrid *tgtGrid; 118 119 GridSearchBins srcBins; 120 GridSearchBins tgtBins; 121 122 GridPointSearch gps; 123 GridCellSearch gcs; 124 }; 125 126 struct // RemapType 127 #ifdef WARN_UNUSED 128 [[gnu::warn_unused]] 129 #endif 130 RemapType 131 { 132 int nused; 133 int gridID; 134 size_t gridsize; 135 size_t nmiss; 136 RemapGrid src_grid; 137 RemapGrid tgt_grid; 138 RemapVars vars; 139 RemapSearch search; 140 }; 141 142 #define REMAP_WRITE_REMAP 2 143 #define REMAP_MAX_ITER 3 144 #define REMAP_NUM_SRCH_BINS 4 145 #define REMAP_GENWEIGHTS 5 146 147 void remap_set_threshhold(double threshhold); 148 void remap_set_int(int remapvar, int value); 149 150 void remap_init_grids(RemapMethod mapType, bool lextrapolate, int gridID1, RemapGrid &src_grid, int gridID2, RemapGrid &tgt_grid); 151 152 void remap_grid_init(RemapGrid &grid); 153 void remap_grid_free(RemapGrid &grid); 154 void remap_grid_alloc(RemapMethod mapType, RemapGrid &grid); 155 void remap_search_init(RemapMethod mapType, RemapSearch &search, RemapGrid &src_grid, RemapGrid &tgt_grid); 156 void remap_search_free(RemapSearch &search); 157 158 void remap_search_points(RemapSearch &rsearch, const LonLatPoint &llp, knnWeightsType &knnWeights); 159 int remap_search_square(RemapSearch &rsearch, const LonLatPoint &llp, size_t (&src_add)[4], double (&src_lats)[4], double (&src_lons)[4]); 160 size_t remap_search_cells(RemapSearch &rsearch, bool isReg2dCell, grid_cell &gridCell, Varray<size_t> &srchAddr); 161 162 void remap_bilinear_weights(RemapSearch &rsearch, RemapVars &rv); 163 void remap_bicubic_weights(RemapSearch &rsearch, RemapVars &rv); 164 void remap_distwgt_weights(size_t numNeighbors, RemapSearch &rsearch, RemapVars &rv); 165 void remap_conserv_weights(RemapSearch &rsearch, RemapVars &rv); 166 void remap_conserv_weights_scrip(RemapSearch &rsearch, RemapVars &rv); 167 168 void remap_bilinear(RemapSearch &rsearch, const Field &field1, Field &field2); 169 void remap_bicubic(RemapSearch &rsearch, const Field &field1, Field &field2); 170 void remap_dist_wgt(size_t numNeighbors, RemapSearch &rsearch, const Field &field1, Field &field2); 171 void remap_conserv(NormOpt normOpt, RemapSearch &rsearch, const Field &field1, Field &field2); 172 173 void remap_stat(int remapOrder, RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv, const Field &field1, const Field &field2); 174 175 template <typename T> 176 void remap_gradients(RemapGrid &grid, const Varray<short> &mask, const Varray<T> &array, RemapGradients &gradients); 177 void remap_gradients(RemapGrid &grid, const Field &field, RemapGradients &gradients); 178 179 void sort_add(size_t num_links, size_t num_wts, size_t *add1, size_t *add2, double *weights); 180 void sort_iter(size_t num_links, size_t num_wts, size_t *add1, size_t *add2, double *weights, int parent); 181 182 void remap_write_data_scrip(const char *interp_file, RemapMethod mapType, SubmapType submapType, int numNeighbors, int remapOrder, 183 RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv); 184 void remap_read_data_scrip(const char *interp_file, int gridID1, int gridID2, RemapMethod &mapType, SubmapType &submapType, 185 int &numNeighbors, int &remapOrder, RemapGrid &src_grid, RemapGrid &tgt_grid, RemapVars &rv); 186 187 void calc_lat_bins(GridSearchBins &searchBins); 188 size_t get_srch_cells(size_t tgt_cell_addr, GridSearchBins &tgtBins, GridSearchBins &srcBins, float *tgt_cell_bound_box, 189 Varray<size_t> &srch_add); 190 191 int grid_search_square_reg_2d_NN (size_t nx, size_t ny, size_t *nbr_add, double *nbr_dist, double plat, double plon, 192 const Varray<double> &src_center_lat, const Varray<double> &src_center_lon); 193 194 int grid_search_square_reg_2d(RemapGrid *src_grid, size_t (&src_add)[4], double (&src_lats)[4], double (&src_lons)[4], double plat, double plon); 195 196 bool point_in_quad(bool isCyclic, size_t nx, size_t ny, size_t i, size_t j, size_t adds[4], double lons[4], double lats[4], 197 double plon, double plat, const double *centerLon, const double *centerLat); 198 199 int grid_search_square_curv_2d_scrip(RemapGrid *src_grid, size_t (&src_add)[4], double (&src_lats)[4], double (&src_lons)[4], double plat, double plon, 200 GridSearchBins &srcBins); 201 202 bool remap_find_weights(const LonLatPoint &llpoint, const double (&src_lons)[4], const double (&src_lats)[4], double *ig, double *jg); 203 int rect_grid_search(size_t &ii, size_t &jj, double x, double y, size_t nxm, size_t nym, const Varray<double> &xm, const Varray<double> &ym); 204 205 LonLatPoint remapgrid_get_lonlat(RemapGrid *grid, size_t cell_add); 206 207 void remap_check_area(size_t grid_size, double *cell_area, const char *name); 208 209 #endif /* REMAP_H */ 210