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_lambert_conformal{
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 lambert_conformal */
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_lambert_conformal;
70 
71 extern grib_nearest_class* grib_nearest_class_gen;
72 
73 static grib_nearest_class _grib_nearest_class_lambert_conformal = {
74     &grib_nearest_class_gen,                         /* super                     */
75     "lambert_conformal",                         /* name                      */
76     sizeof(grib_nearest_lambert_conformal),      /* 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_lambert_conformal = &_grib_nearest_class_lambert_conformal;
85 
86 
init_class(grib_nearest_class * c)87 static void init_class(grib_nearest_class* c)
88 {
89 }
90 /* END_CLASS_IMP */
91 
init(grib_nearest * nearest,grib_handle * h,grib_arguments * args)92 static int init(grib_nearest* nearest,grib_handle* h,grib_arguments* args)
93 {
94     grib_nearest_lambert_conformal* self = (grib_nearest_lambert_conformal*) nearest;
95     self->Ni  = grib_arguments_get_name(h,args,self->cargs++);
96     self->Nj  = grib_arguments_get_name(h,args,self->cargs++);
97     self->i=(int*)grib_context_malloc(h->context,2*sizeof(int));
98     self->j=(int*)grib_context_malloc(h->context,2*sizeof(int));
99     return 0;
100 }
101 
compare_doubles(const void * a,const void * b,int ascending)102 static int compare_doubles( const void* a, const void* b, int ascending )
103 {
104     /* ascending is a boolean: 0 or 1 */
105     double* arg1 = (double*) a;
106     double* arg2 = (double*) b;
107     if (ascending) {
108         if( *arg1 < *arg2 ) return -1;		/*Smaller values come before larger ones*/
109     } else {
110         if( *arg1 > *arg2 ) return -1;		/*Larger values come before smaller ones*/
111     }
112     if( *arg1 == *arg2 ) return 0;
113     else return 1;
114 }
115 
compare_doubles_ascending(const void * a,const void * b)116 static int compare_doubles_ascending( const void* a, const void* b )
117 {
118     return compare_doubles(a,b,1);
119 }
120 
121 typedef struct PointStore {
122     double m_lat;
123     double m_lon;
124     double m_dist;
125     double m_value;
126     int    m_index;
127 } PointStore ;
128 
129 /* Comparison function to sort points by distance */
compare_points(const void * a,const void * b)130 int compare_points(const void* a, const void* b)
131 {
132     PointStore *pA = (PointStore*)a;
133     PointStore *pB = (PointStore*)b;
134 
135     if (pA->m_dist < pB->m_dist) return -1;
136     if (pA->m_dist > pB->m_dist) return 1;
137     return 0;
138 }
139 
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)140 static int find(grib_nearest* nearest, grib_handle* h,
141         double inlat, double inlon,unsigned long flags,
142         double* outlats,double* outlons,
143         double *values,double *distances,int* indexes, size_t *len)
144 {
145     grib_nearest_lambert_conformal* self = (grib_nearest_lambert_conformal*) nearest;
146     int ret=0, i=0;
147     size_t nvalues=0;
148     long iradius;
149     double radius;
150     grib_iterator* iter=NULL;
151     double lat=0,lon=0;
152 
153     /* array of candidates for nearest neighbours */
154     PointStore* neighbours = NULL;
155 
156     while (inlon<0) inlon+=360;
157     while (inlon>360) inlon-=360;
158 
159     if( (ret =  grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS)
160         return ret;
161     nearest->values_count = nvalues;
162 
163     if (grib_is_missing(h,self->radius,&ret)) {
164         grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->radius);
165         return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
166     }
167 
168     if( (ret =  grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS)
169         return ret;
170     radius=((double)iradius)/1000.0;
171 
172     neighbours = (PointStore*)grib_context_malloc( nearest->context, nvalues*sizeof(PointStore) );
173     for(i=0; i<nvalues; ++i) {
174         neighbours[i].m_dist = 1e10; /* set all distances to large number to begin with */
175         neighbours[i].m_lat=0;
176         neighbours[i].m_lon=0;
177         neighbours[i].m_value=0;
178         neighbours[i].m_index=0;
179     }
180 
181     /* GRIB_NEAREST_SAME_GRID not yet implemented */
182     {
183         double the_value = 0;
184         double min_dist = 1e10;
185         size_t the_index = 0;
186         int ilat=0, ilon=0;
187         int idx_upper=0, idx_lower=0;
188         double lat1=0, lat2=0; /* inlat will be between these */
189         double dist=0;
190         const double LAT_DELTA = 10.0; /* in degrees */
191 
192         if (grib_is_missing(h,self->Ni,&ret)) {
193             grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Ni);
194             return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
195         }
196 
197         if (grib_is_missing(h,self->Nj,&ret)) {
198             grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Nj);
199             return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
200         }
201 
202         self->lons_count=nvalues;  /* Maybe overestimate but safe */
203         self->lats_count=nvalues;
204 
205         if (self->lats) grib_context_free(nearest->context,self->lats);
206         self->lats=(double*)grib_context_malloc( nearest->context, nvalues* sizeof(double));
207         if (!self->lats) return GRIB_OUT_OF_MEMORY;
208 
209         if (self->lons) grib_context_free(nearest->context,self->lons);
210         self->lons=(double*)grib_context_malloc( nearest->context, nvalues*sizeof(double));
211         if (!self->lons) return GRIB_OUT_OF_MEMORY;
212 
213         iter=grib_iterator_new(h,0,&ret);
214         if (ret) return ret;
215         /* First pass: collect all latitudes and longitudes */
216         while(grib_iterator_next(iter,&lat,&lon,&the_value))
217         {
218             ++the_index;
219             Assert(ilat < self->lats_count);
220             Assert(ilon < self->lons_count);
221             self->lats[ilat++]=lat;
222             self->lons[ilon++]=lon;
223         }
224 
225         /* See between which 2 latitudes our point lies */
226         qsort(self->lats, nvalues, sizeof(double), &compare_doubles_ascending);
227         grib_binary_search(self->lats, self->lats_count-1, inlat, &idx_upper, &idx_lower);
228         lat2 = self->lats[idx_upper];
229         lat1 = self->lats[idx_lower];
230         Assert(lat1<=lat2);
231 
232         /* Second pass: Iterate again and collect candidate neighbours */
233         grib_iterator_reset(iter);
234         the_index=0;
235         i = 0;
236         while(grib_iterator_next(iter,&lat,&lon,&the_value))
237         {
238             if (lat > lat2+LAT_DELTA || lat < lat1-LAT_DELTA) {
239                 /* Ignore latitudes too far from our point */
240             }
241             else {
242                 dist = grib_nearest_distance(radius, inlon, inlat, lon, lat);
243                 if (dist < min_dist) min_dist = dist;
244                 /*printf("Candidate: lat=%.5f lon=%.5f dist=%f Idx=%ld Val=%f\n",lat,lon,dist,the_index,the_value);*/
245                 /* store this candidate point */
246                 neighbours[i].m_dist  = dist;
247                 neighbours[i].m_index = the_index;
248                 neighbours[i].m_lat   = lat;
249                 neighbours[i].m_lon   = lon;
250                 neighbours[i].m_value = the_value;
251                 i++;
252             }
253             ++the_index;
254         }
255         /* Sort the candidate neighbours in ascending order of distance */
256         /* The first 4 entries will now be the closest 4 neighbours */
257         qsort(neighbours, nvalues, sizeof(PointStore), &compare_points);
258 
259         grib_iterator_delete(iter);
260     }
261     nearest->h=h;
262 
263     /* Sanity check for sorting */
264 #ifdef DEBUG
265     for(i=0; i<nvalues-1; ++i) {
266         Assert( neighbours[i].m_dist <= neighbours[i+1].m_dist);
267     }
268 #endif
269 
270     /* GRIB_NEAREST_SAME_XXX not yet implemented */
271     if (!self->distances) {
272         self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double));
273     }
274     self->distances[0] = neighbours[0].m_dist;
275     self->distances[1] = neighbours[1].m_dist;
276     self->distances[2] = neighbours[2].m_dist;
277     self->distances[3] = neighbours[3].m_dist;
278 
279     for(i=0; i <4; ++i)
280     {
281         distances[i] = neighbours[i].m_dist;
282         outlats[i]   = neighbours[i].m_lat;
283         outlons[i]   = neighbours[i].m_lon;
284         indexes[i]   = neighbours[i].m_index;
285         values[i]    = neighbours[i].m_value;
286         /*printf("(%f,%f)  i=%d  d=%f  v=%f\n",outlats[i],outlons[i],indexes[i],distances[i],values[i]);*/
287     }
288 
289     free(neighbours);
290     return GRIB_SUCCESS;
291 }
292 
destroy(grib_nearest * nearest)293 static int destroy(grib_nearest* nearest)
294 {
295     grib_nearest_lambert_conformal* self = (grib_nearest_lambert_conformal*) nearest;
296     if (self->lats) grib_context_free(nearest->context,self->lats);
297     if (self->lons) grib_context_free(nearest->context,self->lons);
298     if (self->i) grib_context_free(nearest->context,self->i);
299     if (self->j) grib_context_free(nearest->context,self->j);
300     if (self->k) grib_context_free(nearest->context,self->k);
301     if (self->distances) grib_context_free(nearest->context,self->distances);
302     return GRIB_SUCCESS;
303 }
304