Merge pull request opengribs#269 from did-g/fix_period_unit

fix Accumulation time range use the right unit
alcir · Jan 3, 2021 · b58041c · b58041c
2 parents 9e3e41b + b5c2a3d
commit b58041c
Show file tree

Hide file tree

Showing 4 changed files with 27 additions and 12 deletions.
diff --git a/src/Grib2Record.cpp b/src/Grib2Record.cpp
@@ -401,7 +401,7 @@ void Grib2Record::analyseProductDefinitionTemplate (gribfield  *gfld)
         if (unit_of_time_range(gfld->ipdtmpl[25]) >= 3600) {
             periodP1 = gfld->ipdtmpl[8] *unit_of_time_range(gfld->ipdtmpl[25])/3600;
             periodP2 = periodP1 + gfld->ipdtmpl[26] *unit_of_time_range(gfld->ipdtmpl[25])/3600; // time_length
-            periodsec = 3600;
+            resosec = 3600;
         }
         else {
             DBG("Can't determine forecast date");

diff --git a/src/GribReader.cpp b/src/GribReader.cpp
@@ -569,9 +569,11 @@ void  GribReader::computeAccumulationRecords (DataCode dtc)
 		if ( !rec )
 			continue;
 
-		// XXX double check reference date and timerange 
 		if (prev != nullptr ) {
-			// printf("rec %d, prev %d %d %d\n", rec->getPeriodP1(), prev->getPeriodP1(), prev->getPeriodP2(), prev->getTimeRange());
+#if 0
+			printf("rec %d, prev %d %d timerange %d per sec %d\n", rec->getPeriodP1(), prev->getPeriodP1(), prev->getPeriodP2(), prev->getTimeRange(),
+			 prev->getResoSec());
+#endif
 			if (prev->getPeriodP1() == rec->getPeriodP1()) {
 				if (rec->getTimeRange() == 4) {
 					// accumulation 
@@ -588,7 +590,7 @@ void  GribReader::computeAccumulationRecords (DataCode dtc)
 			}
 			if (p2 > p1 && rec->getTimeRange() == 4) {
 				//p1 and p2 units is getPeriodSec second convert to hour
-				prev->multiplyAllData( 3600.0/((p2 -p1)* prev->getPeriodSec()) );
+				prev->multiplyAllData( 3600.0/((p2 -p1)* prev->getResoSec()) );
 			}
 		}
 		prev = rec;
@@ -597,7 +599,7 @@ void  GribReader::computeAccumulationRecords (DataCode dtc)
 	}
 	if (prev != nullptr && p2 > p1 && prev->getTimeRange() == 4 ) {
 	    // the last one
-        prev->multiplyAllData( 3600.0/((p2 -p1) * prev->getPeriodSec()));
+        prev->multiplyAllData( 3600.0/((p2 -p1) * prev->getResoSec()));
 	}
 }
 

diff --git a/src/GribRecord.cpp b/src/GribRecord.cpp
@@ -1019,6 +1019,7 @@ bool GribRecord::readGribSection1_PDS(ZUFILE* file) {
     periodP2  = data1[19];
     timeRange = data1[20];
     periodsec = periodSeconds(data1[17],data1[18],data1[19],timeRange);
+    resosec   = resoSecond(data1[17]);
     curDate = UTC_mktime(refyear,refmonth,refday,refhour,refminute,periodsec);
 
     int decim;
@@ -1372,7 +1373,7 @@ zuint GribRecord::makeInt3(zuchar a, zuchar b, zuchar c) {
 zuint GribRecord::makeInt2(zuchar b, zuchar c) {
     return ((zuint)b<<8)+(zuint)c;
 }
-//----------------------------------------------
+//-----------------------------[A-----------------
 zuint GribRecord::readPackedBits(const zuchar *buf, zuint first, zuint nbBits)
 {
     zuint oct = first / 8;
@@ -1396,9 +1397,10 @@ void  GribRecord::setRecordCurrentDate (time_t t)
 	zuint minute = date->tm_min;
 	sprintf(strCurDate, "%04d-%02d-%02d %02d:%02d", year,month,day,hour,minute);
 }
+
 //----------------------------------------------
-zuint GribRecord::periodSeconds(zuchar unit,zuchar P1,zuchar P2,zuchar range) {
-    zuint res, dur;
+zuint GribRecord::resoSecond(zuchar unit) const {
+    zuint res;
     switch (unit) {
         case 0: //	Minute
             res = 60; break;
@@ -1420,9 +1422,18 @@ zuint GribRecord::periodSeconds(zuchar unit,zuchar P1,zuchar P2,zuchar range) {
         case 6: //	Normal (30 years)
         case 7: //	Century (100 years)
         default:
-            erreur("id=%d: unknown time unit in PDS b18=%d",id,unit);
             res = 0;
-            ok = false;
+    }
+    return res;
+}
+//----------------------------------------------
+zuint GribRecord::periodSeconds(zuchar unit,zuchar P1,zuchar P2,zuchar range) {
+    zuint res, dur;
+    res = resoSecond(unit);
+    if (res == 0) {
+        erreur("id=%d: unknown time unit in PDS b18=%d",id,unit);
+        res = 0;
+        ok = false;
     }
     debug("id=%d: PDS (time range) b21=%d P1=%d P2=%d",id,range,P1,P2);
     dur = 0;
@@ -1447,7 +1458,6 @@ zuint GribRecord::periodSeconds(zuchar unit,zuchar P1,zuchar P2,zuchar range) {
     return res*dur;
 }
 
-
 //===============================================================================================
 data_t GribRecord::getInterpolatedValue (double lon, double lat, bool interpolate) const
 {

diff --git a/src/GribRecord.h b/src/GribRecord.h
@@ -84,6 +84,7 @@ class GribRecord : public RegularGridRecord
         int    getPeriodP2() const  { return periodP2; }
         zuchar getTimeRange() const { return timeRange; }
         zuint  getPeriodSec() const { return periodsec; }
+        zuint  getResoSec() const { return resosec; }
 
         // Nombre de points de la grille
         int     getNi() const override { return Ni; }
@@ -177,7 +178,8 @@ class GribRecord : public RegularGridRecord
         zuint  refyear, refmonth, refday, refhour, refminute;
         zuchar periodP1{0}, periodP2{0};
         zuchar timeRange{255};
-        zuint  periodsec{0}; // period in seconds
+        zuint  resosec{0};   // period resolution in second
+        zuint  periodsec{0}; // seconds from reference time
         time_t refDate;      // Reference date
         time_t curDate;      // Current date
         double  decimalFactorD;
@@ -239,6 +241,7 @@ class GribRecord : public RegularGridRecord
         zuint  makeInt2(zuchar b, zuchar c);
 
         inline bool   hasValueInBitBMS (int i, int j) const;
+        zuint  resoSecond(zuchar unit) const;
 		zuint  periodSeconds(zuchar unit, zuchar P1, zuchar P2, zuchar range);
 
 		void   checkOrientation ();