grib_record.h

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/>.
***********************************************************************/
 
#ifndef GRIBRECORD_H
#define GRIBRECORD_H
 
#include <iostream>
#include <cmath>
 
#define DEBUG_INFO false
#define DEBUG_ERROR true
#define grib_debug(format, ...)             \
  {                                         \
    if (DEBUG_INFO) {                       \
      fprintf(stderr, format, __VA_ARGS__); \
      fprintf(stderr, "\n");                \
    }                                       \
  }
#define erreur(format, ...)                 \
  {                                         \
    if (DEBUG_ERROR) {                      \
      fprintf(stderr, "Grib ERROR: ");      \
      fprintf(stderr, format, __VA_ARGS__); \
      fprintf(stderr, "\n");                \
    }                                       \
  }
 
#define zuint unsigned int
#define zuchar unsigned char
 
#define GRIB_NOTDEF (-999999999)
 
//--------------------------------------------------------
// dataTypes    Cf function translateDataType()
//--------------------------------------------------------
#define GRB_PRESSURE 2   /* Pa     */
#define GRB_GEOPOT_HGT 7 /* gpm    */
#define GRB_TEMP 11      /* K      */
#define GRB_TPOT 13      /* K      */
#define GRB_TMAX 15      /* K      */
#define GRB_TMIN 16      /* K      */
#define GRB_DEWPOINT 17  /* K      */
 
#define GRB_WIND_DIR 31   /* Deg. Wind Direction */
#define GRB_WIND_SPEED 32 /* m/s  Wind Speed     */
#define GRB_WIND_VX 33    /* m/s U  */
#define GRB_WIND_VY 34    /* m/s V  */
 
#define GRB_CUR_DIR 47   /* Deg. Direction of current  */
#define GRB_CUR_SPEED 48 /* m/s Speed of current       */
#define GRB_UOGRD 49     /*"u-component of current", "m/s" */
#define GRB_VOGRD 50     /*"v-component of current", "m/s" */
 
#define GRB_HUMID_SPEC 51  /* kg/kg  */
#define GRB_HUMID_REL 52   /* %      */
#define GRB_PRECIP_RATE 59 /* l/m2/s */
#define GRB_PRECIP_TOT 61  /* l/m2   */
#define GRB_SNOW_DEPTH 66  /* m      */
#define GRB_CLOUD_TOT 71   /* %      */
#define GRB_HTSGW 100      /* m      */
#define GRB_WTMP 80        /* "Water Temperature", "K" */
#define GRB_COMP_REFL 212  /* dBZ */
 
#define GRB_WVDIR 101
#define GRB_WVHGT 102
#define GRB_WVPER 103
#define GRB_SWDIR 104
#define GRB_SWELL 105
#define GRB_SWPER 106
#define GRB_DIRPW 107
#define GRB_PERPW 108
#define GRB_DIRSW 109
#define GRB_PERSW 110
#define GRB_PER 209
#define GRB_DIR 210
 
#define GRB_CRAIN 140         /* "Categorical rain", "yes=1;no=0" */
#define GRB_FRZRAIN_CATEG 141 /* 1=yes 0=no */
#define GRB_SNOW_CATEG 143    /* 1=yes 0=no */
#define GRB_CAPE 157          /* J/kg   */
 
#define GRB_TSEC \
  171 /* "Seconds prior to initial reference time (defined in bytes 18-20)" */
#define GRB_WIND_GUST 180    /* m/s "wind gust */
#define GRB_WIND_GUST_VX 181 /* m/s */
#define GRB_WIND_GUST_VY 182 /* m/s */
 
#define GRB_USCT 190 /* Scatterometer estimated U Wind, NCEP Center 7  */
#define GRB_VSCT 191 /* Scatterometer estimated V Wind, NCEP Center 7  */
 
#define GRB_WIND_XY2D 250    /* private : GRB_WIND_VX+GRB_WIND_VX */
#define GRB_DIFF_TEMPDEW 251 /* private : GRB_TEMP-GRB_DEWPOINT */
 
//--------------------------------------------------------
// Levels types (altitude reference)
//--------------------------------------------------------
#define LV_GND_SURF 1
#define LV_ISOTHERM0 4
#define LV_ISOBARIC 100
#define LV_MSL 102
#define LV_ABOV_MSL 103
#define LV_ABOV_GND 105
#define LV_SIGMA 107
#define LV_ATMOS_ENT 10
#define LV_ATMOS_ALL 200
 
//---------------------------------------------------------
enum DataCenterModel {
  NOAA_GFS,
  NOAA_NCEP_WW3,
  NOAA_NCEP_SST,
  NOAA_RTOFS,
  FNMOC_WW3_GLB,
  FNMOC_WW3_MED,
  NORWAY_METNO,
  ECMWF_ERA5,
  KNMI_HIRLAM,
  KNMI_HARMONIE_AROME,
  OTHER_DATA_CENTER
};
 
//----------------------------------------------
class GribCode {
public:
  static zuint makeCode(zuchar dataType, zuchar levelType, zuint levelValue) {
    return ((levelValue & 0xFFFF) << 16) + ((levelType & 0xFF) << 8) + dataType;
  }
  static zuchar getDataType(zuint code) { return code & 0xFF; }
  static zuchar getLevelType(zuint code) { return (code >> 8) & 0xFF; }
  static zuint getLevelValue(zuint code) { return (code >> 16) & 0xFFFF; }
};
 
class GribRecord {
public:
  GribRecord(const GribRecord &rec);
  GribRecord() { is_filled = false; }
 
  virtual ~GribRecord();
 
  static GribRecord *InterpolatedRecord(const GribRecord &rec1,
                                        const GribRecord &rec2, double d,
                                        bool dir = false);
  static GribRecord *Interpolated2DRecord(GribRecord *&rety,
                                          const GribRecord &rec1x,
                                          const GribRecord &rec1y,
                                          const GribRecord &rec2x,
                                          const GribRecord &rec2y, double d);
 
  static GribRecord *MagnitudeRecord(const GribRecord &rec1,
                                     const GribRecord &rec2);
 
  static void Polar2UV(GribRecord *pDIR, GribRecord *pSPEED);
 
  void multiplyAllData(double k);
  void Substract(const GribRecord &rec, bool positive = true);
  void Average(const GribRecord &rec);
 
  [[nodiscard]] bool IsOk() const { return ok; };
  [[nodiscard]] bool IsDataKnown() const { return known_data; };
  [[nodiscard]] bool IsEof() const { return eof; };
  [[nodiscard]] bool isDuplicated() const { return is_duplicated; };
  [[nodiscard]] zuchar GetDataType() const { return data_type; }
  void SetDataType(zuchar t);
 
  [[nodiscard]] zuchar GetLevelType() const { return level_type; }
  [[nodiscard]] zuint GetLevelValue() const { return level_value; }
  [[nodiscard]] zuint GetDataCenterModel() const { return data_center_model; }
  //-----------------------------------------
 
  [[nodiscard]] zuchar GetIdCenter() const { return id_center; }
  [[nodiscard]] zuchar GetIdModel() const { return id_model; }
  [[nodiscard]] zuchar GetIdGrid() const { return id_grid; }
 
  //-----------------------------------------
  [[nodiscard]] std::string GetKey() const { return data_key; }
  static std::string MakeKey(int dataType, int levelType, int levelValue);
 
  //-----------------------------------------
  [[nodiscard]] int GetPeriodP1() const { return static_cast<int>(period_p1); }
  [[nodiscard]] int GetPeriodP2() const { return static_cast<int>(period_p2); }
  [[nodiscard]] zuint GetPeriodSec() const { return periodsec; }
  [[nodiscard]] zuchar GetTimeRange() const { return time_range; }
 
  // Number of points in the grid
  [[nodiscard]] int GetNi() const { return static_cast<int>(Ni); }
  [[nodiscard]] int GetNj() const { return static_cast<int>(Nj); }
  [[nodiscard]] double GetDi() const { return Di; }
  [[nodiscard]] double GetDj() const { return Dj; }
 
  [[nodiscard]] double GetValue(int i, int j) const { return data[j * Ni + i]; }
 
  void SetValue(zuint i, zuint j, double v) {
    if (i < Ni && j < Nj) data[j * Ni + i] = v;
  }
 
  [[nodiscard]] double GetInterpolatedValue(double px, double py,
                                            bool numericalInterpolation = true,
                                            bool dir = false) const;
 
  static bool GetInterpolatedValues(double &M, double &A, const GribRecord *GRX,
                                    const GribRecord *GRY, double px, double py,
                                    bool numericalInterpolation = true);
 
  [[nodiscard]] inline double GetX(int i) const { return Lo1 + i * Di; }
  [[nodiscard]] inline double GetY(int j) const { return La1 + j * Dj; }
  void getXY(int i, int j, double *x, double *y) const {
    *x = GetX(i);
    *y = GetY(j);
  };
 
  [[nodiscard]] double GetLatMin() const { return lat_min; }
  [[nodiscard]] double GetLonMin() const { return lon_min; }
  [[nodiscard]] double GetLatMax() const { return lat_max; }
  [[nodiscard]] double GetLonMax() const { return lon_max; }
 
  // Is there a value at a particular grid point ?
  [[nodiscard]] inline bool HasValue(int i, int j) const;
  // Is there a value that is not GRIB_NOTDEF ?
  [[nodiscard]] inline bool IsDefined(int i, int j) const {
    return HasValue(i, j) && GetValue(i, j) != GRIB_NOTDEF;
  }
 
  // Reference date Date (file creation date)
  [[nodiscard]] time_t GetRecordRefDate() const { return ref_date; }
  [[nodiscard]] const char *GetStrRecordRefDate() const { return str_ref_date; }
 
  // Date courante des prévisions
  [[nodiscard]] time_t GetRecordCurrentDate() const { return cur_date; }
  [[nodiscard]] const char *GetStrRecordCurDate() const { return str_cur_date; }
  void SetRecordCurrentDate(time_t t);
  void Print() const;
  [[nodiscard]] bool IsFilled() const { return is_filled; }
  void SetFilled(bool val = true) { is_filled = val; }
 
private:
  // Is a point within the extent of the grid?
  [[nodiscard]] inline bool IsPointInMap(double x, double y) const;
  [[nodiscard]] inline bool IsXInMap(double x) const;
  [[nodiscard]] inline bool IsYInMap(double y) const;
 
protected:
  // private:
  static bool 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);
 
  int id;
  bool ok;
  bool known_data;
  bool wave_data;
  bool is_duplicated;
  bool eof;
  std::string data_key;
  char str_ref_date[32];
  char str_cur_date[32];
  int data_center_model;
  bool is_filled;
 
  //---------------------------------------------
  // SECTION 0: THE INDICATOR SECTION (IS)
  //---------------------------------------------
  zuchar edition_number;
 
  // SECTION 1: THE PRODUCT DEFINITION SECTION (PDS)
  zuchar id_center;
  zuchar id_model;
  zuchar id_grid;
  zuchar data_type;  // octet 9 = parameters and units
  zuchar level_type;
  zuint level_value;
 
  bool hasBMS;
  zuint refyear, refmonth, refday, refhour, refminute;
  // zuchar period_p1, period_p2;
  zuint period_p1, period_p2;
  zuchar time_range;
  zuint periodsec;
  time_t ref_date;
  time_t cur_date;
  // SECTION 2: THE GRID DESCRIPTION SECTION (GDS)
  zuchar NV, PV;
  zuchar grid_type;
  zuint Ni, Nj;
  double La1, Lo1;  
  double La2, Lo2;  
  double lat_min, lon_min, lat_max, lon_max;
  double Di, Dj;
  zuchar resol_flags, scan_flags;
  bool has_di_dj;
  bool is_earth_spheric;
  bool is_ueast_vnorth;
  bool is_scan_i_positive;
  bool is_scan_j_positive;
  bool is_adjacent_i;
  // SECTION 3: BIT MAP SECTION (BMS)
  zuint bms_size;
  zuchar *bms_bits;
  // SECTION 4: BINARY DATA SECTION (BDS)
  double *data;
  // SECTION 5: END SECTION (ES)
 
  time_t MakeDate(zuint year, zuint month, zuint day, zuint hour, zuint min,
                  zuint sec);
 
  //        void   print();
};
 
//==========================================================================
inline bool GribRecord::HasValue(int i, int j) const {
  // is data present in BMS ?
  if (!hasBMS) {
    return true;
  }
  int bit;
  if (is_adjacent_i) {
    bit = static_cast<int>(j * Ni + i);
  } else {
    bit = static_cast<int>(i * Nj + j);
  }
  zuchar c = bms_bits[bit / 8];
  zuchar m = static_cast<zuchar>(128) >> (bit % 8);
  return (m & c) != 0;
}
 
//-----------------------------------------------------------------
inline bool GribRecord::IsPointInMap(double x, double y) const {
  return IsXInMap(x) && IsYInMap(y);
  /*    if (Dj < 0)
          return x>=Lo1 && y<=La1 && x<=Lo1+(Ni-1)*Di && y>=La1+(Nj-1)*Dj;
      else
          return x>=Lo1 && y>=La1 && x<=Lo1+(Ni-1)*Di && y<=La1+(Nj-1)*Dj;*/
}
//-----------------------------------------------------------------
inline bool GribRecord::IsXInMap(double x) const {
  //    return x>=Lo1 && x<=Lo1+(Ni-1)*Di;
  // printf ("%f %f %f\n", Lo1, Lo2, x);
  if (Di > 0) {
    double maxLo = Lo2;
    if (Lo2 + Di >= 360) /* grib that covers the whole world */
      maxLo += Di;
    return x >= Lo1 && x <= maxLo;
  } else {
    double maxLo = Lo1;
    if (Lo2 + Di >= 360) /* grib that covers the whole world */
      maxLo += Di;
    return x >= Lo2 && x <= maxLo;
  }
}
//-----------------------------------------------------------------
inline bool GribRecord::IsYInMap(double y) const {
  if (Dj < 0)
    return y <= La1 && y >= La2;
  else
    return y >= La1 && y <= La2;
}
 
#endif