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