1 /*
2  * Copyright 2005-2018 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 #include "grib_api_internal.h"
12 
13 /*
14    This is used by make_class.pl
15 
16    START_CLASS_DEF
17    CLASS      = nearest
18    SUPER      = grib_nearest_class_gen
19    IMPLEMENTS = init
20    IMPLEMENTS = find;destroy
21    MEMBERS    = double* lats
22    MEMBERS    = int  lats_count
23    MEMBERS    = double* lons
24    MEMBERS    = int  lons_count
25    MEMBERS    = double* distances
26    MEMBERS    = int* k
27    MEMBERS    = int* i
28    MEMBERS    = int* j
29    MEMBERS    = const char* Ni
30    MEMBERS    = const char* Nj
31    END_CLASS_DEF
32 
33  */
34 
35 /* START_CLASS_IMP */
36 
37 /*
38 
39 Don't edit anything between START_CLASS_IMP and END_CLASS_IMP
40 Instead edit values between START_CLASS_DEF and END_CLASS_DEF
41 or edit "nearest.class" and rerun ./make_class.pl
42 
43 */
44 
45 
46 static void init_class              (grib_nearest_class*);
47 
48 static int init               (grib_nearest* nearest,grib_handle* h,grib_arguments* args);
49 static int find(grib_nearest* nearest, grib_handle* h,double inlat, double inlon, unsigned long flags, double* outlats,double* outlons, double *values,double *distances, int *indexes,size_t *len);
50 static int destroy            (grib_nearest* nearest);
51 
52 typedef struct grib_nearest_regular{
53   grib_nearest nearest;
54 /* Members defined in gen */
55 	const char* values_key;
56 	const char* radius;
57 	int cargs;
58 /* Members defined in regular */
59 	double* lats;
60 	int  lats_count;
61 	double* lons;
62 	int  lons_count;
63 	double* distances;
64 	int* k;
65 	int* i;
66 	int* j;
67 	const char* Ni;
68 	const char* Nj;
69 } grib_nearest_regular;
70 
71 extern grib_nearest_class* grib_nearest_class_gen;
72 
73 static grib_nearest_class _grib_nearest_class_regular = {
74     &grib_nearest_class_gen,                         /* super                     */
75     "regular",                         /* name                      */
76     sizeof(grib_nearest_regular),      /* size of instance          */
77     0,                              /* inited */
78     &init_class,                    /* init_class */
79     &init,                          /* constructor               */
80     &destroy,                       /* destructor                */
81     &find,                          /* find nearest              */
82 };
83 
84 grib_nearest_class* grib_nearest_class_regular = &_grib_nearest_class_regular;
85 
86 
init_class(grib_nearest_class * c)87 static void init_class(grib_nearest_class* c)
88 {
89 }
90 /* END_CLASS_IMP */
91 
92 
init(grib_nearest * nearest,grib_handle * h,grib_arguments * args)93 static int init(grib_nearest* nearest,grib_handle* h,grib_arguments* args)
94 {
95     grib_nearest_regular* self = (grib_nearest_regular*) nearest;
96     self->Ni  = grib_arguments_get_name(h,args,self->cargs++);
97     self->Nj  = grib_arguments_get_name(h,args,self->cargs++);
98     self->i=(int*)grib_context_malloc(h->context,2*sizeof(int));
99     self->j=(int*)grib_context_malloc(h->context,2*sizeof(int));
100     return 0;
101 }
102 
103 #if 0
104 static int find(grib_nearest* nearest, grib_handle* h,
105         double inlat, double inlon,unsigned long flags,
106         double* outlats,double* outlons,
107         double *values,double *distances,int* indexes, size_t *len) {
108     grib_nearest_regular* self = (grib_nearest_regular*) nearest;
109     int ret=0,kk=0,ii=0,jj=0;
110     size_t nvalues=0;
111 
112     long iradius;
113     double radius;
114 
115     if( (ret =  grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS)
116         return ret;
117     radius=((double)iradius)/1000.0;
118 
119     if (!nearest->h || (flags & GRIB_NEAREST_SAME_DATA)==0 || nearest->h!=h) {
120         grib_iterator* iter=NULL;
121         double lat=0,lon=0;
122 
123         if( (ret =  grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS)
124             return ret;
125         nearest->values_count = nvalues;
126         if (nearest->values) grib_context_free(nearest->context,nearest->values);
127         nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double));
128         if (!nearest->values) return GRIB_OUT_OF_MEMORY;
129 
130         ret=grib_get_double_array_internal( h,self->values_key,
131                 nearest->values,&(nearest->values_count));
132         if (ret!=GRIB_SUCCESS) grib_context_log(nearest->context,GRIB_LOG_ERROR,
133                 "nearest: unable to get values array");
134 
135         if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID)==0) {
136             double dummy=0;
137             double olat=1.e10, olon=1.e10;
138             int ilat=0,ilon=0;
139             long n=0;
140 
141             if( (ret =  grib_get_long(h,self->Ni,&n))!= GRIB_SUCCESS)
142                 return ret;
143             self->lons_count=n;
144 
145             if( (ret =  grib_get_long(h,self->Nj,&n))!= GRIB_SUCCESS)
146                 return ret;
147             self->lats_count=n;
148 
149             if (self->lats) grib_context_free(nearest->context,self->lats);
150             self->lats=grib_context_malloc( nearest->context,
151                     self->lats_count* sizeof(double));
152             if (!self->lats) return GRIB_OUT_OF_MEMORY;
153 
154             if (self->lons) grib_context_free(nearest->context,self->lons);
155             self->lons=grib_context_malloc( nearest->context,
156                     self->lons_count*sizeof(double));
157             if (!self->lons) return GRIB_OUT_OF_MEMORY;
158 
159             iter=grib_iterator_new(h,0,&ret);
160             if (ret) {
161                 grib_context_log(nearest->context,GRIB_LOG_ERROR,"unable to create iterator");
162                 return ret;
163             }
164             while(grib_iterator_next(iter,&lat,&lon,&dummy)) {
165                 if (olat != lat) {
166                     self->lats[ilat++]=lat;
167                     olat=lat;
168                 }
169                 if (ilon<self->lons_count && olon != lon) {
170                     self->lons[ilon++]=lon;
171                     olon=lon;
172                 }
173             }
174             grib_iterator_delete(iter);
175 
176         }
177         nearest->h=h;
178 
179     }
180 
181     if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT)==0
182             || (flags & GRIB_NEAREST_SAME_GRID)==0) {
183 
184         grib_binary_search(self->lats,self->lats_count-1,inlat,
185                 &(self->j[0]),&(self->j[1]));
186         grib_binary_search(self->lons,self->lons_count-1,inlon,
187                 &(self->i[0]),&(self->i[1]));
188         if (!self->distances)
189             self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double));
190         if (!self->k)
191             self->k=(int*)grib_context_malloc( nearest->context,4*sizeof(int));
192         kk=0;
193         for (ii=0;ii<2;ii++) {
194             for (jj=0;jj<2;jj++) {
195                 self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]-1;
196                 self->distances[kk]=grib_nearest_distance(radius,inlon,inlat,
197                         self->lons[self->i[ii]],self->lats[self->j[jj]]);
198                 kk++;
199             }
200         }
201     }
202 
203     kk=0;
204     for (ii=0;ii<2;ii++) {
205         for (jj=0;jj<2;jj++) {
206             distances[kk]=self->distances[kk];
207             outlats[kk]=self->lats[self->j[jj]];
208             outlons[kk]=self->lons[self->i[ii]];
209             values[kk]=nearest->values[self->k[kk]];
210             indexes[kk]=self->k[kk];
211             kk++;
212         }
213     }
214 
215     return GRIB_SUCCESS;
216 }
217 #else
is_rotated_grid(grib_handle * h)218 static int is_rotated_grid(grib_handle* h)
219 {
220     long is_rotated = 0;
221     int err = grib_get_long(h, "is_rotated_grid", &is_rotated);
222     if (!err && is_rotated ) return 1;
223     return 0;
224 }
225 
find(grib_nearest * nearest,grib_handle * h,double inlat,double inlon,unsigned long flags,double * outlats,double * outlons,double * values,double * distances,int * indexes,size_t * len)226 static int find(grib_nearest* nearest, grib_handle* h,
227         double inlat, double inlon,unsigned long flags,
228         double* outlats,double* outlons,
229         double *values,double *distances,int* indexes, size_t *len)
230 {
231     grib_nearest_regular* self = (grib_nearest_regular*) nearest;
232     int ret=0,kk=0,ii=0,jj=0;
233     size_t nvalues=0;
234     long iradius;
235     double radius;
236 
237     grib_iterator* iter=NULL;
238     double lat=0,lon=0;
239 
240     while (inlon<0) inlon+=360;
241     while (inlon>360) inlon-=360;
242 
243     if( (ret =  grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS)
244         return ret;
245     nearest->values_count = nvalues;
246 
247     if (grib_is_missing(h,self->radius,&ret)) {
248         grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->radius);
249         return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
250     }
251 
252     if( (ret =  grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS)
253         return ret;
254     radius=((double)iradius)/1000.0;
255 
256     if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID)==0) {
257         double dummy=0;
258         double olat=1.e10, olon=1.e10;
259         int ilat=0,ilon=0;
260         long n=0;
261 
262         if (grib_is_missing(h,self->Ni,&ret)) {
263             grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Ni);
264             return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
265         }
266 
267         if (grib_is_missing(h,self->Nj,&ret)) {
268             grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Nj);
269             return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
270         }
271 
272         /* Support for rotated grids not yet implemented */
273         if (is_rotated_grid(h)) {
274             grib_context_log(h->context,GRIB_LOG_ERROR,
275                     "Nearest neighbour functionality is not supported for rotated grids.");
276             return GRIB_NOT_IMPLEMENTED;
277         }
278 
279         if ((ret =  grib_get_long(h,self->Ni,&n))!= GRIB_SUCCESS)
280             return ret;
281         self->lons_count=n;
282 
283         if ((ret =  grib_get_long(h,self->Nj,&n))!= GRIB_SUCCESS)
284             return ret;
285         self->lats_count=n;
286 
287         if (self->lats) grib_context_free(nearest->context,self->lats);
288         self->lats=(double*)grib_context_malloc( nearest->context,
289                 self->lats_count* sizeof(double));
290         if (!self->lats) return GRIB_OUT_OF_MEMORY;
291 
292         if (self->lons) grib_context_free(nearest->context,self->lons);
293         self->lons=(double*)grib_context_malloc( nearest->context,
294                 self->lons_count*sizeof(double));
295         if (!self->lons) return GRIB_OUT_OF_MEMORY;
296 
297         iter=grib_iterator_new(h,0,&ret);
298         if (ret != GRIB_SUCCESS) {
299             grib_context_log(h->context,GRIB_LOG_ERROR,"Unable to create lat/lon iterator");
300             return ret;
301         }
302         while(grib_iterator_next(iter,&lat,&lon,&dummy)) {
303             if (olat != lat) {
304                 Assert( ilat < self->lats_count );
305                 self->lats[ilat++]=lat;
306                 olat=lat;
307             }
308             if (ilon<self->lons_count && olon != lon) {
309                 self->lons[ilon++]=lon ;
310                 olon=lon;
311             }
312         }
313         grib_iterator_delete(iter);
314 
315     }
316     nearest->h=h;
317 
318     if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT)==0
319             || (flags & GRIB_NEAREST_SAME_GRID)==0) {
320         int nearest_lons_found=0;
321 
322         if (self->lats[self->lats_count-1] > self->lats[0]) {
323             if (inlat<self->lats[0] || inlat>self->lats[self->lats_count-1])
324                 return GRIB_OUT_OF_AREA;
325         } else {
326             if (inlat > self->lats[0] || inlat < self->lats[self->lats_count-1])
327                 return GRIB_OUT_OF_AREA;
328         }
329 
330         if (self->lons[self->lons_count-1] > self->lons[0]) {
331             if (inlon<self->lons[0] || inlon>self->lons[self->lons_count-1]) {
332                 /* try to scale*/
333                 if (inlon>0) inlon-=360;
334                 else inlon+=360;
335 
336                 if (inlon<self->lons[0] || inlon>self->lons[self->lons_count-1]) {
337                     if (	self->lons[0]+360-self->lons[self->lons_count-1]<=
338                             self->lons[1]-self->lons[0]) {
339                         /*it's a global field in longitude*/
340                         self->i[0]=0;
341                         self->i[1]=self->lons_count-1;
342                         nearest_lons_found=1;
343                     } else
344                         return GRIB_OUT_OF_AREA;
345                 }
346             }
347         } else {
348             if (inlon>self->lons[0] || inlon<self->lons[self->lons_count-1]) {
349                 /* try to scale*/
350                 if (inlon>0) inlon-=360;
351                 else inlon+=360;
352                 if (self->lons[0]-self->lons[self->lons_count-1]-360 <=
353                         self->lons[0]-self->lons[1]) {
354                     /*it's a global field in longitude*/
355                     self->i[0]=0;
356                     self->i[1]=self->lons_count-1;
357                     nearest_lons_found=1;
358                 } else if (inlon>self->lons[0] || inlon<self->lons[self->lons_count-1])
359                     return GRIB_OUT_OF_AREA;
360             }
361         }
362 
363         grib_binary_search(self->lats,self->lats_count-1,inlat,
364                 &(self->j[0]),&(self->j[1]));
365 
366         if (!nearest_lons_found)
367             grib_binary_search(self->lons,self->lons_count-1,inlon,
368                     &(self->i[0]),&(self->i[1]));
369 
370         if (!self->distances)
371             self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double));
372         if (!self->k)
373             self->k=(int*)grib_context_malloc( nearest->context,4*sizeof(int));
374         kk=0;
375         for (jj=0;jj<2;jj++) {
376             for (ii=0;ii<2;ii++) {
377                 self->k[kk]=self->i[ii]+self->lons_count*self->j[jj];
378                 self->distances[kk]=grib_nearest_distance(radius,inlon,inlat,
379                         self->lons[self->i[ii]],self->lats[self->j[jj]]);
380                 kk++;
381             }
382         }
383     }
384 
385     kk=0;
386 
387     /*
388      * Brute force algorithm:
389      * First unpack all the values into an array. Then when we need the 4 points
390      * we just index into this array so no need to call grib_get_double_element_internal
391      *
392      *   if (nearest->values) grib_context_free(nearest->context,nearest->values);
393      *   nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double));
394      *   if (!nearest->values) return GRIB_OUT_OF_MEMORY;
395      *   ret = grib_get_double_array(h, self->values_key, nearest->values ,&nvalues);
396      *   if (ret) return ret;
397      */
398 
399     for (jj=0;jj<2;jj++) {
400         for (ii=0;ii<2;ii++) {
401             distances[kk]=self->distances[kk];
402             outlats[kk]=self->lats[self->j[jj]];
403             outlons[kk]=self->lons[self->i[ii]];
404             grib_get_double_element_internal(h,self->values_key,self->k[kk],&(values[kk]));
405             /* Using the brute force approach described above */
406             /* Assert(self->k[kk] < nvalues); */
407             /* values[kk]=nearest->values[self->k[kk]]; */
408 
409             indexes[kk]=self->k[kk];
410             kk++;
411         }
412     }
413 
414     return GRIB_SUCCESS;
415 }
416 #endif
417 
destroy(grib_nearest * nearest)418 static int destroy(grib_nearest* nearest) {
419     grib_nearest_regular* self = (grib_nearest_regular*) nearest;
420     if (self->lats) grib_context_free(nearest->context,self->lats);
421     if (self->lons) grib_context_free(nearest->context,self->lons);
422     if (self->i) grib_context_free(nearest->context,self->i);
423     if (self->j) grib_context_free(nearest->context,self->j);
424     if (self->k) grib_context_free(nearest->context,self->k);
425     if (self->distances) grib_context_free(nearest->context,self->distances);
426     return GRIB_SUCCESS;
427 }
428