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