Updating gfsfile to support cll2xy grid lookups.

cmake/Dependencies.cmake

 TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
    SUBPACKAGES_DIRS_CLASSIFICATIONS_OPTREQS
      bug  radixbug      SS OPTIONAL
-     io   radixio       SS OPTIONAL
      math radixmath     SS OPTIONAL
+     io   radixio       SS OPTIONAL
      geometry radixgeometry SS OPTIONAL
      plot radixplot     SS OPTIONAL
      widgets radixwidgets SS OPTIONAL

radixio/cmake/Dependencies.cmake

 TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
-LIB_REQUIRED_PACKAGES
+LIB_REQUIRED_PACKAGES radixmath
 LIB_OPTIONAL_PACKAGES
 TEST_REQUIRED_PACKAGES testframework
 TEST_OPTIONAL_PACKAGES

radixio/gfsfile.cc

@@ -2,6 +2,7 @@
 
 #include "radixio/eafstream.hh"
 #include "radixbug/bug.hh"
+#include "radixmath/util.hh"
 
 namespace radix
 {
@@ -194,8 +195,8 @@ GFSFile::GFSFile(std::string file)
                << " " << mLabel.hour;
             mEndTime = ss.str();
         }
-//        radix_tagged_line("Profile time: " << mProfiles[mProfiles.size()-1]
-//                << mLabel.toString());
+        //        radix_tagged_line("Profile time: " << mProfiles[mProfiles.size()-1]
+        //                << mLabel.toString());
     }
     rstream->close();
     delete rstream;
@@ -230,11 +231,72 @@ std::pair<float, float> GFSFile::gbl2xy(float clat
     radix_tagged_line("\tcomputed x =" << result.first);
     return result;
 }
+
+std::pair<float, float> GFSFile::cll2xy(float clat, float clon) const
+{
+    std::pair<float, float> result;
+    float almost1 = .99999;
+    float gamma = std::sin(mHeader.tang_lat*PI_BELOW_180);
+    float dlat = clat;
+    float strcmp2 = cspanf(mHeader.ref_lon, -180, 180);
+    float dlong = cspanf(clon - strcmp2, -180, 180);
+    dlong = dlong * PI_ON_180;
+    float gdlong = gamma * dlong;
+    float csdgam = 0.0, sndgam = 0.0;
+    if(std::abs(gdlong) < 0.01)
+    {
+        // For gamma small or zero. avoids round-off error or division
+        // by zero in the case of mercator or near-mercator projections.
+        gdlong = gdlong * gdlong;
+        sndgam = dlong * (1.-1./6. * gdlong *
+                          (1.-1./20. * gdlong *
+                           (1.-1./42. * gdlong)));
+        csdgam = dlong * dlong * .5 *
+                (1.-1./12. * gdlong *
+                 (1.-1./30. * gdlong *
+                  (1.-1./56. * gdlong)));
+    } else
+    {
+        sndgam = std::sin(gdlong)/gamma;
+        csdgam = (1.-std::cos(gdlong))/gamma/gamma;
+    }
+    float slat = std::sin(dlat*PI_ON_180);
+    if((slat >= almost1) || (slat <= -almost1))
+    {
+        result.first = 0.0;
+        result.second = 1./gamma;
+    }
+    float mercy = .5 * std::log((1.+slat)/(1.-slat));
+    float gmercy = gamma * mercy;
+    float rhog1 = 0;
+    if(std::abs(gmercy) < .001)
+    {
+        // For gamma small or zero. avoids round-off error or division
+        // by zero in the case of mercator or near-mercator projections.
+        rhog1 = mercy * ( 1. -.5 * gmercy *
+                          (1.-1./3. * gmercy *
+                           (1.-1./4. * gmercy)));
+    } else
+    {
+        rhog1 = (1. - std::exp(-gmercy)) / gamma;
+    }
+    result.first = (1.-gamma*rhog1)*sndgam;
+    result.second = rhog1 + (1.-gamma*rhog1)*gamma*csdgam;
+
+    return result;
+}
 std::pair<int,int> GFSFile::nearestPoint(float lat, float lon) const
 {
-    std::pair<float,float> point = gbl2xy(lat, lon
-                                          , mHeader.sync_lat, mHeader.ref_lat
-                                          , mHeader.sync_lon, mHeader.ref_lon);
+    std::pair<float,float> point;
+    if(mHeader.size == 0)
+    {
+        point = gbl2xy(lat, lon
+                       , mHeader.sync_lat, mHeader.ref_lat
+                       , mHeader.sync_lon, mHeader.ref_lon);
+    } else
+    {
+        point = cll2xy(lat, lon);
+    }
     std::pair<int,int> ipoint;
     ipoint.first =  (int)std::round(point.first);
     ipoint.second = (int)std::round(point.second);
@@ -268,8 +330,8 @@ std::vector<std::vector<float>> GFSFile::query(float lat
     // assume class was correctly initialized
     // get the grid points for the lon, lat in the met file
     std::pair<int,int> point = nearestPoint(lat, lon); //gbl2xy(lat, lon
-                                          //, mHeader.sync_lat, mHeader.ref_lat
-                                          //, mHeader.sync_lon, mHeader.ref_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);

radixio/gfsfile.hh

@@ -170,7 +170,7 @@ protected:
     Label mLabel;
     Header mHeader;
     /**
-    *  \brief gbl2xy converts lat&lon to xy coordinates for grid lookup
+    *  \brief gbl2xy converts global lat&lon to xy coordinates for grid lookup
     */
     std::pair<float,float> gbl2xy(float clat
                                   , float clon
@@ -179,6 +179,14 @@ protected:
                                   , float sync_lon
                                   , float ref_lon) const;
 
+    /**
+     * @brief cll2xy converts geographic coordinates lat&lon to xy coordinates for grid lookup
+     * @param clat real latitude
+     * @param clon real longitude
+     * @return std::pair<float,float>
+     */
+    std::pair<float,float> cll2xy(float clat, float clon) const;
+
 public:
     /**
      * @brief GFSFile Interface around an on disk GFS file

radixmath/constants.hh

@@ -26,6 +26,7 @@ namespace radix
     const Real minReal = std::numeric_limits<Real>::min();
     const Real PI = 2 * std::acos( Real( 0 ) );//3.1415926535897932384;
     const Real TWO_PI = PI * 2;//6.2831853071795864769;
+    const Real PI_BELOW_180 = 180 / PI;
     const Real PI_ON_180 = PI / 180;//0.0174532925199432957;
     const Real invPI = 1 / PI;//0.3183098861837906715;
     const Real invTWO_PI = 1 / TWO_PI;//0.1591549430918953358;

radixmath/tests/tstUtil.cc

@@ -30,6 +30,21 @@ TEST(radix, cspanf)
     lat = -180;
     v = cspanf(lat, -90, 90);
     EXPECT_FLOAT_EQ(v, 0);
+    lat = -180.5;
+    v = cspanf(lat, -90, 90);
+    EXPECT_FLOAT_EQ(v, -0.5);
+    lat = 5.5;
+    v = cspanf(lat, -90, 90);
+    EXPECT_FLOAT_EQ(v, 5.5);
+    lat = -5.5;
+    v = cspanf(lat, -90, 90);
+    EXPECT_FLOAT_EQ(v, -5.5);
+    lat = -90;
+    v = cspanf(lat, -90, 90);
+    EXPECT_FLOAT_EQ(v, 90);
+    lat = 90;
+    v = cspanf(lat, -90, 90);
+    EXPECT_FLOAT_EQ(v, 90);
 }
 
 TEST(radix, greatCircleDistance){