Merge branch 'met_update' into 'master'

Met update

See merge request !23

radixio/gfsfile.cc

@@ -88,7 +88,7 @@ GFSFile::GFSFile(std::string file)
   std::string header = rstream->readString(108);
   // initialize
   mLabel.expand(label);
-  mHeader.expand(header);
+  mHeader.expand(header, mLabel);
   mHeader.latlon = false;
   mHeader.global = false;
   mHeader.gbldat = false;
@@ -109,7 +109,7 @@ GFSFile::GFSFile(std::string file)
     label            = recl.substr(0, 50);
     header           = recl.substr(50);
     mLabel.expand(label);
-    mHeader.expand(header);
+    mHeader.expand(header, mLabel);
     radix_tagged_line("Found grid: " << mHeader.toString());
   }
   //
@@ -137,29 +137,33 @@ GFSFile::GFSFile(std::string file)
     stcm1p(mHeader.sync_xp, mHeader.sync_yp, mHeader.sync_lat, mHeader.sync_lon,
            mHeader.ref_lat, mHeader.ref_lon, mHeader.size, mHeader.orient);
   }
-  int kol  = 108;
-  int nrec = nndx;
-  int nlvl = mHeader.nz;
+  int kol     = 108;
+  int nrec    = nndx;
+  size_t nlvl = size_t(mHeader.nz);
+  radix_line("Number of levels (" << nlvl << ")");
   mNumVarb.clear();
   mNumVarb.resize(nlvl);
   mVarbId.clear();
   mVarbId.resize(nlvl);
   mHeight.clear();
   mHeight.resize(nlvl);
-  std::vector<std::vector<int>> chk_sum(mHeader.nz);
-  for (int l = 0; l < nlvl; ++l)
+  std::vector<std::vector<int>> chk_sum(size_t(mHeader.nz));
+  for (size_t l = 0; l < nlvl; ++l)
   {
-    mHeight[l]  = (float)std::atof(header.substr(kol, 6).c_str());
+    mHeight[l]  = float(std::atof(header.substr(kol, 6).c_str()));
     mNumVarb[l] = std::atoi(header.substr(kol + 6, 2).c_str());
+    radix_line("Height(" << mHeight[l] << ") NumVarb(" << mNumVarb[l] << ")");
 
     kol += 8;
 
-    mVarbId[l].resize(mNumVarb[l]);
-    chk_sum[l].resize(mNumVarb[l]);
+    mVarbId[l].resize(size_t(mNumVarb[l]));
+    chk_sum[l].resize(size_t(mNumVarb[l]));
     for (int k = 0; k < mNumVarb[l]; ++k)
     {
       mVarbId[l][k] = header.substr(kol, 4);
       chk_sum[l][k] = std::atoi(header.substr(kol + 4, 3).c_str());
+      radix_line("\tVarbId(" << mVarbId[l][k] << ") chk_sum(" << chk_sum[l][k]
+                             << ")");
 
       kol += 8;
       nrec++;
@@ -629,6 +633,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
     matchingIndex.push_back(i);
   }
   float sfcp = 1013.0f;
+  float msle = 1013.15f;
   float sfct = 0.0f;
   int lp     = 0;
   std::vector<std::vector<float>> vdata(mvar);
@@ -698,6 +703,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
     // of the vertical coordinate system (mHeight = signma*scaling)
     if (varb.compare("PRSS") == 0) sfcp = vdata[kvar][ll];
     if (varb.compare("SHGT") == 0) sfct = vdata[kvar][ll];
+    if (varb.compare("MSLE") == 0) msle = vdata[kvar][ll];
 
     // load the winds for subsequent rotation to true
     if (varb.compare("UWND") == 0 || varb.compare("U10M") == 0)
@@ -723,7 +729,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
     int nvar = mNumVarb[ll];
     // default vertical motion units in mb/s
     for (size_t nn = 0; nn < mUnits.size(); ++nn)
-      if (0 == mVarb[ll].compare("WWND")) mUnits[nn] = "mb/h";
+      if (0 == mVarb[nn].compare("WWND")) mUnits[nn] = "mb/h";
 
     if (mHeader.z_flag == 1)
     {
@@ -797,7 +803,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
     if (hIt != columns.end())
     {
       results[ll][hIt - columns.begin()] =
-          hpaToAltitude(plevel) - surfaceAltitude;
+          hpaToAltitude(plevel, temp, msle) - surfaceAltitude;
     }
     // check for time
     auto timeIt = std::find(columns.begin(), columns.end(), "TIME");
@@ -835,7 +841,7 @@ std::vector<std::vector<float>> GFSFile::query(float lat, float lon, int month,
       if (varb.compare("RH2M") == 0) varb = "RELH";
       if (varb.compare("PRSS") == 0)
       {
-        surfaceAltitude = hpaToAltitude(vdata[kk][ll]) - 2.f;
+        surfaceAltitude = hpaToAltitude(vdata[kk][ll], temp, msle) - 2.f;
       }
       auto it = std::find(columns.begin(), columns.end(), varb);
       if (it != columns.end())

radixio/gfsfile.hh

@@ -62,7 +62,7 @@ class RADIX_PUBLIC GFSFile
     int hour;
     int ic;
     int il;
-    int cgrid;
+    std::string cgrid;
     std::string kvar;
     int nexp;
     float prec;
@@ -75,11 +75,11 @@ class RADIX_PUBLIC GFSFile
       hour  = std::atoi(val.substr(6, 2).c_str());
       ic    = std::atoi(val.substr(8, 2).c_str());
       il    = std::atoi(val.substr(10, 2).c_str());
-      cgrid = std::atoi(val.substr(12, 2).c_str());
+      cgrid = val.substr(12, 2).c_str();
       kvar  = val.substr(14, 4);
       nexp  = std::atoi(val.substr(18, 4).c_str());
-      prec  = (float)std::atof(val.substr(22, 14).c_str());
-      var1  = (float)std::atof(val.substr(36, 14).c_str());
+      prec  = float(std::atof(val.substr(22, 14).c_str()));
+      var1  = float(std::atof(val.substr(36, 14).c_str()));
     }
     float totalSeconds() { return total_seconds(year, month, day, hour); }
     std::string toString()
@@ -124,7 +124,7 @@ class RADIX_PUBLIC GFSFile
     int nz;
     int z_flag;
     int lenh;
-    void expand(std::string val)
+    void expand(std::string val, const Label &label)
     {
       model_id = val.substr(0, 4);
       icx      = std::atoi(val.substr(4, 3).c_str());
@@ -146,6 +146,15 @@ class RADIX_PUBLIC GFSFile
       nz       = std::atoi(val.substr(99, 3).c_str());
       z_flag   = std::atoi(val.substr(102, 2).c_str());
       lenh     = std::atoi(val.substr(104, 4).c_str());
+
+      int knx = ord(label.cgrid[0]);
+      int kny = ord(label.cgrid[1]);
+      // Check for the grid domain extending beyond 3 digits
+      if (knx >= 64 || kny >= 64)
+      {
+        nx = (knx - 64) * 1000 + nx;
+        ny = (kny - 64) * 1000 + ny;
+      }
     }
     std::string toString()
     {

radixmath/tests/tstUtil.cc

@@ -8,9 +8,9 @@ using namespace radix;
 TEST(radix, hpaToAltitude)
 {
   Real altitude = hpaToAltitude(835);
-  EXPECT_NEAR(1601.6289, altitude, 1e-1);
+  EXPECT_NEAR(1601.9551358843219, altitude, 1e-1);
   altitude = hpaToAltitude(67.7);
-  EXPECT_NEAR(17830.553, altitude, 1e-1);
+  EXPECT_NEAR(17835.433862090933, altitude, 1e-1);
 }
 
 TEST(radix, cspanf)

radixmath/util.cc

@@ -156,9 +156,9 @@ Real cspanf(Real value, Real begin, Real end)
   }
 }
 
-Real hpaToAltitude(Real hpa)
+Real hpaToAltitude(Real hpa, Real temperature, Real msle)
 {
-  return 0.3048 * (1 - std::pow(hpa / 1013.25, 0.190284)) * 145366.45;
+  return ((1 - std::pow(hpa / msle, 0.19022256039)) * temperature) / 0.0065;
 }
 
 /** Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1)
@@ -227,4 +227,5 @@ double exponentialIntegral(double d)
 
   return numer / denom / (d * std::exp(d));
 }
+
 }  // namespace radix

radixmath/util.hh

@@ -20,9 +20,12 @@ namespace radix
 /**
  * @brief hpaToAltitude converts hectorpascals or millibars to altitude (meters)
  * @param hpa hectorpascals in millibars
+ * @param temperature in kelvin (default=288.15)
+ * @param msle mean sea level pressure (default=1013.25)
  * @return altitude in meters
  */
-Real RADIX_PUBLIC hpaToAltitude(Real hpa);
+Real RADIX_PUBLIC hpaToAltitude(Real hpa, Real temperature = 288.15,
+                                Real msle = 1013.25);
 /**
  * @brief cspanf returns a value in the interval (begin,end]
  * This function is used to reduce periodic variables to a standard range.