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
20 /******************************************
21 Elément de base d'un fichier GRIB
22 ******************************************/
23
24 #ifndef GRIBRECORD_H
25 #define GRIBRECORD_H
26
27 #include <iostream>
28 #include <cmath>
29
30 #define DEBUG_INFO false
31 #define DEBUG_ERROR true
32 #define grib_debug(format, ...) {if(DEBUG_INFO) {fprintf(stderr,format,__VA_ARGS__);fprintf(stderr,"\n");}}
33 #define erreur(format, ...) {if(DEBUG_ERROR) {fprintf(stderr,"Grib ERROR: ");fprintf(stderr,format,__VA_ARGS__);fprintf(stderr,"\n");}}
34
35 #define zuint unsigned int
36 #define zuchar unsigned char
37
38 #define GRIB_NOTDEF -999999999
39
40 //--------------------------------------------------------
41 // dataTypes Cf function translateDataType()
42 //--------------------------------------------------------
43 #define GRB_PRESSURE 2 /* Pa */
44 #define GRB_GEOPOT_HGT 7 /* gpm */
45 #define GRB_TEMP 11 /* K */
46 #define GRB_TPOT 13 /* K */
47 #define GRB_TMAX 15 /* K */
48 #define GRB_TMIN 16 /* K */
49 #define GRB_DEWPOINT 17 /* K */
50
51 #define GRB_WIND_DIR 31 /* Deg. Wind Direction */
52 #define GRB_WIND_SPEED 32 /* m/s Wind Speed */
53 #define GRB_WIND_VX 33 /* m/s U */
54 #define GRB_WIND_VY 34 /* m/s V */
55
56 #define GRB_CUR_DIR 47 /* Deg. Direction of current */
57 #define GRB_CUR_SPEED 48 /* m/s Speed of current */
58 #define GRB_UOGRD 49 /*"u-component of current", "m/s" */
59 #define GRB_VOGRD 50 /*"v-component of current", "m/s" */
60
61 #define GRB_HUMID_SPEC 51 /* kg/kg */
62 #define GRB_HUMID_REL 52 /* % */
63 #define GRB_PRECIP_RATE 59 /* l/m2/s */
64 #define GRB_PRECIP_TOT 61 /* l/m2 */
65 #define GRB_SNOW_DEPTH 66 /* m */
66 #define GRB_CLOUD_TOT 71 /* % */
67 #define GRB_HTSGW 100 /* m */
68 #define GRB_WTMP 80 /* "Water Temperature", "K" */
69 #define GRB_COMP_REFL 212 /* dBZ */
70
71 #define GRB_WVDIR 101
72 #define GRB_WVHGT 102
73 #define GRB_WVPER 103
74 #define GRB_SWDIR 104
75 #define GRB_SWELL 105
76 #define GRB_SWPER 106
77 #define GRB_DIRPW 107
78 #define GRB_PERPW 108
79 #define GRB_DIRSW 109
80 #define GRB_PERSW 110
81 #define GRB_PER 209
82 #define GRB_DIR 210
83
84 #define GRB_CRAIN 140 /* "Categorical rain", "yes=1;no=0" */
85 #define GRB_FRZRAIN_CATEG 141 /* 1=yes 0=no */
86 #define GRB_SNOW_CATEG 143 /* 1=yes 0=no */
87 #define GRB_CAPE 157 /* J/kg */
88
89 #define GRB_TSEC 171 /* "Seconds prior to initial reference time (defined in bytes 18-20)" */
90 #define GRB_WIND_GUST 180 /* m/s "wind gust */
91 #define GRB_WIND_GUST_VX 181 /* m/s */
92 #define GRB_WIND_GUST_VY 182 /* m/s */
93
94 #define GRB_USCT 190 /* Scatterometer estimated U Wind, NCEP Center 7 */
95 #define GRB_VSCT 191 /* Scatterometer estimated V Wind, NCEP Center 7 */
96
97 #define GRB_WIND_XY2D 250 /* private : GRB_WIND_VX+GRB_WIND_VX */
98 #define GRB_DIFF_TEMPDEW 251 /* private : GRB_TEMP-GRB_DEWPOINT */
99
100 //--------------------------------------------------------
101 // Levels types (altitude reference)
102 //--------------------------------------------------------
103 #define LV_GND_SURF 1
104 #define LV_ISOTHERM0 4
105 #define LV_ISOBARIC 100
106 #define LV_MSL 102
107 #define LV_ABOV_MSL 103
108 #define LV_ABOV_GND 105
109 #define LV_SIGMA 107
110 #define LV_ATMOS_ENT 10
111 #define LV_ATMOS_ALL 200
112 //---------------------------------------------------------
113 enum DataCenterModel {
114 NOAA_GFS,
115 NOAA_NCEP_WW3,
116 NOAA_NCEP_SST,
117 NOAA_RTOFS,
118 FNMOC_WW3_GLB,
119 FNMOC_WW3_MED,
120 NORWAY_METNO,
121 ECMWF_ERA5,
122 KNMI_HIRLAM,
123 KNMI_HARMONIE_AROME,
124 OTHER_DATA_CENTER
125 };
126
127 //----------------------------------------------
128 class GribCode
129 {
130 public:
makeCode(zuchar dataType,zuchar levelType,zuint levelValue)131 static zuint makeCode (zuchar dataType, zuchar levelType, zuint levelValue) {
132 return ((levelValue&0xFFFF)<<16)+((levelType&0xFF)<<8)+dataType;
133 }
getDataType(zuint code)134 static zuchar getDataType (zuint code) {
135 return code&0xFF;
136 }
getLevelType(zuint code)137 static zuchar getLevelType (zuint code) {
138 return (code>>8)&0xFF;
139 }
getLevelValue(zuint code)140 static zuint getLevelValue (zuint code) {
141 return (code>>16)&0xFFFF;
142 }
143 };
144
145 //----------------------------------------------
146 class GribRecord
147 {
148 public:
149 GribRecord(const GribRecord &rec);
GribRecord()150 GribRecord() { m_bfilled = false;}
151
152 virtual ~GribRecord();
153
154
155 static GribRecord *InterpolatedRecord(const GribRecord &rec1, const GribRecord &rec2, double d, bool dir=false);
156 static GribRecord *Interpolated2DRecord(GribRecord *&rety,
157 const GribRecord &rec1x, const GribRecord &rec1y,
158 const GribRecord &rec2x, const GribRecord &rec2y, double d);
159
160 static GribRecord *MagnitudeRecord(const GribRecord &rec1, const GribRecord &rec2);
161
162 static void Polar2UV(GribRecord *pDIR, GribRecord *pSPEED);
163
164 void multiplyAllData(double k);
165 void Substract(const GribRecord &rec, bool positive=true);
166 void Average(const GribRecord &rec);
167
isOk()168 bool isOk() const {return ok;};
isDataKnown()169 bool isDataKnown() const {return knownData;};
isEof()170 bool isEof() const {return eof;};
isDuplicated()171 bool isDuplicated() const {return IsDuplicated;};
172 //-----------------------------------------
getDataType()173 zuchar getDataType() const { return dataType; }
174 void setDataType(const zuchar t);
175
getLevelType()176 zuchar getLevelType() const { return levelType; }
getLevelValue()177 zuint getLevelValue() const { return levelValue; }
getDataCenterModel()178 zuint getDataCenterModel() const { return dataCenterModel; }
179 //-----------------------------------------
180
getIdCenter()181 zuchar getIdCenter() const { return idCenter; }
getIdModel()182 zuchar getIdModel() const { return idModel; }
getIdGrid()183 zuchar getIdGrid() const { return idGrid; }
184
185 //-----------------------------------------
getKey()186 std::string getKey() const { return dataKey; }
187 static std::string makeKey(int dataType,int levelType,int levelValue);
188
189 //-----------------------------------------
getPeriodP1()190 int getPeriodP1() const { return periodP1; }
getPeriodP2()191 int getPeriodP2() const { return periodP2; }
getPeriodSec()192 zuint getPeriodSec() const { return periodsec; }
getTimeRange()193 zuchar getTimeRange() const { return timeRange; }
194
195 // Number of points in the grid
getNi()196 int getNi() const { return Ni; }
getNj()197 int getNj() const { return Nj; }
getDi()198 double getDi() const { return Di; }
getDj()199 double getDj() const { return Dj; }
200
201 // Value at one point of the grid
getValue(int i,int j)202 double getValue(int i, int j) const { return data[j*Ni+i];}
203
setValue(zuint i,zuint j,double v)204 void setValue(zuint i, zuint j, double v)
205 { if (i<Ni && j<Nj)
206 data[j*Ni+i] = v; }
207
208 // Value for one point interpolated
209 double getInterpolatedValue(double px, double py, bool numericalInterpolation=true, bool dir=false) const;
210
211 // Value for polar interpolation of vectors
212 static bool getInterpolatedValues(double &M, double &A,
213 const GribRecord *GRX, const GribRecord *GRY,
214 double px, double py, bool numericalInterpolation=true);
215
216 // coordiantes of grid point
getX(int i)217 inline double getX(int i) const { return Lo1+i*Di;}
getY(int j)218 inline double getY(int j) const { return La1+j*Dj;}
getXY(int i,int j,double * x,double * y)219 void getXY(int i, int j, double *x, double *y) const { *x = getX(i); *y = getY(j);};
220
getLatMin()221 double getLatMin() const { return latMin;}
getLonMin()222 double getLonMin() const { return lonMin;}
getLatMax()223 double getLatMax() const { return latMax;}
getLonMax()224 double getLonMax() const { return lonMax;}
225
226 // Is there a value at a particular grid point ?
227 inline bool hasValue(int i, int j) const;
228 // Is there a value that is not GRIB_NOTDEF ?
isDefined(int i,int j)229 inline bool isDefined(int i, int j) const
230 { return hasValue(i, j) && getValue(i, j) != GRIB_NOTDEF; }
231
232 // Reference date Date (file creation date)
getRecordRefDate()233 time_t getRecordRefDate () const { return refDate; }
getStrRecordRefDate()234 const char* getStrRecordRefDate () const { return strRefDate; }
235
236 // Date courante des prévisions
getRecordCurrentDate()237 time_t getRecordCurrentDate () const { return curDate; }
getStrRecordCurDate()238 const char* getStrRecordCurDate () const { return strCurDate; }
239 void setRecordCurrentDate (time_t t);
240 void print();
isFilled()241 bool isFilled(){ return m_bfilled; }
242 void setFilled(bool val=true){ m_bfilled = val;}
243
244 private:
245 // Is a point within the extent of the grid?
246 inline bool isPointInMap(double x, double y) const;
247 inline bool isXInMap(double x) const;
248 inline bool isYInMap(double y) const;
249
250 protected:
251 //private:
252 static bool GetInterpolatedParameters
253 (const GribRecord &rec1, const GribRecord &rec2,
254 double &La1, double &Lo1, double &La2, double &Lo2, double &Di, double &Dj,
255 int &im1, int &jm1, int &im2, int &jm2,
256 int &Ni, int &Nj, int &rec1offi, int &rec1offj, int &rec2offi, int &rec2offj );
257
258 int id; // unique identifiant
259 bool ok; // valid?
260 bool knownData; // type de donnée connu
261 bool waveData;
262 bool IsDuplicated;
263 bool eof;
264 std::string dataKey;
265 char strRefDate [32];
266 char strCurDate [32];
267 int dataCenterModel;
268 bool m_bfilled;
269
270 //---------------------------------------------
271 // SECTION 0: THE INDICATOR SECTION (IS)
272 //---------------------------------------------
273 zuchar editionNumber;
274
275 // SECTION 1: THE PRODUCT DEFINITION SECTION (PDS)
276 zuchar idCenter;
277 zuchar idModel;
278 zuchar idGrid;
279 zuchar dataType; // octet 9 = parameters and units
280 zuchar levelType;
281 zuint levelValue;
282
283 bool hasBMS;
284 zuint refyear, refmonth, refday, refhour, refminute;
285 //zuchar periodP1, periodP2;
286 zuint periodP1, periodP2;
287 zuchar timeRange;
288 zuint periodsec; // period in seconds
289 time_t refDate; // C reference date
290 time_t curDate; // C current date
291 // SECTION 2: THE GRID DESCRIPTION SECTION (GDS)
292 zuchar NV, PV;
293 zuchar gridType;
294 zuint Ni, Nj;
295 double La1, Lo1, La2, Lo2;
296 double latMin, lonMin, latMax, lonMax;
297 double Di, Dj;
298 zuchar resolFlags, scanFlags;
299 bool hasDiDj;
300 bool isEarthSpheric;
301 bool isUeastVnorth;
302 bool isScanIpositive;
303 bool isScanJpositive;
304 bool isAdjacentI;
305 // SECTION 3: BIT MAP SECTION (BMS)
306 zuint BMSsize;
307 zuchar *BMSbits;
308 // SECTION 4: BINARY DATA SECTION (BDS)
309 double *data;
310 // SECTION 5: END SECTION (ES)
311
312 time_t makeDate(zuint year,zuint month,zuint day,zuint hour,zuint min,zuint sec);
313
314 // void print();
315 };
316
317 //==========================================================================
hasValue(int i,int j)318 inline bool GribRecord::hasValue(int i, int j) const
319 {
320 // is data present in BMS ?
321 if (!hasBMS) {
322 return true;
323 }
324 int bit;
325 if (isAdjacentI) {
326 bit = j*Ni + i;
327 }
328 else {
329 bit = i*Nj + j;
330 }
331 zuchar c = BMSbits[bit/8];
332 zuchar m = (zuchar)128 >> (bit % 8);
333 return (m & c) != 0;
334 }
335
336 //-----------------------------------------------------------------
isPointInMap(double x,double y)337 inline bool GribRecord::isPointInMap(double x, double y) const
338 {
339 return isXInMap(x) && isYInMap(y);
340 /* if (Dj < 0)
341 return x>=Lo1 && y<=La1 && x<=Lo1+(Ni-1)*Di && y>=La1+(Nj-1)*Dj;
342 else
343 return x>=Lo1 && y>=La1 && x<=Lo1+(Ni-1)*Di && y<=La1+(Nj-1)*Dj;*/
344 }
345 //-----------------------------------------------------------------
isXInMap(double x)346 inline bool GribRecord::isXInMap(double x) const
347 {
348 // return x>=Lo1 && x<=Lo1+(Ni-1)*Di;
349 //printf ("%f %f %f\n", Lo1, Lo2, x);
350 if (Di > 0) {
351 double maxLo = Lo2;
352 if(Lo2+Di >= 360) /* grib that covers the whole world */
353 maxLo += Di;
354 return x>=Lo1 && x<=maxLo;
355 } else {
356 double maxLo = Lo1;
357 if(Lo2+Di >= 360) /* grib that covers the whole world */
358 maxLo += Di;
359 return x>=Lo2 && x<=maxLo;
360 }
361 }
362 //-----------------------------------------------------------------
isYInMap(double y)363 inline bool GribRecord::isYInMap(double y) const
364 {
365 if (Dj < 0)
366 return y<=La1 && y>=La2;
367 else
368 return y>=La1 && y<=La2;
369 }
370
371 #endif
372
373
374
375