WIP ARL format read. Need to develop a smaller file for inclusion in tests

radixio/arldatastream.cc

@@ -2,6 +2,8 @@
 
 #include "radixbug/bug.hh"
 #include "radixcore/stringfunctions.hh"
+#include "radixio/eafstream.hh"
+#include "radixmath/util.hh"
 
 #include <fstream>
 namespace radix
@@ -10,15 +12,16 @@ class ARLDataStream::PImpl
 {
  public:
   std::string file;
-  std::shared_ptr<std::fstream> stream;
+  std::shared_ptr<radix::eafstream> stream;
+  int recordSize = 0;
 };
 
 void ARLDataStream::expand(const std::string& val, ARLIndexHeader& index)
 {
-  if (val.size() < 50)
+  if (val.size() < indexHeaderLength)
   {
     throw std::runtime_error(
-        "Incorrect size for string expandsion to ARLIndexHeader.");
+        "Incorrect size for string expansion to ARLIndexHeader.");
   }
   index.year  = std::atoi(val.substr(0, 2).c_str());
   index.month = std::atoi(val.substr(2, 2).c_str());
@@ -33,9 +36,9 @@ void ARLDataStream::expand(const std::string& val, ARLIndexHeader& index)
   index.var1  = float(std::atof(val.substr(36, 14).c_str()));
 }
 void ARLDataStream::expand(const std::string& val, const ARLIndexHeader& index,
-                           ARLHeader& header)
+                           ARLRecordHeader& header)
 {
-  if (val.size() < 108)
+  if (val.size() < recordHeaderLength)
   {
     throw std::runtime_error(
         "Incorrect size for string expansion to ARLHeader.");
@@ -74,7 +77,111 @@ void ARLDataStream::expand(const std::string& val, const ARLIndexHeader& index,
 ARLDataStream::ARLDataStream(const std::string& file)
     : p(new PImpl(), [](PImpl* impl) { delete impl; })
 {
-  p->file = file;
+  p->file   = file;
+  p->stream = std::make_shared<radix::eafstream>(
+      p->file.c_str(), std::ios::in | std::ios::binary);
+  if (!p->stream->is_open())
+  {
+    radix_line("Failed to open ARL file: " << file);
+  }
+}
+
+bool ARLDataStream::read_index_header(ARLIndexHeader& iheader)
+{
+  std::string headerString = p->stream->readString(indexHeaderLength);
+  expand(headerString, iheader);
+
+  radix_line("Read index header: year = " << iheader.year);
+  radix_line("                  month = " << iheader.month);
+  radix_line("                    day = " << iheader.day);
+  radix_line("                   hour = " << iheader.hour);
+  radix_line("               variable = " << iheader.kvar);
+  radix_line("       scaling exponent = " << iheader.nexp);
+  radix_line("          initial value = " << iheader.var1);
+  radix_line(headerString);
+  return true;
+}
+
+bool ARLDataStream::read_next_index_header(ARLIndexHeader& iheader)
+{
+  int bytesToSkip = roundUpInt(p->stream->bytesRead(), p->recordSize) -
+                    p->stream->bytesRead();
+
+  p->stream->skipBytes(bytesToSkip);
+
+  return read_index_header(iheader);
+}
+
+bool ARLDataStream::read_record_header(const ARLIndexHeader& iheader,
+                                       ARLRecordHeader& rheader)
+{
+  std::string headerString = p->stream->readString(recordHeaderLength);
+  expand(headerString, iheader, rheader);
+
+  // Calculate and save record size
+  p->recordSize = (rheader.nx * rheader.ny) + indexHeaderLength;
+
+  radix_line("  Read record header: model = " << rheader.model_id);
+  radix_line("                         nx = " << rheader.nx);
+  radix_line("                         ny = " << rheader.ny);
+  radix_line("                         nz = " << rheader.nz);
+  radix_line("        Size of each record = " << p->recordSize);
+  radix_line(headerString);
+  return true;
+}
+
+bool ARLDataStream::read_record(const ARLIndexHeader& iheader,
+                                const ARLRecordHeader& rheader,
+                                std::vector<std::vector<double> >& record)
+{
+  // Set up the vector size
+  record.resize(rheader.nx);
+  for (std::vector<double>& vec : record)
+  {
+    vec.resize(rheader.ny);
+  }
+
+  // Calculate the scaling factor
+  double scaleFactor = pow(2.0, 7.0 - (double)iheader.nexp), lastValue = 0.0;
+
+  radix_line("    Reading record data:");
+  for (int y = 0; y < rheader.ny; ++y)
+  {
+    for (int x = 0; x < rheader.nx; ++x)
+    {
+      // Read the raw value
+      unsigned char ch = (unsigned char)p->stream->readChar();
+      int packedValue  = ch - 127;
+      //      int packedValue = p->stream->readChar() - 127;
+
+      // Calculate the unpacked value
+      double unpackedValue = 0.0;
+      // Get the correct 'last value' if we are at a 0 index
+      if (x == 0)
+      {
+        if (y == 0)
+        {
+          lastValue = iheader.var1;
+        }
+        else
+        {
+          lastValue = record[x][y - 1];
+        }
+      }
+
+      unpackedValue = (packedValue / scaleFactor) + lastValue;
+      record[x][y]  = unpackedValue;
+      lastValue     = unpackedValue;
+
+      if ((x < 2 && y < 2) || ((rheader.nx - x < 3) && (rheader.ny - y < 3)))
+      {
+        radix_line("      [" << x << "," << y << "]: packed = " << packedValue
+                             << ", unpacked = " << unpackedValue);
+      }
+    }
+  }
+
+  return true;
 }
 
 /**

radixio/arldatastream.hh

@@ -16,7 +16,7 @@ namespace radix
 //
 // Forward declare ARLRecord and ARLHeader
 struct ARLIndexHeader;
-struct ARLHeader;
+struct ARLRecordHeader;
 class ARLRecord;
 
 class RADIX_PUBLIC ARLDataStream
@@ -27,6 +27,8 @@ class RADIX_PUBLIC ARLDataStream
   std::unique_ptr<PImpl, void (*)(PImpl*)> p;
 
  private:
+  static const int indexHeaderLength = 50, recordHeaderLength = 108;
+
   /**
    * @brief expand Expands string of 50 characters to an ARLIndexHeader
    * @param val String for expansion
@@ -42,7 +44,7 @@ class RADIX_PUBLIC ARLDataStream
    * @param header ARLHeader
    */
   static void expand(const std::string& val, const ARLIndexHeader& index,
-                     ARLHeader& header);
+                     ARLRecordHeader& header);
 
  public:
   ARLDataStream() = delete;
@@ -54,22 +56,32 @@ class RADIX_PUBLIC ARLDataStream
    * @return
    */
   bool read_index_header(ARLIndexHeader& iheader);
+  /**
+   * @brief next_index_header Sets the stream position to that of the next index
+   * header and reads that header
+   * @return
+   */
+  bool read_next_index_header(ARLIndexHeader& iheader);
+
+  bool read_record_header(const ARLIndexHeader& iheader,
+                          ARLRecordHeader& rheader);
 
-  bool read_record_header(ARLHeader& rheader);
   /**
    * @brief read_record Reads a single record from the stream
    * @param ARLRecord& record
    *
    * @return bool on if inventory was populated
    */
-  bool read_record(ARLRecord& record);
+  bool read_record(const ARLIndexHeader& iheader,
+                   const ARLRecordHeader& rheader,
+                   std::vector<std::vector<double> >& record);
 
   // TODO
   // bool write_record(const ARLRecord& record);
 
 };  // class ARLDataStream
 
-struct RADIX_PUBLIC ARLHeader
+struct RADIX_PUBLIC ARLRecordHeader
 {
  public:
   bool prime;
@@ -128,4 +140,6 @@ class RADIX_PUBLIC ARLRecord
 };  // ARLRecord
 }  // namespace radix
 
+//#include "arldatastream.i.hh"
+
 #endif /** RADIX_RADIXIO_ARLDATASTREAM_HH_ */

radixio/eafstream.cc

@@ -620,6 +620,25 @@ void eafstream::writeString(const std::string &var)
   mBytesWritten += var.size();
 }
 
+char eafstream::readChar()
+{
+  char ch;
+  this->read(&ch, 1);
+
+  if (this->bad() || this->eof())
+  {
+    std::stringstream ss;
+    ss << "Error, cannot read the string from file"
+       << ((this->bad() != false) ? "- this->bad()" : "- this->eof()");
+    std::cerr << ss.str() << std::endl;
+    throw std::logic_error(ss.str());
+  }
+
+  mBytesRead++;
+
+  return ch;
+}
+
 void eafstream::writeString(const std::string &var, size_t length, char filler)
 {
   std::string copy(var);

radixio/eafstream.hh

@@ -102,6 +102,8 @@ class RADIX_PUBLIC eafstream : public std::fstream
   void writeString(const std::string &var);
   void writeString(const std::string &var, size_t length, char filler = ' ');
 
+  char readChar();
+
   /**
    * \brief Returns the number of bytes read since the last header marker
    * \return int

radixio/tests/CMakeLists.txt

@@ -7,3 +7,4 @@ ADD_GOOGLE_TEST(tstF71Stream.cc NP 1)
 ADD_GOOGLE_TEST(tstDecayDb.cc NP 1)
 ADD_GOOGLE_TEST(tstHysplitCDump.cc NP 1)
 ADD_GOOGLE_TEST(tstSpectrum.cc NP 1)
+ADD_GOOGLE_TEST(tstArlDataStream.cc NP 1)

radixio/tests/tstArlDataStream.cc

+#include "gtest/gtest.h"
+
+#include <vector>
+
+#include "radixio/arldatastream.hh"
+
+using namespace radix;
+
+TEST(RadixIO, ReadArlStream)
+{
+  std::string testFile = "/home/ohp/data/met/RP195111.gbl";
+  ARLDataStream testStream(testFile);
+
+  // Read the initial index header
+  ARLIndexHeader testTopIndexHeader;
+  bool success = testStream.read_index_header(testTopIndexHeader);
+  ASSERT_TRUE(success);
+  // Test the content
+  EXPECT_EQ(51, testTopIndexHeader.year);
+  EXPECT_EQ(11, testTopIndexHeader.month);
+  EXPECT_EQ(1, testTopIndexHeader.day);
+  EXPECT_EQ(0, testTopIndexHeader.hour);
+  EXPECT_STREQ("INDX", testTopIndexHeader.kvar.c_str());
+
+  ARLRecordHeader testRecordHeader;
+  success = testStream.read_record_header(testTopIndexHeader, testRecordHeader);
+  ASSERT_TRUE(success);
+  // Test the content
+  EXPECT_EQ(144, testRecordHeader.nx);
+  EXPECT_EQ(73, testRecordHeader.ny);
+  EXPECT_EQ(18, testRecordHeader.nz);
+
+  // Read the data starting from the next index header
+  ARLIndexHeader testIndexHeader;
+  success = testStream.read_next_index_header(testIndexHeader);
+  ASSERT_TRUE(success);
+  // Test the content
+  EXPECT_EQ(51, testIndexHeader.year);
+  EXPECT_EQ(11, testIndexHeader.month);
+  EXPECT_EQ(1, testIndexHeader.day);
+  EXPECT_EQ(0, testIndexHeader.hour);
+  EXPECT_STREQ("PRSS", testIndexHeader.kvar.c_str());
+  // Read data
+  std::vector<std::vector<double>> testRecord;
+  success =
+      testStream.read_record(testIndexHeader, testRecordHeader, testRecord);
+}

radixmath/tests/tstUtil.cc

@@ -675,3 +675,38 @@ TEST(radix, SpecificHumidityPoint25)
       specificHumidityToRelativeHumidity(specificHumidity, temp3, pressure3);
   ASSERT_NEAR(expectRelHum3, relHum, expectRelHum3 * tolerance);
 }
+
+TEST(radix, RoundUpInt)
+{
+  int testIntValue = 17, testIntMultiple = 15;
+  EXPECT_EQ(30, roundUpInt(testIntValue, testIntMultiple));
+
+  testIntValue = 3, testIntMultiple = 3;
+  EXPECT_EQ(3, roundUpInt(testIntValue, testIntMultiple));
+
+  testIntValue = 17, testIntMultiple = 0;
+  EXPECT_EQ(17, roundUpInt(testIntValue, testIntMultiple));
+
+  testIntValue = 17, testIntMultiple = 250;
+  EXPECT_EQ(250, roundUpInt(testIntValue, testIntMultiple));
+
+  testIntValue = 0, testIntMultiple = 12;
+  EXPECT_EQ(0, roundUpInt(testIntValue, testIntMultiple));
+}
+TEST(radix, RoundDownInt)
+{
+  int testIntValue = 17, testIntMultiple = 15;
+  EXPECT_EQ(15, roundDownInt(testIntValue, testIntMultiple));
+
+  testIntValue = 3, testIntMultiple = 3;
+  EXPECT_EQ(3, roundDownInt(testIntValue, testIntMultiple));
+
+  testIntValue = 17, testIntMultiple = 0;
+  EXPECT_EQ(17, roundDownInt(testIntValue, testIntMultiple));
+
+  testIntValue = 502, testIntMultiple = 16;
+  EXPECT_EQ(496, roundDownInt(testIntValue, testIntMultiple));
+
+  testIntValue = 0, testIntMultiple = 12;
+  EXPECT_EQ(0, roundDownInt(testIntValue, testIntMultiple));
+}

radixmath/util.hh

@@ -260,6 +260,29 @@ double RADIX_PUBLIC gammaRayAttenuationInAir(double energy,
  * @return double - the exponential integral of the given value
  */
 double RADIX_PUBLIC exponentialIntegral(double d);
+
+inline int RADIX_PUBLIC roundUpInt(int value, int multiple)
+{
+  if (multiple == 0)
+  {
+    return value;
+  }
+
+  int remainder = value % multiple;
+  return (remainder == 0 ? value : value + multiple - remainder);
+}
+
+inline int RADIX_PUBLIC roundDownInt(int value, int multiple)
+{
+  if (multiple == 0)
+  {
+    return value;
+  }
+
+  int remainder = value % multiple;
+  return (remainder == 0 ? value : value - remainder);
+}
+
 }  // namespace radix
 
 #include "radixmath/util.i.hh"