Cdump utils

.gitlab-ci.yml

@@ -145,11 +145,7 @@ linux_analysis:
         -D ExperimentalBuild
         -D ExperimentalTest
         -D ExperimentalMemCheck
-    - gcovr --filter=".*/radix.*"
-        --xml-pretty
-        --print-summary
-        -o coverage.xml
-        -r $CI_PROJECT_DIR
+    - cmake --build . --target gcovr_xml
   coverage: /^\s*lines:\s*\d+.\d+\%/
   allow_failure: true
 

CMakeLists.txt

@@ -38,6 +38,31 @@ IF (NOT TRIBITS_PROCESSING_PACKAGE)
   IF(ENABLE_PYTHON_WRAPPERS AND NOT BUILD_SHARED_LIBS)
      MESSAGE(FATAL_ERROR "ENABLE_PYTHON_WRAPPERS=ON requires BUILD_SHARED_LIBS=ON")
   ENDIF()
+  ##---------------------------------------------------------------------------##
+  ## Code coverage testing
+  ##---------------------------------------------------------------------------##
+  if(${PROJECT_NAME}_ENABLE_COVERAGE_TESTING)
+    message(STATUS "Enabling coverage testing")
+    find_program(gcovr_program
+            NAMES gcovr)
+    if(NOT gcovr_program)
+      MESSAGE(STATUS "Could not find gcovr. Build targets disabled.")
+    else()
+      message(STATUS "gcovr found '${govr_program}'")
+      set(gcovr_filter "\".*/${PROJECT_NAME}.*\"")
+      message(STATUS "Adding gcovr_xml target")
+      add_custom_target(gcovr_xml
+              COMMAND ${gcovr_program}
+              --filter=${gcovr_filter} --xml-pretty --print-summary -o coverage.xml
+              -r "${PROJECT_SOURCE_DIR}")
+              #gcovr --filter=".*/src.*" --xml-pretty --print-summary -o coverage.xml -r $CI_PROJECT_DIR
+      message(STATUS "Adding gcovr_html target")
+      add_custom_target(gcovr_html
+              COMMAND ${gcovr_program}
+              --filter=${gcovr_filter} --html -o coverage.html
+              -r "${PROJECT_SOURCE_DIR}")
+    endif()
+  endif()
   #
   # For windows with BUILD_SHARED_LIBS we must use CMAKE_RUNTIME_OUTPUT_DIRECTORY
   # to place all *dll and *exe in the same directory so unit tests will work

radixio/hysplitcdump.cc

 
 #include "radixio/hysplitcdump.hh"
 
+#include "radixcore/stringfunctions.hh"
+#include "radixmath/util.hh"
+
 namespace radix
 {
+std::string HysplitCDump::fortran_fill(const std::string &value, size_t len)
+{
+  std::string v;
+  v.resize(len);
+  size_t vlen = value.length();
+  for (size_t i = 0; i < len; ++i) v[i] = '_';
+  for (size_t i = 0; i < vlen && i < len; ++i)
+  {
+    v[i] = value[i];
+  }
+  return v;
+}
+
 std::string HysplitCDump::id() const { return mId; }
 
 void HysplitCDump::setId(const std::string &id) { mId = id; }
@@ -392,4 +408,173 @@ HysplitPoint HysplitCDump::point(int timeIndex, int pollutantIndex,
   return p;
 }
 
+std::string HysplitCDump::encode_mass_name(const std::string &name)
+{
+  char label[5];
+  std::string v = fortran_fill(name, 4);
+  snprintf(label, sizeof(label), "C%s", v.c_str());
+  label[4] = '\0';  // null terminate
+  return label;
+}
+
+std::string HysplitCDump::encode_dilution_factor_name(const std::string &name)
+{
+  char label[5];
+  std::string v = fortran_fill(name, 4);
+  snprintf(label, sizeof(label), "D%s", v.c_str());
+  label[4] = '\0';  // null terminate
+  return label;
+}
+
+std::string HysplitCDump::encode_eqf_name(int masschain)
+{
+  char label[5];
+  snprintf(label, sizeof(label), "E%03d", masschain);
+  label[4] = '\0';  // null terminate
+  return label;
+}
+
+std::string HysplitCDump::encode_eqf_ratio_name(int masschainA, int masschainB)
+{
+  char label[5];
+  sprintf(label, "R%c%c_", radix::itoa(masschainA), radix::itoa(masschainB));
+  label[4] = '\0';  // null terminate
+  return label;
+}
+
+std::string HysplitCDump::encode_exposure_name(HysplitUnit unit)
+{
+  char label[5];
+  sprintf(label, "X%c__", radix::itoa(int(unit)));
+  label[4] = '\0';  // null terminate
+  return label;
+}
+
+std::string HysplitCDump::encode_activity_name(HysplitUnit unit)
+{
+  char label[5];
+  label[0] = 'A';
+  if (HysplitUnit::CI == unit)
+  {
+    label[1] = 'C';
+    label[2] = 'i';
+    label[3] = '_';  // Fortran binary doesn't like whitespace
+  }
+  else if (HysplitUnit::MICROCI == unit)
+  {
+    label[1] = 'u';
+    label[2] = 'C';
+    label[3] = 'i';
+  }
+  else
+  {
+    radix_not_implemented("HYSPLIT activity unit not implemented.");
+  }
+  label[4] = '\0';  // null terminate
+  return label;
+}
+
+std::string HysplitCDump::decode_name(const std::string &name)
+{
+  if (name.size() == 0) return std::string();
+  char key = name[0];
+  std::string result;
+  switch (key)
+  {
+    case '%':
+      // Probability of exeedance
+      result = "Probability of exceedance";
+      break;
+    case 'C':
+      // mass concentration
+      result = "Mass (g)";
+      break;
+    case 'D':
+      result = "Mass Dilution Factor";
+      break;
+    case 'A':
+      result = std::string("Activity (");
+      if ('C' == name[1])
+      {
+        result = result + std::string(name.substr(1, 2));
+      }
+      else if ('u' == name[1])
+      {
+        result = result + name.substr(1);
+      }
+      else
+      {
+        radix_not_implemented("Decode of activity unit not implemented.");
+      }
+      result = result + ")";
+      break;
+    case 'X':  // exposure
+      result = std::string("Exposure Rate");
+      if (name.size() > 1)
+      {  // convert unit
+
+        HysplitUnit unit = static_cast<HysplitUnit>(radix::ordinal(name[1]));
+        if (HysplitUnit::REM_HR == unit)
+        {
+          result = result + std::string(" (rem/hr)");
+        }
+        else if (HysplitUnit::R_HR == unit)
+        {
+          result = result + std::string(" (R/hr)");
+        }
+        else
+        {
+          radix_not_implemented(
+              "Decode of exposure rate unit not implemented.");
+        }
+      }
+      break;
+    case 'E':  // equivalent fissions
+      result = "Equivalent Fissions (";
+      if (name.size() > 1)
+      {
+        // capture the mass chain
+        int masschain = std::atoi(name.substr(1).c_str());
+        result        = result + std::to_string(masschain);
+      }
+      result = result + std::string(")");
+      break;
+    case 'R':  // equivalent fission ratio
+      result = "Equivalent Fission Ratio (";
+      if (name.size() > 1)
+      {
+        // capture the mass chain
+        int masschain = radix::ordinal(name[1]);
+        result        = result + std::to_string(masschain);
+      }
+      result = result + std::string("/");
+      if (name.size() > 2)
+      {
+        // capture the mass chain
+        int masschain = radix::ordinal(name[2]);
+        result        = result + std::to_string(masschain);
+      }
+      result = result + std::string(")");
+      break;
+    default:
+      result = name;
+  }
+  return result;
+}
+
+double HysplitCDump::cell_area(short yi) const
+{
+  double area = 0.;
+
+  double lat1;
+  double cLat = cornerLat();
+  double dLat = deltaLat();
+  double dLon = deltaLon();
+  // determine cell center latitude
+  lat1                = cLat + dLat * double(yi);
+  double cos_lat_dgpr = std::cos(radix::toRadians(lat1));
+  area = dLat * dLon * cos_lat_dgpr * std::pow(111198.41730306273, 2.);
+  return area;
+}
+
 }  // namespace radix

radixio/hysplitcdump.hh

@@ -33,6 +33,17 @@ struct RADIX_PUBLIC HysplitTime
   }
 };  // struct Time
 
+enum class HysplitUnit
+{
+  MASS  = -1,
+  DMASS = 0,
+  R_HR  = 1,
+  REM_HR,
+  EQF,
+  CI,
+  MICROCI
+};
+
 struct RADIX_PUBLIC HysplitPoint
 {
   short xi;
@@ -46,6 +57,7 @@ class RADIX_PUBLIC HysplitCDump
   typedef std::shared_ptr<HysplitCDump> SP;
 
  private:
+  static std::string fortran_fill(const std::string &value, size_t len);
   // identification
   std::string mId;
   // initial time
@@ -150,6 +162,62 @@ class RADIX_PUBLIC HysplitCDump
              float &concentration) const;
   // for python bindings
   HysplitPoint point(int ti, int pi, int li, int nzi) const;
+  /**
+   * @brief encode_mass_name encode mass pollutant into 4 characters
+   * Prepends decode character of 'C' and first 3 characters of name
+   * @param name
+   * @return
+   */
+  static std::string encode_mass_name(const std::string &name);
+  /**
+   * @brief encode_dilution_factor_name encode mass dilution factor into 4
+   * characters Prepends decode character of 'D' and first 3 characters of name
+   * @param name
+   * @return
+   */
+  static std::string encode_dilution_factor_name(const std::string &name);
+  /**
+   * @brief encode_eqf_name encode equivalient fission pollutant into 4
+   * characters Prepends decode character of 'E' and an ascii encoding of
+   * masschain
+   * @param name
+   * @return
+   */
+  static std::string encode_eqf_name(int masschain);
+  /**
+   * @brief encode_eqf_ratio_name encode equivalient fission ratio pollutant
+   * into 4 characters Prepends decode character of 'R' and an ascii encoding of
+   * masschains
+   * @param name
+   * @return
+   */
+  static std::string encode_eqf_ratio_name(int masschainA, int masschainB);
+  /**
+   * @brief encode_exposure_name encode mass dilution factor into 4
+   * characters Prepends decode character of 'X' and integer to ascii of unit
+   * @return
+   */
+  static std::string encode_exposure_name(HysplitUnit unit);
+  /**
+   * @brief encode_activity_name encode mass dilution factor into 4
+   * characters Prepends decode character of 'A' and unit
+   * @return
+   */
+  static std::string encode_activity_name(HysplitUnit unit);
+  /**
+   * @brief decode_name Decode pollutant names into verbose names
+   * @param name
+   * @return
+   */
+  static std::string decode_name(const std::string &name);
+
+  /**
+   * @brief cell_area Calculate the area for a specific grid cell latitude
+   * index
+   * @param yi latitude index
+   * @return grid cell area in meters^2
+   */
+  double cell_area(short yi) const;
 };  // class HysplitCDump
 /**
  * @class HysplitCDump

radixio/tests/test_radixio.py

@@ -212,6 +212,83 @@ class TestRadixIOBindings(unittest.TestCase):
                         self.assertEqual(p1.yi, p2.yi)
                         self.assertEqual(p1.concentration, p2.concentration)
 
+    def test_hysplit_cdump_encoding(self):
+        self.assertEqual("Cmas", radixio.HysplitCDump.encode_mass_name("mass"))
+        self.assertEqual("Cnam", radixio.HysplitCDump.encode_mass_name("name"))
+        self.assertEqual("Cm__", radixio.HysplitCDump.encode_mass_name("m"))
+        self.assertEqual("Dmas", radixio.HysplitCDump.encode_dilution_factor_name("mass"))
+        self.assertEqual("Dm__", radixio.HysplitCDump.encode_dilution_factor_name("m"))
+        self.assertEqual("E099", radixio.HysplitCDump.encode_eqf_name(99))
+        self.assertEqual("E140", radixio.HysplitCDump.encode_eqf_name(140))
+        # following three are not aligning because of unicode encoding I think
+        #self.assertEqual(r'''Rc\udc8c_''', radixio.HysplitCDump.encode_eqf_ratio_name(99, 140))
+        #self.assertEqual(r'''X\x1__''', radixio.HysplitCDump.encode_exposure_name(radixio.HysplitUnit_R_HR))
+        #self.assertEqual(r'''X\x2__''', radixio.HysplitCDump.encode_exposure_name(radixio.HysplitUnit_REM_HR))
+        self.assertEqual("ACi_", radixio.HysplitCDump.encode_activity_name(radixio.HysplitUnit_CI))
+        self.assertEqual("AuCi", radixio.HysplitCDump.encode_activity_name(radixio.HysplitUnit_MICROCI))
+
+    def test_hysplit_cdump_cell_area_north(self):
+        print("Testing hysplit cdump cell area northern hemisphere...")
+        blessed = [
+            12365088010.706079, 12177234539.621647, 11619381959.028811,
+            10708480337.301855, 9472206959.2784901, 7948125365.965446,
+            6182544005.3530407, 4229109173.6562023, 2147174999.7503204,
+            215800541.59562227,  539501.3700674492]
+
+        cvf = radixio.HysplitCDump()
+        cvf.setCornerLat(0.)
+        cvf.setDeltaLat(1.)
+        cvf.setDeltaLon(1.)
+
+        increment = 10
+        values = []
+        for yi in range(0, 90, increment):
+          area = cvf.cell_area(yi)
+          values.append(area)
+
+        area = cvf.cell_area(89)
+        values.append(area)
+        # change to a 0.05 degree delta
+        cvf.setCornerLat(89.)
+        cvf.setDeltaLat(0.05)
+        cvf.setDeltaLon(0.05)
+        area = cvf.cell_area(0)
+        values.append(area)
+        self.assertEqual(len(blessed), len(values))
+        for i in range(len(blessed)):
+          self.assertEqual(blessed[i], values[i])
+
+    def test_hysplit_cdump_cell_area_south(self):
+        print("Testing hysplit cdump cell area southern hemisphere...")
+        blessed = [
+             539501.3700674492,  215800541.59562227,  2147174999.7503204,
+             4229109173.6562023, 6182544005.3530407, 7948125365.965446,
+             9472206959.2784901, 10708480337.301855, 11619381959.028811,
+             12177234539.621647, 12365088010.706079]
+
+        values = []
+        cvf = radixio.HysplitCDump()
+         # change to a 0.05 degree delta
+        cvf.setCornerLat(-89.)
+        cvf.setDeltaLat(0.05)
+        cvf.setDeltaLon(0.05)
+        values.append(cvf.cell_area(0))
+
+        cvf.setCornerLat(-89.)
+        cvf.setDeltaLat(1)
+        cvf.setDeltaLon(1)
+
+        values.append(cvf.cell_area(0))
+        cvf.setCornerLat(-80)
+        increment = 10
+        for yi in range(0,90,increment):
+          area = cvf.cell_area(yi)
+          values.append(area)
+
+        self.assertEqual(len(blessed), len(values))
+        for i in range(len(blessed)):
+          self.assertEqual(blessed[i], values[i])
+
     @unittest.skip("origen.rev04.mpdkxgam.data not in the repo")
     def test_emitdb_xgam(self):
         print("Testing XGamDb...")

radixio/tests/tstHysplitCDump.cc

 
+#include <numeric>
+
 #include "gtest/gtest.h"
 #include "radixio/hysplitcdump.hh"
 
 using namespace radix;
+TEST(RadioIO, GridCellAreaPositive)
+{
+  std::vector<double> blessed{
+      12365088010.706079, 12177234539.621647, 11619381959.028811,
+      10708480337.301855, 9472206959.2784901, 7948125366.00,
+      6182544005.3530407, 4229109173.6562023, 2147174999.7503204,
+      215800541.60,       539501.37};
+  typedef HysplitCDump CVF;
+  CVF::SP cvf = std::make_shared<CVF>();
+  cvf->setCornerLat(0.f);
+  cvf->setDeltaLat(1.f);
+  cvf->setDeltaLon(1.f);
+
+  short increment = 10;
+  std::vector<double> values;
+  for (short yi = 0; yi < 90; yi += increment)
+  {
+    double area = cvf->cell_area(yi);
+    values.push_back(area);
+    std::cout << "Latitude \t" << (cvf->cornerLat() + (yi * cvf->deltaLat()))
+              << "\t area \t" << std::setprecision(10) << area << std::endl;
+  }
+  double area = cvf->cell_area(89);  // cLat + 89 * deltaLat
+  std::cout << "Latitude \t" << 89 << "\t area \t" << std::setprecision(10)
+            << area << std::endl;
+  values.push_back(area);
+  // change to a 0.05 degree delta
+  cvf->setCornerLat(89.f);
+  cvf->setDeltaLat(0.05f);
+  cvf->setDeltaLon(0.05f);
+  area = cvf->cell_area(0);
+  std::cout << "Latitude \t" << 89 << "\t area \t" << std::setprecision(10)
+            << area << std::endl;
+  values.push_back(area);
+  ASSERT_EQ(blessed.size(), values.size());
+  for (size_t i = 0; i < blessed.size(); ++i)
+  {
+    EXPECT_NEAR(blessed[i], values[i], 1e-1);
+  }
+  double total_blessed = std::accumulate(blessed.begin(), blessed.end(), 0.);
+  double total_value   = std::accumulate(values.begin(), values.end(), 0.);
+  EXPECT_NEAR(total_blessed, total_value, 1e-1);
+
+  std::cout << "Calculated/Blessed = " << (100. * (total_value / total_blessed))
+            << "%" << std::endl;
+}
+TEST(APTool, GridCellAreaNegative)
+{
+  std::vector<double> blessed{
+      539501.3700674492,  215800541.59562227, 2147174999.7503204,
+      4229109173.6562023, 6182544005.3530407, 7948125365.965446,
+      9472206959.2784901, 10708480337.301855, 11619381959.028811,
+      12177234539.621647, 12365088010.706079};
+  typedef HysplitCDump CVF;
+  CVF::SP cvf = std::make_shared<CVF>();
+  std::vector<double> values;
+  // change to a 0.05 degree delta
+  cvf->setCornerLat(-89.f);
+  cvf->setDeltaLat(0.05f);
+  cvf->setDeltaLon(0.05f);
+  double area = cvf->cell_area(0);
+  std::cout << "Latitude \t" << cvf->cornerLat() << "\t area \t"
+            << std::setprecision(10) << area << std::endl;
+  values.push_back(area);
 
+  cvf->setCornerLat(-89.);
+  cvf->setDeltaLat(1.f);
+  cvf->setDeltaLon(1.f);
+
+  area = cvf->cell_area(0);
+  values.push_back(area);
+  std::cout << "Latitude \t" << (cvf->cornerLat() + (0 * cvf->deltaLat()))
+            << "\t area \t" << std::setprecision(10) << area << std::endl;
+  cvf->setCornerLat(-80);
+  short increment = 10;
+  for (short yi = 0; yi < 90; yi += increment)
+  {
+    area = cvf->cell_area(yi);
+    values.push_back(area);
+    std::cout << "Latitude \t" << (cvf->cornerLat() + (yi * cvf->deltaLat()))
+              << "\t area \t" << std::setprecision(10) << area << std::endl;
+  }
+  ASSERT_EQ(blessed.size(), values.size());
+  for (size_t i = 0; i < blessed.size(); ++i)
+  {
+    EXPECT_NEAR(blessed[i], values[i], 1e-1);
+  }
+
+  double total_blessed = std::accumulate(blessed.begin(), blessed.end(), 0.);
+  double total_value   = std::accumulate(values.begin(), values.end(), 0.);
+  EXPECT_NEAR(total_blessed, total_value, 1e-1);
+  std::cout << "Calculated/Blessed = " << (100. * (total_value / total_blessed))
+            << "%" << std::endl;
+}
+TEST(RadixIO, HysplitCDumpEncoding)
+{
+  EXPECT_EQ("Cmas", HysplitCDump::encode_mass_name("mass"));
+  EXPECT_EQ("Cnam", HysplitCDump::encode_mass_name("name"));
+  EXPECT_EQ("Cm__", HysplitCDump::encode_mass_name("m"));
+  EXPECT_EQ("Dmas", HysplitCDump::encode_dilution_factor_name("mass"));
+  EXPECT_EQ("Dm__", HysplitCDump::encode_dilution_factor_name("m"));
+  EXPECT_EQ("E099", HysplitCDump::encode_eqf_name(99));
+  EXPECT_EQ("E140", HysplitCDump::encode_eqf_name(140));
+  EXPECT_EQ("Rc\x8C_", HysplitCDump::encode_eqf_ratio_name(99, 140));
+  EXPECT_EQ("X\x1__", HysplitCDump::encode_exposure_name(HysplitUnit::R_HR));
+  EXPECT_EQ("X\x2__", HysplitCDump::encode_exposure_name(HysplitUnit::REM_HR));
+  EXPECT_EQ("ACi_", HysplitCDump::encode_activity_name(HysplitUnit::CI));
+  EXPECT_EQ("AuCi", HysplitCDump::encode_activity_name(HysplitUnit::MICROCI));
+}
 TEST(RadixIO, HysplitCDump)
 {
   std::string file("cdump.bin");