Adding debug capabilities to gfsfile.

radixio/gfsfile.cc

 #include "radixio/gfsfile.hh"
 
 #include "radixio/eafstream.hh"
+#include "radixbug/bug.hh"
 
 namespace radix
 {
@@ -93,6 +94,7 @@ std::vector<float> GFSFile::mFact = { 1.0f, 1.0f, 1000.f, 1000.f, 1000.f, 60000.
 GFSFile::GFSFile(std::string file)
     : mFile(file)
 {
+    radix_tagged_line("GFSFile(" << file << ")" );
     // initialize data structure
     eafstream * rstream = new eafstream(file.c_str(), std::ifstream::in | std::ifstream::binary);
     std::string label = rstream->readString(50);
@@ -117,6 +119,7 @@ GFSFile::GFSFile(std::string file)
         header = recl.substr(50);
         mLabel.expand(label);
         mHeader.expand(header);
+        radix_tagged_line("Found grid: " << mHeader.toString());
     }
     int kol = 108;
     int nrec = nndx;
@@ -161,7 +164,6 @@ GFSFile::GFSFile(std::string file)
         }
         label = recl.substr(0,50);
         mLabel.expand(label);
-        header = recl.substr(50);
         mRecordTimes.push_back(mLabel.totalSeconds());
         if(!first_date_loaded)
         {
@@ -192,6 +194,8 @@ GFSFile::GFSFile(std::string file)
                << " " << mLabel.hour;
             mEndTime = ss.str();
         }
+//        radix_tagged_line("Profile time: " << mProfiles[mProfiles.size()-1]
+//                << mLabel.toString());
     }
     rstream->close();
     delete rstream;
@@ -203,11 +207,18 @@ std::pair<float, float> GFSFile::gbl2xy(float clat
                                         , float sync_lon
                                         , float ref_lon) const
 {
+    radix_tagged_line("gbl2xy("
+                      << clat << ","
+                      << clon << ","
+                      << sync_lat << ","
+                      << sync_lon << ","
+                      << ref_lon << ")");
     float tlat = clat;
     std::pair<float,float> result;
     if(tlat > 90.0f) tlat = 180.0f-tlat;
     if(tlat < -90.0f) tlat = -180.0f-tlat;
     result.second = 1.0f+(tlat-sync_lat)/ref_lat;
+    radix_tagged_line("\tcomputed y =" << result.second);
 
     float tlon = clon;
     // default PRIME section
@@ -216,6 +227,7 @@ std::pair<float, float> GFSFile::gbl2xy(float clat
     tlon = tlon-sync_lon;
     if(tlon < 0.0f) tlon = tlon+360.0f;
     result.first = 1.0f+tlon/ref_lon;
+    radix_tagged_line("\tcomputed x =" << result.first);
     return result;
 }
 std::pair<int,int> GFSFile::nearestPoint(float lat, float lon) const
@@ -255,11 +267,18 @@ std::vector<std::vector<float>> GFSFile::query(float lat
     float searchTime = total_seconds(year, month, day, hour);
     // assume class was correctly initialized
     // get the grid points for the lon, lat in the met file
-    std::pair<float,float> point = gbl2xy(lat, lon
-                                          , mHeader.sync_lat, mHeader.ref_lat
-                                          , mHeader.sync_lon, mHeader.ref_lon);
-    int x =  (int)std::round(point.first);
-    int y = (int)std::round(point.second);
+    std::pair<int,int> point = nearestPoint(lat, lon); //gbl2xy(lat, lon
+                                          //, mHeader.sync_lat, mHeader.ref_lat
+                                          //, mHeader.sync_lon, mHeader.ref_lon);
+
+    int x = point.first;//(int)std::round(point.first);
+    int y = point.second;//(int)std::round(point.second);
+    if(x < 0 || x >= mHeader.nx
+            || y < 0 || y >= mHeader.ny)
+    {
+        std::cerr << "Selected location is outside of file boundary." << std::endl;
+        return std::vector<std::vector<float>>();
+    }
     float minDelta = 99999999.f;
     size_t minIndex = 999999999;
     for(size_t i = 0; i < mRecordTimes.size(); ++i)

radixio/gfsfile.hh

@@ -39,9 +39,9 @@ std::vector<std::vector<float>> pakinp(const std::string& cvar
 class GFSFile {
 
 private:
-	static std::vector<std::string> mVarb;
-	static std::vector<std::string> mUnits;
-	static std::vector<float> mFact;
+    static std::vector<std::string> mVarb;
+    static std::vector<std::string> mUnits;
+    static std::vector<float> mFact;
     std::string mFile;
     // length of record
     int mLrec;
@@ -83,6 +83,20 @@ protected:
         float totalSeconds() {
             return total_seconds(year, month, day, hour);
         }
+        std::string toString()
+        {
+            std::stringstream ss;
+            ss << "year=" << year << " "
+               << "month=" << month << " "
+               << "day=" << day << " "
+               << "ic=" << ic << " "
+               << "il=" << il << " "
+               << "cgrid=" << cgrid << " "
+               << "kvar=" << kvar << " "
+               << "nexp=" << nexp << " "
+               << "var1=" << var1;
+            return ss.str();
+        }
     }; // class Label
     class Header {
     public:
@@ -128,6 +142,29 @@ protected:
             z_flag = std::atoi(val.substr(102,2).c_str());
             lenh = std::atoi(val.substr(104,4).c_str());
         }
+        std::string toString()
+        {
+            std::stringstream ss;
+            ss << "icx=" << icx << " "
+               << "mn=" << mn << " "
+               << "pole_lat=" << pole_lat << " "
+               << "pole_lon=" << pole_lon << " "
+               << "ref_lon=" << ref_lon << " "
+               << "ref_lat=" << ref_lat << " "
+               << "size=" << size << " "
+               << "orient=" << orient << " "
+               << "sync_xp=" << sync_xp << " "
+               << "sync_yp=" << sync_yp << " "
+               << "sync_lat=" << sync_lat << " "
+               << "sync_lon=" << sync_lon << " "
+               << "dummy=" << dummy << " "
+               << "nx=" << nx << " "
+               << "ny=" << ny << " "
+               << "nz=" << nz << " "
+               << "z_flag=" << z_flag << " "
+               << "lenh=" << lenh;
+            return ss.str();
+        }
     }; // class Header
 
     Label mLabel;