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 "remap.h"
9
10 static void
nbr_add_and_dist(int & search_result,size_t jj,size_t ii,size_t nx,double & dist_min,double distance,size_t * nbr_add,double * nbr_dist)11 nbr_add_and_dist(int &search_result, size_t jj, size_t ii, size_t nx, double &dist_min, double distance, size_t *nbr_add,
12 double *nbr_dist)
13 {
14 if (distance < dist_min)
15 {
16 const auto srch_add = jj * nx + ii;
17 for (unsigned n = 0; n < 4; ++n)
18 {
19 if (distance < nbr_dist[n])
20 {
21 for (unsigned i = 3; i > n; --i) nbr_add[i] = nbr_add[i - 1];
22 for (unsigned i = 3; i > n; --i) nbr_dist[i] = nbr_dist[i - 1];
23 search_result = -1;
24 nbr_add[n] = srch_add;
25 nbr_dist[n] = distance;
26 dist_min = nbr_dist[3];
27 break;
28 }
29 }
30 }
31 }
32
33 int
grid_search_square_reg_2d_NN(size_t nx,size_t ny,size_t * nbr_add,double * nbr_dist,double plat,double plon,const Varray<double> & src_center_lat,const Varray<double> & src_center_lon)34 grid_search_square_reg_2d_NN(size_t nx, size_t ny, size_t *nbr_add, double *nbr_dist, double plat, double plon,
35 const Varray<double> &src_center_lat, const Varray<double> &src_center_lon)
36 {
37 int search_result = 0;
38
39 const auto coslat_dst = std::cos(plat);
40 const auto sinlat_dst = std::sin(plat);
41 const auto coslon_dst = std::cos(plon);
42 const auto sinlon_dst = std::sin(plon);
43
44 double dist_min = DBL_MAX;
45 for (unsigned n = 0; n < 4; ++n) nbr_dist[n] = DBL_MAX;
46
47 size_t jjf = 0, jjl = ny - 1;
48 if (plon >= src_center_lon[0] && plon <= src_center_lon[nx - 1])
49 {
50 if (src_center_lat[0] < src_center_lat[ny - 1])
51 {
52 if (plat <= src_center_lat[0])
53 jjl = (ny == 1) ? 0 : 1;
54 else
55 jjf = (ny == 1) ? 0 : ny - 2;
56 }
57 else
58 {
59 if (plat >= src_center_lat[0])
60 jjl = (ny == 1) ? 0 : 1;
61 else
62 jjf = (ny == 1) ? 0 : ny - 2;
63 }
64 }
65
66 std::vector<double> sincoslon(nx);
67
68 for (size_t ii = 0; ii < nx; ++ii)
69 sincoslon[ii] = coslon_dst * std::cos(src_center_lon[ii]) + sinlon_dst * std::sin(src_center_lon[ii]);
70
71 for (size_t jj = jjf; jj <= jjl; ++jj)
72 {
73 const auto coslat = coslat_dst * std::cos(src_center_lat[jj]);
74 const auto sinlat = sinlat_dst * std::sin(src_center_lat[jj]);
75
76 const auto jjskip = (jj > 1 && jj < (ny - 2));
77
78 if (jjskip)
79 {
80 size_t ii = 0;
81 auto distance = std::acos(coslat * sincoslon[ii] + sinlat);
82 nbr_add_and_dist(search_result, jj, ii, nx, dist_min, distance, nbr_add, nbr_dist);
83 ii = nx - 1;
84 distance = std::acos(coslat * sincoslon[ii] + sinlat);
85 nbr_add_and_dist(search_result, jj, ii, nx, dist_min, distance, nbr_add, nbr_dist);
86 }
87 else
88 {
89 for (size_t ii = 0; ii < nx; ++ii)
90 {
91 const auto distance = std::acos(coslat * sincoslon[ii] + sinlat);
92 nbr_add_and_dist(search_result, jj, ii, nx, dist_min, distance, nbr_add, nbr_dist);
93 }
94 }
95 }
96
97 for (unsigned n = 0; n < 4; ++n) nbr_dist[n] = 1.0 / (nbr_dist[n] + TINY);
98 double distance = 0.0;
99 for (unsigned n = 0; n < 4; ++n) distance += nbr_dist[n];
100 for (unsigned n = 0; n < 4; ++n) nbr_dist[n] /= distance;
101
102 return search_result;
103 }
104
105 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)106 grid_search_square_reg_2d(RemapGrid *src_grid, size_t (&src_add)[4], double (&src_lats)[4], double (&src_lons)[4], double plat,
107 double plon)
108 {
109 /*
110 Input variables:
111
112 plat : latitude of the search point
113 plon : longitude of the search point
114
115 Output variables:
116
117 src_add[4] : address of each corner point enclosing P
118 src_lats[4] : latitudes of the four corner points
119 src_lons[4] : longitudes of the four corner points
120 */
121 int search_result = 0;
122 const auto &src_center_lat = src_grid->reg2d_center_lat;
123 const auto &src_center_lon = src_grid->reg2d_center_lon;
124
125 for (unsigned n = 0; n < 4; ++n) src_add[n] = 0;
126
127 const auto nx = src_grid->dims[0];
128 const auto ny = src_grid->dims[1];
129
130 const auto nxm = src_grid->is_cyclic ? nx + 1 : nx;
131
132 if (plon < src_center_lon[0]) plon += PI2;
133 if (plon > src_center_lon[nxm - 1]) plon -= PI2;
134
135 size_t ii, jj;
136 int lfound = rect_grid_search(ii, jj, plon, plat, nxm, ny, src_center_lon, src_center_lat);
137 if (lfound)
138 {
139 auto iix = ii;
140 if (src_grid->is_cyclic && iix == (nxm - 1)) iix = 0;
141 src_add[0] = (jj - 1) * nx + (ii - 1);
142 src_add[1] = (jj - 1) * nx + (iix);
143 src_add[2] = (jj) *nx + (iix);
144 src_add[3] = (jj) *nx + (ii - 1);
145
146 src_lons[0] = src_center_lon[ii - 1];
147 src_lons[1] = src_center_lon[iix];
148 // For consistency, we must make sure all lons are in same 2pi interval
149 if (src_lons[0] > PI2) src_lons[0] -= PI2;
150 if (src_lons[0] < 0) src_lons[0] += PI2;
151 if (src_lons[1] > PI2) src_lons[1] -= PI2;
152 if (src_lons[1] < 0) src_lons[1] += PI2;
153 src_lons[2] = src_lons[1];
154 src_lons[3] = src_lons[0];
155
156 src_lats[0] = src_center_lat[jj - 1];
157 src_lats[1] = src_lats[0];
158 src_lats[2] = src_center_lat[jj];
159 src_lats[3] = src_lats[2];
160
161 search_result = 1;
162
163 return search_result;
164 }
165
166 /*
167 If no cell found, point is likely either in a box that straddles either pole
168 or is outside the grid. Fall back to a distance-weighted average of the four
169 closest points. Go ahead and compute weights here, but store in src_lats and
170 return -add to prevent the parent routine from computing bilinear weights.
171 */
172 if (!src_grid->lextrapolate) return search_result;
173
174 search_result = grid_search_square_reg_2d_NN(nx, ny, src_add, src_lats, plat, plon, src_center_lat, src_center_lon);
175
176 return search_result;
177 }
178