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