Refactored radixsnd2arl to use radixmath/interpolate

radixio/examples/radixsnd2arl.cc

@@ -4,6 +4,7 @@
  */
 
 #include <algorithm>
+#include <cmath>
 #include <ctime>
 #include <iostream>
 #include <string>
@@ -12,6 +13,7 @@
 #include "radixcommand/commandline.hh"
 #include "radixcore/stringfunctions.hh"
 #include "radixmath/constants.hh"
+#include "radixmath/interpolate.hh"
 #include "radixmath/util.hh"
 
 #include "radixio/arldatastream.hh"
@@ -37,137 +39,6 @@ void addHour(struct tm *time, int hours)
  * definitely a better phrase for this)
  * @param valuesToInterpolate
  */
-bool interpolateValues(const std::vector<float> &pressureValues,
-                       std::vector<float> &valuesToInterpolate,
-                       bool circular = false)
-{
-  float missingValue = -9999.f;
-
-  // Ensure pressure and other vector are same length for interpolation
-  if (pressureValues.size() != valuesToInterpolate.size())
-  {
-    std::cout << "Error! Pressure vector is not same size ("
-              << pressureValues.size() << ") as vector to interpolate ("
-              << valuesToInterpolate.size() << ") - can't use to interpolate"
-              << std::endl;
-    return false;
-  }
-
-  // Ensure we have no missing values in the pressure dataset
-  for (float pressure : pressureValues)
-  {
-    if (pressure == -9999.f)
-    {
-      std::cout << "Error! Pressure vector contains " << missingValue
-                << " values - can't use to interpolate" << std::endl;
-      std::cout << "  Pressure vector:";
-      for (float p : pressureValues)
-      {
-        std::cout << " " << p << ";";
-      }
-      std::cout << std::endl;
-      return false;
-    }
-  }
-
-  // Write out initial values
-  std::cout << "Interpolating " << valuesToInterpolate.size()
-            << " values; initial:" << std::endl
-            << "  ";
-  for (float f : valuesToInterpolate)
-  {
-    std::cout << f << " ";
-  }
-  std::cout << std::endl;
-
-  // Loop through vector first to find first/last non-missing indices
-  size_t firstIndex = valuesToInterpolate.size(), lastIndex = 0;
-  bool foundFirst = false;
-  for (size_t i = 0; i < valuesToInterpolate.size(); ++i)
-  {
-    if (valuesToInterpolate[i] != missingValue)
-    {
-      if (!foundFirst)
-      {
-        firstIndex = i;
-        foundFirst = true;
-      }
-      lastIndex = i;
-    }
-  }
-  // Fill in values before first and after last indices
-  std::cout << "  Found first (" << firstIndex << ") and last (" << lastIndex
-            << ") non-missing indices; filling in values before & after..."
-            << std::endl;
-  for (size_t i = 0; i < firstIndex; ++i)
-  {
-    valuesToInterpolate[i] = valuesToInterpolate[firstIndex];
-  }
-  for (size_t i = lastIndex + 1; i < valuesToInterpolate.size(); ++i)
-  {
-    valuesToInterpolate[i] = valuesToInterpolate[lastIndex];
-  }
-
-  // Fill in missing values in central parts of vector
-  std::cout << "  Filling in missing data in central part of vector..."
-            << std::endl;
-  for (size_t i = firstIndex + 1; i < lastIndex; ++i)
-  {
-    // Search for a missing value
-    size_t lastGood = i - 1;
-    if (valuesToInterpolate[i] == missingValue)
-    {
-      // Get the next good value
-      while (valuesToInterpolate[i] == missingValue)
-      {
-        i++;
-      }
-      size_t nextGood = i;
-      // Interpolate between the two good values
-      for (size_t j = lastGood + 1; j < nextGood; ++j)
-      {
-        float lastGoodValue = valuesToInterpolate[lastGood],
-              nextGoodValue = valuesToInterpolate[nextGood];
-
-        if (circular && fabs(lastGoodValue - nextGoodValue) > 180.0)
-        {
-          std::cout << "    Circular interpolation with distance > 180 "
-                       "degrees: performing correction"
-                    << std::endl;
-          if (lastGoodValue < nextGoodValue)
-          {
-            lastGoodValue += 360.0;
-          }
-          else
-          {
-            nextGoodValue += 360.0;
-          }
-        }
-
-        valuesToInterpolate[j] =
-            lastGoodValue +
-            ((nextGoodValue - lastGoodValue) *
-             ((pressureValues[j] - pressureValues[lastGood]) /
-              (pressureValues[nextGood] - pressureValues[lastGood])));
-
-        if (circular)
-        {
-          valuesToInterpolate[j] = fmod(valuesToInterpolate[j], 360.0);
-        }
-      }
-    }
-  }
-
-  // Write out final values
-  std::cout << "Interpolation complete; final:" << std::endl << "  ";
-  for (float f : valuesToInterpolate)
-  {
-    std::cout << f << " ";
-  }
-  std::cout << std::endl;
-
-  return true;
-}
 
 int main(int argc, char **argv)
 {
@@ -521,23 +392,23 @@ int main(int argc, char **argv)
   // Interpolate the data to remove missing values
   std::cout << "Interpolating values to remove missing data..." << std::endl;
   std::cout << "Temperature:" << std::endl;
-  interpolateValues(inputPressures, inputTemps);
+  inputTemps = interpolateValues(inputPressures, inputTemps);
   if (usingRelHum)
   {
     std::cout << "Relative humidity:" << std::endl;
-    interpolateValues(inputPressures, inputRelHums);
+    inputRelHums = interpolateValues(inputPressures, inputRelHums);
   }
   if (usingDewPt)
   {
     std::cout << "Dew point:" << std::endl;
-    interpolateValues(inputPressures, inputDewPts);
+    inputDewPts = interpolateValues(inputPressures, inputDewPts);
   }
   std::cout << "Wind speed:" << std::endl;
-  interpolateValues(inputPressures, inputWSpds);
+  inputWSpds = interpolateValues(inputPressures, inputWSpds);
   std::cout << "Wind direction:" << std::endl;
-  interpolateValues(inputPressures, inputWDirs, true);
+  inputWDirs = interpolateValues(inputPressures, inputWDirs, true);
   std::cout << "Height:" << std::endl;
-  interpolateValues(inputPressures, inputHeights);
+  inputHeights = interpolateValues(inputPressures, inputHeights);
   std::cout << "Interpolation complete." << std::endl;
 
   // Convert the data into ARL format