1 /*
2 * (C) Copyright 2005- ECMWF.
3 *
4 * This software is licensed under the terms of the Apache Licence Version 2.0
5 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6 *
7 * In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
8 * virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
9 */
10
11 /**
12 * Author: Enrico Fucile
13 * date: 31/10/2007
14 *
15 */
16
17 #include "grib_api_internal.h"
18
19 /* Note: The 'values' argument can be NULL in which case the data section will not be decoded
20 * See ECC-499
21 */
grib_nearest_find(grib_nearest * nearest,const grib_handle * ch,double inlat,double inlon,unsigned long flags,double * outlats,double * outlons,double * values,double * distances,int * indexes,size_t * len)22 int grib_nearest_find(
23 grib_nearest* nearest, const grib_handle* ch,
24 double inlat, double inlon,
25 unsigned long flags,
26 double* outlats, double* outlons,
27 double* values, double* distances, int* indexes, size_t* len)
28 {
29 grib_handle* h = (grib_handle*)ch;
30 grib_nearest_class* c = NULL;
31 if (!nearest)
32 return GRIB_INVALID_ARGUMENT;
33 c = nearest->cclass;
34 Assert(flags <= (GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA | GRIB_NEAREST_SAME_POINT));
35
36 while (c) {
37 grib_nearest_class* s = c->super ? *(c->super) : NULL;
38 if (c->find) {
39 int ret = c->find(nearest, h, inlat, inlon, flags, outlats, outlons, values, distances, indexes, len);
40 if (ret != GRIB_SUCCESS) {
41 if (inlon > 0)
42 inlon -= 360;
43 else
44 inlon += 360;
45 ret = c->find(nearest, h, inlat, inlon, flags, outlats, outlons, values, distances, indexes, len);
46 }
47 return ret;
48 }
49 c = s;
50 }
51 Assert(0);
52 return 0;
53 }
54
55 /* For this one, ALL init are called */
init_nearest(grib_nearest_class * c,grib_nearest * i,grib_handle * h,grib_arguments * args)56 static int init_nearest(grib_nearest_class* c, grib_nearest* i, grib_handle* h, grib_arguments* args)
57 {
58 if (c) {
59 int ret = GRIB_SUCCESS;
60 grib_nearest_class* s = c->super ? *(c->super) : NULL;
61 if (!c->inited) {
62 if (c->init_class)
63 c->init_class(c);
64 c->inited = 1;
65 }
66 if (s)
67 ret = init_nearest(s, i, h, args);
68
69 if (ret != GRIB_SUCCESS)
70 return ret;
71
72 if (c->init)
73 return c->init(i, h, args);
74 }
75 return GRIB_INTERNAL_ERROR;
76 }
77
grib_nearest_init(grib_nearest * i,grib_handle * h,grib_arguments * args)78 int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args)
79 {
80 return init_nearest(i->cclass, i, h, args);
81 }
82
83 /* For this one, ALL destroy are called */
84
grib_nearest_delete(grib_nearest * i)85 int grib_nearest_delete(grib_nearest* i)
86 {
87 grib_nearest_class* c = NULL;
88 if (!i)
89 return GRIB_INVALID_ARGUMENT;
90 c = i->cclass;
91 while (c) {
92 grib_nearest_class* s = c->super ? *(c->super) : NULL;
93 if (c->destroy)
94 c->destroy(i);
95 c = s;
96 }
97 return 0;
98 }
99
grib_binary_search(double xx[],const unsigned long n,double x,int * ju,int * jl)100 void grib_binary_search(double xx[], const unsigned long n, double x,
101 int* ju, int* jl)
102 {
103 size_t jm = 0;
104 int ascending = 0;
105 *jl = 0;
106 *ju = n;
107 ascending = (xx[n] >= xx[0]);
108 while (*ju - *jl > 1) {
109 jm = (*ju + *jl) >> 1;
110 if ((x >= xx[jm]) == ascending)
111 *jl = jm;
112 else
113 *ju = jm;
114 }
115 }
116
grib_nearest_find_multiple(const grib_handle * h,int is_lsm,const double * inlats,const double * inlons,long npoints,double * outlats,double * outlons,double * values,double * distances,int * indexes)117 int grib_nearest_find_multiple(
118 const grib_handle* h, int is_lsm,
119 const double* inlats, const double* inlons, long npoints,
120 double* outlats, double* outlons,
121 double* values, double* distances, int* indexes)
122 {
123 grib_nearest* nearest = 0;
124 double* pdistances = distances;
125 double* poutlats = outlats;
126 double* poutlons = outlons;
127 double* pvalues = values;
128 int* pindexes = indexes;
129 int idx = 0, ii = 0;
130 double max, min;
131 double qdistances[4] = {0,};
132 double qoutlats[4] = {0,};
133 double qoutlons[4] = {0,};
134 double qvalues[4] = {0,};
135 double* rvalues = NULL;
136 int qindexes[4] = {0,};
137 int ret = 0;
138 long i = 0;
139 size_t len = 4;
140 int flags = GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA;
141
142 if (values)
143 rvalues = qvalues;
144
145 nearest = grib_nearest_new(h, &ret);
146 if (ret != GRIB_SUCCESS)
147 return ret;
148
149 if (is_lsm) {
150 int noland = 1;
151 /* ECC-499: In land-sea mask mode, 'values' cannot be NULL because we need to query whether >= 0.5 */
152 Assert(values);
153 for (i = 0; i < npoints; i++) {
154 ret = grib_nearest_find(nearest, h, inlats[i], inlons[i], flags, qoutlats, qoutlons,
155 qvalues, qdistances, qindexes, &len);
156 max = qdistances[0];
157 for (ii = 0; ii < 4; ii++) {
158 if (max < qdistances[ii]) {
159 max = qdistances[ii];
160 idx = ii;
161 }
162 if (qvalues[ii] >= 0.5)
163 noland = 0;
164 }
165 min = max;
166 for (ii = 0; ii < 4; ii++) {
167 if ((min >= qdistances[ii]) && (noland || (qvalues[ii] >= 0.5))) {
168 min = qdistances[ii];
169 idx = ii;
170 }
171 }
172 *poutlats = qoutlats[idx];
173 poutlats++;
174 *poutlons = qoutlons[idx];
175 poutlons++;
176 *pvalues = qvalues[idx];
177 pvalues++;
178 *pdistances = qdistances[idx];
179 pdistances++;
180 *pindexes = qindexes[idx];
181 pindexes++;
182 }
183 }
184 else {
185 /* ECC-499: 'values' can be NULL */
186 for (i = 0; i < npoints; i++) {
187 ret = grib_nearest_find(nearest, h, inlats[i], inlons[i], flags, qoutlats, qoutlons,
188 rvalues, qdistances, qindexes, &len);
189 min = qdistances[0];
190 for (ii = 0; ii < 4; ii++) {
191 if ((min >= qdistances[ii])) {
192 min = qdistances[ii];
193 idx = ii;
194 }
195 }
196 *poutlats = qoutlats[idx];
197 poutlats++;
198 *poutlons = qoutlons[idx];
199 poutlons++;
200 if (values) {
201 *pvalues = qvalues[idx];
202 pvalues++;
203 }
204 *pdistances = qdistances[idx];
205 pdistances++;
206 *pindexes = qindexes[idx];
207 pindexes++;
208 }
209 }
210
211 grib_nearest_delete(nearest);
212
213 return ret;
214 }
215
216
217 /* Generic implementation of nearest for Lambert, Polar stereo, Mercator etc */
compare_doubles(const void * a,const void * b,int ascending)218 static int compare_doubles(const void* a, const void* b, int ascending)
219 {
220 /* ascending is a boolean: 0 or 1 */
221 double* arg1 = (double*)a;
222 double* arg2 = (double*)b;
223 if (ascending) {
224 if (*arg1 < *arg2)
225 return -1; /*Smaller values come before larger ones*/
226 }
227 else {
228 if (*arg1 > *arg2)
229 return -1; /*Larger values come before smaller ones*/
230 }
231 if (*arg1 == *arg2)
232 return 0;
233 else
234 return 1;
235 }
236
compare_doubles_ascending(const void * a,const void * b)237 static int compare_doubles_ascending(const void* a, const void* b)
238 {
239 return compare_doubles(a, b, 1);
240 }
241
242 typedef struct PointStore
243 {
244 double m_lat;
245 double m_lon;
246 double m_dist;
247 double m_value;
248 int m_index;
249 } PointStore;
250
251 /* Comparison function to sort points by distance */
compare_points(const void * a,const void * b)252 static int compare_points(const void* a, const void* b)
253 {
254 PointStore* pA = (PointStore*)a;
255 PointStore* pB = (PointStore*)b;
256
257 if (pA->m_dist < pB->m_dist) return -1;
258 if (pA->m_dist > pB->m_dist) return 1;
259 return 0;
260 }
261
grib_nearest_find_generic(grib_nearest * nearest,grib_handle * h,double inlat,double inlon,unsigned long flags,const char * values_keyname,const char * radius_keyname,const char * Ni_keyname,const char * Nj_keyname,double ** out_lats,int * out_lats_count,double ** out_lons,int * out_lons_count,double ** out_distances,double * outlats,double * outlons,double * values,double * distances,int * indexes,size_t * len)262 int grib_nearest_find_generic(
263 grib_nearest* nearest, grib_handle* h,
264 double inlat, double inlon, unsigned long flags,
265
266 const char* values_keyname,
267 const char* radius_keyname,
268 const char* Ni_keyname,
269 const char* Nj_keyname,
270 double** out_lats,
271 int* out_lats_count,
272 double** out_lons,
273 int* out_lons_count,
274 double** out_distances,
275
276 double* outlats, double* outlons,
277 double* values, double* distances, int* indexes, size_t* len)
278 {
279 int ret = 0, i = 0;
280 size_t nvalues = 0, nneighbours = 0;
281 double radiusInMetres, radiusInKm;
282 grib_iterator* iter = NULL;
283 double lat = 0, lon = 0;
284
285 /* array of candidates for nearest neighbours */
286 PointStore* neighbours = NULL;
287
288 inlon = normalise_longitude_in_degrees(inlon);
289
290 if ((ret = grib_get_size(h, values_keyname, &nvalues)) != GRIB_SUCCESS)
291 return ret;
292 nearest->values_count = nvalues;
293
294 /* We need the radius to calculate the nearest distance. For an oblate earth
295 approximate this using the average of the semimajor and semiminor axes */
296 if ((ret = grib_get_double(h, radius_keyname, &radiusInMetres)) == GRIB_SUCCESS &&
297 !grib_is_missing(h, radius_keyname, &ret)) {
298 radiusInKm = radiusInMetres / 1000.0;
299 }
300 else {
301 double minor = 0, major = 0;
302 if ((ret = grib_get_double_internal(h, "earthMinorAxisInMetres", &minor)) != GRIB_SUCCESS) return ret;
303 if ((ret = grib_get_double_internal(h, "earthMajorAxisInMetres", &major)) != GRIB_SUCCESS) return ret;
304 if (grib_is_missing(h, "earthMinorAxisInMetres", &ret)) return GRIB_GEOCALCULUS_PROBLEM;
305 if (grib_is_missing(h, "earthMajorAxisInMetres", &ret)) return GRIB_GEOCALCULUS_PROBLEM;
306 radiusInMetres = (major + minor) / 2;
307 radiusInKm = radiusInMetres / 1000.0;
308 }
309
310 neighbours = (PointStore*)grib_context_malloc(nearest->context, nvalues * sizeof(PointStore));
311 for (i = 0; i < nvalues; ++i) {
312 neighbours[i].m_dist = 1e10; /* set all distances to large number to begin with */
313 neighbours[i].m_lat = 0;
314 neighbours[i].m_lon = 0;
315 neighbours[i].m_value = 0;
316 neighbours[i].m_index = 0;
317 }
318
319 /* GRIB_NEAREST_SAME_GRID not yet implemented */
320 {
321 double the_value = 0;
322 double min_dist = 1e10;
323 size_t the_index = 0;
324 int ilat = 0, ilon = 0;
325 int idx_upper = 0, idx_lower = 0;
326 double lat1 = 0, lat2 = 0; /* inlat will be between these */
327 const double LAT_DELTA = 10.0; /* in degrees */
328
329 if (grib_is_missing(h, Ni_keyname, &ret)) {
330 grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Ni_keyname);
331 return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
332 }
333
334 if (grib_is_missing(h, Nj_keyname, &ret)) {
335 grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Nj_keyname);
336 return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
337 }
338
339 *out_lons_count = nvalues; /* Maybe overestimate but safe */
340 *out_lats_count = nvalues;
341
342 if (*out_lats)
343 grib_context_free(nearest->context, *out_lats);
344 *out_lats = (double*)grib_context_malloc(nearest->context, nvalues * sizeof(double));
345 if (!*out_lats)
346 return GRIB_OUT_OF_MEMORY;
347
348 if (*out_lons)
349 grib_context_free(nearest->context, *out_lons);
350 *out_lons = (double*)grib_context_malloc(nearest->context, nvalues * sizeof(double));
351 if (!*out_lons)
352 return GRIB_OUT_OF_MEMORY;
353
354 iter = grib_iterator_new(h, 0, &ret);
355 if (ret)
356 return ret;
357 /* First pass: collect all latitudes and longitudes */
358 while (grib_iterator_next(iter, &lat, &lon, &the_value)) {
359 ++the_index;
360 Assert(ilat < *out_lats_count);
361 Assert(ilon < *out_lons_count);
362 (*out_lats)[ilat++] = lat;
363 (*out_lons)[ilon++] = lon;
364 }
365
366 /* See between which 2 latitudes our point lies */
367 qsort(*out_lats, nvalues, sizeof(double), &compare_doubles_ascending);
368 grib_binary_search(*out_lats, *out_lats_count - 1, inlat, &idx_upper, &idx_lower);
369 lat2 = (*out_lats)[idx_upper];
370 lat1 = (*out_lats)[idx_lower];
371 Assert(lat1 <= lat2);
372
373 /* Second pass: Iterate again and collect candidate neighbours */
374 grib_iterator_reset(iter);
375 the_index = 0;
376 i = 0;
377 while (grib_iterator_next(iter, &lat, &lon, &the_value)) {
378 if (lat > lat2 + LAT_DELTA || lat < lat1 - LAT_DELTA) {
379 /* Ignore latitudes too far from our point */
380 }
381 else {
382 double dist = geographic_distance_spherical(radiusInKm, inlon, inlat, lon, lat);
383 if (dist < min_dist)
384 min_dist = dist;
385 /*printf("Candidate: lat=%.5f lon=%.5f dist=%f Idx=%ld Val=%f\n",lat,lon,dist,the_index,the_value);*/
386 /* store this candidate point */
387 neighbours[i].m_dist = dist;
388 neighbours[i].m_index = the_index;
389 neighbours[i].m_lat = lat;
390 neighbours[i].m_lon = lon;
391 neighbours[i].m_value = the_value;
392 i++;
393 }
394 ++the_index;
395 }
396 nneighbours = i;
397 /* Sort the candidate neighbours in ascending order of distance */
398 /* The first 4 entries will now be the closest 4 neighbours */
399 qsort(neighbours, nneighbours, sizeof(PointStore), &compare_points);
400
401 grib_iterator_delete(iter);
402 }
403 nearest->h = h;
404
405 /* Sanity check for sorting */
406 #ifdef DEBUG
407 for (i = 0; i < nneighbours - 1; ++i) {
408 Assert(neighbours[i].m_dist <= neighbours[i + 1].m_dist);
409 }
410 #endif
411
412 /* GRIB_NEAREST_SAME_XXX not yet implemented */
413 if (!*out_distances) {
414 *out_distances = (double*)grib_context_malloc(nearest->context, 4 * sizeof(double));
415 }
416 (*out_distances)[0] = neighbours[0].m_dist;
417 (*out_distances)[1] = neighbours[1].m_dist;
418 (*out_distances)[2] = neighbours[2].m_dist;
419 (*out_distances)[3] = neighbours[3].m_dist;
420
421 for (i = 0; i < 4; ++i) {
422 distances[i] = neighbours[i].m_dist;
423 outlats[i] = neighbours[i].m_lat;
424 outlons[i] = neighbours[i].m_lon;
425 indexes[i] = neighbours[i].m_index;
426 values[i] = neighbours[i].m_value;
427 /*printf("(%f,%f) i=%d d=%f v=%f\n",outlats[i],outlons[i],indexes[i],distances[i],values[i]);*/
428 }
429
430 free(neighbours);
431 return GRIB_SUCCESS;
432 }
433