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 #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 {
54     grib_nearest nearest;
55     /* Members defined in gen */
56     const char* values_key;
57     const char* radius;
58     int cargs;
59     /* Members defined in regular */
60     double* lats;
61     int lats_count;
62     double* lons;
63     int lons_count;
64     double* distances;
65     int* k;
66     int* i;
67     int* j;
68     const char* Ni;
69     const char* Nj;
70 } grib_nearest_regular;
71 
72 extern grib_nearest_class* grib_nearest_class_gen;
73 
74 static grib_nearest_class _grib_nearest_class_regular = {
75     &grib_nearest_class_gen,      /* super                     */
76     "regular",                    /* name                      */
77     sizeof(grib_nearest_regular), /* size of instance          */
78     0,                            /* inited */
79     &init_class,                  /* init_class */
80     &init,                        /* constructor               */
81     &destroy,                     /* destructor                */
82     &find,                        /* find nearest              */
83 };
84 
85 grib_nearest_class* grib_nearest_class_regular = &_grib_nearest_class_regular;
86 
87 
init_class(grib_nearest_class * c)88 static void init_class(grib_nearest_class* c)
89 {
90 }
91 /* END_CLASS_IMP */
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     if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT)==0
181             || (flags & GRIB_NEAREST_SAME_GRID)==0) {
182 
183         grib_binary_search(self->lats,self->lats_count-1,inlat,
184                 &(self->j[0]),&(self->j[1]));
185         grib_binary_search(self->lons,self->lons_count-1,inlon,
186                 &(self->i[0]),&(self->i[1]));
187         if (!self->distances)
188             self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double));
189         if (!self->k)
190             self->k=(int*)grib_context_malloc( nearest->context,4*sizeof(int));
191         kk=0;
192         for (ii=0;ii<2;ii++) {
193             for (jj=0;jj<2;jj++) {
194                 self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]-1;
195                 self->distances[kk]=geographic_distance_spherical(radius,inlon,inlat,
196                         self->lons[self->i[ii]],self->lats[self->j[jj]]);
197                 kk++;
198             }
199         }
200     }
201 
202     kk=0;
203     for (ii=0;ii<2;ii++) {
204         for (jj=0;jj<2;jj++) {
205             distances[kk]=self->distances[kk];
206             outlats[kk]=self->lats[self->j[jj]];
207             outlons[kk]=self->lons[self->i[ii]];
208             values[kk]=nearest->values[self->k[kk]];
209             indexes[kk]=self->k[kk];
210             kk++;
211         }
212     }
213 
214     return GRIB_SUCCESS;
215 }
216 #else
is_rotated_grid(grib_handle * h)217 static int is_rotated_grid(grib_handle* h)
218 {
219     long is_rotated = 0;
220     int err         = grib_get_long(h, "isRotatedGrid", &is_rotated);
221     if (!err && is_rotated)
222         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     const int is_rotated   = is_rotated_grid(h);
240     double angleOfRotation = 0, southPoleLat = 0, southPoleLon = 0;
241 
242     while (inlon < 0)
243         inlon += 360;
244     while (inlon > 360)
245         inlon -= 360;
246 
247     if ((ret = grib_get_size(h, self->values_key, &nvalues)) != GRIB_SUCCESS)
248         return ret;
249     nearest->values_count = nvalues;
250 
251     if (grib_is_missing(h, self->radius, &ret)) {
252         grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->radius);
253         return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
254     }
255 
256     if ((ret = grib_get_long(h, self->radius, &iradius)) != GRIB_SUCCESS)
257         return ret;
258     radius = ((double)iradius) / 1000.0;
259 
260     /* Compute lat/lon info, create iterator etc if it's the 1st time or different grid.
261      * This is for performance: if the grid has not changed, we only do this once
262      * and reuse for other messages */
263     if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID) == 0) {
264         double dummy = 0;
265         double olat = 1.e10, olon = 1.e10;
266         int ilat = 0, ilon = 0;
267         long n = 0;
268 
269         if (grib_is_missing(h, self->Ni, &ret)) {
270             grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->Ni);
271             return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
272         }
273 
274         if (grib_is_missing(h, self->Nj, &ret)) {
275             grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->Nj);
276             return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
277         }
278 
279         /* ECC-600: Support for rotated grids
280          * First:   rotate the input point
281          * Then:    run the lat/lon iterator over the rotated grid (disableUnrotate)
282          * Finally: unrotate the resulting point
283          */
284         if (is_rotated) {
285             double new_lat = 0, new_lon = 0;
286             ret = grib_get_double_internal(h, "angleOfRotation", &angleOfRotation);
287             if (ret)
288                 return ret;
289             ret = grib_get_double_internal(h, "latitudeOfSouthernPoleInDegrees", &southPoleLat);
290             if (ret)
291                 return ret;
292             ret = grib_get_double_internal(h, "longitudeOfSouthernPoleInDegrees", &southPoleLon);
293             if (ret)
294                 return ret;
295             ret = grib_set_long(h, "iteratorDisableUnrotate", 1);
296             if (ret)
297                 return ret;
298             /* Rotate the inlat, inlon */
299             rotate(inlat, inlon, angleOfRotation, southPoleLat, southPoleLon, &new_lat, &new_lon);
300             inlat = new_lat;
301             inlon = new_lon;
302             /*if(h->context->debug) printf("nearest find: rotated grid: new point=(%g,%g)\n",new_lat,new_lon);*/
303         }
304 
305         if ((ret = grib_get_long(h, self->Ni, &n)) != GRIB_SUCCESS)
306             return ret;
307         self->lons_count = n;
308 
309         if ((ret = grib_get_long(h, self->Nj, &n)) != GRIB_SUCCESS)
310             return ret;
311         self->lats_count = n;
312 
313         if (self->lats)
314             grib_context_free(nearest->context, self->lats);
315         self->lats = (double*)grib_context_malloc(nearest->context,
316                                                   self->lats_count * sizeof(double));
317         if (!self->lats)
318             return GRIB_OUT_OF_MEMORY;
319 
320         if (self->lons)
321             grib_context_free(nearest->context, self->lons);
322         self->lons = (double*)grib_context_malloc(nearest->context,
323                                                   self->lons_count * sizeof(double));
324         if (!self->lons)
325             return GRIB_OUT_OF_MEMORY;
326 
327         iter = grib_iterator_new(h, 0, &ret);
328         if (ret != GRIB_SUCCESS) {
329             grib_context_log(h->context, GRIB_LOG_ERROR, "Unable to create lat/lon iterator");
330             return ret;
331         }
332         while (grib_iterator_next(iter, &lat, &lon, &dummy)) {
333             if (olat != lat) {
334                 Assert(ilat < self->lats_count);
335                 self->lats[ilat++] = lat;
336                 olat               = lat;
337             }
338             if (ilon < self->lons_count && olon != lon) {
339                 self->lons[ilon++] = lon;
340                 olon               = lon;
341             }
342         }
343         grib_iterator_delete(iter);
344     }
345     nearest->h = h;
346 
347     /* Compute distances if it's the 1st time or different point or different grid.
348      * This is for performance: if the grid and the input point have not changed
349      * we only do this once and reuse for other messages */
350     if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT) == 0 || (flags & GRIB_NEAREST_SAME_GRID) == 0) {
351         int nearest_lons_found = 0;
352 
353         if (self->lats[self->lats_count - 1] > self->lats[0]) {
354             if (inlat < self->lats[0] || inlat > self->lats[self->lats_count - 1])
355                 return GRIB_OUT_OF_AREA;
356         }
357         else {
358             if (inlat > self->lats[0] || inlat < self->lats[self->lats_count - 1])
359                 return GRIB_OUT_OF_AREA;
360         }
361 
362         if (self->lons[self->lons_count - 1] > self->lons[0]) {
363             if (inlon < self->lons[0] || inlon > self->lons[self->lons_count - 1]) {
364                 /* try to scale*/
365                 if (inlon > 0)
366                     inlon -= 360;
367                 else
368                     inlon += 360;
369 
370                 if (inlon < self->lons[0] || inlon > self->lons[self->lons_count - 1]) {
371                     if (self->lons[0] + 360 - self->lons[self->lons_count - 1] <=
372                         self->lons[1] - self->lons[0]) {
373                         /*it's a global field in longitude*/
374                         self->i[0]         = 0;
375                         self->i[1]         = self->lons_count - 1;
376                         nearest_lons_found = 1;
377                     }
378                     else
379                         return GRIB_OUT_OF_AREA;
380                 }
381             }
382         }
383         else {
384             if (inlon > self->lons[0] || inlon < self->lons[self->lons_count - 1]) {
385                 /* try to scale*/
386                 if (inlon > 0)
387                     inlon -= 360;
388                 else
389                     inlon += 360;
390                 if (self->lons[0] - self->lons[self->lons_count - 1] - 360 <=
391                     self->lons[0] - self->lons[1]) {
392                     /*it's a global field in longitude*/
393                     self->i[0]         = 0;
394                     self->i[1]         = self->lons_count - 1;
395                     nearest_lons_found = 1;
396                 }
397                 else if (inlon > self->lons[0] || inlon < self->lons[self->lons_count - 1])
398                     return GRIB_OUT_OF_AREA;
399             }
400         }
401 
402         grib_binary_search(self->lats, self->lats_count - 1, inlat,
403                            &(self->j[0]), &(self->j[1]));
404 
405         if (!nearest_lons_found)
406             grib_binary_search(self->lons, self->lons_count - 1, inlon,
407                                &(self->i[0]), &(self->i[1]));
408 
409         if (!self->distances)
410             self->distances = (double*)grib_context_malloc(nearest->context, 4 * sizeof(double));
411         if (!self->k)
412             self->k = (int*)grib_context_malloc(nearest->context, 4 * sizeof(int));
413         kk = 0;
414         for (jj = 0; jj < 2; jj++) {
415             for (ii = 0; ii < 2; ii++) {
416                 self->k[kk]         = self->i[ii] + self->lons_count * self->j[jj];
417                 self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
418                                                             self->lons[self->i[ii]], self->lats[self->j[jj]]);
419                 kk++;
420             }
421         }
422     }
423 
424     kk = 0;
425 
426     /*
427      * Brute force algorithm:
428      * First unpack all the values into an array. Then when we need the 4 points
429      * we just index into this array so no need to call grib_get_double_element_internal
430      *
431      *   if (nearest->values) grib_context_free(nearest->context,nearest->values);
432      *   nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double));
433      *   if (!nearest->values) return GRIB_OUT_OF_MEMORY;
434      *   ret = grib_get_double_array(h, self->values_key, nearest->values ,&nvalues);
435      *   if (ret) return ret;
436      */
437 
438     for (jj = 0; jj < 2; jj++) {
439         for (ii = 0; ii < 2; ii++) {
440             distances[kk] = self->distances[kk];
441             outlats[kk]   = self->lats[self->j[jj]];
442             outlons[kk]   = self->lons[self->i[ii]];
443             if (is_rotated) {
444                 /* Unrotate resulting lat/lon */
445                 double new_lat = 0, new_lon = 0;
446                 unrotate(outlats[kk], outlons[kk], angleOfRotation, southPoleLat, southPoleLon, &new_lat, &new_lon);
447                 outlats[kk] = new_lat;
448                 outlons[kk] = new_lon;
449             }
450             if (values) { /* ECC-499 */
451                 grib_get_double_element_internal(h, self->values_key, self->k[kk], &(values[kk]));
452             }
453             /* Using the brute force approach described above */
454             /* Assert(self->k[kk] < nvalues); */
455             /* values[kk]=nearest->values[self->k[kk]]; */
456 
457             indexes[kk] = self->k[kk];
458             kk++;
459         }
460     }
461 
462     return GRIB_SUCCESS;
463 }
464 #endif
465 
destroy(grib_nearest * nearest)466 static int destroy(grib_nearest* nearest)
467 {
468     grib_nearest_regular* self = (grib_nearest_regular*)nearest;
469     if (self->lats)
470         grib_context_free(nearest->context, self->lats);
471     if (self->lons)
472         grib_context_free(nearest->context, self->lons);
473     if (self->i)
474         grib_context_free(nearest->context, self->i);
475     if (self->j)
476         grib_context_free(nearest->context, self->j);
477     if (self->k)
478         grib_context_free(nearest->context, self->k);
479     if (self->distances)
480         grib_context_free(nearest->context, self->distances);
481     return GRIB_SUCCESS;
482 }
483