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
8 #include <cdi.h>
9
10 #include "cdo_math.h"
11 #include "remap.h"
12 #include <mpim_grid.h>
13
14 static size_t
fill_src_add(bool is_cyclic,long nx,long ny,long ii,long jj,long k,size_t * psrc_add)15 fill_src_add(bool is_cyclic, long nx, long ny, long ii, long jj, long k, size_t *psrc_add)
16 {
17 k /= 2;
18
19 auto j0 = jj - k;
20 auto jn = jj + k;
21 auto i0 = ii - k;
22 auto in = ii + k;
23 if (j0 < 0) j0 = 0;
24 if (jn >= ny) jn = ny - 1;
25 if ((in - i0) > nx)
26 {
27 i0 = 0;
28 in = nx - 1;
29 }
30
31 size_t num_add = 0;
32
33 for (long j = j0; j <= jn; ++j)
34 for (long i = i0; i <= in; ++i)
35 {
36 auto ix = i;
37 if (is_cyclic && ix < 0) ix += nx;
38 if (is_cyclic && ix >= nx) ix -= nx;
39 if (ix >= 0 && ix < nx && j >= 0 && j < ny) psrc_add[num_add++] = j * nx + ix;
40 }
41
42 return num_add;
43 }
44
45 static void
store_distance(GridPointSearch & gps,double plon,double plat,knnWeightsType & knnWeights,size_t nx,size_t num_add,size_t * psrc_add)46 store_distance(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights, size_t nx, size_t num_add,
47 size_t *psrc_add)
48 {
49 const auto &coslon = gps.coslon;
50 const auto &sinlon = gps.sinlon;
51 const auto &coslat = gps.coslat;
52 const auto &sinlat = gps.sinlat;
53
54 double xyz[3], query_pt[3];
55 gcLLtoXYZ(plon, plat, query_pt);
56 const auto sqrSearchRadius = cdo::sqr(gps.searchRadius);
57
58 for (size_t na = 0; na < num_add; ++na)
59 {
60 const auto nadd = psrc_add[na];
61 const auto iy = nadd / nx;
62 const auto ix = nadd - iy * nx;
63
64 xyz[0] = coslat[iy] * coslon[ix];
65 xyz[1] = coslat[iy] * sinlon[ix];
66 xyz[2] = sinlat[iy];
67 // Find distance to this point
68 const double sqrDist = (float) squareDistance(query_pt, xyz);
69 if (sqrDist <= sqrSearchRadius)
70 {
71 // Store the address and distance if this is one of the smallest so far
72 knnWeights.storeDistance(nadd, std::sqrt(sqrDist));
73 }
74 }
75
76 knnWeights.checkDistance();
77 }
78
79 // This routine finds the closest numNeighbor points to a search point and computes a distance to each of the neighbors
80 static void
gridSearchPointReg2d(GridPointSearch & gps,double plon,double plat,knnWeightsType & knnWeights)81 gridSearchPointReg2d(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights)
82 {
83 /*
84 Input variables:
85
86 plat : latitude of the search point
87 plon : longitude of the search point
88
89 Output variables:
90
91 knnWeights.m_addr[numNeighbors] : address of each of the closest points
92 knnWeights.m_dist[numNeighbors] : distance to each of the closest points
93 */
94 const auto numNeighbors = knnWeights.maxNeighbors();
95 auto &nbr_add = knnWeights.m_addr;
96 auto &nbr_dist = knnWeights.m_dist;
97
98 // Initialize distance and address arrays
99 knnWeights.initAddr();
100 knnWeights.initDist();
101
102 const auto &src_center_lon = gps.reg2d_center_lon;
103 const auto &src_center_lat = gps.reg2d_center_lat;
104
105 long nx = gps.dims[0];
106 long ny = gps.dims[1];
107 size_t nxm = gps.is_cyclic ? nx + 1 : nx;
108
109 if (plon < src_center_lon[0]) plon += PI2;
110 if (plon > src_center_lon[nxm - 1]) plon -= PI2;
111
112 size_t ii, jj;
113 auto lfound = rect_grid_search(ii, jj, plon, plat, nxm, ny, src_center_lon, src_center_lat);
114 if (lfound)
115 {
116 if (gps.is_cyclic && ii == (nxm - 1)) ii = 0;
117
118 constexpr size_t MAX_SEARCH_CELLS = 25;
119 size_t src_add[MAX_SEARCH_CELLS];
120 size_t *psrc_add = src_add;
121
122 size_t k;
123 for (k = 3; k < 10000; k += 2)
124 if (numNeighbors <= (size_t)(k - 2) * (k - 2)) break;
125
126 std::vector<size_t> src_add_tmp;
127 if ((k * k) > MAX_SEARCH_CELLS)
128 {
129 src_add_tmp.resize(k * k);
130 psrc_add = src_add_tmp.data();
131 }
132
133 const auto num_add = fill_src_add(gps.is_cyclic, nx, ny, ii, jj, k, psrc_add);
134
135 store_distance(gps, plon, plat, knnWeights, nx, num_add, psrc_add);
136 }
137 else if (gps.extrapolate)
138 {
139 int search_result = 0;
140
141 if (numNeighbors < 4)
142 {
143 size_t nbr_add4[4];
144 double nbr_dist4[4];
145 for (size_t n = 0; n < numNeighbors; ++n) nbr_add4[n] = SIZE_MAX;
146 search_result = grid_search_square_reg_2d_NN(nx, ny, nbr_add4, nbr_dist4, plat, plon, src_center_lat, src_center_lon);
147 if (search_result < 0)
148 {
149 for (size_t n = 0; n < numNeighbors; ++n) nbr_add[n] = nbr_add4[n];
150 for (size_t n = 0; n < numNeighbors; ++n) nbr_dist[n] = nbr_dist4[n];
151 }
152 }
153 else
154 {
155 search_result
156 = grid_search_square_reg_2d_NN(nx, ny, nbr_add.data(), nbr_dist.data(), plat, plon, src_center_lat, src_center_lon);
157 }
158
159 if (search_result >= 0)
160 for (size_t n = 0; n < numNeighbors; ++n) nbr_add[n] = SIZE_MAX;
161 }
162 }
163
164 void
grid_search_point(GridPointSearch & gps,double plon,double plat,knnWeightsType & knnWeights)165 grid_search_point(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights)
166 {
167 /*
168 Input variables:
169
170 plat : latitude of the search point
171 plon : longitude of the search point
172
173 Output variables:
174
175 knnWeights.m_addr[numNeighbors] : address of each of the closest points
176 knnWeights.m_dist[numNeighbors] : distance to each of the closest points
177 */
178 auto numNeighbors = knnWeights.maxNeighbors();
179
180 // check some more points if distance is the same use the smaller index (nadd)
181 auto ndist = (numNeighbors > 8) ? numNeighbors + 8 : numNeighbors * 2;
182 if (ndist > gps.n) ndist = gps.n;
183
184 if (knnWeights.m_tmpaddr.empty()) knnWeights.m_tmpaddr.resize(ndist);
185 if (knnWeights.m_tmpdist.empty()) knnWeights.m_tmpdist.resize(ndist);
186 auto &adds = knnWeights.m_tmpaddr;
187 auto &dist = knnWeights.m_tmpdist;
188
189 size_t nadds = 0;
190 if (numNeighbors == 1)
191 nadds = grid_point_search_nearest(gps, plon, plat, adds.data(), dist.data());
192 else
193 nadds = grid_point_search_qnearest(gps, plon, plat, ndist, adds.data(), dist.data());
194
195 ndist = nadds;
196 if (ndist < numNeighbors) numNeighbors = ndist;
197
198 // Initialize distance and address arrays
199 knnWeights.initAddr();
200 knnWeights.initDist();
201 for (size_t i = 0; i < ndist; ++i) knnWeights.storeDistance(adds[i], dist[i], numNeighbors);
202
203 knnWeights.checkDistance();
204 }
205
206 void
grid_search_point_smooth(GridPointSearch & gps,double plon,double plat,knnWeightsType & knnWeights)207 grid_search_point_smooth(GridPointSearch &gps, double plon, double plat, knnWeightsType &knnWeights)
208 {
209 /*
210 Input variables:
211
212 plat : latitude of the search point
213 plon : longitude of the search point
214
215 Output variables:
216
217 knnWeights.m_addr[numNeighbors] : address of each of the closest points
218 knnWeights.m_dist[numNeighbors] : distance to each of the closest points
219 */
220 auto numNeighbors = knnWeights.maxNeighbors();
221 const auto checkDistance = (numNeighbors <= 32);
222
223 // check some more points if distance is the same use the smaller index (nadd)
224 auto ndist = checkDistance ? ((numNeighbors > 8) ? numNeighbors + 8 : numNeighbors * 2) : numNeighbors;
225 if (ndist > gps.n) ndist = gps.n;
226
227 if (knnWeights.m_tmpaddr.empty()) knnWeights.m_tmpaddr.resize(ndist);
228 if (knnWeights.m_tmpdist.empty()) knnWeights.m_tmpdist.resize(ndist);
229 auto &adds = knnWeights.m_tmpaddr;
230 auto &dist = knnWeights.m_tmpdist;
231
232 size_t nadds = 0;
233 if (numNeighbors == 1)
234 nadds = grid_point_search_nearest(gps, plon, plat, adds.data(), dist.data());
235 else
236 nadds = grid_point_search_qnearest(gps, plon, plat, ndist, adds.data(), dist.data());
237
238 ndist = nadds;
239
240 if (checkDistance)
241 {
242 if (ndist < numNeighbors) numNeighbors = ndist;
243
244 // Initialize distance and address arrays
245 knnWeights.initAddr(numNeighbors);
246 knnWeights.initDist(numNeighbors);
247 for (size_t i = 0; i < ndist; ++i) knnWeights.storeDistance(adds[i], dist[i], numNeighbors);
248 }
249 else
250 {
251 knnWeights.m_numNeighbors = ndist;
252 for (size_t i = 0; i < ndist; ++i) knnWeights.m_addr[i] = adds[i];
253 for (size_t i = 0; i < ndist; ++i) knnWeights.m_dist[i] = dist[i];
254 }
255
256 knnWeights.checkDistance();
257 }
258
259 void
remap_search_points(RemapSearch & rsearch,const LonLatPoint & llp,knnWeightsType & knnWeights)260 remap_search_points(RemapSearch &rsearch, const LonLatPoint &llp, knnWeightsType &knnWeights)
261 {
262 if (rsearch.srcGrid->type == RemapGridType::Reg2D)
263 gridSearchPointReg2d(rsearch.gps, llp.lon, llp.lat, knnWeights);
264 else
265 grid_search_point(rsearch.gps, llp.lon, llp.lat, knnWeights);
266 }
267
268 static int
gridSearchSquareCurv2d(GridPointSearch & gps,RemapGrid * rgrid,size_t (& src_add)[4],double (& src_lats)[4],double (& src_lons)[4],double plat,double plon)269 gridSearchSquareCurv2d(GridPointSearch &gps, RemapGrid *rgrid, size_t (&src_add)[4], double (&src_lats)[4], double (&src_lons)[4],
270 double plat, double plon)
271 {
272 /*
273 Input variables:
274
275 plat : latitude of the search point
276 plon : longitude of the search point
277
278 Output variables:
279
280 src_add[4] : address of each corner point enclosing P
281 src_lats[4] : latitudes of the four corner points
282 src_lons[4] : longitudes of the four corner points
283 */
284 int search_result = 0;
285
286 for (unsigned n = 0; n < 4; ++n) src_add[n] = 0;
287
288 double dist;
289 size_t addr;
290 size_t nadds = grid_point_search_nearest(gps, plon, plat, &addr, &dist);
291 if (nadds > 0)
292 {
293 const auto nx = rgrid->dims[0];
294 const auto ny = rgrid->dims[1];
295
296 for (unsigned k = 0; k < 4; ++k)
297 {
298 // Determine neighbor addresses
299 auto j = addr / nx;
300 auto i = addr - j * nx;
301 if (k == 0 || k == 2) i = (i > 0) ? i - 1 : rgrid->is_cyclic ? nx - 1 : 0;
302 if (k == 0 || k == 1) j = (j > 0) ? j - 1 : 0;
303 if (point_in_quad(rgrid->is_cyclic, nx, ny, i, j, src_add, src_lons, src_lats, plon, plat, rgrid->cell_center_lon.data(),
304 rgrid->cell_center_lat.data()))
305 {
306 search_result = 1;
307 return search_result;
308 }
309 }
310 }
311
312 /*
313 If no cell found, point is likely either in a box that straddles either pole or is outside the grid.
314 Fall back to a distance-weighted average of the four closest points. Go ahead and compute weights here,
315 but store in src_lats and return -add to prevent the parent routine from computing bilinear weights.
316 */
317 if (!rgrid->lextrapolate) return search_result;
318
319 size_t ndist = 4;
320 nadds = grid_point_search_qnearest(gps, plon, plat, ndist, src_add, src_lats);
321 if (nadds == 4)
322 {
323 for (unsigned n = 0; n < 4; ++n) src_lats[n] = 1.0 / (src_lats[n] + TINY);
324 double distance = 0.0;
325 for (unsigned n = 0; n < 4; ++n) distance += src_lats[n];
326 for (unsigned n = 0; n < 4; ++n) src_lats[n] /= distance;
327 search_result = -1;
328 }
329
330 return search_result;
331 }
332
333 int
remap_search_square(RemapSearch & rsearch,const LonLatPoint & llp,size_t (& src_add)[4],double (& src_lats)[4],double (& src_lons)[4])334 remap_search_square(RemapSearch &rsearch, const LonLatPoint &llp, size_t (&src_add)[4], double (&src_lats)[4],
335 double (&src_lons)[4])
336 {
337 if (rsearch.srcGrid->type == RemapGridType::Reg2D)
338 return grid_search_square_reg_2d(rsearch.srcGrid, src_add, src_lats, src_lons, llp.lat, llp.lon);
339 else if (rsearch.gps.in_use)
340 return gridSearchSquareCurv2d(rsearch.gps, rsearch.srcGrid, src_add, src_lats, src_lons, llp.lat, llp.lon);
341 else
342 return grid_search_square_curv_2d_scrip(rsearch.srcGrid, src_add, src_lats, src_lons, llp.lat, llp.lon, rsearch.srcBins);
343 }
344