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