1 /**********************************************************************
2 zyGrib: meteorological GRIB file viewer
3 Copyright (C) 2008 - Jacques Zaninetti - http://www.zygrib.org
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 ***********************************************************************/
18
19 #include "wx/wxprec.h"
20
21 #ifndef WX_PRECOMP
22 #include "wx/wx.h"
23 #endif //precompiled headers
24
25 //#include "dychart.h" // for some compile time fixups
26 //#include "cutil.h"
27 #include <stdlib.h>
28
29 //#include <QDateTime>
30
31 #include "GribRecord.h"
32
33 // interpolate two angles in range +- 180 or +-PI, with resulting angle in the same range
interp_angle(double a0,double a1,double d,double p)34 static double interp_angle(double a0, double a1, double d, double p)
35 {
36 if(a0 - a1 > p) a0 -= 2*p;
37 else if(a1 - a0 > p) a1 -= 2*p;
38 double a = (1-d)*a0 + d*a1;
39 if(a < ( p == 180. ? 0. : -p ) ) a += 2*p;
40 return a;
41 }
42
43 //-------------------------------------------------------------------------------
print()44 void GribRecord::print()
45 {
46 printf("%d: idCenter=%d idModel=%d idGrid=%d dataType=%d levelType=%d levelValue=%d hr=%f\n",
47 id, idCenter, idModel, idGrid, dataType, levelType,levelValue,
48 (curDate-refDate)/3600.0
49 );
50 }
51
52 //-------------------------------------------------------------------------------
53 // Constructeur de recopie
54 //-------------------------------------------------------------------------------
GribRecord(const GribRecord & rec)55 GribRecord::GribRecord(const GribRecord &rec)
56 {
57 *this = rec;
58 IsDuplicated = true;
59 // recopie les champs de bits
60 if (rec.data != NULL) {
61 int size = rec.Ni*rec.Nj;
62 this->data = new double[size];
63 for (int i=0; i<size; i++)
64 this->data[i] = rec.data[i];
65 }
66 if (rec.BMSbits != NULL) {
67 int size = rec.BMSsize;
68 this->BMSbits = new zuchar[size];
69 for (int i=0; i<size; i++)
70 this->BMSbits[i] = rec.BMSbits[i];
71 }
72 }
73
GetInterpolatedParameters(const GribRecord & rec1,const GribRecord & rec2,double & La1,double & Lo1,double & La2,double & Lo2,double & Di,double & Dj,int & im1,int & jm1,int & im2,int & jm2,int & Ni,int & Nj,int & rec1offi,int & rec1offj,int & rec2offi,int & rec2offj)74 bool GribRecord::GetInterpolatedParameters
75 (const GribRecord &rec1, const GribRecord &rec2,
76 double &La1, double &Lo1, double &La2, double &Lo2, double &Di, double &Dj,
77 int &im1, int &jm1, int &im2, int &jm2,
78 int &Ni, int &Nj, int &rec1offi, int &rec1offj, int &rec2offi, int &rec2offj )
79 {
80 if(!rec1.isOk() || !rec2.isOk())
81 return false;
82
83 /* make sure Dj both have same sign */
84 if(rec1.getDj() * rec2.getDj() <= 0)
85 return false;
86
87 Di = wxMax(rec1.getDi(), rec2.getDi());
88 Dj = rec1.getDj() > 0 ?
89 wxMax(rec1.getDj(), rec2.getDj()):
90 wxMin(rec1.getDj(), rec2.getDj());
91
92 /* get overlapping region */
93 if(Dj > 0)
94 La1 = wxMax(rec1.La1, rec2.La1), La2 = wxMin(rec1.La2, rec2.La2);
95 else
96 La1 = wxMin(rec1.La1, rec2.La1), La2 = wxMax(rec1.La2, rec2.La2);
97
98 Lo1 = wxMax(rec1.Lo1, rec2.Lo1), Lo2 = wxMin(rec1.Lo2, rec2.Lo2);
99
100 // align gribs on integer boundaries
101 int i, j;
102 // shut up compiler warning 'may be used uninitialized'
103 // rec2.Dj / rec1.Dj > 0
104 // XXX Is it true for rec2.Di / rec1.Di ?
105 double rec1offdi = 0, rec2offdi = 0;
106 double rec1offdj = 0., rec2offdj = 0.;
107
108 double iiters = rec2.Di / rec1.Di;
109 if(iiters < 1) {
110 iiters = 1/iiters;
111 im1 = 1, im2 = iiters;
112 } else
113 im1 = iiters, im2 = 1;
114
115 for(i=0; i<iiters; i++) {
116 rec1offdi = (Lo1 - rec1.Lo1)/rec1.Di;
117 rec2offdi = (Lo1 - rec2.Lo1)/rec2.Di;
118 if(rec1offdi == floor(rec1offdi) && rec2offdi == floor(rec2offdi))
119 break;
120
121 Lo1 += wxMin(rec1.Di, rec2.Di);
122 }
123 if(i == iiters) // failed to align, would need spacial interpolation to work
124 return false;
125
126 double jiters = rec2.Dj / rec1.Dj;
127 if(jiters < 1) {
128 jiters = 1/jiters;
129 jm1 = 1, jm2 = jiters;
130 } else
131 jm1 = jiters, jm2 = 1;
132
133 for(j=0; j<jiters; j++) {
134 rec1offdj = (La1 - rec1.La1)/rec1.Dj;
135 rec2offdj = (La1 - rec2.La1)/rec2.Dj;
136 if(rec1offdj == floor(rec1offdj) && rec2offdj == floor(rec2offdj))
137 break;
138
139 La1 += Dj < 0 ?
140 wxMax(rec1.getDj(), rec2.getDj()):
141 wxMin(rec1.getDj(), rec2.getDj());
142 }
143 if(j == jiters) // failed to align
144 return false;
145
146 /* no overlap */
147 if(La1*Dj > La2*Dj || Lo1 > Lo2)
148 return false;
149
150 /* compute integer sizes for data array */
151 Ni = (Lo2-Lo1)/Di + 1, Nj = (La2-La1)/Dj + 1;
152
153 /* back-compute final La2 and Lo2 to fit this integer boundary */
154 Lo2 = Lo1 + (Ni-1)*Di, La2 = La1 + (Nj-1)*Dj;
155
156 rec1offi = rec1offdi, rec2offi = rec2offdi;
157 rec1offj = rec1offdj, rec2offj = rec2offdj;
158
159 if (!rec1.data || !rec2.data)
160 return false;
161
162 return true;
163 }
164
165 //-------------------------------------------------------------------------------
166 // Constructeur de interpolate
167 //-------------------------------------------------------------------------------
InterpolatedRecord(const GribRecord & rec1,const GribRecord & rec2,double d,bool dir)168 GribRecord * GribRecord::InterpolatedRecord(const GribRecord &rec1, const GribRecord &rec2, double d, bool dir)
169 {
170 double La1, Lo1, La2, Lo2, Di, Dj;
171 int im1, jm1, im2, jm2;
172 int Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj;
173 if(!GetInterpolatedParameters(rec1, rec2, La1, Lo1, La2, Lo2, Di, Dj,
174 im1, jm1, im2, jm2,
175 Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj))
176 return NULL;
177
178 // recopie les champs de bits
179 int size = Ni*Nj;
180 double *data = new double[size];
181
182 zuchar *BMSbits = NULL;
183 if (rec1.BMSbits != NULL && rec2.BMSbits != NULL)
184 BMSbits = new zuchar[(Ni*Nj-1)/8+1]();
185
186 for (int i=0; i<Ni; i++)
187 for (int j=0; j<Nj; j++) {
188 int in=j*Ni+i;
189 int i1 = (j*jm1+rec1offj)*rec1.Ni + i*im1+rec1offi;
190 int i2 = (j*jm2+rec2offj)*rec2.Ni + i*im2+rec2offi;
191 double data1 = rec1.data[i1], data2 = rec2.data[i2];
192 if(data1 == GRIB_NOTDEF || data2 == GRIB_NOTDEF)
193 data[in] = GRIB_NOTDEF;
194 else {
195 if( !dir )
196 data[in] = (1-d)*data1 + d*data2;
197 else
198 data[in] = interp_angle(data1, data2, d, 180.);
199 }
200
201 if(BMSbits) {
202 int b1 = rec1.BMSbits[i1>>3] & 1<<(i1&7);
203 int b2 = rec2.BMSbits[i2>>3] & 1<<(i2&7);
204 if(b1 && b2)
205 BMSbits[in>>3] |= 1<<(in&7);
206 else
207 BMSbits[in>>3] &= ~(1<<(in&7));
208 }
209 }
210
211 /* should maybe update strCurDate ? */
212
213 GribRecord *ret = new GribRecord;
214 *ret = rec1;
215
216 ret->Di = Di, ret->Dj = Dj;
217 ret->Ni = Ni, ret->Nj = Nj;
218
219 ret->La1 = La1, ret->La2 = La2;
220 ret->Lo1 = Lo1, ret->Lo2 = Lo2;
221
222 ret->data = data;
223 ret->BMSbits = BMSbits;
224
225 ret->latMin = wxMin(La1, La2), ret->latMax = wxMax(La1, La2);
226 ret->lonMin = Lo1, ret->lonMax = Lo2;
227
228 ret->m_bfilled = false;
229
230 return ret;
231 }
232
233 /* for interpolation for x and y records, we must do them together because otherwise
234 we end up with a vector interpolation which is not what we want.. instead we want
235 to interpolate from the polar magnitude, and angles */
Interpolated2DRecord(GribRecord * & rety,const GribRecord & rec1x,const GribRecord & rec1y,const GribRecord & rec2x,const GribRecord & rec2y,double d)236 GribRecord *GribRecord::Interpolated2DRecord(GribRecord *&rety,
237 const GribRecord &rec1x, const GribRecord &rec1y,
238 const GribRecord &rec2x, const GribRecord &rec2y, double d)
239 {
240 double La1, Lo1, La2, Lo2, Di, Dj;
241 int im1, jm1, im2, jm2;
242 int Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj;
243
244 rety = 0;
245 if(!GetInterpolatedParameters(rec1x, rec2x, La1, Lo1, La2, Lo2, Di, Dj,
246 im1, jm1, im2, jm2,
247 Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj))
248 return NULL;
249
250
251 if(!rec1y.data || !rec2y.data || !rec1y.isOk() || !rec2y.isOk() ||
252 rec1x.Di != rec1y.Di ||rec1x.Dj != rec1y.Dj ||
253 rec2x.Di != rec2y.Di ||rec2x.Dj != rec2y.Dj ||
254 rec1x.Ni != rec1y.Ni ||rec1x.Nj != rec1y.Nj ||
255 rec2x.Ni != rec2y.Ni ||rec2x.Nj != rec2y.Nj)
256 {
257 // could also make sure lat and lon min/max are the same...
258 // copy first
259 rety = new GribRecord(rec1y);
260
261 return new GribRecord(rec1x);
262 }
263 // recopie les champs de bits
264 int size = Ni*Nj;
265 double *datax = new double[size], *datay = new double[size];
266 for (int i=0; i<Ni; i++) {
267 for (int j=0; j<Nj; j++) {
268 int in=j*Ni+i;
269 int i1 = (j*jm1+rec1offj)*rec1x.Ni + i*im1+rec1offi;
270 int i2 = (j*jm2+rec2offj)*rec2x.Ni + i*im2+rec2offi;
271 double data1x = rec1x.data[i1], data1y = rec1y.data[i1];
272 double data2x = rec2x.data[i2], data2y = rec2y.data[i2];
273 if(data1x == GRIB_NOTDEF || data1y == GRIB_NOTDEF ||
274 data2x == GRIB_NOTDEF || data2y == GRIB_NOTDEF) {
275 datax[in] = GRIB_NOTDEF;
276 datay[in] = GRIB_NOTDEF;
277 } else {
278 double data1m = sqrt(pow(data1x, 2) + pow(data1y, 2));
279 double data2m = sqrt(pow(data2x, 2) + pow(data2y, 2));
280 double datam = (1-d)*data1m + d*data2m;
281
282 double data1a = atan2(data1y, data1x);
283 double data2a = atan2(data2y, data2x);
284 if(data1a - data2a > M_PI) data1a -= 2*M_PI;
285 else if(data2a - data1a > M_PI) data2a -= 2*M_PI;
286 double dataa = (1-d)*data1a + d*data2a;
287
288 datax[in] = datam*cos(dataa);
289 datay[in] = datam*sin(dataa);
290 }
291 }
292 }
293
294 /* should maybe update strCurDate ? */
295
296 GribRecord *ret = new GribRecord;
297
298 *ret = rec1x;
299
300 ret->Di = Di, ret->Dj = Dj;
301 ret->Ni = Ni, ret->Nj = Nj;
302
303 ret->La1 = La1, ret->La2 = La2;
304 ret->Lo1 = Lo1, ret->Lo2 = Lo2;
305
306 ret->data = datax;
307 ret->BMSbits = NULL;
308 ret->hasBMS = false; // I don't think wind or current ever use BMS correct?
309
310 ret->latMin = wxMin(La1, La2), ret->latMax = wxMax(La1, La2);
311 ret->lonMin = Lo1, ret->lonMax = Lo2;
312
313 rety = new GribRecord;
314 *rety = *ret;
315 rety->dataType = rec1y.dataType;
316 rety->data = datay;
317 rety->BMSbits = NULL;
318 rety->hasBMS = false;
319
320 return ret;
321 }
322
MagnitudeRecord(const GribRecord & rec1,const GribRecord & rec2)323 GribRecord *GribRecord::MagnitudeRecord(const GribRecord &rec1, const GribRecord &rec2)
324 {
325 GribRecord *rec = new GribRecord(rec1);
326
327 /* generate a record which is the combined magnitude of two records */
328 if (rec1.data && rec2.data && rec1.Ni == rec2.Ni && rec1.Nj == rec2.Nj) {
329 int size = rec1.Ni*rec1.Nj;
330 for (int i=0; i<size; i++)
331 if(rec1.data[i] == GRIB_NOTDEF || rec2.data[i] == GRIB_NOTDEF)
332 rec->data[i] = GRIB_NOTDEF;
333 else
334 rec->data[i] = sqrt(pow(rec1.data[i], 2) + pow(rec2.data[i], 2));
335 } else
336 rec->ok=false;
337
338 if (rec1.BMSbits != NULL && rec2.BMSbits != NULL) {
339 if(rec1.BMSsize == rec2.BMSsize) {
340 int size = rec1.BMSsize;
341 for (int i=0; i<size; i++)
342 rec->BMSbits[i] = rec1.BMSbits[i] & rec2.BMSbits[i];
343 } else
344 rec->ok = false;
345 }
346
347 return rec;
348 }
349
Polar2UV(GribRecord * pDIR,GribRecord * pSPEED)350 void GribRecord::Polar2UV(GribRecord *pDIR, GribRecord *pSPEED)
351 {
352 if (pDIR->data && pSPEED->data && pDIR->Ni == pSPEED->Ni && pDIR->Nj == pSPEED->Nj) {
353 int size = pDIR->Ni*pDIR->Nj;
354 for (int i=0; i<size; i++) {
355 if(pDIR->data[i] != GRIB_NOTDEF && pSPEED->data[i] != GRIB_NOTDEF) {
356 double dir = pDIR->data[i];
357 double speed = pSPEED->data[i];
358 pDIR->data[i] = -speed * sin ( dir *M_PI/180.);
359 pSPEED->data[i] = -speed * cos ( dir *M_PI/180.);
360 }
361 }
362 if (pDIR->dataType == GRB_WIND_DIR) {
363 pDIR->dataType = GRB_WIND_VX;
364 pSPEED->dataType = GRB_WIND_VY;
365 }
366 else {
367 pDIR->dataType = GRB_UOGRD;
368 pSPEED->dataType = GRB_VOGRD;
369 }
370 }
371 }
372
Substract(const GribRecord & rec,bool pos)373 void GribRecord::Substract(const GribRecord &rec, bool pos)
374 {
375 // for now only substract records of same size
376 if (rec.data == 0 || !rec.isOk())
377 return;
378
379 if (data == 0 || !isOk())
380 return;
381
382 if (Ni != rec.Ni || Nj != rec.Nj)
383 return;
384
385 zuint size = Ni *Nj;
386 for (zuint i=0; i<size; i++) {
387 if (rec.data[i] == GRIB_NOTDEF)
388 continue;
389 if (data[i] == GRIB_NOTDEF) {
390 data[i] = -rec.data[i];
391 if (BMSbits != 0) {
392 if (BMSsize > i) {
393 BMSbits[i >>3] |= 1 << (i&7);
394 }
395 }
396 }
397 else
398 data[i] -= rec.data[i];
399 if (data[i] < 0. && pos) {
400 // data type should be positive...
401 data[i] = 0.;
402 }
403 }
404 }
405
406 //------------------------------------------------------------------------------
Average(const GribRecord & rec)407 void GribRecord::Average(const GribRecord &rec)
408 {
409
410 // for now only average records of same size
411 // this : 6-12
412 // rec : 6-9
413 // compute average 9-12
414 //
415 // this : 0-12
416 // rec : 0-11
417 // compute average 11-12
418
419 if (rec.data == 0 || !rec.isOk())
420 return;
421
422 if (data == 0 || !isOk())
423 return;
424
425 if (Ni != rec.Ni || Nj != rec.Nj)
426 return;
427
428 if (getPeriodP1() != rec.getPeriodP1())
429 return;
430
431 double d2 = getPeriodP2() - getPeriodP1();
432 double d1 = rec.getPeriodP2() - rec.getPeriodP1();
433
434 if (d2 <= d1)
435 return;
436
437 zuint size = Ni *Nj;
438 double diff = d2 -d1;
439 for (zuint i=0; i<size; i++) {
440 if (rec.data[i] == GRIB_NOTDEF)
441 continue;
442 if (data[i] == GRIB_NOTDEF)
443 continue;
444
445 data[i] = (data[i]*d2 -rec.data[i]*d1)/diff;
446 }
447 }
448
449 //-------------------------------------------------------------------------------
setDataType(const zuchar t)450 void GribRecord::setDataType(const zuchar t)
451 {
452 dataType = t;
453 dataKey = makeKey(dataType, levelType, levelValue);
454 }
455 //------------------------------------------------------------------------------
makeKey(int dataType,int levelType,int levelValue)456 std::string GribRecord::makeKey(int dataType,int levelType,int levelValue)
457 { // Make data type key sample:'11-100-850'
458 // char ktmp[32];
459 // wxSnprintf((wxChar *)ktmp, 32, "%d-%d-%d", dataType, levelType, levelValue);
460 // return std::string(ktmp);
461
462 wxString k;
463 k.Printf(_T("%d-%d-%d"), dataType, levelType, levelValue);
464 return std::string(k.mb_str());
465
466 }
467 //-----------------------------------------
~GribRecord()468 GribRecord::~GribRecord()
469 {
470 if (data) {
471 delete [] data;
472 data = NULL;
473 }
474 if (BMSbits) {
475 delete [] BMSbits;
476 BMSbits = NULL;
477 }
478
479 //if (dataType==GRB_TEMP) printf("record destroyed %s %d\n", dataKey.mb_str(), (int)curDate/3600);
480 }
481
482 //-------------------------------------------------------------------------------
multiplyAllData(double k)483 void GribRecord::multiplyAllData(double k)
484 {
485 if (data == 0 || !isOk())
486 return;
487
488 for (zuint j=0; j<Nj; j++) {
489 for (zuint i=0; i<Ni; i++)
490 {
491 if (isDefined(i,j)) {
492 data[j*Ni+i] *= k;
493 }
494 }
495 }
496 }
497
498 //----------------------------------------------
setRecordCurrentDate(time_t t)499 void GribRecord::setRecordCurrentDate (time_t t)
500 {
501 curDate = t;
502
503 struct tm *date = gmtime(&t);
504
505 zuint year = date->tm_year+1900;
506 zuint month = date->tm_mon+1;
507 zuint day = date->tm_mday;
508 zuint hour = date->tm_hour;
509 zuint minute = date->tm_min;
510 sprintf(strCurDate, "%04d-%02d-%02d %02d:%02d", year,month,day,hour,minute);
511 }
512
513 //----------------------------------------------
isleapyear(zuint y)514 static bool isleapyear( zuint y) { return ( ( y %4 ==0) && ( y % 100 != 0) ) || ( y %400 ==0); }
515
makeDate(zuint year,zuint month,zuint day,zuint hour,zuint min,zuint sec)516 time_t GribRecord::makeDate(
517 zuint year,zuint month,zuint day,zuint hour,zuint min,zuint sec) {
518 if (year<1970 || year>2200 || month<1 || month>12 || day<1)
519 return -1;
520 time_t r = 0;
521
522 // TODO : optimize (precomputed data)
523 for (zuint y=1970; y<year; y++) {
524 r += 365*24*3600;
525 if (isleapyear(y))
526 r += 24*3600;
527 }
528 for (zuint m=1; m<month; m++) {
529 if (m==2) {
530 r += 28*24*3600;
531 if (isleapyear(year))
532 r += 24*3600;
533 }
534 else if (m==1||m==3||m==5||m==7||m==8||m==10||m==12) {
535 r += 31*24*3600;
536 }
537 else {
538 r += 30*24*3600;
539 }
540 }
541 r += (day-1)*24*3600;
542 r += hour*3600;
543 r += min*60;
544 r += sec;
545 return r;
546 }
547
548
549 //===============================================================================================
550
getInterpolatedValue(double px,double py,bool numericalInterpolation,bool dir) const551 double GribRecord::getInterpolatedValue(double px, double py, bool numericalInterpolation, bool dir) const
552 {
553 if (!ok || Di==0 || Dj==0)
554 return GRIB_NOTDEF;
555
556 if (!isPointInMap(px,py)) {
557 px += 360.0; // tour du monde à droite ?
558 if (!isPointInMap(px,py)) {
559 px -= 2*360.0; // tour du monde à gauche ?
560 if (!isPointInMap(px,py)) {
561 return GRIB_NOTDEF;
562 }
563 }
564 }
565 double pi, pj; // coord. in grid unit
566 pi = (px-Lo1)/Di;
567 pj = (py-La1)/Dj;
568
569 // 00 10 point is in a square
570 // 01 11
571 int i0 = (int) pi; // point 00
572 int j0 = (int) pj;
573
574 unsigned int i1 = pi+1, j1 = pj+1;
575
576 if(i1 >= Ni)
577 i1 = i0;
578
579 if(j1 >= Nj)
580 j1 = j0;
581
582 // distances to 00
583 double dx = pi-i0;
584 double dy = pj-j0;
585
586 if (! numericalInterpolation)
587 {
588 if (dx >= 0.5)
589 i0 = i1;
590 if (dy >= 0.5)
591 j0 = j1;
592
593 return getValue(i0, j0);
594 }
595
596 // bool h00,h01,h10,h11;
597 // int nbval = 0; // how many values in grid ?
598 // if ((h00=isDefined(i0, j0)))
599 // nbval ++;
600 // if ((h10=isDefined(i1, j0)))
601 // nbval ++;
602 // if ((h01=isDefined(i0, j1)))
603 // nbval ++;
604 // if ((h11=isDefined(i1, j1)))
605 // nbval ++;
606
607 int nbval = 0; // how many values in grid ?
608 if (getValue(i0, j0) != GRIB_NOTDEF)
609 nbval ++;
610 if (getValue(i1, j0) != GRIB_NOTDEF)
611 nbval ++;
612 if (getValue(i0, j1) != GRIB_NOTDEF)
613 nbval ++;
614 if (getValue(i1, j1) != GRIB_NOTDEF)
615 nbval ++;
616
617
618
619
620
621 if (nbval < 3)
622 return GRIB_NOTDEF;
623
624 dx = (3.0 - 2.0*dx)*dx*dx; // pseudo hermite interpolation
625 dy = (3.0 - 2.0*dy)*dy*dy;
626
627 double xa, xb, xc, kx, ky;
628 // Triangle :
629 // xa xb
630 // xc
631 // kx = distance(xa,x)
632 // ky = distance(xa,y)
633 if (nbval == 4)
634 {
635 double x00 = getValue(i0, j0);
636 double x01 = getValue(i0, j1);
637 double x10 = getValue(i1, j0);
638 double x11 = getValue(i1, j1);
639 if( !dir ) {
640 double x1 = (1.0-dx)*x00 + dx*x10;
641 double x2 = (1.0-dx)*x01 + dx*x11;
642 return (1.0-dy)*x1 + dy*x2;
643 } else {
644 double x1 = interp_angle(x00, x01, dx, 180.);
645 double x2 = interp_angle(x10, x11, dx, 180.);
646 return interp_angle(x1, x2, dy, 180.);
647 }
648 }
649
650 //interpolation with only three points is too hazardous for angles
651 if( dir ) return GRIB_NOTDEF;
652
653 // here nbval==3, check the corner without data
654 if (getValue(i0, j0) == GRIB_NOTDEF) {
655 //printf("! h00 %f %f\n", dx,dy);
656 xa = getValue(i1, j1); // A = point 11
657 xb = getValue(i0, j1); // B = point 01
658 xc = getValue(i1, j0); // C = point 10
659 kx = 1-dx;
660 ky = 1-dy;
661 }
662 else if (getValue(i0, j1) == GRIB_NOTDEF) {
663 //printf("! h01 %f %f\n", dx,dy);
664 xa = getValue(i1, j0); // A = point 10
665 xb = getValue(i1, j1); // B = point 11
666 xc = getValue(i0, j0); // C = point 00
667 kx = dy;
668 ky = 1-dx;
669 }
670 else if (getValue(i1, j0) == GRIB_NOTDEF) {
671 //printf("! h10 %f %f\n", dx,dy);
672 xa = getValue(i0, j1); // A = point 01
673 xb = getValue(i0, j0); // B = point 00
674 xc = getValue(i1, j1); // C = point 11
675 kx = 1-dy;
676 ky = dx;
677 }
678 else {
679 //printf("! h11 %f %f\n", dx,dy);
680 xa = getValue(i0, j0); // A = point 00
681 xb = getValue(i1, j0); // B = point 10
682 xc = getValue(i0, j1); // C = point 01
683 kx = dx;
684 ky = dy;
685 }
686
687 double k = kx + ky;
688 if (k<0 || k>1)
689 return GRIB_NOTDEF;
690
691 if (k == 0)
692 return xa;
693
694 // axes interpolation
695 double vx = k*xb + (1-k)*xa;
696 double vy = k*xc + (1-k)*xa;
697 // diagonal interpolation
698 double k2 = kx / k;
699 return k2*vx + (1-k2)*vy;
700 }
701
getInterpolatedValues(double & M,double & A,const GribRecord * GRX,const GribRecord * GRY,double px,double py,bool numericalInterpolation)702 bool GribRecord::getInterpolatedValues(double &M, double &A,
703 const GribRecord *GRX, const GribRecord *GRY,
704 double px, double py, bool numericalInterpolation)
705 {
706 if(!GRX || !GRY)
707 return false;
708
709 if (!GRX->ok || !GRY->ok || GRX->Di==0 || GRX->Dj==0)
710 return false;
711
712 if (!GRX->isPointInMap(px,py) || !GRY->isPointInMap(px,py)) {
713 px += 360.0; // tour du monde à droite ?
714 if (!GRX->isPointInMap(px,py) || !GRY->isPointInMap(px,py)) {
715 px -= 2*360.0; // tour du monde à gauche ?
716 if (!GRX->isPointInMap(px,py) || !GRY->isPointInMap(px,py)) {
717 return false;
718 }
719 }
720 }
721 double pi, pj; // coord. in grid unit
722 pi = (px-GRX->Lo1)/GRX->Di;
723 pj = (py-GRX->La1)/GRX->Dj;
724
725 // 00 10 point is in a square
726 // 01 11
727 int i0 = (int) pi; // point 00
728 int j0 = (int) pj;
729
730 unsigned int i1 = pi+1, j1 = pj+1;
731 if(i1 >= GRX->Ni)
732 i1 = i0;
733
734 if(j1 >= GRX->Nj)
735 j1 = j0;
736
737 // distances to 00
738 double dx = pi-i0;
739 double dy = pj-j0;
740
741 if (! numericalInterpolation)
742 {
743 double vx, vy;
744 if (dx >= 0.5)
745 i0 = i1;
746 if (dy >= 0.5)
747 j0 = j1;
748
749 vx = GRX->getValue(i0, j0);
750 vy = GRY->getValue(i0, j0);
751 if (vx == GRIB_NOTDEF || vy == GRIB_NOTDEF)
752 return false;
753
754 M = sqrt(vx*vx + vy*vy);
755 A = atan2(-vx, -vy) * 180 / M_PI;
756 return true;
757 }
758
759 // bool h00,h01,h10,h11;
760 // int nbval = 0; // how many values in grid ?
761 // if ((h00=GRX->isDefined(i0, j0) && GRX->isDefined(i0, j0)))
762 // nbval ++;
763 // if ((h10=GRX->isDefined(i1, j0) && GRY->isDefined(i1, j0)))
764 // nbval ++;
765 // if ((h01=GRX->isDefined(i0, j1) && GRY->isDefined(i0, j1)))
766 // nbval ++;
767 // if ((h11=GRX->isDefined(i1, j1) && GRY->isDefined(i1, j1)))
768 // nbval ++;
769
770 int nbval = 0; // how many values in grid ?
771 if (GRY->getValue(i0, j0) != GRIB_NOTDEF)
772 nbval ++;
773 if (GRY->getValue(i1, j0) != GRIB_NOTDEF)
774 nbval ++;
775 if (GRY->getValue(i0, j1) != GRIB_NOTDEF)
776 nbval ++;
777 if (GRY->getValue(i1, j1) != GRIB_NOTDEF)
778 nbval ++;
779
780 if (nbval <= 3)
781 return false;
782
783 nbval = 0; // how many values in grid ?
784 if (GRX->getValue(i0, j0) != GRIB_NOTDEF)
785 nbval ++;
786 if (GRX->getValue(i1, j0) != GRIB_NOTDEF)
787 nbval ++;
788 if (GRX->getValue(i0, j1) != GRIB_NOTDEF)
789 nbval ++;
790 if (GRX->getValue(i1, j1) != GRIB_NOTDEF)
791 nbval ++;
792
793 if (nbval <= 3)
794 return false;
795
796 dx = (3.0 - 2.0*dx)*dx*dx; // pseudo hermite interpolation
797 dy = (3.0 - 2.0*dy)*dy*dy;
798
799 // Triangle :
800 // xa xb
801 // xc
802 // kx = distance(xa,x)
803 // ky = distance(xa,y)
804 if (nbval == 4)
805 {
806 double x00x = GRX->getValue(i0, j0), x00y = GRY->getValue(i0, j0);
807 double x00m = sqrt(x00x*x00x + x00y*x00y), x00a = atan2(x00x, x00y);
808
809 double x01x = GRX->getValue(i0, j1), x01y = GRY->getValue(i0, j1);
810 double x01m = sqrt(x01x*x01x + x01y*x01y), x01a = atan2(x01x, x01y);
811
812 double x10x = GRX->getValue(i1, j0), x10y = GRY->getValue(i1, j0);
813 double x10m = sqrt(x10x*x10x + x10y*x10y), x10a = atan2(x10x, x10y);
814
815 double x11x = GRX->getValue(i1, j1), x11y = GRY->getValue(i1, j1);
816 double x11m = sqrt(x11x*x11x + x11y*x11y), x11a = atan2(x11x, x11y);
817
818 double x0m = (1-dx)*x00m + dx*x10m, x0a = interp_angle(x00a, x10a, dx, M_PI);
819
820 double x1m = (1-dx)*x01m + dx*x11m, x1a = interp_angle(x01a, x11a, dx, M_PI);
821
822 M = (1-dy)*x0m + dy*x1m;
823 A = interp_angle(x0a, x1a, dy, M_PI);
824 A *= 180 / M_PI; // degrees
825 A += 180;
826
827 return true;
828 }
829
830 return false; // TODO: make this work in the cases of only 3 points
831 #if 0
832 double xa, xb, xc, kx, ky;
833 // here nbval==3, check the corner without data
834 if (!h00) {
835 //printf("! h00 %f %f\n", dx,dy);
836 xa = getValue(i1, j1); // A = point 11
837 xb = getValue(i0, j1); // B = point 01
838 xc = getValue(i1, j0); // C = point 10
839 kx = 1-dx;
840 ky = 1-dy;
841 }
842 else if (!h01) {
843 //printf("! h01 %f %f\n", dx,dy);
844 xa = getValue(i1, j0); // A = point 10
845 xb = getValue(i1, j1); // B = point 11
846 xc = getValue(i0, j0); // C = point 00
847 kx = dy;
848 ky = 1-dx;
849 }
850 else if (!h10) {
851 //printf("! h10 %f %f\n", dx,dy);
852 xa = getValue(i0, j1); // A = point 01
853 xb = getValue(i0, j0); // B = point 00
854 xc = getValue(i1, j1); // C = point 11
855 kx = 1-dy;
856 ky = dx;
857 }
858 else {
859 //printf("! h11 %f %f\n", dx,dy);
860 xa = getValue(i0, j0); // A = point 00
861 xb = getValue(i1, j0); // B = point 10
862 xc = getValue(i0, j1); // C = point 01
863 kx = dx;
864 ky = dy;
865 }
866 }
867 double k = kx + ky;
868 if (k<0 || k>1) {
869 val = GRIB_NOTDEF;
870 }
871 else if (k == 0) {
872 val = xa;
873 }
874 else {
875 // axes interpolation
876 double vx = k*xb + (1-k)*xa;
877 double vy = k*xc + (1-k)*xa;
878 // diagonal interpolation
879 double k2 = kx / k;
880 val = k2*vx + (1-k2)*vy;
881 }
882 return val;
883 #endif
884 }
885