GribRecord.cpp

Go to the documentation of this file.

/**********************************************************************
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 <http://www.gnu.org/licenses/>.
***********************************************************************/
#include "wx/wxprec.h"
 
#ifndef WX_PRECOMP
#include "wx/wx.h"
#endif  // precompiled headers
 
#include <stdlib.h>
 
// #include <QDateTime>
 
#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 != nullptr) {
    int size = rec.Ni * rec.Nj;
    this->data = new double[size];
    for (int i = 0; i < size; i++) this->data[i] = rec.data[i];
  }
  if (rec.BMSbits != nullptr) {
    int size = rec.BMSsize;
    this->BMSbits = new zuchar[size];
    for (int i = 0; i < size; i++) this->BMSbits[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 < iiters; i++) {
    rec1offdi = (Lo1 - rec1.Lo1) / rec1.Di;
    rec2offdi = (Lo1 - rec2.Lo1) / rec2.Di;
    if (rec1offdi == floor(rec1offdi) && rec2offdi == floor(rec2offdi)) break;
 
    Lo1 += wxMin(rec1.Di, rec2.Di);
  }
  if (i == iiters)  // failed to align, would need spacial interpolation to work
    return false;
 
  double jiters = rec2.Dj / rec1.Dj;
  if (jiters < 1) {
    jiters = 1 / jiters;
    jm1 = 1, jm2 = jiters;
  } else
    jm1 = jiters, jm2 = 1;
 
  for (j = 0; j < jiters; j++) {
    rec1offdj = (La1 - rec1.La1) / rec1.Dj;
    rec2offdj = (La1 - rec2.La1) / rec2.Dj;
    if (rec1offdj == floor(rec1offdj) && rec2offdj == floor(rec2offdj)) break;
 
    La1 += Dj < 0 ? wxMax(rec1.getDj(), rec2.getDj())
                  : wxMin(rec1.getDj(), rec2.getDj());
  }
  if (j == jiters)  // failed to align
    return false;
 
  /* no overlap */
  if (La1 * Dj > 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 nullptr;
 
  // recopie les champs de bits
  int size = Ni * Nj;
  double *data = new double[size];
 
  zuchar *BMSbits = nullptr;
  if (rec1.BMSbits != nullptr && rec2.BMSbits != nullptr)
    BMSbits = new zuchar[(Ni * Nj - 1) / 8 + 1]();
 
  for (int i = 0; i < Ni; i++)
    for (int j = 0; j < Nj; j++) {
      int in = j * Ni + i;
      int i1 = (j * jm1 + rec1offj) * rec1.Ni + i * im1 + rec1offi;
      int i2 = (j * jm2 + rec2offj) * rec2.Ni + i * im2 + rec2offi;
      double data1 = rec1.data[i1], data2 = rec2.data[i2];
      if (data1 == GRIB_NOTDEF || data2 == GRIB_NOTDEF)
        data[in] = GRIB_NOTDEF;
      else {
        if (!dir)
          data[in] = (1 - d) * data1 + d * data2;
        else
          data[in] = interp_angle(data1, data2, d, 180.);
      }
 
      if (BMSbits) {
        int b1 = rec1.BMSbits[i1 >> 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 nullptr;
 
  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 < Ni; i++) {
    for (int j = 0; j < Nj; j++) {
      int in = j * Ni + i;
      int i1 = (j * jm1 + rec1offj) * rec1x.Ni + i * im1 + rec1offi;
      int i2 = (j * jm2 + rec2offj) * rec2x.Ni + i * im2 + rec2offi;
      double data1x = rec1x.data[i1], data1y = rec1y.data[i1];
      double data2x = rec2x.data[i2], data2y = rec2y.data[i2];
      if (data1x == GRIB_NOTDEF || data1y == GRIB_NOTDEF ||
          data2x == GRIB_NOTDEF || data2y == GRIB_NOTDEF) {
        datax[in] = GRIB_NOTDEF;
        datay[in] = GRIB_NOTDEF;
      } else {
        double data1m = sqrt(pow(data1x, 2) + pow(data1y, 2));
        double data2m = sqrt(pow(data2x, 2) + pow(data2y, 2));
        double datam = (1 - d) * data1m + d * data2m;
 
        double data1a = atan2(data1y, data1x);
        double data2a = atan2(data2y, data2x);
        if (data1a - data2a > 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 = nullptr;
  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 = nullptr;
  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; i < size; i++)
      if (rec1.data[i] == GRIB_NOTDEF || rec2.data[i] == GRIB_NOTDEF)
        rec->data[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 != nullptr && rec2.BMSbits != nullptr) {
    if (rec1.BMSsize == rec2.BMSsize) {
      int size = rec1.BMSsize;
      for (int i = 0; i < size; i++)
        rec->BMSbits[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; i < size; i++) {
      if (pDIR->data[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 < size; i++) {
    if (rec.data[i] == GRIB_NOTDEF) continue;
    if (data[i] == GRIB_NOTDEF) {
      data[i] = -rec.data[i];
      if (BMSbits != 0) {
        if (BMSsize > 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; i < size; i++) {
    if (rec.data[i] == GRIB_NOTDEF) continue;
    if (data[i] == GRIB_NOTDEF) continue;
 
    data[i] = (data[i] * d2 - rec.data[i] * d1) / diff;
  }
}
 
//-------------------------------------------------------------------------------
void GribRecord::setDataType(const zuchar t) {
  dataType = t;
  dataKey = makeKey(dataType, levelType, levelValue);
}
//------------------------------------------------------------------------------
std::string GribRecord::makeKey(
    int dataType, int levelType,
    int levelValue) {  // Make data type key  sample:'11-100-850'
                       //  char ktmp[32];
  //  wxSnprintf((wxChar *)ktmp, 32, "%d-%d-%d", dataType, levelType,
  //  levelValue); return std::string(ktmp);
 
  wxString k;
  k.Printf(_T("%d-%d-%d"), dataType, levelType, levelValue);
  return std::string(k.mb_str());
}
//-----------------------------------------
GribRecord::~GribRecord() {
  if (data) {
    delete[] data;
    data = nullptr;
  }
  if (BMSbits) {
    delete[] BMSbits;
    BMSbits = nullptr;
  }
 
  // if (dataType==GRB_TEMP) printf("record destroyed %s   %d\n",
  // dataKey.mb_str(), (int)curDate/3600);
}
 
//-------------------------------------------------------------------------------
void GribRecord::multiplyAllData(double k) {
  if (data == 0 || !isOk()) return;
 
  for (zuint j = 0; j < Nj; j++) {
    for (zuint i = 0; i < Ni; i++) {
      if (isDefined(i, j)) {
        data[j * Ni + i] *= k;
      }
    }
  }
}
 
//----------------------------------------------
void GribRecord::setRecordCurrentDate(time_t t) {
  curDate = t;
 
  struct tm *date = gmtime(&t);
 
  zuint year = date->tm_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 < year; y++) {
    r += 365 * 24 * 3600;
    if (isleapyear(y)) r += 24 * 3600;
  }
  for (zuint m = 1; m < month; m++) {
    if (m == 2) {
      r += 28 * 24 * 3600;
      if (isleapyear(year)) r += 24 * 3600;
    } else if (m == 1 || m == 3 || m == 5 || m == 7 || m == 8 || m == 10 ||
               m == 12) {
      r += 31 * 24 * 3600;
    } else {
      r += 30 * 24 * 3600;
    }
  }
  r += (day - 1) * 24 * 3600;
  r += hour * 3600;
  r += min * 60;
  r += sec;
  return r;
}
 
//===============================================================================================
 
double GribRecord::getInterpolatedValue(double px, double py,
                                        bool numericalInterpolation,
                                        bool dir) const {
  if (!ok || Di == 0 || Dj == 0) return GRIB_NOTDEF;
 
  if (!isPointInMap(px, py)) {
    px += 360.0;  // tour du monde à droite ?
    if (!isPointInMap(px, py)) {
      px -= 2 * 360.0;  // tour du monde à gauche ?
      if (!isPointInMap(px, py)) {
        return GRIB_NOTDEF;
      }
    }
  }
  double pi, pj;  // coord. in grid unit
  pi = (px - Lo1) / Di;
  pj = (py - La1) / 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 >= 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
}