/**********************************************************************
zyGrib: meteorological GRIB file viewer
Copyright (C) 2008 - Jacques Zaninetti - http://www.zygrib.org
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
***********************************************************************/
#include "wx/wxprec.h"
#ifndef WX_PRECOMP
#include "wx/wx.h"
#endif //precompiled headers
//#include "dychart.h" // for some compile time fixups
//#include "cutil.h"
#include
//#include
#include "GribRecord.h"
// interpolate two angles in range +- 180 or +-PI, with resulting angle in the same range
static double interp_angle(double a0, double a1, double d, double p)
{
if(a0 - a1 > p) a0 -= 2*p;
else if(a1 - a0 > p) a1 -= 2*p;
double a = (1-d)*a0 + d*a1;
if(a < ( p == 180. ? 0. : -p ) ) a += 2*p;
return a;
}
//-------------------------------------------------------------------------------
void GribRecord::print()
{
printf("%d: idCenter=%d idModel=%d idGrid=%d dataType=%d levelType=%d levelValue=%d hr=%f\n",
id, idCenter, idModel, idGrid, dataType, levelType,levelValue,
(curDate-refDate)/3600.0
);
}
//-------------------------------------------------------------------------------
// Constructeur de recopie
//-------------------------------------------------------------------------------
GribRecord::GribRecord(const GribRecord &rec)
{
*this = rec;
IsDuplicated = true;
// recopie les champs de bits
if (rec.data != NULL) {
int size = rec.Ni*rec.Nj;
this->data = new double[size];
for (int i=0; idata[i] = rec.data[i];
}
if (rec.BMSbits != NULL) {
int size = rec.BMSsize;
this->BMSbits = new zuchar[size];
for (int i=0; iBMSbits[i] = rec.BMSbits[i];
}
}
bool GribRecord::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 )
{
if(!rec1.isOk() || !rec2.isOk())
return false;
/* make sure Dj both have same sign */
if(rec1.getDj() * rec2.getDj() <= 0)
return false;
Di = wxMax(rec1.getDi(), rec2.getDi());
Dj = rec1.getDj() > 0 ?
wxMax(rec1.getDj(), rec2.getDj()):
wxMin(rec1.getDj(), rec2.getDj());
/* get overlapping region */
if(Dj > 0)
La1 = wxMax(rec1.La1, rec2.La1), La2 = wxMin(rec1.La2, rec2.La2);
else
La1 = wxMin(rec1.La1, rec2.La1), La2 = wxMax(rec1.La2, rec2.La2);
Lo1 = wxMax(rec1.Lo1, rec2.Lo1), Lo2 = wxMin(rec1.Lo2, rec2.Lo2);
// align gribs on integer boundaries
int i, j;
// shut up compiler warning 'may be used uninitialized'
// rec2.Dj / rec1.Dj > 0
// XXX Is it true for rec2.Di / rec1.Di ?
double rec1offdi = 0, rec2offdi = 0;
double rec1offdj = 0., rec2offdj = 0.;
double iiters = rec2.Di / rec1.Di;
if(iiters < 1) {
iiters = 1/iiters;
im1 = 1, im2 = iiters;
} else
im1 = iiters, im2 = 1;
for(i=0; i La2*Dj || Lo1 > Lo2)
return false;
/* compute integer sizes for data array */
Ni = (Lo2-Lo1)/Di + 1, Nj = (La2-La1)/Dj + 1;
/* back-compute final La2 and Lo2 to fit this integer boundary */
Lo2 = Lo1 + (Ni-1)*Di, La2 = La1 + (Nj-1)*Dj;
rec1offi = rec1offdi, rec2offi = rec2offdi;
rec1offj = rec1offdj, rec2offj = rec2offdj;
if (!rec1.data || !rec2.data)
return false;
return true;
}
//-------------------------------------------------------------------------------
// Constructeur de interpolate
//-------------------------------------------------------------------------------
GribRecord * GribRecord::InterpolatedRecord(const GribRecord &rec1, const GribRecord &rec2, double d, bool dir)
{
double La1, Lo1, La2, Lo2, Di, Dj;
int im1, jm1, im2, jm2;
int Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj;
if(!GetInterpolatedParameters(rec1, rec2, La1, Lo1, La2, Lo2, Di, Dj,
im1, jm1, im2, jm2,
Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj))
return NULL;
// recopie les champs de bits
int size = Ni*Nj;
double *data = new double[size];
zuchar *BMSbits = NULL;
if (rec1.BMSbits != NULL && rec2.BMSbits != NULL)
BMSbits = new zuchar[(Ni*Nj-1)/8+1]();
for (int i=0; i>3] & 1<<(i1&7);
int b2 = rec2.BMSbits[i2>>3] & 1<<(i2&7);
if(b1 && b2)
BMSbits[in>>3] |= 1<<(in&7);
else
BMSbits[in>>3] &= ~(1<<(in&7));
}
}
/* should maybe update strCurDate ? */
GribRecord *ret = new GribRecord;
*ret = rec1;
ret->Di = Di, ret->Dj = Dj;
ret->Ni = Ni, ret->Nj = Nj;
ret->La1 = La1, ret->La2 = La2;
ret->Lo1 = Lo1, ret->Lo2 = Lo2;
ret->data = data;
ret->BMSbits = BMSbits;
ret->latMin = wxMin(La1, La2), ret->latMax = wxMax(La1, La2);
ret->lonMin = Lo1, ret->lonMax = Lo2;
ret->m_bfilled = false;
return ret;
}
/* for interpolation for x and y records, we must do them together because otherwise
we end up with a vector interpolation which is not what we want.. instead we want
to interpolate from the polar magnitude, and angles */
GribRecord *GribRecord::Interpolated2DRecord(GribRecord *&rety,
const GribRecord &rec1x, const GribRecord &rec1y,
const GribRecord &rec2x, const GribRecord &rec2y, double d)
{
double La1, Lo1, La2, Lo2, Di, Dj;
int im1, jm1, im2, jm2;
int Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj;
rety = 0;
if(!GetInterpolatedParameters(rec1x, rec2x, La1, Lo1, La2, Lo2, Di, Dj,
im1, jm1, im2, jm2,
Ni, Nj, rec1offi, rec1offj, rec2offi, rec2offj))
return NULL;
if(!rec1y.data || !rec2y.data || !rec1y.isOk() || !rec2y.isOk() ||
rec1x.Di != rec1y.Di ||rec1x.Dj != rec1y.Dj ||
rec2x.Di != rec2y.Di ||rec2x.Dj != rec2y.Dj ||
rec1x.Ni != rec1y.Ni ||rec1x.Nj != rec1y.Nj ||
rec2x.Ni != rec2y.Ni ||rec2x.Nj != rec2y.Nj)
{
// could also make sure lat and lon min/max are the same...
// copy first
rety = new GribRecord(rec1y);
return new GribRecord(rec1x);
}
// recopie les champs de bits
int size = Ni*Nj;
double *datax = new double[size], *datay = new double[size];
for (int i=0; i M_PI) data1a -= 2*M_PI;
else if(data2a - data1a > M_PI) data2a -= 2*M_PI;
double dataa = (1-d)*data1a + d*data2a;
datax[in] = datam*cos(dataa);
datay[in] = datam*sin(dataa);
}
}
}
/* should maybe update strCurDate ? */
GribRecord *ret = new GribRecord;
*ret = rec1x;
ret->Di = Di, ret->Dj = Dj;
ret->Ni = Ni, ret->Nj = Nj;
ret->La1 = La1, ret->La2 = La2;
ret->Lo1 = Lo1, ret->Lo2 = Lo2;
ret->data = datax;
ret->BMSbits = NULL;
ret->hasBMS = false; // I don't think wind or current ever use BMS correct?
ret->latMin = wxMin(La1, La2), ret->latMax = wxMax(La1, La2);
ret->lonMin = Lo1, ret->lonMax = Lo2;
rety = new GribRecord;
*rety = *ret;
rety->dataType = rec1y.dataType;
rety->data = datay;
rety->BMSbits = NULL;
rety->hasBMS = false;
return ret;
}
GribRecord *GribRecord::MagnitudeRecord(const GribRecord &rec1, const GribRecord &rec2)
{
GribRecord *rec = new GribRecord(rec1);
/* generate a record which is the combined magnitude of two records */
if (rec1.data && rec2.data && rec1.Ni == rec2.Ni && rec1.Nj == rec2.Nj) {
int size = rec1.Ni*rec1.Nj;
for (int i=0; idata[i] = GRIB_NOTDEF;
else
rec->data[i] = sqrt(pow(rec1.data[i], 2) + pow(rec2.data[i], 2));
} else
rec->ok=false;
if (rec1.BMSbits != NULL && rec2.BMSbits != NULL) {
if(rec1.BMSsize == rec2.BMSsize) {
int size = rec1.BMSsize;
for (int i=0; iBMSbits[i] = rec1.BMSbits[i] & rec2.BMSbits[i];
} else
rec->ok = false;
}
return rec;
}
void GribRecord::Polar2UV(GribRecord *pDIR, GribRecord *pSPEED)
{
if (pDIR->data && pSPEED->data && pDIR->Ni == pSPEED->Ni && pDIR->Nj == pSPEED->Nj) {
int size = pDIR->Ni*pDIR->Nj;
for (int i=0; idata[i] != GRIB_NOTDEF && pSPEED->data[i] != GRIB_NOTDEF) {
double dir = pDIR->data[i];
double speed = pSPEED->data[i];
pDIR->data[i] = -speed * sin ( dir *M_PI/180.);
pSPEED->data[i] = -speed * cos ( dir *M_PI/180.);
}
}
if (pDIR->dataType == GRB_WIND_DIR) {
pDIR->dataType = GRB_WIND_VX;
pSPEED->dataType = GRB_WIND_VY;
}
else {
pDIR->dataType = GRB_UOGRD;
pSPEED->dataType = GRB_VOGRD;
}
}
}
void GribRecord::Substract(const GribRecord &rec, bool pos)
{
// for now only substract records of same size
if (rec.data == 0 || !rec.isOk())
return;
if (data == 0 || !isOk())
return;
if (Ni != rec.Ni || Nj != rec.Nj)
return;
zuint size = Ni *Nj;
for (zuint i=0; i i) {
BMSbits[i >>3] |= 1 << (i&7);
}
}
}
else
data[i] -= rec.data[i];
if (data[i] < 0. && pos) {
// data type should be positive...
data[i] = 0.;
}
}
}
//------------------------------------------------------------------------------
void GribRecord::Average(const GribRecord &rec)
{
// for now only average records of same size
// this : 6-12
// rec : 6-9
// compute average 9-12
//
// this : 0-12
// rec : 0-11
// compute average 11-12
if (rec.data == 0 || !rec.isOk())
return;
if (data == 0 || !isOk())
return;
if (Ni != rec.Ni || Nj != rec.Nj)
return;
if (getPeriodP1() != rec.getPeriodP1())
return;
double d2 = getPeriodP2() - getPeriodP1();
double d1 = rec.getPeriodP2() - rec.getPeriodP1();
if (d2 <= d1)
return;
zuint size = Ni *Nj;
double diff = d2 -d1;
for (zuint i=0; itm_year+1900;
zuint month = date->tm_mon+1;
zuint day = date->tm_mday;
zuint hour = date->tm_hour;
zuint minute = date->tm_min;
sprintf(strCurDate, "%04d-%02d-%02d %02d:%02d", year,month,day,hour,minute);
}
//----------------------------------------------
static bool isleapyear( zuint y) { return ( ( y %4 ==0) && ( y % 100 != 0) ) || ( y %400 ==0); }
time_t GribRecord::makeDate(
zuint year,zuint month,zuint day,zuint hour,zuint min,zuint sec) {
if (year<1970 || year>2200 || month<1 || month>12 || day<1)
return -1;
time_t r = 0;
// TODO : optimize (precomputed data)
for (zuint y=1970; y= Ni)
i1 = i0;
if(j1 >= Nj)
j1 = j0;
// distances to 00
double dx = pi-i0;
double dy = pj-j0;
if (! numericalInterpolation)
{
if (dx >= 0.5)
i0 = i1;
if (dy >= 0.5)
j0 = j1;
return getValue(i0, j0);
}
// bool h00,h01,h10,h11;
// int nbval = 0; // how many values in grid ?
// if ((h00=isDefined(i0, j0)))
// nbval ++;
// if ((h10=isDefined(i1, j0)))
// nbval ++;
// if ((h01=isDefined(i0, j1)))
// nbval ++;
// if ((h11=isDefined(i1, j1)))
// nbval ++;
int nbval = 0; // how many values in grid ?
if (getValue(i0, j0) != GRIB_NOTDEF)
nbval ++;
if (getValue(i1, j0) != GRIB_NOTDEF)
nbval ++;
if (getValue(i0, j1) != GRIB_NOTDEF)
nbval ++;
if (getValue(i1, j1) != GRIB_NOTDEF)
nbval ++;
if (nbval < 3)
return GRIB_NOTDEF;
dx = (3.0 - 2.0*dx)*dx*dx; // pseudo hermite interpolation
dy = (3.0 - 2.0*dy)*dy*dy;
double xa, xb, xc, kx, ky;
// Triangle :
// xa xb
// xc
// kx = distance(xa,x)
// ky = distance(xa,y)
if (nbval == 4)
{
double x00 = getValue(i0, j0);
double x01 = getValue(i0, j1);
double x10 = getValue(i1, j0);
double x11 = getValue(i1, j1);
if( !dir ) {
double x1 = (1.0-dx)*x00 + dx*x10;
double x2 = (1.0-dx)*x01 + dx*x11;
return (1.0-dy)*x1 + dy*x2;
} else {
double x1 = interp_angle(x00, x01, dx, 180.);
double x2 = interp_angle(x10, x11, dx, 180.);
return interp_angle(x1, x2, dy, 180.);
}
}
//interpolation with only three points is too hazardous for angles
if( dir ) return GRIB_NOTDEF;
// here nbval==3, check the corner without data
if (getValue(i0, j0) == GRIB_NOTDEF) {
//printf("! h00 %f %f\n", dx,dy);
xa = getValue(i1, j1); // A = point 11
xb = getValue(i0, j1); // B = point 01
xc = getValue(i1, j0); // C = point 10
kx = 1-dx;
ky = 1-dy;
}
else if (getValue(i0, j1) == GRIB_NOTDEF) {
//printf("! h01 %f %f\n", dx,dy);
xa = getValue(i1, j0); // A = point 10
xb = getValue(i1, j1); // B = point 11
xc = getValue(i0, j0); // C = point 00
kx = dy;
ky = 1-dx;
}
else if (getValue(i1, j0) == GRIB_NOTDEF) {
//printf("! h10 %f %f\n", dx,dy);
xa = getValue(i0, j1); // A = point 01
xb = getValue(i0, j0); // B = point 00
xc = getValue(i1, j1); // C = point 11
kx = 1-dy;
ky = dx;
}
else {
//printf("! h11 %f %f\n", dx,dy);
xa = getValue(i0, j0); // A = point 00
xb = getValue(i1, j0); // B = point 10
xc = getValue(i0, j1); // C = point 01
kx = dx;
ky = dy;
}
double k = kx + ky;
if (k<0 || k>1)
return GRIB_NOTDEF;
if (k == 0)
return xa;
// axes interpolation
double vx = k*xb + (1-k)*xa;
double vy = k*xc + (1-k)*xa;
// diagonal interpolation
double k2 = kx / k;
return k2*vx + (1-k2)*vy;
}
bool GribRecord::getInterpolatedValues(double &M, double &A,
const GribRecord *GRX, const GribRecord *GRY,
double px, double py, bool numericalInterpolation)
{
if(!GRX || !GRY)
return false;
if (!GRX->ok || !GRY->ok || GRX->Di==0 || GRX->Dj==0)
return false;
if (!GRX->isPointInMap(px,py) || !GRY->isPointInMap(px,py)) {
px += 360.0; // tour du monde à droite ?
if (!GRX->isPointInMap(px,py) || !GRY->isPointInMap(px,py)) {
px -= 2*360.0; // tour du monde à gauche ?
if (!GRX->isPointInMap(px,py) || !GRY->isPointInMap(px,py)) {
return false;
}
}
}
double pi, pj; // coord. in grid unit
pi = (px-GRX->Lo1)/GRX->Di;
pj = (py-GRX->La1)/GRX->Dj;
// 00 10 point is in a square
// 01 11
int i0 = (int) pi; // point 00
int j0 = (int) pj;
unsigned int i1 = pi+1, j1 = pj+1;
if(i1 >= GRX->Ni)
i1 = i0;
if(j1 >= GRX->Nj)
j1 = j0;
// distances to 00
double dx = pi-i0;
double dy = pj-j0;
if (! numericalInterpolation)
{
double vx, vy;
if (dx >= 0.5)
i0 = i1;
if (dy >= 0.5)
j0 = j1;
vx = GRX->getValue(i0, j0);
vy = GRY->getValue(i0, j0);
if (vx == GRIB_NOTDEF || vy == GRIB_NOTDEF)
return false;
M = sqrt(vx*vx + vy*vy);
A = atan2(-vx, -vy) * 180 / M_PI;
return true;
}
// bool h00,h01,h10,h11;
// int nbval = 0; // how many values in grid ?
// if ((h00=GRX->isDefined(i0, j0) && GRX->isDefined(i0, j0)))
// nbval ++;
// if ((h10=GRX->isDefined(i1, j0) && GRY->isDefined(i1, j0)))
// nbval ++;
// if ((h01=GRX->isDefined(i0, j1) && GRY->isDefined(i0, j1)))
// nbval ++;
// if ((h11=GRX->isDefined(i1, j1) && GRY->isDefined(i1, j1)))
// nbval ++;
int nbval = 0; // how many values in grid ?
if (GRY->getValue(i0, j0) != GRIB_NOTDEF)
nbval ++;
if (GRY->getValue(i1, j0) != GRIB_NOTDEF)
nbval ++;
if (GRY->getValue(i0, j1) != GRIB_NOTDEF)
nbval ++;
if (GRY->getValue(i1, j1) != GRIB_NOTDEF)
nbval ++;
if (nbval <= 3)
return false;
nbval = 0; // how many values in grid ?
if (GRX->getValue(i0, j0) != GRIB_NOTDEF)
nbval ++;
if (GRX->getValue(i1, j0) != GRIB_NOTDEF)
nbval ++;
if (GRX->getValue(i0, j1) != GRIB_NOTDEF)
nbval ++;
if (GRX->getValue(i1, j1) != GRIB_NOTDEF)
nbval ++;
if (nbval <= 3)
return false;
dx = (3.0 - 2.0*dx)*dx*dx; // pseudo hermite interpolation
dy = (3.0 - 2.0*dy)*dy*dy;
// Triangle :
// xa xb
// xc
// kx = distance(xa,x)
// ky = distance(xa,y)
if (nbval == 4)
{
double x00x = GRX->getValue(i0, j0), x00y = GRY->getValue(i0, j0);
double x00m = sqrt(x00x*x00x + x00y*x00y), x00a = atan2(x00x, x00y);
double x01x = GRX->getValue(i0, j1), x01y = GRY->getValue(i0, j1);
double x01m = sqrt(x01x*x01x + x01y*x01y), x01a = atan2(x01x, x01y);
double x10x = GRX->getValue(i1, j0), x10y = GRY->getValue(i1, j0);
double x10m = sqrt(x10x*x10x + x10y*x10y), x10a = atan2(x10x, x10y);
double x11x = GRX->getValue(i1, j1), x11y = GRY->getValue(i1, j1);
double x11m = sqrt(x11x*x11x + x11y*x11y), x11a = atan2(x11x, x11y);
double x0m = (1-dx)*x00m + dx*x10m, x0a = interp_angle(x00a, x10a, dx, M_PI);
double x1m = (1-dx)*x01m + dx*x11m, x1a = interp_angle(x01a, x11a, dx, M_PI);
M = (1-dy)*x0m + dy*x1m;
A = interp_angle(x0a, x1a, dy, M_PI);
A *= 180 / M_PI; // degrees
A += 180;
return true;
}
return false; // TODO: make this work in the cases of only 3 points
#if 0
double xa, xb, xc, kx, ky;
// here nbval==3, check the corner without data
if (!h00) {
//printf("! h00 %f %f\n", dx,dy);
xa = getValue(i1, j1); // A = point 11
xb = getValue(i0, j1); // B = point 01
xc = getValue(i1, j0); // C = point 10
kx = 1-dx;
ky = 1-dy;
}
else if (!h01) {
//printf("! h01 %f %f\n", dx,dy);
xa = getValue(i1, j0); // A = point 10
xb = getValue(i1, j1); // B = point 11
xc = getValue(i0, j0); // C = point 00
kx = dy;
ky = 1-dx;
}
else if (!h10) {
//printf("! h10 %f %f\n", dx,dy);
xa = getValue(i0, j1); // A = point 01
xb = getValue(i0, j0); // B = point 00
xc = getValue(i1, j1); // C = point 11
kx = 1-dy;
ky = dx;
}
else {
//printf("! h11 %f %f\n", dx,dy);
xa = getValue(i0, j0); // A = point 00
xb = getValue(i1, j0); // B = point 10
xc = getValue(i0, j1); // C = point 01
kx = dx;
ky = dy;
}
}
double k = kx + ky;
if (k<0 || k>1) {
val = GRIB_NOTDEF;
}
else if (k == 0) {
val = xa;
}
else {
// axes interpolation
double vx = k*xb + (1-k)*xa;
double vy = k*xc + (1-k)*xa;
// diagonal interpolation
double k2 = kx / k;
val = k2*vx + (1-k2)*vy;
}
return val;
#endif
}